Orientation Dynamics of Sedimenting Anisotropic Particles in Turbulence

TL;DR

This study investigates how small anisotropic particles orient and settle in turbulence, revealing non-Gaussian orientation distributions and the influence of particle shape and flow parameters on settling behavior.

Contribution

It combines direct numerical simulations with theory to analyze the orientation dynamics and settling velocities of spheroids in turbulence, highlighting non-monotonic effects and deviations at finite Stokes numbers.

Findings

Orientation distributions are non-Gaussian and localized near broadside-on orientation.

Settling velocities depend non-monotonically on spheroid aspect ratio.

Preferential sweeping causes settling velocities to exceed orientational averages at higher Stokes numbers.

Abstract

We examine the dynamics of small anisotropic particles (spheroids) sedimenting through homogeneous isotropic turbulence using direct numerical simulations and theory. The gravity-induced inertial torque acting on sub-Kolmogorov spheroids leads to pronouncedly non-Gaussian orientation distributions localized about the broadside-on(to gravity) orientation. Orientation distributions and average settling velocities are obtained over a wide range of spheroid aspect ratios, Stokes and Froude numbers. Orientational moments from the simulations compare well with analytical predictions in the inertialess rapid-settling limit, with both exhibiting a non-monotonic dependence on spheroid aspect ratio. Deviations arise at Stokes numbers of order unity due to a spatially inhomogeneous particle concentration field resulting from a preferential sweeping effect; as a consequence, the time-averaged particle settling velocities exceed the orientationally averaged estimates.

Tables

Table 1. Table 1: Parameters for DNS runs in Figures 1,2 and 4 in the main manuscript

Parameters in pumir_18	$\kappa$	$F{r}_{\eta }=U_s/u_{\eta }$	$S{t}_{\eta }$	$R{e}_s$	$R{e}_{\dot{\gamma }}$	$L\left(\mu m\right)$	$L/\eta$
$R_{\lambda }=200$	0.05	4.2	0.4	0.8	0.036	107	0.22
	0.05	3.1	0.29	0.5	0.026	92	0.18
	0.05	1.05	0.1	0.1	0.01	54	0.1
$ϵ=0.51m^2{s}^{-3}$	0.02	3.1	0.296	0.8	0.066	145	0.29
	0.02	2.3	0.216	0.5	0.048	124	0.25
	0.02	0.78	0.074	0.1	0.016	73	0.15
$\nu =1.1\times {10}^{-5}{m}^2{s}^{-1}$	0.01	2.5	0.24	0.8	0.1	183	0.37
	0.01	1.8	0.17	0.5	0.076	157	0.32
	0.01	0.6	0.06	0.1	0.026	92	0.18
$R_{\lambda }=150$	0.05	3.7	0.4	0.8	0.04	107	0.22
	0.05	2.7	0.29	0.5	0.03	92	0.18
	0.05	0.9	0.1	0.1	0.01	54	0.1
$ϵ=0.0246m^2{s}^{-3}$	0.02	2.7	0.296	0.8	0.08	145	0.29
	0.02	2.0	0.216	0.5	0.06	124	0.25
	0.02	0.7	0.074	0.1	0.02	73	0.15
$\nu =1.1\times {10}^{-5}{m}^2{s}^{-1}$	0.01	2.1	0.24	0.8	0.14	183	0.37
	0.01	1.6	0.17	0.5	0.1	157	0.32
	0.01	0.5	0.06	0.1	0.03	92	0.18
$R_{\lambda }=96$	0.05	7.4	0.1	0.8	0.01	107	0.1
	0.05	5.4	0.074	0.5	0.008	92	0.09
	0.05	1.8	0.025	0.1	0.003	54	0.05
$ϵ=0.001562m^2{s}^{-3}$	0.02	5.4	0.075	0.8	0.02	145	0.15
	0.02	4.0	0.05	0.5	0.016	124	0.13
	0.02	1.4	0.019	0.1	0.005	73	0.07
$\nu =1.1\times {10}^{-5}{m}^2{s}^{-1}$	0.01	4.3	0.06	0.8	0.03	183	0.19
	0.01	3.2	0.04	0.5	0.02	157	0.16
	0.01	1.1	0.015	0.1	0.008	92	0.09
$R_{\lambda }=47$	0.05	13.4	0.025	0.8	0.00354	107	0.05
	0.05	9.8	0.019	0.5	0.0026	92	0.04
	0.05	3.4	0.006	0.1	0.00088	54	0.03
$ϵ=9.8\times {10}^{-5}{m}^2{s}^{-3}$	0.02	10	0.019	0.8	0.006	145	0.07
	0.02	7.3	0.014	0.5	0.0047	124	0.06
	0.02	2.5	0.0047	0.1	0.0016	73	0.04
$\nu =1.1\times {10}^{-5}{m}^2{s}^{-1}$	0.01	7.9	0.015	0.8	0.01	183	0.09
	0.01	5.8	0.011	0.5	0.0075	157	0.08
	0.01	2	0.0037	0.1	0.0026	92	0.05
	0.1	17	0.032	0.8	0.0022	85	0.043
	0.1	12.4	0.023	0.5	0.0016	73	0.037
	0.1	4.2	0.008	0.1	0.00056	43	0.02

Table 2. Table 2: Parameters for DNS runs in Figure 3 in the main manuscript

Parameters in pumir_18	$\kappa$	$F{r}_{\eta }=U_s/u_{\eta }$	$S{t}_{\eta }$	$R{e}_s$	$R{e}_{\dot{\gamma }}$	$L\left(\mu m\right)$	$L/\eta$
$R_{\lambda }=150$	0.999	5.5	0.64	0.36	0.004	33	0.07
	0.91	5.1	0.6	0.34	0.004	33	0.07
$ϵ=0.0246m^2{s}^{-3}$	0.67	3.9	0.46	0.26	0.004	33	0.07
	0.53	3.2	0.37	0.21	0.004	33	0.07
$\nu =1.1\times {10}^{-5}{m}^2{s}^{-1}$	0.1	0.64	0.07	0.04	0.004	33	0.07
	0.05	0.32	0.04	0.02	0.004	33	0.07
$R_{\lambda }=96$	0.999	11.9	0.16	0.4	0.001	33	0.03
	0.91	11	0.15	0.37	0.001	33	0.03
$ϵ=0.001562m^2{s}^{-3}$	0.67	8.5	0.11	0.28	0.001	33	0.03
	0.53	6.9	0.09	0.23	0.001	33	0.03
$\nu =1.1\times {10}^{-5}{m}^2{s}^{-1}$	0.1	1.4	0.019	0.05	0.001	33	0.03
	0.05	0.7	0.0095	0.02	0.001	33	0.03
$R_{\lambda }=47$	0.999	31.6	0.06	0.63	0.0004	39	0.02
	0.91	29.3	0.05	0.59	0.0004	39	0.02
$ϵ=9.8\times {10}^{-5}{m}^2{s}^{-3}$	0.67	22.5	0.04	0.45	0.0004	39	0.02
	0.53	18.3	0.03	0.36	0.0004	39	0.02
$\nu =1.1\times {10}^{-5}{m}^2{s}^{-1}$	0.1	3.7	0.007	0.07	0.0004	39	0.02
	0.05	1.9	0.003	0.04	0.0004	39	0.02

Equations

$$\frac{d{\mathbf{U}}_p}{dt}=\mathbf{g}+\frac{1}{{\tau }_p{X}_A}{\mathbf{M}}_t^{-1}\cdot (\mathbf{u}-{\mathbf{U}}_p),$$

