Importance of fluid inertia for the orientation of spheroids settling in turbulent flow

TL;DR

This paper investigates how fluid inertia influences the orientation of settling spheroids in turbulent flows, revealing that inertia dominates at high Reynolds numbers, affecting particle orientation and drag.

Contribution

It combines effects of fluid inertia and turbulent velocity gradients to determine their relative influence on particle orientation in different flow regimes.

Findings

Fluid inertia dominates at high Reynolds numbers, aligning particles broad side first.

Velocity gradient effects dominate at low Reynolds numbers, aligning particles tip first.

The ratio of torques predicts the dominant orientation mechanism in turbulent flows.

Abstract

How non-spherical particles orient as they settle in a flow has important practical implications in a number of scientific and engineering problems. In a quiescent fluid, a slowly settling particle orients so that it settles with its broad side first. This is an effect of the torque due to convective inertia of the fluid set in motion by the settling particle, which maximises the drag experienced by the particle. Turbulent flows tend to randomise the particle orientation. Recently the settling of non-spherical particles in turbulence was analysed neglecting the effect of convective fluid inertia, but taking into account the effect of the turbulent fluid-velocity gradients on the particle orientation. These studies reached the opposite conclusion, namely that a rod settles preferentially with its tip first, wheras a disk settles with its edge first, therefore minimizing the drag on the particle. Here, we consider both effects, the convective inertial torque as well as the torque due to fluctuating velocity gradients, and ask under which circumstances either one or the other dominate. To this end we estimate the ratio of the magnitudes of the two torques. Our estimates suggest that the fluid-inertia torque prevails in high-Reynolds number flows. In this case non-spherical particles are expected to settle with a maximal drag. But when the Reynolds number is small then the torque due to fluid-velocity gradients may dominate, causing the particle to settle with its broad side first.

Tables

Table 1. (a)

	$ℛ$	$w_o$	${\mathrm{Re}}_p$	$S{v}_L$	$S{v}_{\eta }$	$S$
1	0.6	0.21	0.029	0.17	0.31	1.08
2	1.0	0.36	0.047	0.32	0.57	1.13
3	2.5	0.57	0.074	0.50	0.87	1.15
4	6.3	0.88	0.117	0.87	1.55	1.14
5	11	1.13	0.147	1.12	2.01	1.12

Table 2. (a)

	$ℛ$	$w_o$	${\mathrm{Re}}_p$	$S{v}_L$	$S{v}_{\eta }$	$S$
1	0.6	0.21	0.029	0.17	0.31	1.08
2	1.0	0.36	0.047	0.32	0.57	1.13
3	2.5	0.57	0.074	0.50	0.87	1.15
4	6.3	0.88	0.117	0.87	1.55	1.14
5	11	1.13	0.147	1.12	2.01	1.12

Table 3. (b)

	$ℛ$	$w_o$	${\mathrm{Re}}_p$	$S{v}_L$	$S{v}_{\eta }$	$S$
1	0.6	0.20	0.051	0.20	0.35	1.04
2	1.0	0.37	0.084	0.36	0.64	1.15
3	2.2	0.54	0.123	0.56	1.00	1.14
4	4.8	0.81	0.181	0.87	1.54	1.20
5	8.0	1.03	0.222	1.13	2.02	1.20

Table 4. (a)

	$ℛ$	St	$w_o$	${\mathrm{Re}}_p$	$S{v}_L$	$S{v}_{\eta }$	$S$
1	10	0.01	1.08	0.047	1.12	2.01	1.07
2	10	0.04	1.06	0.091	1.09	1.95	1.10
3	11	0.10	1.13	0.147	1.12	2.01	1.12
4	12	0.40	1.28	0.332	1.25	2.23	1.12

Table 5. (a)

	$ℛ$	St	$w_o$	${\mathrm{Re}}_p$	$S{v}_L$	$S{v}_{\eta }$	$S$
1	10	0.01	1.08	0.047	1.12	2.01	1.07
2	10	0.04	1.06	0.091	1.09	1.95	1.10
3	11	0.10	1.13	0.147	1.12	2.01	1.12
4	12	0.40	1.28	0.332	1.25	2.23	1.12

Table 6. (b)

	$ℛ$	St	$w_o$	${\mathrm{Re}}_p$	$S{v}_L$	$S{v}_{\eta }$	$S$
1	7.0	0.01	0.96	0.065	1.06	1.89	1.15
2	7.5	0.04	0.975	0.129	1.07	1.90	1.18
3	8.0	0.10	1.03	0.222	1.13	2.02	1.20
4	8.0	0.40	1.12	0.439	1.16	2.06	1.14

Equations

$$\frac{d\mathbf{v}}{dt}=\mathbf{g}+\frac{6\pi \mu \tilde{a}}{m_p}{\text{M}}_{St}\cdot (\mathbf{u}-\mathbf{v})+\frac{1}{m_p}\frac{9\pi }{8}{\rho }_f{\tilde{a}}^2|\mathbf{u}-\mathbf{v}|{\text{M}}_I\cdot (\mathbf{u}-\mathbf{v}),$$

$${\text{M}}_{St}=X^A\mathbf{n}\mathbf{n}+Y^A(\mathbb{I}-\mathbf{n}\mathbf{n}),$$

$${\text{M}}_I=[3X^A-(X^A{\cos }^2\alpha +Y^A{\sin }^2\alpha )]X^A\mathbf{n}\mathbf{n}+[3Y^A-(X^A{\cos }^2\alpha +Y^A{\sin }^2\alpha )]Y^A(\mathbb{I}-\mathbf{n}\mathbf{n}),$$