$$\frac{d{\mathbf{\omega }}_p}{dt}+{\mathbf{I}}_p^{-1}\cdot [{\mathbf{\omega }}_p\wedge ({\mathbf{I}}_p\cdot {\mathbf{\omega }}_p)]=K_{sed}{\mathbf{I}}_p^{-1}\cdot [({\mathbf{M}}_t\cdot \hat{\mathbf{g}})\cdot \mathbf{p}({\mathbf{M}}_t\cdot \hat{\mathbf{g}})\wedge \mathbf{p}]$$

$$+8\pi \mu L^3{\mathbf{I}}_p^{-1}\cdot [{\mathbf{M}}_r^{-1}\cdot (\frac{1}{2}\mathbf{\Omega }-{\mathbf{\omega }}_p)-Y_H(\mathbf{E}\cdot \mathbf{p})\wedge \mathbf{p}],$$

$$\dot{\mathbf{p}}=$$

$$+\frac{1}{2}\mathbf{\Omega }\wedge \mathbf{p}+\frac{Y_H}{Y_c}[\mathbf{E}\cdot \mathbf{p}-\mathbf{E}:\mathbf{ppp}].$$

$$\mathbf{p}\cdot (\mathbf{I}-\hat{\mathbf{g}}\hat{\mathbf{g}})\approx \frac{1}{F{r}_{\eta }^2}\frac{8\pi Y_A{Y}_c{\tau }_{\eta }}{f_I(\kappa )X_A}(\mathbf{S}+\frac{Y_H}{Y_c}\mathbf{E})\cdot \mathbf{p}$$

$$⟨1-p_3^2⟩\approx \frac{32{\pi }^2{Y}_A^2{Y}_c^2}{f_I^2(\kappa )X_A^2}\left(\frac{1}{3}+\frac{Y_H^2}{5Y_c^2}\right)\frac{1}{F{r}_{\eta }^4}.$$

$$⟨(1-p_3{)}^2⟩\approx \frac{1}{4}(\frac{8\pi Y_A{Y}_c}{f_I(\kappa )X_A}{)}^4[M_1+B^2{M}_2+B^4{M}_3]\frac{1}{F{r}_{\eta }^8},$$

$$\frac{d{\mathbf{U}}_p}{dt}=\mathbf{g}+\frac{1}{{\tau }_p{X}_A}{\mathbf{M}}_t^{-1}\cdot (\mathbf{u}-{\mathbf{U}}_p),$$

$$\frac{d{\mathbf{\omega }}_p}{dt}+{\mathbf{I}}_p^{-1}\cdot [{\mathbf{\omega }}_p\wedge ({\mathbf{I}}_p\cdot {\mathbf{\omega }}_p)]=K_{sed}{\mathbf{I}}_p^{-1}\cdot [({\mathbf{M}}_t\cdot \hat{\mathbf{g}})\cdot \mathbf{p}({\mathbf{M}}_t\cdot \hat{\mathbf{g}})\wedge \mathbf{p}]$$

$$+8\pi \mu L^3{\mathbf{I}}_p^{-1}\cdot [{\mathbf{M}}_r^{-1}\cdot (\frac{1}{2}\mathbf{\Omega }-{\mathbf{\omega }}_p)-Y_H(\mathbf{E}\cdot \mathbf{p})\wedge \mathbf{p}].$$

$$X_A$$

$$Y_A$$

$$X_c$$

$$Y_c$$

$$X_A$$

$$Y_A$$

$$X_c$$

$$Y_c$$

$${\mathbf{I}}_p$$

$$f_I(\kappa )$$

$$+210(20{\kappa }^6-13{\kappa }^4-12{\kappa }^2+7)(1-{\kappa }^2{)}^{1/2}\kappa {\sin }^{-1}((1-{\kappa }^2{)}^{1/2})$$

$$-105(24{\kappa }^6-55{\kappa }^4+50{\kappa }^2-19){\sin }^{-1}((1-{\kappa }^2{)}^{1/2}{)}^2]/D,$$

$$f_I(\kappa )$$

$$+105(24{\kappa }^2-31)({\kappa }^2-1{)}^{1/2}{\kappa }^2{\tanh }^{-1}(\frac{({\kappa }^2-1{)}^{1/2}}{\kappa }{)}^2$$

$$+1995({\kappa }^2-1{)}^{1/2}{\tanh }^{-1}(\frac{({\kappa }^2-1{)}^{1/2}}{\kappa }{)}^2$$

$$+4757({\kappa }^2-1{)}^{1/2}{\kappa }^6-4282({\kappa }^2-1{)}^{1/2}{\kappa }^4$$

$$-210\kappa (20{\kappa }^6-13{\kappa }^4-12{\kappa }^2+7){\tanh }^{-1}(\frac{({\kappa }^2-1{)}^{1/2}}{\kappa })]/E,$$

$$\dot{\mathbf{p}}={\mathbf{\omega }}_p\wedge \mathbf{p},$$

$$K_{sed}(M_{(t)kl}{\hat{g}}_l{p}_k)({ϵ}_{imn}{M}_{(t)mj}{\hat{g}}_j{p}_n)=-8\pi \mu L^3{M}_{(r)ij}^{-1}{\omega }_{pj}^{sed}.$$

$${\omega }_{pj}^{sed}$$

$${\dot{p}}_i^{sed}$$

$${\dot{p}}_i^{jeff}$$

$${\dot{p}}_i$$

$$p_1$$

$$p_2$$

$$p_3$$

$$⟨p_1^2⟩\approx \frac{32{\pi }^2{Y}_A^2{Y}_c^2}{f_I^2(\kappa )X_A^2}(⟨S_{ji}{S}_{km}⟩+\frac{Y_H^2}{Y_c^2}⟨E_{ji}{E}_{km}⟩)p_j{p}_k{\delta }_{i1}{\delta }_{m1}\frac{1}{F{r}_{\eta }^4},$$

Full text

Suspended inertial anisotropic particles show up in a variety of scenarios ranging from pollen dispersion to soot emission. Prominent examples in nature include ice crystals suspended in high-altitude Cirrus clouds which are a crucial element in the planetary greenhouse effect Liou ; Pandit2015 . The radiative properties of such clouds depend sensitively on the orientation distribution of ice crystals baran . The latter come in a variety of pristine shapes with sizes ranging from tens to thousands of microns review , smaller than the typical Kolmogorov scales, about a millimeter, for atmospheric turbulence. Therefore, a first step towards understanding Cirrus cloud radiation is to examine how sub-Kolmogorov anisotropic particles orient themselves while sedimenting in a turbulent flow.

The critical role of turbulence in gravitational settling has been investigated in-depth only for inertial spherical particles Yau2000 ; monchaux2012 ; mehlig2016 . In this simpler scenario, relevant to the dynamics of water droplets in warm clouds, for instance, we now have a detailed understanding of the role of turbulence in enhancing single-particle sedimentation maxey_87 ; maxey_93 ; BecRay2014 as well as collision Collins2016a ; Collins2016b ; Wang1 ; Wang2 ; Bodenschatz ; James and coalescence BecPRE rates which control raindrop formation falkovich2002 ; shaw2003 .

The effect of inertia for anisotropic particles is far more involved owing to additional rotational degrees of freedom Voth2017 . Most earlier studies ignore either inertia Pumir-NJP ; VincenziPRE (the suspended particles acting as probes for the turbulent velocity-gradient tensor Voth2017 ; Meneveau2011 ) or gravity RoyPRE ; Gupta . Experiments have also largely focussed on neutrally buoyant anisotropic tracers in turbulence Voth1 ; Voth2 ; Voth3 . Thus, gravitational settling of heavy anisotropic particles, beyond simple laminar flows under Stokesian conditions Maxey1 ; Maxey2 , remains largely unexplored Voth2017 . Recent efforts address the issue of how such particles sediment in non-trivial flows pumirmehlig2017 ; siewert1 ; siewert2 ; pumir , but the effect of gravity on rotational dynamics is not accounted for, leading to orientation distributions that are far from being representative. There exist efforts analyzing the motion of anisotropic particles in turbulent channel flow, the object of interest often being the particle deposition rate onto walls Zhao2006 ; Fan1995 ; Zhao2014 ; Mortensen2008 ; Coletti2019 . Gravity is omitted in most of these efforts; those that do include gravity again neglect its role in the rotational dynamics Mclaughlin2001 . In this work, using direct numerical simulations (DNSs) and theory, we characterize the distribution of particle orientations in suspensions of spheroids sedimenting in an ambient homogeneous isotropic turbulent field. Rigorously accounting for the effects of gravity on both the particle translational and rotational degrees of freedom, we find, in contrast to earlier efforts pumirmehlig2017 ; siewert1 ; siewert2 ; pumir , that the orientation distributions always peak at the broadside-on (to gravity) orientation. Further, although the particle settling velocities equal the orientationally averaged estimates in the rapid-settling limit, they consistently exceed the latter when effects of particle inertia become significant.

We perform direct numerical simulations of non-interacting spheroids sedimenting through homogeneous isotropic turbulence with a mass loading assumed small enough for carrier-fluid turbulence to remain unaffected(a one-way coupled framework)SM . The fluid velocity and pressure fields satisfy the incompressible Navier-Stokes equations for a fluid with density $\rho_{f}$ and kinematic viscosity $\nu$. Turbulence is maintained in a statistically stationary homogeneous isotropic state via injection of energy at the lowest wavenumbers ($1\leq k_{f}\leq 2$) samriddhi1 . The simulations are pseudospectral in space and involve a second-order Adams-Bashforth scheme for time marching. A spatial resolution of $512^{3}$ collocation points is used, with the choice of four different kinematic viscosities corresponding to Taylor-scale Reynolds numbers, $R_{\lambda}=u_{rms}^{2}\sqrt{15/\epsilon\nu}$, of $47$, $96$, $150$ and $200$ ($u_{rms}$ is the root-mean-square velocity and $\epsilon=2\nu\langle\bm{E}\!\!:\!\!\bm{E}\rangle$ is the averaged dissipation rate). For each $R_{\lambda}$, we follow the motion of $100000$ oblate (prolate) spheroids, with aspect ratios ($\kappa$) ranging from $0.1$ to $0.01$ ($10$ to $100$); here, $\kappa=a/b$, $a$ and $b$ being the semi-axis lengths along and orthogonal to the spheroid symmetry axis $\bm{p}$. The particles are initialized at random positions with their translational velocities set equal to fluid values and angular velocities set equal to those of anisotropic tracers Jeffery1922 at their locations. The initial orientations, as characterized by normalized quaternions evans77 , are uniformly distributed over the unit sphere. The simulations are run for $5-6$ integral-scale eddy turnover times, sufficient to attain a statistical steady state.

The equations governing the particle dynamics are:

[TABLE]

where $\bm{U}_{p}$ and $\bm{\omega}_{p}$ are the translational and angular velocities of the particles, ${\bm{g}}$ is the gravitational acceleration ($\hat{\bm{g}}$ being the corresponding unit vector), $L$ is the largest particle dimension and $\tau_{p}$ is the particle relaxation time(see SM ). $\bm{I}_{p}$ in Eq.2 is the moment of inertia tensor, while ${\bm{M}}_{t}$ and ${\bm{M}}_{r}$ denote the Stokesian translational and rotational mobility tensors for the spheroid, with $\bm{M}_{t(r)}=X_{A(C)}^{-1}(\kappa)\bm{pp}+Y_{A(C)}^{-1}(\kappa)(\bm{I}-\bm{pp})$, the principal resistance coefficients ($X_{A}-Y_{C}$) being well known functions of $\kappa$ Kimkarrila . The large particle-to-fluid density ratio ($\rho_{p}/\rho_{f}$), relevant to the atmospheric scenario, implies the neglect of Basset and added mass forces in Eq.1. The particle Reynolds numbers based on both the Kolmogorov shear rate ($\dot{\gamma}_{\eta}=(\epsilon/\nu)^{\frac{1}{2}}$) and the nominal slip velocity ($U_{s}=\tau_{p}g$) are assumed small ($Re_{\dot{\gamma}_{\eta}}=\dot{\gamma}_{\eta}L^{2}/\nu,Re_{s}=U_{s}L/\nu<1$), so particles are acted on, at leading order, by the sum of the gravitational force and quasi-steady Stokes drag proportional to the slip velocity footnote1 ; see SM . Since sub-Kolmogorov spheroids experience turbulence as a fluctuating linear flow, the Jeffery relation Kimkarrila ; Jeffery1922 is used for the turbulent torque in Eq.2 with the ratio $Y_{H}/Y_{C}=(\kappa^{2}-1)/(\kappa^{2}+1)$ being the Bretherton constant $B$ Bretherton1962 . Eq.2 includes, in addition, the gravity-induced torque acting to orient an anisotropic particle, sedimenting in a quiescent fluid at small but finite $Re_{s}$, broadside-on to gravity cox1965 ; khayatcox89 ; navaneeth2015 ; an expression for this torque was obtained in navaneeth2015 . The superposition of the gravity and shear-induced torques in Eq.2 has been used footnote2 earlier to determine the orientation dynamics of particles sedimenting through simple shear flow subkoch2005 ; subkoch2006 . The quantity $T\!R=\frac{K_{sed}}{\mu L^{3}\dot{\gamma}_{\eta}}\sim Fr_{\eta}^{2}f_{I}(\kappa)$ characterizes the relative magnitudes of these torques in Eq.2, where $K_{sed}=Re_{s}\mu U_{s}L^{2}f_{I}(\kappa)X_{A}^{2}$, with the aspect-ratio dependent function, $f_{I}(\kappa)$, having been obtained in navaneeth2015 , and $Fr_{\eta}=\tau_{p}g/u_{\eta}$ being the Froude number based on the Kolmogorov velocity scale (${\bm{u}}_{\eta}=(\nu\epsilon)^{\frac{1}{4}})$. In Eq.1 and Eq.2, $\bm{u}$, $\bm{\Omega}$ and $\bm{E}$ denote the undisturbed turbulent velocity, vorticity and rate-of-strain fields interpolated at the particle positions SM .