$$\frac{d\hat{\mathbf{\omega }}}{dt}=\left(\begin{array}{c}{\hat{T}}_{xx}/{\hat{I}}_{xx}\\ {\hat{T}}_{yy}/{\hat{I}}_{yy}\\ {\hat{T}}_{zz}/{\hat{I}}_{zz}\end{array}\right)+\left(\begin{array}{c}{\hat{\omega }}_y{\hat{\omega }}_z\frac{{\hat{I}}_{yy}-{\hat{I}}_{zz}}{{\hat{I}}_{xx}}\\ {\hat{\omega }}_z{\hat{\omega }}_x\frac{{\hat{I}}_{zz}-{\hat{I}}_{xx}}{{\hat{I}}_{yy}}\\ {\hat{\omega }}_x{\hat{\omega }}_y\frac{{\hat{I}}_{xx}-{\hat{I}}_{yy}}{{\hat{I}}_{zz}}\end{array}\right).$$

$${\hat{\mathbf{T}}}_{St}=\frac{16}{3}\pi \mu a^3\beta \left(\begin{array}{ccc}\frac{1+{\beta }^2}{{\alpha }_0+{\beta }^2{\gamma }_0}& 0& 0\\ 0& \frac{1+{\beta }^2}{{\alpha }_0+{\beta }^2{\gamma }_0}& 0\\ 0& 0& \frac{1}{{\alpha }_0}\end{array}\right)\cdot \left(\begin{array}{c}\frac{1-{\beta }^2}{1+{\beta }^2}{\hat{S}}_{zy}+({\hat{\Omega }}_{zy}-{\hat{\omega }}_x)\\ \frac{{\beta }^2-1}{1+{\beta }^2}{\hat{S}}_{xz}+({\hat{\Omega }}_{xz}-{\hat{\omega }}_y)\\ {\hat{\Omega }}_{yx}-{\hat{\omega }}_z\end{array}\right),$$

$${\hat{\mathbf{T}}}_I={\mathbf{R}}^{\mathbf{rot}}{\mathbf{T}}_I,\mathrm{with}{\mathbf{T}}_I=-{\rho }_f|\mathbf{u}-\mathbf{v}{|}^2{\tilde{a}}^3{F}_{\beta }\left(\frac{{\mathbf{u}}_s}{|{\mathbf{u}}_s|}\cdot \mathbf{n}\right)\left(\frac{{\mathbf{u}}_s}{|{\mathbf{u}}_s|}\times \mathbf{n}\right),$$

$${\hat{T}}_{I,x}/{\hat{T}}_{St,x}\sim {\hat{T}}_{I,y}/{\hat{T}}_{St,y}\sim \frac{5u_s^2}{8\nu s\log \beta },$$

$${\hat{T}}_{I,z}/{\hat{T}}_{St,z}\sim \frac{5u_s^2{\beta }^2}{16\nu s(\log \beta {)}^2},$$

$${\hat{T}}_I/{\hat{T}}_{St}\sim \frac{3}{32}{F}_0\frac{u_s^2}{\nu s},$$

$$|{\mathbf{T}}_{St}|/|{\mathbf{T}}_I|\sim \mathcal{R},\mathrm{with}\mathcal{R}\equiv \frac{|\mathbf{u}-\mathbf{v}{|}^2}{\nu |\mathbf{\Omega }-\mathbf{\omega }|},$$

$$\mathscr{R}\sim\Big{(}\frac{W_{s}}{U_{0}}\Big{)}^{2}Re_{f}^{1/2}\,.$$

$$\mathscr{R}\sim\Big{(}\frac{W_{s}}{U_{0}}\Big{)}^{2}Re_{f}.$$

$$X^A=\frac{8e^3}{3\left(-2e+(1+e^2)\ln \left(\frac{1+e}{1-e}\right)\right)},Y^A=\frac{16e^3}{3\left(2e+(-1+3e^2)\ln \left(\frac{1+e}{1-e}\right)\right)}.$$

$$X^A=\frac{4e^3}{3\left((2e^2-1){\cot }^{-1}\left(\frac{\sqrt{1-e^2}}{e}\right)+e\sqrt{1-e^2}\right)},Y^A=\frac{8e^3}{3\left((2e^2+1){\cot }^{-1}\left(\frac{\sqrt{1-e^2}}{e}\right)-e\sqrt{1-e^2}\right)}.$$

$${\alpha }_0=\frac{{\beta }^2}{{\beta }^2-1}-\frac{\text{arccosh}(\beta )\beta }{({\beta }^2-1{)}^{3/2}},{\gamma }_0=\frac{-2}{{\beta }^2-1}+\frac{2\text{arccosh}(\beta )\beta }{({\beta }^2-1{)}^{3/2}},$$

$$\begin{array}{rl}{F}_{\beta }=& \frac{-\pi e^2(420e+2240e^3+4249e^5-2152e^7)}{315((e^2+1){\tanh }^{-1}e-e{)}^2((1-3e^2){\tanh }^{-1}e-e)}\\ & +\frac{\pi e^2(420+3360e^2+1890e^4-1470e^6){\tanh }^{-1}e}{315((e^2+1){\tanh }^{-1}e-e{)}^2((1-3e^2){\tanh }^{-1}e-e)}\\ & -\frac{\pi e^2(1260e-1995e^3+2730e^5-1995e^7)({\tanh }^{-1}e{)}^2}{315((e^2+1){\tanh }^{-1}e-e{)}^2((1-3e^2){\tanh }^{-1}e-e)}.\end{array}$$

$${\alpha }_0=\frac{{\beta }^2}{{\beta }^2-1}+\frac{\arccos (\beta )\beta }{(1-{\beta }^2{)}^{3/2}},{\gamma }_0=\frac{2}{1-{\beta }^2}-\frac{2\arccos (\beta )\beta }{(1-{\beta }^2{)}^{3/2}},$$

$$\begin{array}{rl}{F}_{\beta }=& \frac{\pi e^3\sqrt{1-e^2}(-420+3500e^2-9989e^4+4757e^6)}{315\sqrt{1-e^2}(-e\sqrt{1-e^2}+(1+2e^2){\sin }^{-1}e)(e\sqrt{1-e^2}+(2e^2-1){\sin }^{-1}e{)}^2}\\ & +\frac{210\pi e^2(2-24e^2+69e^4-67e^6+20e^8){\sin }^{-1}e}{315\sqrt{1-e^2}(-e\sqrt{1-e^2}+(1+2e^2){\sin }^{-1}e)(e\sqrt{1-e^2}+(2e^2-1){\sin }^{-1}e{)}^2}\\ & +\frac{105\pi e^3(12-17e^2+24e^4)({\sin }^{-1}e{)}^2}{315(-e\sqrt{1-e^2}+(1+2e^2){\sin }^{-1}e)(e\sqrt{1-e^2}+(2e^2-1){\sin }^{-1}e{)}^2}.\end{array}$$

Full text

I Introduction

The settling of non-spherical particles in a flow is of importance in several scientific and engineering problems (Bagheri & Bonadonna, 2016). One example is given by deep cumulus clouds. In such clouds, the temperature falls well below $0^{o}C$, inducing the formation of small ice crystals, which play a very significant role in the formation of precipitation (Pruppacher & Klett, 1997; Chen & Lamb, 1994). The orientation of such small ice crystals also determines how electromagnetic radiation is reflected from clouds (Yang et al., 2015). A second example is the sedimentation of organic and anorganic matter in the turbulent ocean (Kiørboe, 2001; Ruiz et al., 2004). The dynamics of motile micro-organisms in turbulence (Durham et al., 2013; Zhan et al., 2014; Gustavsson et al., 2016; Borgnino et al., 2019), usually slightly heavier than water, has important consequences for the population dynamics and the ecology of the ocean system.

More generally, understanding the angular dynamics of non-spherical particles in turbulent flows is a challenging fundamental problem (Voth & Soldati, 2017). When the particles are heavier than the surrounding fluid, the problem is significantly complicated. Even in the simplest case of spherical particles suspended in a turbulent flow, a detailed understanding of their settling properties is only beginning to emerge (Maxey, 1983; Nielsen, 1993; Good et al., 2014; Bec et al., 2014; Gustavsson et al., 2014b). Clearly, understanding the angular dynamics is essential to describe the settling of non-spherical particles.

An exact theoretical description of the problem requires the solution of the full Navier-Stokes equations, imposing no-slip boundary conditions at the surface of the solid (Naso & Prosperetti, 2010; Homann & Bec, 2010; Candelier et al., 2019). An alternative, much more tractable approach consists in using simplified descriptions, based on known solutions of the Navier-Stokes equations. For example, it is frequently assumed that the particle dynamics can be described by the Stokes approximation, neglecting the acceleration of the surrounding fluid set in motion by the particle. Thus one assumes that the hydrodynamic force on the particle is simply Stokes force, and that the torque is Jeffery’s torque (Jeffery, 1922). The resulting set of equations has been used to study the orientation distribution of non-spherical particles in turbulence (Marchioli et al., 2010; Gustavsson et al., 2014a; Siewert et al., 2014a, b; Gustavsson et al., 2017; Jucha et al., 2018; Naso et al., 2018). This model predicts that the orientation of non-spherical particles settling in turbulence is biased, despite the fact that the turbulence is isotropic. It was found that rods tend to settle tip first, and disks are biased towards settling edge first, i.e. that the orientation of the spheroids preferentially minimises drag.

This model, however, neglects the effect of fluid inertia. How convective fluid inertia affects the orientation of particles settling in a fluid has been studied mostly for a fluid at rest. Whereas the corrections to the translational motion, in particular for a sphere, are very well understood for small particle Reynolds numbers (Happel & Brenner, 1983), the effect of the convective fluid inertia upon the angular dynamics in the same regime has been less studied. Analytic expressions were proposed by Cox (1965) for nearly spherical particles, by Khayat & Cox (1989) for slender bodies, and more recently by Dabade et al. (2015) for spheroids of arbitrary aspect ratios.

Much less is known when particles settle through fluid in motion. Lopez & Guazzelli (2017) measured the angular distribution of slender rods settling in a two-dimensional (2D) steady vortex flow. In most cases they found that the distributions are strongly peaked so that the rods tend to settle edge first, as they would in a quiescent fluid. This indicates that the effect of convective fluid inertia dominates the angular dynamics. Lopez & Guazzelli (2017) performed numerical simulations of a model taking into account convective fluid inertia, generally in good agreement with their experimental results. Kramel (2017) analysed the angular dynamics of slender rods settling in turbulence at different Reynolds numbers, and showed that rods settle with their broad side first. The observed distribution of orientation is very narrow, the more so as the settling number, $Sv$ (a dimensionless measure of the settling velocity) is large. Klett (1995) arrived at the same conclusion, although with a different prediction for the orientation distribution. The study of Gustavsson et al. (2019) provides a consistent description of the orientation statistics.