Apart from $R_{\lambda}$, $\kappa$ and $Fr_{\eta}$, the dynamics as governed by Eq.1 and Eq.2, on length scales of the order of the Kolmogorov scale ($l_{\eta}=(\nu^{3}/\epsilon)^{\frac{1}{4}}$) or smaller, is a function of the Kolmogorov Stokes number ($St_{\eta}=\tau_{p}/\tau_{\eta}$ with $\tau_{\eta}=\dot{\gamma}_{\eta}^{-1}$ the Kolmogorov time scale). Using parameters characteristic of the atmospheric scenario, including ice crystal sizes and turbulence dissipation rates from pumir , the simulations reported here correspond to $St_{\eta}\in(0.0037,0.4)$ and $Fr_{\eta}\in(0.5,17)$. For a given $R_{\lambda}$, the dynamics of the thinnest (disk-like) spheroids corresponds to the smallest Stokes and Froude numbers. The torque ratio, $T\!R$ ranges from $1-800$ for all ice crystal sizes and turbulence intensities considered here. Thus, the gravity-induced torque is expected to be dominant for typical ice clouds. This is borne out in Fig. 1 which shows the distribution of orientations $\cos\theta_{0}=|\hat{\bm{g}}\cdot{\bm{p}}|$ (since ${\bm{p}}$ and $-{\bm{p}}$ correspond to the same spheroid orientation, we take the modulus), obtained from our DNSs for (a) $R_{\lambda}=150$ and (b) $R_{\lambda}=47$. For each $R_{\lambda}$, we show results for oblate spheroids of different aspect ratios (see legend), both with($Re_{s}\!=\!0.8$) and without($Re_{s}\!=\!0$) the gravity-induced torque. The gravity-induced torque causes the distributions to be sharply localized about the broadside-on orientation ($\cos\theta_{0}=1$), especially for the smaller $R_{\lambda}$. In contrast, as emphasized in the insets of Fig. 1, neglect of this torque leads to distributions peaked at the longside-on orientation ($\cos\theta_{0}=0$ for oblate spheroids), although this maximum is quite shallow, consistent with earlier studies siewert1 ; siewert2 ; pumir . The continuous curves in Fig. 1 are a guide to the eye; the comparison with a Gaussian in Fig. 1(a) nevertheless conveys the pronouncedly non-Gaussian character of the distributions for $Re_{s}\!=\!0.8$.

Analytical progress is possible in the rapid-settling limit(henceforth, RST or ‘Rapid-Settling Theory’), $l_{\eta}/U_{s}\ll\tau_{\eta}$ or $Fr_{\eta}\gg 1$, when a particle settles through a Kolmogorov eddy much faster than the eddy decorrelates anu2019 (see SM ). Further, assuming $St_{\eta}\ll 1$, and neglecting the angular acceleration in Eq.2, the rate of change of spheroid orientation, $\dot{\bm{p}}={\bm{\omega}}_{p}\wedge\bm{p}$, is given by:

[TABLE]

As already seen, the torque ratio $T\!R\sim f_{I}(\kappa)Fr_{\eta}^{2}$ with $f_{I}(\kappa)\sim\mathcal{O}(1)$ for oblate spheroids. For large $Fr_{\eta}$, the weak turbulent shear only leads to small fluctuations about the broadside-on orientation. For such orientations, with $\hat{\bm{g}}\!=\!\bm{1}_{3}$, one has $\bm{p}\cdot\hat{\bm{g}}=p_{3}\approx 1$ and $p_{1,2}\ll 1$. Furthermore, the rotation rate of the nearly broadside-on spheroid, in any plane containing $\hat{\bm{g}}$, is asymptotically small since the gravity-induced torque vanishes for the broadside-on orientation. Thus, there is a near-balance between the $1$ and $2$ components of the turbulent and gravity-induced torques at leading order, the terms proportional to $\dot{p}_{1,2}$ in (S14) being $\mathcal{O}(Fr_{\eta}^{-2})$ smaller. This gives

[TABLE]

for the projection of the spheroid axis in the plane transverse to gravity; here $\bm{S}=\frac{1}{2}\bm{\epsilon}\cdot\bm{\Omega}$ is the vorticity tensor and $\bm{\epsilon}$ being the Levi-Civita symbol. The components $p_{1,2}$ transverse to gravity are linear functionals of the turbulent velocity gradient tensor. Turbulent velocity gradients are dominated by the smallest (Kolmogorov) scales, and are pronouncedly non-Gaussianfootnote4 ; hence the orientation distributions, in the rapid-settling limit, are non-Gaussian(characterized below via the second and fourth moments) despite the localization about the broadside-on orientation.

Since $p_{3}=\cos\theta_{0}\approx 1-\frac{\theta_{0}^{2}}{2}$ for $\theta_{0}\ll 1$, $\langle 1-p_{3}^{2}\rangle=\langle p_{1}^{2}+p_{2}^{2}\rangle\approx\langle\theta_{0}^{2}\rangle$ corresponds to the variance of the orientation distribution about the broadside-on orientation. With $p_{1,2}$ linear in ${\bm{E}}$ and ${\bm{S}}$, calculating $\langle 1-p_{3}^{2}\rangle$ requires the variance of the turbulent rate of strain and vorticity tensors over a particle settling trajectory. For $St_{\eta}\ll 1,Fr_{\eta}\gg 1$, one expects no preferential sampling and the average along a settling trajectory, $\langle\cdot\rangle$, above may be replaced by the usual fluid ensemble average Batchelor1953 . For homogeneous isotropic turbulence, the ensemble averages are: $\langle E_{ij}E_{kl}\rangle=\frac{\dot{\gamma}_{\eta}^{2}}{20}(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}-\frac{2}{3}\delta_{ij}\delta_{kl}),$ $\langle S_{ij}S_{kl}\rangle=\frac{\dot{\gamma}_{\eta}^{2}}{12}(\delta_{ik}\delta_{jl}-\delta_{il}\delta_{jk}),$ and $\langle S_{ij}E_{kl}\rangle=0$ BrunkKoch1997 ; PopeBook . Using these SM , one finds:

[TABLE]

Fig. 2(a) compares the DNS results for $\langle 1-p_{3}^{2}\rangle$ to Eq.5 and demonstrates the good agreement for large $Fr_{\eta}$, with deviations arising for $Fr_{\eta}$ of order unity and smaller, in which case $\langle 1-p_{3}^{2}\rangle$ approaches a plateau.

A more sensitive measure of the orientation distributions is $\langle(1-p_{3})^{2}\rangle$. For a distribution localized about the broadside-on orientation, $\langle(1-p_{3})^{2}\rangle\propto\langle\theta_{0}^{4}\rangle$, and is therefore a measure of the fourth moment. Proceeding along lines sketched above, $\langle(1-p_{3})^{2}\rangle\approx\frac{1}{4}\langle(p_{1}^{2}+p_{2}^{2})^{2}\rangle$ with $p_{1,2}$ as given above, and the calculation involves the fourth moment of the turbulent velocity gradient tensor SM . One obtains:

[TABLE]

where $M_{1}=\frac{3G_{1}}{2}+\frac{32G_{2}}{15}-\frac{162G_{3}}{15}$, $M_{2}=-3G_{1}+\frac{8G_{2}}{15}+\frac{54G_{3}}{5}$ and $M_{3}=\frac{3G_{1}}{2}$, with $G_{1}\!=\!\langle(\tau_{\eta}\partial{u_{1}}/\partial{x_{1}})^{4}\rangle$, $G_{2}\!=\!\langle(\tau_{\eta}\partial{u_{1}}/\partial{x_{2}})^{4}\rangle$ and $G_{3}\!=\!\langle\tau_{\eta}^{4}(\partial{u_{1}}/\partial{x_{1}})^{2}(\partial{u_{1}}/\partial{x_{2}})^{2}\rangle$ being the independent (non-dimensional) scalar components involving the fourth moment of the velocity gradient. Unlike the second moment, the pre-factor multiplying $Fr_{\eta}^{-8}$ is both a function of $\kappa$ and $R_{\lambda}$, the latter dependence arising from dissipation-range intermittency referred to above. Fig.  2(b) compares Eq.S24 with DNS results, the pattern of agreement being similar to that of the second moment above footnote3 . Since $\langle 1-p_{3}^{2}\rangle=\langle\theta_{o}^{2}\rangle$ and $\langle(1-p_{3})^{2}\rangle=\frac{1}{4}\langle\theta_{o}^{4}\rangle$ for large $Fr_{\eta}$, the ratio $\frac{4\langle(1-p_{3})^{2}\rangle}{3\langle 1-p_{3}^{2}\rangle^{2}}$, which is independent of $Fr_{\eta}$, characterizes the departure from Gaussianity. This ratio, which is unity for a Gaussian, is plotted as an inset in Fig. 2(b) for $\kappa\rightarrow 0$ (a flat disk); it is well above unity and increases with increasing $R_{\lambda}$. One therefore expects orientation distributions in the atmospheric case, with $R_{\lambda}$’s one to two orders of magnitude higher than those in our simulations  shaw2003 , to have similar variances but be significantly more intermittent.