The aim of the present work is to determine the conditions under which it may be justified to neglect fluid inertia for the angular dynamics of spheroids settling in turbulence. For this, we use the same model as Klett (1995) and Lopez & Guazzelli (2017): we assume that the torque on the particle is simply the sum of Jeffery’s torque and the convective fluid-inertia torque for a particle moving slowly and steadily in a fluid at rest (Dabade et al., 2015). Our model also includes the first fluid-inertia correction due to slip in the centre-of-mass dynamics (Brenner, 1961). We estimate in Section III the ratio $\mathscr{R}$ between the magnitudes of the convective fluid-inertia torque, $\mathbf{T}_{I}$, and Jeffery’s torque, $\mathbf{T}_{St}$. The torque induced by fluid inertia therefore dominates when $\mathscr{R}\gg 1$, so that a spheroid settles with its broad side first (rods with their long edge first, discs broadside down). When $\mathscr{R}\ll 1$, by contrast, we expect a spheroidal particle to settle with its narrow edge first.

In the case of a laminar flow, with a small Reynolds number $Re_{f}=U_{0}L/\nu$ ($U_{0}$ and $L$ are respectively the characteristic velocity and lengthscale of the flow, our estimate of $\mathscr{R}$ leads to conclusions consistent with those of Lopez & Guazzelli (2017), who studied the settling of rods in the overdamped limit. In the small flow Reynolds number case, the value of $\mathscr{R}$ can be either large or small, depending on the ratio between the settling velocity and the flow velocity. As a result, the two possible orientation biases (particles settling with their narrow edge first ($\mathscr{R}\ll 1$), or with their broad edge first ($\mathscr{R}\gg 1$) can be observed.

In turbulent flows, we derive a different estimate of $\mathscr{R}$. This estimate implies that the convective torque always dominates when the settling velocity is larger than the typical flow velocity. As the orientation bias obtained when ignoring completely fluid inertia requires the settling velocity to be larger than the flow velocity (Gustavsson et al., 2017; Jucha et al., 2018), we reach the conclusion that spheroids settling with their narrow side first, as described in Siewert et al. (2014a, b); Gustavsson et al. (2017); Jucha et al. (2018), cannot be observed in turbulent flows.

These predictions are in good agreement with results of numerical simulations presented below, in Section IV.

II Equations of motion

In an incompressible fluid with density $\rho_{f}$ and kinematic viscosity $\nu=\mu/\rho_{f}$, we consider solid spheroids of density $\rho_{p}$ and of semi-axes $a$ and $a\beta$. With our convention, oblate (prolate) spheroids correspond to $\beta<1$ ($\beta>1$). The semi-major axis of the particle, $\tilde{a}$, is equal to $a$ ($a\beta$) for oblate (prolate) particles, and the mass of the particle is $m_{p}=\frac{4}{3}\pi a^{3}\beta\rho_{p}$. We denote by $\mathbf{n}$ a unit vector aligned with the symmetry axis of the particle, $\mathbf{n}$ (Fig. 1). The motion of solid particles is entirely characterized by the velocity of their centers of mass, $\mathbf{v}$, and their angular velocity, ${\boldsymbol{\omega}}$. In the following, we restrict ourselves to particles heavier than the fluid, $\rho_{p}\gg\rho_{f}$.

We restrict ourselves to particles very small compared to the size of the smallest eddies $\eta$, the Kolmogorov length: $\tilde{a}<\eta$. To describe the dynamics of the settling particle we use the same model as Klett (1995) and Lopez & Guazzelli (2017). This model can be rigorously justified only for a particle settling in a homogeneous flow at small particle Reynolds number $Re_{p}\equiv W_{s}\tilde{a}/\nu$ (Cox, 1965; Khayat & Cox, 1989; Happel & Brenner, 1983; Dabade et al., 2015), and neglecting the effects of shear and unsteadiness (Candelier et al., 2016). But the model was shown to qualitatively describe the angular dynamics of rods settling in a two-dimensional vortex flow (Lopez & Guazzelli, 2017). In this model, the equation of motion describing the evolution of the particle velocity $\mathbf{v}(t)$ in the laboratory frame reads:

[TABLE]

where ${\bf g}$ is gravity and $\mathbf{u}=\mathbf{u}(\mathbf{x}(t),t)$ represents the undisturbed velocity of the fluid at the point where the particle is located, $\mathbf{x}(t)$, at time $t$. Aside from gravitational settling, the force acting on the particle results from the slip velocity, $\mathbf{u}_{s}=\mathbf{v}-\mathbf{u}$. Together with the equation $d\mathbf{x}/dt=\mathbf{v}$, Eq. (1) determines the position of the particle center of mass. The second term on the right-hand side (r.h.s.) of Eq. (1) expresses the Stokes drag force, appropriate for a spheroidal particle. The expression of this force involves the anisotropic drag tensor

[TABLE]

where ${\bf n}$ is the unit vector defined above (see Fig. 1), and $X^{A}$ and $Y^{A}$ are dimensionless coefficients which only depend on the shape of the spheroid through the parameter $\beta$. Their expressions can be found in the Appendix. The third term on the r.h.s. of Eq.(1) represents the inertial contribution to the force. The dimensionless tensor $\textrm{\bf M}_{I}$ takes the form:

[TABLE]

where $\alpha$ is the angle between the slip velocity ${\bf u}_{s}={\bf v}-{\bf u}$ and ${\bf n}$ ($\alpha=\cos^{-1}({\bf u}_{s}\cdot{\bf n}/|{\bf u}_{s}|)$).

The angular equation of motion reads:

[TABLE]

We expressed here the equations of motion in the particle frame $(\hat{x},\hat{y},\hat{z})$ with axes $\hat{\bf e}_{x}$, $\hat{\bf e}_{y}$, $\hat{\bf e}_{z}$ attached to the particle. Hereafter, all the quantities expressed in this frame carry a hat. The relation between a tensor expressed in the reference frame of the particle, $\hat{\mathbf{K}}$, and in the reference frame of the laboratory, $\mathbf{K}$ is given by: $\hat{\mathbf{K}}=\mathbf{R^{rot}}\,\mathbf{K}\,\mathbf{R^{rot}}^{-1}$, where $\mathbf{R^{rot}}$ is the rotation matrix between the two frames. We note that the equation of evolution of $\mathbf{R^{rot}}$ is simply $dR_{ij}^{rot}/dt=\varepsilon_{ikl}\Omega_{k}R_{lj}^{rot}$, where $\varepsilon_{ijk}$ is the Levi-Civita symbol. In our simulations, rotations were represented in terms of quaternions (Jucha et al., 2018; Naso et al., 2018; Almondo et al., 2018). The tensor $\hat{\mathbf{I}}$ is the moment-of-inertia tensor, which is diagonal in the reference frame of the particle $(\hat{x},\hat{y},\hat{z})$: $\hat{\mathbf{I}}=(4\pi\rho_{p}a^{5}\beta/15)\,{\rm diag}(1+\beta^{2},1+\beta^{2},2)$. The second term on the r.h.s. of Eq. (4) represents the torque induced by the non-Galilean nature of the particle frame (Landau & Lifschitz, 1976). The first term on the r.h.s. of Eq. (4) results from the hydrodynamic torque $\hat{\mathbf{T}}$ acting on the particle.

Following Lopez & Guazzelli (2017), we express $\hat{\mathbf{T}}$ as the sum of Jeffery’s torque (Jeffery, 1922), $\hat{\mathbf{T}}_{St}$, and of the contribution due to fluid inertia (Cox, 1965; Khayat & Cox, 1989; Dabade et al., 2015), $\hat{\mathbf{T}}_{I}$ for a particle moving steadily in a homogeneous flow: $\hat{\mathbf{T}}=\hat{\mathbf{T}}_{St}+\hat{\mathbf{T}}_{I}$. Jeffery’s torque is given by (Jeffery, 1922):

[TABLE]

where $\hat{S}_{ij}$ and $\hat{\Omega}_{i}$ are the components of the strain and vorticity, expressed in the reference frame of the particle. The quantities $\alpha_{0}$ and $\gamma_{0}$ are dimensionless, and depend only on $\beta$, see Appendix. The torque induced by convective fluid inertia reads, derived in (Dabade et al., 2015), reads (to order $Re_{p}$):

[TABLE]

where $F_{\beta}$ is a dimensionless shape factory which depends only on $\beta$. Its expression is given in the Appendix. We have neglected here possible fluid-inertia corrections due to shear. This is justified provided that the Oseen length, ${\ell}_{O}\approx\nu/W_{s}$, is much smaller than the Saffman length, ${\ell}_{S}\approx(\nu/s)^{1/2}$, where $s$ is the order of magnitude of the velocity gradients in the flow (Candelier et al., 2019). We must therefore require that ${\ell}_{S}/{\ell}_{O}\gg 1$, an assumption we reconsider in the following Section.

III Relative importance of Jeffery’s torque and fluid-inertia torque

To estimate the relative importance of the inertial effects on the dynamics of the particles orientation, we calculate the ratio $\mathscr{R}\equiv|\hat{\bf T}_{I}|/|\hat{\bf T}_{St}|$. The particles orientational dynamics is therefore expected to be correctly described by the Jeffery torque when $\mathscr{R}\ll 1$. We begin by looking at the case of very slim, prolate spheroids, with $\beta\gg 1$. In this limit, one obtains:

[TABLE]

where we recall that $\mathbf{u}_{s}$ is the slip velocity, and where $s$ is the inverse of a characteristic time scale representing the angular slip velocity $|\boldsymbol{\Omega}-\boldsymbol{\omega}|$ or the local strain $|\bf{S}|$. To derive the estimate for the ratio between $\hat{T}_{I}$ and $\hat{T}_{St}$ in Eq. (8), we used the fact that $F_{\beta}$ is no larger than $1$ for fibers with $\beta\gg 1$, see Fig.6. The terms responsible for the change in orientation are the $\hat{x}$ and $\hat{y}$ components in Eq. (8). For the present purpose, we therefore consider only these components. As we are interested only in orders of magnitude, we will approximate $5/8$ by $1$ in the following, and the $\log(\beta)$ factor, which in practice is close to $1$ if $5\leq\beta\leq 100$, will be omitted.

In the opposite limit of very thin disks, $\beta\ll 1$, one obtains:

[TABLE]

where $u_{s}$ and $s$ are defined as before, and $F_{0}=\lim_{\beta\rightarrow 0}F_{\beta}=38/9-17216/(945\pi^{2})$ (Dabade et al., 2015), so $3F_{0}/32\sim 1$. The factor $F_{\beta}$ has once again been omitted since its absolute value is always smaller than $2.5$, as shown in Fig. 6. Equation (9) for disks therefore only differs from Eq. (8) for rods by the $1/\log(\beta)$ factor, which in practice does not vary much, as mentioned in the previous paragraph.

To summarise, for very flat disks ($\beta\ll 1$) and for thin rods ($5\leq\beta\leq 100$, not too thin so that the value of $\log\beta$ remains moderate), the ratio $\hat{T}_{I}/\hat{T}_{St}$ is of order $\mathscr{R}$, defined by

[TABLE]

where $\mathbf{v}$ is the particle velocity, ${\boldsymbol{\omega}}$ its angular velocity, $\mathbf{u}$ is the fluid velocity, ${\boldsymbol{\Omega}}$ the vorticity and $\nu$ is the viscosity of the fluid. The differences $|\mathbf{u}-\mathbf{v}|$ and $|{\boldsymbol{\Omega}}-{\boldsymbol{\omega}}|$ are the translational and rotational slip velocities. The convective fluid inertia torque dominates when $\mathscr{R}\gg 1$, in this case a spheroid settles with maximal drag. When $\mathscr{R}\ll 1$, by contrast, we expect a spheroid to settle with minimal drag. To provide an estimate for $\mathscr{R}$, we approximate the slip velocity as $u_{s}\sim W_{s}$, where $W_{s}\approx g\tau_{p}$ is the settling velocity of the particle. Here $\tau_{p}$ is an estimate of its response time, defined as $a^{2}\log\beta(\rho_{p}/\rho_{f})/(3\nu)$ for prolate bodies and $\pi a^{2}\beta(\rho_{p}/\rho_{f})/(8\nu)$ for oblate ones (strictly speaking, the particle linear response time depends on its orientation, as stated in Eq. (1)). To estimate the inverse characteristic time scale $s$, we need to distinguish the laminar and turbulent cases.

III.1 Turbulent flows

Turbulent flows generate large gradients. Standard estimates show that the vorticity root-mean square is of order $|\langle{\boldsymbol{\Omega}}^{2}\rangle^{1/2}|\sim(U_{0}/L)Re_{f}^{1/2}$. Substituting this expression, together with $|\mathbf{u}-\mathbf{v}|\sim W_{s}$ in Eq. (10) leads to:

[TABLE]

Using the Kolmogorov velocity scale (Frisch, 1995), $u_{K}\sim U_{0}Re_{f}^{-1/4}$, we can rewrite this estimate as $\mathscr{R}\sim({W_{s}}/{u_{K}})^{2}$, consistent with the formulation in Gustavsson et al. (2019) where the width of the orientation distribution in the limit $\mathscr{R}\gg 1$ is evaluated. In any case, Eq. (11) shows that in the high $Re_{f}$ regime, the fluid-inertia torque can be neglected (i.e., $\mathscr{R}$ can be small) only if $W_{s}/U_{0}$ is small. In this limit, the orientation distribution is approximately uniform, as shown in (Siewert et al., 2014a, b; Gustavsson et al., 2017; Jucha et al., 2018). This means that the orientation bias analysed in these works (preferential minimal drag) cannot be actually observed at large $Re_{f}$.

The estimate for the velocity gradients, $s\sim 1/\tau_{K}$, leads to the ratio between the Saffman and the Oseen lengths (Candelier et al., 2019) ${\ell}_{S}/{\ell}_{O}\sim\mathscr{R}^{1/2}\sim W_{s}/u_{K}$. This demonstrates that neglecting fluid-inertia effects due to shear, as done here, is consistent when $\mathscr{R}\gg 1$.

III.2 Laminar flows

The analysis developed in the previous subsection applies also to a laminar flow, such as the cellular flow considered by Lopez & Guazzelli (2017). In the case of a laminar flow, at small $Re_{f}$, the velocity gradient can be simply estimated as $s\sim U_{0}/L$, so that:

[TABLE]

We remark that the relative magnitude of the fluid-inertia contributions to the angular dynamics [Appendix B of (Lopez & Guazzelli, 2017)] amounts to $Re_{p}w_{0}/\lambda_{p}$, where $Re_{p}=\beta aW_{s}/\nu$, $w_{0}=W_{s}/U_{0}$ and $\lambda_{p}=\beta a/L$. As a result, our expression  (12) coincides with theirs. Equation (12) suggests that the distribution of orientation predicted in (Siewert et al., 2014a, b; Gustavsson et al., 2017; Jucha et al., 2018; Naso et al., 2018) – spheroids settling preferentially with a minimal drag – could be observable at small $Re_{f}$.

We notice that estimating $\mathscr{R}$ in the limit $Re_{f}\rightarrow 0$ requires some care. To fully analyze the problem, it is convenient to introduce, in addition to $U_{0}$ and $W_{s}$, the velocity $V_{\nu}=\nu/L$, where $L$ is the characteristic length scale of the problem and $\nu$ the viscosity. The two dimensionless parameters $Re_{f}$ and $\mathscr{R}$ can simply be expressed as: $Re_{f}=(U_{0}/V_{\nu})$ and $\mathscr{R}=(W_{s}/U_{0})\times(W_{s}/V_{\nu})$. The limit $Re_{f}\rightarrow 0$ may be reached in two very different ways, either by letting $U_{0}\rightarrow 0$ or $V_{\nu}\rightarrow\infty$. In the former case, $U_{0}\rightarrow 0$ at fixed $V_{\nu}$, the parameter $\mathscr{R}$ is very large, so the fluid inertia dominates. In contrast, letting $V_{\nu}\rightarrow\infty$ at fixed $U_{0}$ leads to $\mathscr{R}\rightarrow 0$, so the Jeffery torque prevails. A careful estimate of the ratios between $W_{s}$, $U_{0}$ and $V_{\nu}$ is therefore necessary to conclude about the importance of fluid inertia in the $Re_{f}\rightarrow 0$ limit.

To summarize, the main result of this section is our systematic analysis of the ratio between the torque induced by fluid convection and the Jeffery torque, and quantified by $\mathscr{R}$, introduced by Eq. (10). With our definition, the torque induced by convective forces dominates when $\mathscr{R}\gg 1$. A more refined relation between $\mathscr{R}$, the Reynolds number $Re_{f}$, and the ratio between the settling velocity, $W_{s}$ and the flow velocity, $U_{0}$, depends on whether the flow is turbulent, Eq. (11), or laminar, Eq. (12).

IV Results of numerical simulations

In this section, we test numerically the results of the analysis of the previous paragraph. To this end, we study the dynamics of spheroids in different flow models that describe low-$Re_{f}$ flows, as well as high-$Re_{f}$ (turbulent) flows with many spatial scales.

IV.1 Two-dimensional vortex flow

We performed numerical simulations of prolate spheroids ($\beta=20$) settling in a simple two-dimensional (2D) cellular flow mimicking that of Lopez & Guazzelli (2017). The (marginal) probability distribution function of $\theta$, the angle between the axis of symmetry of the particle and gravity, is shown in Fig. 2 for different values of $\mathscr{R}$. In this 2D flow, a random orientation of the particle would result in a uniform distribution of $\theta$. The distributions shown in Fig. 2 are symmetric with respect to the value $\pi/2$, as expected due the symmetry of the objects. At low values of $\mathscr{R}$ (Fig. 2a), the orientation distribution exhibits a peak around $\theta\approx 0$ and $\theta\approx\pi$. Here $\theta$ is the angle between $\mathbf{n}$ and the direction of gravity (Fig. 1). So this corresponds to rods settling with their tips first. As the value of $\mathscr{R}$ increases (Fig. 2b, c), however, the peaks migrate towards $\theta\approx\pi/2$. This corresponds to rods settling with their broad sides first.

IV.2 Three-dimensional turbulence

We investigate turbulent flows a cubic box with periodic boundary conditions in all directions. Our pseudo-spectral (de-alised) code, see e.g. (Jucha et al., 2018), is fully de-aliased. We used a grid resolution $N^{3}=128^{3}$, corresponding to a Reynolds number based on the Taylor microscale $R_{\lambda}=40$. An external force acts on low-wavenumber modes to ensure a constant injection rate of energy, $\varepsilon=100\;cm^{2}.s^{-3}$, and reach stationary homogeneous and isotropic turbulence. The kinematic viscosity $\nu=0.15\;cm^{2}.s^{-1}$ has been adjusted to get a satisfying grid resolution $k_{max}\eta\sim 3.4$, where $k_{max}$ is the highest wave number and $\eta=(\nu^{3}/\varepsilon)^{1/4}$ is the Kolmogorov scale.

Figures 3 and 4 display the results of our simulations. Shown is the distribution of $n_{z}$, the projection on the vertical axis of the unit vector aligned with the particle symmetry axis (see Fig. 1). In 3D flows, such a distribution would be uniform if the spheroids were randomly oriented. The distribution of $n_{z}$ is first shown in Fig. 3 for different values of $\mathscr{R}$, for prolate spheroids ($\beta=5$), panel a, and for oblate spheroids ($\beta=0.02$), panel b. These results were obtained for a constant value of the Stokes number $St=\tau_{p}/\tau_{K}$, where $\tau_{K}$ is the Kolmogorov time scale.

As expected, Fig. 3 shows that for $\mathscr{R}\lesssim 1$ the distribution of angles is almost uniform. When the parameter $\mathscr{R}$ becomes larger than $1$, the fluid-inertia torque dominates the angular dynamics. For prolate particles (Fig. 3(a)) the distribution peaks at $n_{z}\approx 0$, where $\mathbf{n}\perp\mathbf{g}$. This corresponds to particles settling with their long side first. For oblate spheroids, on the other hand (Fig. 3(b)), the distributions peak at $n_{z}=\pm 1$. This corresponds to disks settling with their broad sides first. For increasing values of $\mathscr{R}$ (not shown here), the distributions of $n_{z}$ become more and more peaked around $n_{z}=0$ for prolate spheroids, and $n_{z}=\pm 1$ for oblate particles.

Our simulations results fully support the prediction that when the settling parameter $w_{0}=W_{s}/U_{0}$ is of order unity or larger, then $\mathscr{R}>1$ in a turbulent flow with large $Re_{f}$. In this case the fluid-inertia torque dominates the angular dynamics, leading to a distribution of orientation maximising the drag, at odds with the predictions of models neglecting this torque (Siewert et al., 2014a, b; Gustavsson et al., 2017; Jucha et al., 2018; Naso et al., 2018). A theory for the orientation distribution for $\mathscr{R}>1$, when the fluid-inertia torque dominates, has been derived by Gustavsson et al. (2019) (see also (Kramel, 2017)).

Figure 4 shows that at a fixed value of $\mathscr{R}$, the Stokes number has only a limited effect on the orientation statistics, for the range of parameters considered here ($St\lesssim 1$), consistent with the findings of Gustavsson et al. (2019) (our numerical data exhibit the same behavior for any value of $\mathscr{R}$).

Tables 1 and 2 list the parameters used in the simulations of Fig. 3 and 4, respectively. We defined the settling numbers $Sv_{L}=g\tau_{p}/U_{0}$, $Sv_{\eta}=g\tau_{p}/u_{K}$, and velocity ratio $w_{0}=W_{s}/U_{0}$. The particle Reynolds number, $Re_{p}$, is defined as $W_{s}\tilde{a}/\nu$. The parameter $S=|{\bf u}-{\bf v}|/W_{s}$ is the ratio between the slip velocity and the settling one. Its value is always close to one, thereby confirming the fact that $|{\bf u}-{\bf v}|$ is, in the present setup, always well approximated by $W_{s}$.

The previous results were obtained for a moderate value of $R_{\lambda}=40$. As shown in Fig. 5, similar results are obtained for $R_{\lambda}=80$. Finally, we stress the fact that, for these very dense particles, similar results can also be obtained by using much simplified models of turbulent flows. In fact, we compared directly the prediction of the Kinematic Simulation (KS) model of turbulence with the results of DNS. Briefly, the KS flow (Fung et al., 1992) is obtained by representing the large range of scales in the flow by a superposition of a few Fourier modes, with a distribution of wavenumber amplitude in $k$-space. Details about our implementation can be found in (Ducasse & Pumir, 2010). As shown in Fig. 5, we found a very good agreement between the DNS and KS results, as it was the case in (Gustavsson et al., 2017). This can be understood by noticing that the velocity gradient tensor acting on a settling particle decorrelates faster than on a tracer particle, and as a result, the particle response is not sensitive to the fine correlations of the velocity gradient.

V Discussion

The analysis presented here, along with the numerical results obtained using very simplified flows, shows that the effect of convective fluid inertia on the orientation of particles settling in a flow becomes important when the parameter $\mathscr{R}$, defined by Eq. (10) is $\gtrsim 1$. This condition implies $(W_{s}/U_{0})^{2}\,Re_{f}^{\alpha}\gtrsim 1$, where $\alpha=1$ in the laminar regime ($Re_{f}<1$), and $\alpha=1/2$ at large Reynolds numbers, when the flow is turbulent.

For a turbulent flow this condition implies that the effect of convective fluid inertia can only be neglected when $W_{s}/U_{0}\ll 1$. In this regime, the orientation distribution is essentially uniform. When the ratio $W_{s}/U_{0}$ becomes of order unity, our study suggests that a particle tends to settle preferentially broad side down (with a maximal drag orientation). This conclusion has been previously reached by considering the restoring torque induced by convective inertia (Klett, 1995), and has also been observed in laboratory experiments (Kramel, 2017). The expression for $\mathscr{R}$ given by Eq. (10) was derived in the limits of infinitely thin ($\beta\gg 1$) or infinitely flat ($\beta\ll 1$) particles, but our numerical results show that this parameter is also relevant for moderate values of $\beta$.

The present analysis was performed by keeping, in addition to the Jeffery torque, only the leading-$Re_{p}$ contribution of the convective fluid-inertia torque for a particle settling slowly and steadily in a fluid at rest. This model has been validated, to some extent, in a recent experimental work (Lopez & Guazzelli, 2017). This opens the possibility to understand quantitatively the orientation distribution of small non-spherical particles settling in turbulent flows (Gustavsson et al., 2019), a question of fundamental importance in cloud physics (Pruppacher & Klett, 1997).

Acknowledgements.

We acknowledge discussions with E. Guazzelli and J. Magnaudet. This work was supported by the IDEXLYON project (Contract ANR-16-IDEX-0005) under University of Lyon auspices. We also aknowledge support from Vetenskapsrådet [grant number 2013-3992], Formas [grant number 2014-585], and from the Knut and Alice Wallenberg Foundation (grant ‘Bottlenecks for particle growth in turbulent aerosols’ 2014.0048). The simulations were performed using the resources provided by PSMN (ENS Lyon).

Appendix A: Expression of the coefficients $X^{A}$, $Y^{A}$, $\alpha_{0}$, $\gamma_{0}$ and $F_{\beta}$

In the following, we introduce the ellipticity $e$, defined for oblate bodies ($\beta<1$) by $e=\sqrt{1-\beta^{2}}$ and for prolate ones ($\beta>1$) by $e=\sqrt{1-(1/\beta)^{2}}$.

Expression of the force We begin by giving the expressions of the coefficients $X^{A}$ and $Y^{A}$ that define the resistance tensor $\textrm{\bf M}_{St}$, determined in Kim & Karrila (1991) and introduced in Eq. (2),(3). For prolate spheroids ($\beta>1$):

[TABLE]

For oblate spheroids ($\beta<1$):

[TABLE]

Expression of the torque The tensor necessary to calculate the Jeffery torque, see Eq. (5), requires the two parameters $\alpha_{0}$ and $\gamma_{0}$ (Jeffery, 1922). Additionally, the expression for the torque due to inertial effects involves a factor $F_{\beta}$ (see Eq. (6)), derived by Dabade et al. (2015) (Eq. (4.1) and (4.2) therein). The expressions of these parameters are given below, and the $\beta$-dependence of $F_{\beta}$ is more specifically shown in Fig. 6.

For prolate spheroids ($\beta>1$):

[TABLE]

[TABLE]

For oblate spheroids ($\beta<1$):

[TABLE]

[TABLE]