In the inset of Fig. 3, we plot orientation distributions as a function of the spheroid aspect ratio, other physical parameters being fixed SM . Interestingly, the localization about the broadside-on orientation first increases as $\kappa$ increases from zero (a flat disk), attains a maximum, before decreasing again as $\kappa$ approaches unity. The non-monotonicity arises because the gravity-induced torque is small for both flat disks (due to the vanishingly small mass of such shapes) and near-spheres (since the torque scales with the square of the small eccentricity). The second moment from the RST framework, Eq.5​, can be rewritten to isolate the $\kappa$-dependence through a change of variable $Fr_{\eta}=Fr_{\eta,sph}*\frac{\kappa}{X_{A}}$, where $Fr_{\eta,sph}=\frac{2\rho_{p}L^{2}g}{9\mu u_{\eta}}$. The resulting $\kappa-$dependence is consistent with the above non-monotonicity; although, within the RST framework, $\langle 1-p_{3}^{2}\rangle\sim\mathcal{O}(\kappa^{-4})$ for $\kappa\rightarrow 0$ and $\langle 1-p_{3}^{2}\rangle\sim\mathcal{O}(\kappa-1)^{-2}$ for $\kappa\rightarrow 1$. Since $\langle 1-p_{3}^{2}\rangle\leq 1$, the divergences above betray a breakdown of the assumption of a localized distribution in the analysis. As shown in Fig. 3, the second moments from our DNS agree with Eq.5 for intermediate values of $\kappa$ (maximum localization of $\cos\theta_{0}$), but plateau in the aforementioned asymptotic limits (corresponding to a uniform distribution of $\cos\theta_{0}$). Overall, the disagreement with theory, expectedly, grows with increasing $R_{\lambda}$.

With increase in the turbulence intensity, $Fr_{\eta}$ decreases while $St_{\eta}$ increases to values of order unity. As already seen in Fig. 2, DNS results depart from RST predictions in this limit. A suspension of spherical particles in a turbulent flow is no longer spatially homogeneous when $Fr_{\eta},St_{\eta}\sim\mathcal{O}(1)$ BecRay2014 ; eatonReview ; Collins2016a ; Collins2016b . Preferential sampling of regions of low vorticity by inertial particles, together with a sweeping effect in presence of gravity, leads to enhanced settling velocities maxey_87 ; maxey_93 ; collins2014 . Fig. 4 shows this to be true for the suspensions of spheroids considered here. For large $Fr_{\eta}$, the time-averaged settling speeds(which scale linearly with $Fr_{\eta}$ on account of being proportional to the acceleration due to gravity) from the DNS agree with the orientational averages for $R_{\lambda}=47$ and $96$ (the $Fr_{\eta}$ required for this agreement increases with increasing $St_{\eta}$). For finite $Fr_{\eta}$ and $St_{\eta}$, the time averages consistently exceed the orientation averaged estimates due to the preferential sweeping effect SM .

In this letter, we have characterized the orientation distributions and settling speeds of spheroids in homogeneous isotropic turbulence. Orientation distributions are localized about the broadside-on (to gravity) orientation, but are pronouncedly non-Gaussian for parameters typical of the atmospheric scenario. This is in contrast to recent studies which neglect the gravity-induced torque, and predict distributions peaked at the longside-on orientation pumirmehlig2017 ; siewert1 ; siewert2 ; pumir . The non-Gaussian distributions found here are also in contrast to earlier analyses reliant on a Gaussian ansatz Klett1994 ; mehligNJP2019 . While the broadside-on peak has been captured in mehligNJP2019 , the simplistic Gaussian ansatz used for the velocity field, and the resulting Gaussian nature of the orientation fluctuations, is incorrect. Furthermore, mehligNJP2019 lacks any discussion on the spatial organization of the particles, and its effect on particle settling speeds. In contrast, we show that the particle concentration field remains homogeneous for $St_{\eta}\ll 1$; for $St_{\eta}\sim\mathcal{O}(1)$, preferential sweeping effects lead to a spatially inhomogeneous concentration and enhanced settling speeds (Fig. 4,SM ). Results for prolate spheroids (not shown) are similar to those discussed above. It would be of interest, in future, to characterize pair-level statistics for anisotropic particles in position-orientation space, as a step towards analyzing ice-water and ice-ice collision efficiencies; the latter thought of as crucial to explaining observed ice-crystal concentrations in mixed-phase clouds and relatively rapid snow-flake formation in ice clouds Khain1 ; Khain2 ; Khain3 .

I The spheroid equations of motion

The equations governing particle translation and rotation are given by:

[TABLE]

The Stokesian translational ($\bm{M}_{t}$) and rotational ($\bm{M}_{r}$) mobility tensors appearing in (S1) and (S2) characterize the viscous force and torque acting on the spheroid. For a spheroid, these are of the form $\bm{M}_{t(r)}=X_{A(C)}^{-1}(\kappa)\bm{pp}+Y_{A(C)}^{-1}(\kappa)(\bm{I}-\bm{pp})$, with the scalar aspect-ratio dependent resistance functions given as follows Kim :

• •

Oblate:

[TABLE]

• •

Prolate:

[TABLE]

Since sub-Kolmogorov spheroids experience turbulence as a fluctuating linear flow, the Jeffery relation Kim ; jeffery has been used for the turbulent shear torque in (S2) with the ratio $Y_{H}/Y_{C}=(\kappa^{2}-1)/(\kappa^{2}+1)$ being the Bretherton constant $B$ Bretherton_1962 . The particle relaxation time, $\tau_{p}$, appearing in equation (S1) is defined as $\tau_{p}=\dfrac{2\rho_{p}b^{2}\kappa}{9\mu X_{A}}$ for oblate($\kappa<1$) and $\tau_{p}=\dfrac{2\rho_{p}a^{2}}{9\mu\kappa^{2}X_{A}}$ for prolate($\kappa>1$) spheroids. The moment of inertia tensor in equation (S2) is defined by:

• •

Oblate:

[TABLE]

• •

Prolate:

[TABLE]

There is an additional and important torque contribution in (S2) due to gravity, arising from the effects of fluid inertia associated with a spheroid settling in an otherwise quiescent fluid. The dominant contribution to this torque is due to inertial forces acting in a region around the spheroid of order its own size, and the torque is therefore proportional to the sedimentation Reynolds number ($Re_{s}$) for $Re_{s}\ll 1$. In other words, the gravity-induced inertial torque emerges as a regular perturbation about the Stokesian limit. Therefore, to $O(Re_{s})$, the functional dependence of the torque on $\hat{\bm{g}}$ and ${\bm{p}}$, viz, the form $(\bm{M}_{t}\!\cdot\!\hat{\bm{g}})\!\cdot\!\bm{p}(\bm{M}_{t}\!\cdot\!\hat{\bm{g}})\!\wedge\!\bm{p}$ in (S2) may be readily inferred using symmetry arguments. The additional aspect ratio dependence, contained in the coefficient $K_{sed}$ in (S2), requires a detailed analysis. This calculation has been done in navaneeth_2015 , using a generalized reciprocal theorem formulation, and one obtains $K_{sed}=Re_{s}\mu U_{s}L^{2}f_{I}(\kappa)X_{A}^{2}$, with the aspect-ratio dependent inertial function, $f_{I}(\kappa)$ specified below:

• •

Oblate:

[TABLE]

where $D=315(\kappa(\kappa^{2}-1)-(1-\kappa^{2})^{1/2}(2\kappa^{2}-3)\sin^{-1}((1-\kappa^{2})^{1/2}))\\ (\kappa(\kappa^{2}-1)+(2\kappa^{2}-1)(1-\kappa^{2})^{1/2}\sin^{-1}((1-\kappa^{2})^{1/2}))^{2}$

• •

Prolate:

[TABLE]

where $E=315\kappa^{3}(\kappa(\kappa^{2}-1)^{1/2}+(1-2\kappa^{2})\tanh^{-1}(\frac{(\kappa^{2}-1)^{1/2}}{\kappa}))^{2}\\ (\kappa(\kappa^{2}-1)^{1/2}+(2\kappa^{2}-3)\tanh^{-1}(\frac{(\kappa^{2}-1)^{1/2}}{\kappa}))$.

On account of inertia being a regular perturbation, use of the generalized reciprocal theorem shows that the inertial correction to the viscous torque in (S2), in the limit $Re_{s},Re_{\dot{\gamma}_{\eta}}\ll 1$ ($Re_{\dot{\gamma}_{\eta}}$ being the Reynolds number based on the Kolmogorov shear rate), may be constructed as a linear superposition subkoch_2005 ; subkoch_2006 of a shear-induced contribution in the absence of gravity ($\sim F_{I}(\kappa)Re_{\dot{\gamma}_{\eta}}\mu L^{3}\dot{\gamma_{\eta}}$, see navaneeth_2016 ; navaneeth_2017 ) and the gravity-induced contribution above that neglects any ambient shear. The ratio of the turbulent shear-induced inertial to the gravity-induced torques then turns out to be $\mathcal{O}(Fr_{\eta}^{-2})(L/l_{\eta})^{2}\ll 1$. Thus, the shear-induced inertial torque may be neglected for the sub-Kolmogorov spheroids examined here, and only the gravity-induced contribution is therefore included in (S2).

It is worth noting that, in contrast to the torque problem, the inertial correction to the Stokes drag in (S1) is not a linear superposition of the gravity and shear contributions since inertial effects enter as a singular perturbation in this case. This is evident from the saffman_1965 and mclaughlin_1991 derivations for the inertial lift on sphere in a simple shear flow. Even in the limit $Re_{s},Re_{\dot{\gamma}_{\eta}}\ll 1$, the scaling and the direction of the inertial force is crucially dependent on the ratio of the two screening length ($\nu/U_{s}$ for sedimentation and $(\nu/\dot{\gamma}_{\eta})^{\frac{1}{2}}$ for shear). Importantly, however, both of these inertial screening lengths are much larger than the size of the (sub-Kolmogorov) spheroid, and the inertial corrections, although non-trivial, are nevertheless small in comparison to the Stokes drag in (S1). Non-linear corrections to the drag become important for particles comparable to the Kolmogorov scale; these corrections are not known in closed form, however, and one then needs fully resolved simulations schneider_2019 .

Based on the above system of equations, we perform direct numerical simulations of non-interacting spheroids sedimenting through homogeneous and isotropic turbulence with a mass loading assumed small enough for carrier-fluid turbulence to remain unaffected; that is to say, a one-way coupled framework. For mass loadings of order unity, one requires a two-way coupled framework horwitz_2016 ; bala_2017 , in which case the fluid velocity needs to be accurately determined at the particle positions in order to estimate the particle forces (which then act to modify the turbulence); this in turn involves rather subtle issues with regard to interpolation schemes horwitz_2016 ; bala_2017 .

II The Rapid-Settling Theory (RST)

II.1 Formulation

We turn our attention to equation (S2), which governs the rotational dynamics of the particles. Our objective is to calculate $\bm{\omega}_{p}$ in the limit where the particles settle rapidly through a Kolmogorov eddy in a time much smaller than the eddy decorrelation time; that is, $l_{\eta}/U_{s}\ll\tau_{\eta}\Rightarrow U_{s}/u_{\eta}$ or $Fr_{\eta}\gg 1$. We also assume that the Stokes number based on the Kolmogorov timescale, $St_{\eta}\ll 1$. The angular acceleration of the particles, in equation (S2), can then be neglected since it is $\mathcal{O}(St_{\eta}/Fr_{\eta}^{2})$ smaller than the gravity-induced torque and $\mathcal{O}(St_{\eta})$ smaller than the turbulent shear-induced torque. The rotation rate of a spheroid is given by:

[TABLE]

The contributions to $\bm{\omega}_{p}$, and hence $\dot{\bm{p}}$, are $\dot{\bm{p}}^{sed}$ from the gravity-induced torque and $\dot{\bm{p}}^{jeff}$ from the turbulent shear-induced torque.

To calculate $\dot{\bm{p}}^{sed}$, we use the Stokesian relation between the gravity-induced torque and the hydrodynamic torque acting on a spheroid rotating in a quiescent fluid, which yields:

[TABLE]

Using $\hat{\bm{g}}=-\bm{1}_{3}$ in the above expression, one obtains:

[TABLE]

Substituting in (S9),

[TABLE]

In the sub-Kolmogorov range, particles see the turbulence as a fluctuating linear flow. In the absence of inertial effects, the particles rotate with an angular velocity given by the Jeffery relation at leading order. This gives:

[TABLE]

where $B$ is the Bretherton’s constant, as defined earlier. Adding the gravity-induced (S12) and the Jeffery (S13) contributions, one obtains:

[TABLE]

For large $Fr_{\eta}$, the dominant gravity-induced torque implies that, for both oblate and prolate spheroids, the weak turbulent shear only leads to small fluctuations about the broadside-on orientation (for oblate spheroids, the broadside equilibrium corresponds to $p_{3}=1,p_{1,2}=0$; for prolate spheroids $p_{3}=0$). The rotation rate of a nearly broadside-on oblate spheroid($p_{1,2}\ll p_{3}\approx 1$), in any plane containing $\hat{\bm{g}}$, is asymptotically small, owing to the vicinity to the aforementioned gravity-induced equilibrium. Thus, in the $1$ and $2$ components of (S14) for oblate spheroids, there is a near-balance between the turbulent and gravity-induced torques at leading order; the terms proportional to $\dot{p}_{1,2}$ are $\mathcal{O}(p_{1,2}\dot{\gamma}_{\eta})$ with the Jeffery contribution being $\mathcal{O}(p_{3}\dot{\gamma}_{\eta})$. Since $p_{3}\approx 1$, and $p_{1,2}$ turn out to be $\mathcal{O}(Fr_{\eta}^{-2})$ (see (S18a) and (S18b) below), the unsteady terms may be neglected with an error of $\mathcal{O}(Fr_{\eta}^{-2})$. The term $\bm{E}:\bm{ppp}$ in (S14) is also ignored since it involves quadratic combinations of $p_{1}$ and $p_{2}$, all of which are asymptotically small in the rapid settling limit, as seen above. With these approximations, one obtains:

[TABLE]

for an oblate spheroid. Similarly, for a nearly broadside-on prolate spheroid ($p_{3}\ll p_{1,2}\approx\mathcal{O}(1)$), the turbulent and gravity-induced torques nearly balance, at leading order, in the $3^{rd}$ component of (S14), and one obtains:

[TABLE]

The moments of interest (with regard to characterizing the orientation distribution) are of the general form $\langle(1-p_{3}^{2})^{n}\rangle$ for an oblate spheroid. Note that the relations (S15) and (S16) above, which express $p_{1,2}$ as linear functionals of the turbulent velocity gradient, imply that the orientational moment of order $n$ requires the $2n^{th}$ moment of the turbulent velocity gradient. In the rapid settling limit, one expects no preferential sampling, and hence, the averages involved in the aforementioned moments, which are along a settling particle trajectory, may be replaced by the usual fluid ensemble averages(batchelor1953theory ). The first two moments $\langle 1-p_{3}^{2}\rangle$ and $\langle(1-p_{3}^{2})^{2}\rangle$, for an oblate spheroid, will be evaluated below. The rapid settling limit corresponds to a small-$\theta_{0}$ approximation ($\theta_{0}$ is the angle between the spheroid axis $\bm{p}$ and $\hat{\bm{g}}$), in which case $\langle 1-p_{3}^{2}\rangle\propto\langle\theta_{0}^{2}\rangle$ and $\langle(1-p_{3}^{2})^{2}\rangle\propto\langle\theta_{0}^{4}\rangle$. Thus, evaluating these will also allow one to characterize the departure of the distribution from Gaussianity.

II.2 The second moments of the turbulent velocity gradient and spheroid orientation

The second moment $\langle 1-p_{3}^{2}\rangle=\langle p_{1}^{2}+p_{2}^{2}\rangle$, and we therefore begin by squaring equations (S15), (S16) and then ensemble averaging. The latter eliminates terms proportional $\langle\bm{SE}\rangle$ from symmetry arguments pertaining to homogeneous isotropic turbulence. This gives:

[TABLE]

For an oblate spheroid, the largest terms in the double contraction in equation S18 are those that involve $p_{3}\,(\approx 1)$, and (S18a) and (S18b) reduce to.

[TABLE]

Evaluating (S19) requires the variance of the turbulent velocity gradient tensor which is given by Pope_Book :

[TABLE]

The expressions for the variance of the rate-of-strain tensor, $\bm{E}=\frac{\bm{\Gamma}+\bm{\Gamma}^{\dagger}}{2}$ and the vorticity tensor, $\bm{S}=\frac{\bm{\Gamma}-\bm{\Gamma}^{\dagger}}{2}$ can be obtained from S20. For example,

[TABLE]

Since all directions in the plane perpendicular to gravity are equivalent, one has $\langle S_{31}S_{31}\rangle=\langle S_{32}S_{32}\rangle$ and $\langle E_{31}E_{31}\rangle=\langle E_{32}E_{32}\rangle$. The explicit expressions for $\langle\bm{EE}\rangle$ and $\langle\bm{SS}\rangle$ have, in fact, already been given in brunk_koch , but we nevertheless calculate them using the full velocity-gradient tensor $\bm{\Gamma}$, since this serves as a prelude to the fourth moment derivation in Section II.4).

Adding equations S19a and S19b and using equation (S20) to calculate the ensemble averages, $\langle 1-p_{3}^{2}\rangle$ for an oblate spheroid turns out to be:

[TABLE]

Proceeding along similar lines, the second moment for a nearly broadside-on prolate spheroid, $\langle p_{3}^{2}\rangle$, is shown to be half of (S22).

II.3 Fourth-moment of the turbulent velocity gradient

In this section, we evaluate the orientational moment $\langle(1-p_{3})^{2}\rangle\approx\frac{1}{4}\langle(1-p_{3}^{2})^{2}\rangle$ which, as indicated above, is proportional to the fourth moment of the orientation distribution when localized about the broadside-on equilibrium. The calculation involves first obtaining the fourth moment of the turbulent velocity gradient, and this involves a rather elaborate effort. We use a graphical approach that allows substantial simplification of the algebra involved. In order to illustrate the approach, we derive expressions for both the third and fourth moments of the velocity gradient tensor, defined below:

• •

Third moment:

[TABLE]

• •

Fourth moment:

[TABLE]

Note that the result for the third moment of the velocity gradient tensor, $\Sigma_{ijkpqr}$, is already known (see page 206 in Pope_Book ), but is nevertheless rederived here for purposes of clarity.

The sixth-order tensor $\Sigma_{ijkpqr}$ and eighth-order tensor $\Pi_{ijkmsrqp}$ must evidently be isotropic, and therefore expressible in terms of tensor products of the Kronecker delta. The usual derivation involves writing down all possible permutations of the Kronecker deltas, although this becomes especially tedious for the fourth and higher moments. To circumvent the algebraic effort involved, we use an alternate method where the aforementioned permutations are represented as ‘graphs’, with terms corresponding to the same graphs being grouped together. Starting off with $\Sigma_{ijkpqr}$, we note that the graph of every term in $\Sigma_{ijkpqr}$ is composed of three lines, each of these connecting a pair of indices in (S23). Thus, each line corresponds to a Kronecker delta tensor in the final permutation sum. There are three types of lines:

• •

A vertical line connects indices belonging to the same partial derivative in $\Sigma_{ijkpqr}$. For example, since the indices ‘i’ and ‘p’ occur in the same partial derivative $\dfrac{\partial{u_{i}}}{\partial{x_{p}}}$, $\delta_{ip}$ is represented as .

• •

A slant line  connects two indices in $\Sigma_{ijkpqr}$ belonging to a $\bm{u}$ and an $\bm{x}$ that correspond to different partial derivatives. For example, the index ‘i’ in $\Sigma_{ijkpqr}$ occurs in $\dfrac{\partial{u_{i}}}{\partial{x_{p}}}$ while ‘q’ occurs in $\dfrac{\partial{u_{j}}}{\partial{x_{q}}}$. So, $\delta_{iq}$ is written as .

• •

A horizontal line  connects indices in $\Sigma_{ijkpqr}$ belonging either to two $\bm{u}$’s or two $\bm{x}$’s; the indices involved obviously correspond to different partial derivatives. For example, the index ‘i’ in $\Sigma_{ijkpqr}$ occurs on $\dfrac{\partial{u_{i}}}{\partial{x_{p}}}$ while ‘m’ occurs on $\dfrac{\partial{u_{k}}}{\partial{x_{r}}}$, and $\delta_{ik}$ is thus written as

The above definitions lead to five distinct graphs for $\Sigma_{ijkpqr}$. Thus, the fifteen terms in the original permutation sum may be divided into five groups with elements in a given group having the same graph. Each of the groups is multiplied by a scalar constant, so the effort reduces to determining five rather than fifteen different constants. The different graphs, along with a representative element corresponding to each one, may be stated as follows: $\delta_{ip}\delta_{jq}\delta_{kr}$: , $\delta_{ip}\delta_{jk}\delta_{qr}$: , $\delta_{ip}\delta_{jr}\delta_{qk}$: , $\delta_{iq}\delta_{pk}\delta_{jr}$: , $\delta_{ij}\delta_{pk}\delta_{qr}$:  (note that the ordering of the lines in a graph does not matter. For instance, $\delta_{ip}\delta_{jk}\delta_{qr}$ can be represented by either

or  ).

In light of the above, $\Sigma_{ijkpqr}$ can be written as:

[TABLE]

In actual tensorial notation, this becomes:

[TABLE]

where the $K^{(n)}$’s are defined as:

[TABLE]

$\Sigma_{ijkpqr}$ is invariant to certain indicial permutations, and the grouping on the right hand side of (S26) is consistent with these invariances. Now, from continuity, one has $\Sigma_{ijkiqr}=0$, which leads to:

[TABLE]

Next, we make use of the homogeneity condition Pope_Book :

[TABLE]

which leads to:

[TABLE]

The four relations between the $c_{i}$’s above imply that there is a single scalar that characterizes $\Sigma_{ijkpqr}$. Based on these four relations, one finds $c_{2}=-\dfrac{4}{3}c_{1},\,c_{3}=-\dfrac{1}{6}c_{1},\,c_{4}=-\dfrac{3}{4}c_{1}$ and $c_{5}=c_{1}$, and the constant $c_{1}$ may in turn be expressed in terms of a particular scalar third moment of the turbulent velocity gradient. A convenient (and standard) choice is the skewness based on the longitudinal velocity gradient $\langle(\tau_{\eta}\frac{\partial{u_{1}}}{\partial{x_{1}}})^{3}\rangle$, which is known to be negative (on account of vortex stretching) for homogeneous isotropic turbulence. Thus, the final expression for the third moment of the turbulent velocity gradient reads as:

[TABLE]

The average, $\langle(\tau_{\eta}\frac{\partial{u_{1}}}{\partial{x_{1}}})^{3}\rangle$, which is dependent on $R_{\lambda}$, may now computed from DNS.

The procedure outlined above is now followed for calculating the fourth moment $\Pi_{ijkmsrqp}$. The 105 possible permutations of the Kronecker deltas may be organized into eight groups with all elements in a group again having the same graph. Thus, the calculation reduces to determining only eight constants, and moreover, not all of these are independent (as already seen above, there are constraints imposed from continuity and homogeneity). Thus, $\Pi_{ijkmsrqp}$ may be written as:

[TABLE]

or, in terms of tensorial notation:

[TABLE]

where the $K^{(n)}$’s are defined as:

[TABLE]

As was the case for the third moment, the expression (S33) is consistent with indicial symmetries of $\Pi_{ijkmsrqp}$. Continuity implies $\Pi_{ijkmirqp}=0$, which leads to:

[TABLE]

Next, the homogeneity condition, $\dfrac{\partial{}}{\partial{x_{i}}}\langle u_{m}\dfrac{\partial{u_{i}}}{\partial{x_{j}}}\dfrac{\partial{u_{j}}}{\partial{x_{k}}}\dfrac{\partial{u_{k}}}{\partial{x_{m}}}\rangle=0$, after some manipulation, leads to the following relation:

[TABLE]

Rather unexpectedly, the above relation turns out to be an identity based on the relations already known above from continuity. In other words, the homogeneity constraint does not lead to new relations between the scalar constants $c_{n}$’s. This implies that, of the original eight, three constants are independent, and evaluating them requires three independent (non-dimensional) scalar combinations of four velocity gradients. These are conveniently chosen as:

[TABLE]

Solving (S35) and (S37) leads to:

[TABLE]

The final expression for the fourth moment of the velocity gradient therefore reads as:

[TABLE]

where the $R_{\lambda}$-dependent values of $G_{1}$, $G_{2}$ and $G_{3}$ are again obtained from DNS.

II.4 Fourth-moment of the spheroid orientation distribution

Having derived the form for the fourth moment of the velocity-gradient tensor in homogeneous isotropic turbulence, we proceed to the derivation of $\langle(1-p_{3})^{2}\rangle$(equation (5) in the main text). In the rapid-settling limit,

[TABLE]

which involves the quantities $\langle p_{1}^{4}\rangle$, $\langle p_{2}^{4}\rangle$ and $\langle p_{1}^{2}p_{2}^{2}\rangle$. We raise equations (S15) and (S16) to the fourth power and also multiply the squares of the two and ensemble average the resulting terms; isotropy arguments are then used to rule out averages of the form $\langle\bm{EEES}\rangle$, for instance. This finally leads to the following expressions:

[TABLE]

keeping in mind that $p_{1,2}\ll p_{3}\approx 1$. To calculate the averages occuring in (S41), for instance, $\langle S_{31}S_{31}S_{31}S_{31}\rangle$, we use $\bm{S}=\frac{\bm{\Gamma}-\bm{\Gamma}^{\dagger}}{2}$ and $\bm{E}=\frac{\bm{\Gamma}+\bm{\Gamma}^{\dagger}}{2}$. This particular ensemble average then takes the form:

[TABLE]

Using  (S39) for the fourth moment, the terms occuring on the right-hand side above can be evaluated. Adopting the same procedure to calculate all the ensemble averages in (S41), we get:

[TABLE]

Hence, we obtain the following expression for the fourth moment of the orientation distribution for an oblate spheroid:

[TABLE]

where $M_{1}=\dfrac{3G_{1}}{2}+\dfrac{32G_{2}}{15}-\dfrac{162G_{3}}{15}$, $M_{2}=-3G_{1}+\dfrac{8G_{2}}{15}+\dfrac{54G_{3}}{5}$ and $M_{3}=\dfrac{3G_{1}}{2}$.

Similarly, the fourth moment for a nearly broadside-on prolate spheroid($p_{3}\ll p_{1,2}\approx 1$) is given by:

[TABLE]

where $M_{1}=\dfrac{9G_{1}}{16}+\dfrac{4G_{2}}{5}-\dfrac{81G_{3}}{20}$, $M_{2}=-\dfrac{9G_{1}}{8}+\dfrac{G_{2}}{5}+\dfrac{81G_{3}}{20}$ and $M_{3}=\dfrac{9G_{1}}{16}$.

III Parameter space for the direct numerical simulations(DNS)

Figures 1, 2 and 4 in the main manuscript showcase DNS data and RST results based on the parameters listed in table 1 while figure 3 in the main manuscript is based on table 2. The physical parameters listed in the table above have been drawn from pumir_18 . Except for $R_{\lambda}$, the other parameters are representative of the atmospheric scenario such as ice-crystal size and aspect ratio (klett_book , Veal_1970 ) and the dissipation rates(siewert_14 ). The $R_{\lambda}$’s are lower than in the atmospheric case and this leads to less intermittent distributions.

IV Results

IV.1 Preferential Concentration

A suspension of spherical particles in a turbulent flow is no longer spatially homogeneous when $Fr_{\eta},St_{\eta}\sim\mathcal{O}(1)$ BecRay_2014 ; eaton_Review . Preferential sampling of regions of low vorticity by inertial particles, together with a sweeping effect in presence of gravity, leads to enhanced settling velocities maxey_1987 ; maxey_1993 ; collins_2014 . Fig. 3 in the manuscript shows this to be true for the suspensions of spheroids considered here. In this figure, it is seen that, for large $Fr_{\eta}$, the time-averaged settling speeds from the DNS agree with the orientational averages for $R_{\lambda}=47$ and $96$ (the $Fr_{\eta}$ required for this agreement increases with increasing $St_{\eta}$). For finite $Fr_{\eta}$ and $St_{\eta}$, however, the time averages consistently exceed the orientation averaged estimates.

Fig. S1 above confirms that the discrepancy between the time and orientation-averaged settling speeds in Fig. 3 (in the main manuscript) is due to the preferential sweeping effect. The insets show instantaneous snapshots of particle positions for (a) $R_{\lambda}=47$ and (b) $150$. The particle concentration field remains spatially homogeneous for $R_{\lambda}=47$, in which case $St_{\eta}=0.025$; while there is clear evidence of clustering for $R_{\lambda}=150$. The spatial inhomogeneity in the particle concentration fields has also been characterized via pair-distribution functions (not shown). The probability distributions for the occurrence of upflow ($u_{3}>0$) and downflow ($u_{3}<0$) along particle trajectories have been shown alongside in Fig.S1. The enhanced sampling of downflow regions for $R_{\lambda}=150$ is evidence of preferential sweeping. Preferential sweeping effects in the anisotropic particle suspensions examined here should not come as a surprise since the particle orientation distributions are localized around the broadside-on orientation, and the variation of settling velocity with orientation is therefore minimal, implying a resemblance to spherical particles.