Clustering and energy spectra in two-dimensional dusty gas turbulence

TL;DR

This study uses DNS to analyze how heavy inertial particles (dust) influence two-dimensional turbulent flows, revealing new clustering behaviors and energy spectrum scaling regimes dependent on particle properties and mass-loading.

Contribution

It demonstrates the dependence of dust clustering and energy spectra on Stokes number and mass-loading, introducing a new scaling regime in the kinetic energy spectra.

Findings

Clustering decreases with increased mass-loading.

A new energy spectrum scaling regime emerges.

Scaling exponent varies with Stokes number and mass-loading.

Abstract

We present Direct Numerical Simulation (DNS) of heavy inertial particles (dust) immersed in two-dimensional turbulent flow (gas). The dust are modeled as mono-dispersed heavy particles capable of modifying the flow through two-way coupling. By varying the Stokes number (St) and the mass-loading parameter $({\phi}_{\rm m})$, we study the clustering phenomenon and the gas phase kinetic energy spectra. We find that the dust-dust correlation dimension $(d_2)$ also depends on ${\phi}_{\rm m}$. In particular, clustering decreases as mass-loading $({\phi}_{\rm m})$, is increased. In the kinetic energy spectra of gas we show: (i) emergence of a new scaling regime, (ii) the scaling exponent in this regime is not unique but rather a function of both St and $({\phi}_{\rm m} )$. Using a scale-by-scale enstrophy budget analysis we show in the new scaling regime, viscous dissipation due to the gas balances back-reaction from the dust.

Equations

$$\frac{d\mathbf{X}(t)}{dt}$$

$$\frac{d\mathbf{V}(t)}{dt}$$

$$D_t\omega (\mathbf{x},t)=\nu {\nabla }^2\omega (\mathbf{x},t)-\alpha \omega (\mathbf{x},t)+f(\mathbf{x},t)+\nabla \times {\mathbf{F}}^{\mathrm{d}\to \mathrm{g}}(\mathbf{x},t)\text{/}.$$

$${\mathbf{F}}^{\mathrm{d}\to \mathrm{g}}(\mathbf{x},t)=\sum _{\mathrm{i}=1}^{N_{\mathrm{p}}}\frac{m}{{\tau }_{\mathrm{p}}{\rho }_{\mathrm{g}}}\left[{\mathbf{V}}_{\mathrm{i}}-\mathbf{u}(\mathbf{x},t)\right]{\delta }^2(\mathbf{x}-{\mathbf{X}}_{\mathrm{i}}),$$

$$\mathbf{u}(\mathbf{X},t)=\sum _{\mathbf{x}}\mathbf{u}(\mathbf{x},t){\delta }_{\mathrm{h}}^2(\mathbf{x}-\mathbf{X})h^2\text{/},$$

$$\delta(x-X)=\begin{dcases}\frac{1}{4h}\bigg{\{}1+\cos\left[\frac{\pi(x-X)}{2h}\right]\bigg{\}},\ \ \ \ |x-X|\leq 2h,\\ \ 0\ \ \ \ \ \text{otherwise}.\end{dcases}$$

$$N(r)\equiv ⟨\frac{2}{N_{\mathrm{p}}(N_{\mathrm{p}}-1)}\sum _{i<j}\Theta (r-\mid {\mathbf{X}}_{\mathrm{i}}-{\mathbf{X}}_{\mathrm{j}}\mid )⟩.$$

$$E(k)\equiv \frac{1}{2}⟨\sum _{k-1/2\le m<k+1/2}\mid {\mathbf{u}}_{\mathbf{m}}{\mid }^2⟩\text{/},$$

$$\Pi (k)=\mathcal{D}(k)-\alpha \Omega (k)+\mathcal{F}(k)+\mathcal{R}(k).$$

$$\Omega (k)\equiv ⟨\sum _{m\le k}\mid {\omega }_m{\mid }^2⟩$$

$$\Pi (k)\equiv ⟨\sum _{m\le k}{\omega }_{\mathbf{m}}(\mathbf{u}\cdot \nabla \mathbf{\omega }{)}_{-\mathbf{m}}⟩$$

$$\mathcal{D}(k)\equiv -\nu ⟨\sum _{m\le k}{m}^2\mid {\omega }_{\mathbf{m}}{\mid }^2⟩\text{/},$$

$$\mathcal{F}(k)\equiv ⟨\sum _{m\le k}{\omega }_{\mathbf{m}}{f}_{-\mathbf{m}}⟩$$

$$\mathcal{R}(k)\equiv ⟨\sum _{m\le k}{\omega }_{\mathbf{m}}(\nabla \times {\mathbf{F}}^{\mathrm{d}\to \mathrm{g}}{)}_{-\mathbf{m}}⟩$$

Full text

1 Introduction

In nature, turbulent flows often include small particles embedded within the flow, typical examples are (a) proto-planetary disks (gas and dust) Armitage (2010), (b) clouds (air and water droplets) Pruppacher and Klett (2010) and (c) aeolian processes (wind and sand) Kok et al. (2012). Analytical, numerical and experimental studies of such multiphase flows have flourished in the last decade (see e.g. Refs. Zaichik et al. (2008); Toschi and Bodenschatz (2009); Pumir and Wilkinson (2016); Gustavsson and Mehlig (2016) for a review). For notational convenience, in the rest of this paper, we shall call the solvent phase “gas” and the solute phase “dust”. Often the simplest model used to study such multiphase flows assume that the dust is a collection of heavy, inertial particles (HIPs) which do not alter the gas flow. The equation of motion of the dust particles are

[TABLE]

where ${\bm{X}}$ is the position, ${\bm{V}}$ is the velocity of a dust particle, and ${\bm{u}}$ is the velocity of the gas at a point ${\bm{X}}$. For incompressible flows, in addition to the Reynolds number, an additional dimensionless number appears, the Stokes number $\mbox{St}\equiv\tau_{\rm p}/\tau_{\rm f}$ where $\tau_{\rm f}$ is a characteristic timescale of the flow of gas. If the size of dust grains are comparable to or larger than the dissipative scales of the flow then the simple approximation encoded in Eq. (1) is not valid any more. Furthermore, Eq. (1) is a reasonable model of reality if the number density of the dust grains is so small that both the dust-dust interaction and the back-reaction from the dust phase to the gas phase can be ignored. In this paper, we study the consequences of relaxing this last assumption. One of our motivations is the recent realization that the dust in astrophysical plasma cannot be treated merely as a passive component. In particular, the inclusion of the back-reaction allows for novel instabilities, e.g. the streaming instability (Youdin and Goodman, 2005; Johansen and Youdin, 2007), to manifest itself.

In the absence of dust, the turbulence in the gas phase has been extensively studied Frisch (1996); Pope (2000); Davidson (2004). The pioneering work of Kolmogorov Kolmogorov (1941) has established that in three-dimensions the (angle-integrated) energy spectrum of the gas shows power-law behavior $E(k)\sim k^{-5/3}$ within the inertial range followed by the dissipation range where the energy spectrum shows exponential decay 111Experiments and recent numerical simulations have demonstrated that the Kolmogorov picture is not complete, but must include corrections due to intermittency. The intermittency corrections to the energy spectrum is small and is ignored in this paper.. More importantly, the inertial range spectral exponent is universal, i.e. it does not depend on the Reynolds number and the mechanism of turbulence generation. Does the presence of dust modifies this energy spectrum? Obviously, in general, the answer depends on the number, size, and shape of the dust grains. In this paper, we study this question using direct numerical simulations (DNS) of the dusty gas flow.

A recent paper Gualtieri et al. (2017a) has suggested that in the presence of dust a new power-law behavior can emerge where $E(k)\sim k^{-4}$ in three dimensions. Is this exponent universal, in the sense that, is it independent of the Stokes number and the dust concentration? It is difficult to provide an answer to this question because an accurate determination of the exponent requires obtaining clean scaling of the energy spectrum over at least a decade. This is a formidable task in three dimensions but is a much simpler proposition in two-dimensions. Hence to understand the universality (or lack thereof) we study this problem in two dimensions.

In two-dimensional gas turbulence Boffetta and Ecke (2012), forced at large scales (small $k$) and in presence of air-drag friction ($\alpha$), the energy spectrum is universal with respect to the Reynolds number but does depend on the air-drag-friction coefficient Nam et al. (2000); Boffetta et al. (2005). The scaling exponent and its non-universality can be understood as an effect of the loss of enstrophy due air-drag-friction Verma (2012). In our simulations we choose an $\alpha$ such that in the absence of dust $E(k)\sim k^{-3.9}$. We then perform extensive simulations of the dusty-gas flow by varying both the Stokes number and the mass-loading-parameter ($\phi_{\rm m}$, ratio of the total mass of the dust to the total fluid mass).

The rest of the paper is organized in the following manner. In section (2) we present our model and describe how it is implemented numerically. The result section (3) is divided in three subsections. First subsection (3.1) is devoted to discussion of the pair distribution function of dust particles where we show that increasing mass-loading parameter reduces the clustering of dust. In the remaining subsections (3.2) and (3.3), we present energy spectra and scale-by-scale enstrophy budget respectively for the gas phase. We show that indeed in the presence dust-gas coupling, a new scaling range emerges in the gas kinetic energy spectra. Furthermore, using a scale-by-scale enstrophy budget analysis we show that the new scaling appears due to a balance between the injection (from the dust to the gas) and viscous dissipation. Our main result is that the scaling exponent is not universal but depends on both St and the mass-loading parameter $\phi_{\rm m}$. Finally, in section (4) we conclude the paper.

2 Model and Numerical Method

The dust is modeled as a system of mono-dispersed spherical particles governed by Eq. (1). Gas is modeled in the Eulerian-framework where the equation for the scalar vorticity field $\omega(\bm{x},t)\equiv\nabla\times\bm{u}(\bm{x},t)$ is

[TABLE]

Here ${\bm{u}}({\bm{x}},t)$ is the incompressible velocity field, $D_{t}=\partial_{t}+\bm{u\cdot}\nabla$ is the material derivative, $\nu$ is the viscosity, $\alpha$ is the Ekman drag coefficient, and $f({\bm{x}},t)=-f_{\rm 0}k_{\rm f}\cos(k_{\rm f}y)$ is the Kolmogorov forcing with amplitude $f_{\rm 0}$ and at wave-number $k_{\rm f}$. The force exerted by the dust particles on the gas is

[TABLE]

where $m$ is the mass of a dust particle and $N_{\rm p}$ is the total number of particles. We use a pseudo-spectral method Canuto et al. (1988); Perlekar et al. (2011) to numerically integrate Eq. (2) in a periodic square box with each side of length $L=2\pi$. The simulation domain is spatially discretized using $N^{2}$ collocation points. For time evolution we employ a second-order Runge-Kutta scheme Press et al. (1992).

In this Eulerian-Lagrangian framework, the position of a dust grain does not, in general, coincides with the Eulerian grids. The gas velocity at the position of a dust in Eq. (1b) is obtained as

[TABLE]

where $h=L/N$, and the $\delta^{2}_{h}(\cdot)$ is a numerical realization of the two-dimensional delta function on grids of linear dimension $h$. We use the following prescription Peskin (2002):

[TABLE]

The same prescription, Eq. (5), is also used to discretize the delta function in Eq. (3). We initialize our simulation with $N_{\rm p}$ randomly placed dust particles. Similar to aerosols in clouds Shaw (2003), we assume $\rho_{\rm d}/\rho_{\rm g}\sim 10^{3}$. The vorticity is initialized as $\omega(\bm{x},0)=-f_{\rm 0}k_{\rm f}\nu\left[\cos(k_{\rm f}x)+\cos(k_{\rm f}y)\right]$.

3 Results

We study the dust-gas turbulence by varying the mass-loading parameter $\phi_{\rm m}\equiv{N_{\rm p}m/(\rho_{\rm g}L^{2})}$ and the Stokes number St.

3.1 Vorticity and Clustering

In Fig. (1) we show the pseudo-color plot of the vorticity field in absence of dust ($\phi_{\rm m}=0$) and at high mass-loading $\phi_{\rm m}=1$. We observe that the in the latter small-scale vortices form in the regions where particle cluster. Our observation is consistent with the earlier study of two-dimensional dusty-gas turbulence Bec et al. (2017). We quantify the clustering by using the cumulative pair distribution function

[TABLE]

Here $\Theta$ is the Heaviside function and the angular brackets denote averaging over different stationary-state turbulent configurations. In Fig. (2a), we plot $N(r)$ versus $r$ for fixed $\phi_{\rm m}=1$ and with different St. In the limit $r{\to 0}$, $N(r)\sim r^{d_{\rm 2}}$, where $d_{\rm 2}$ is the correlation dimension Grassberger and Procaccia (1983). We obtain $d_{\rm 2}$ by performing a least square fit in the range $1<r\eta<10$.

In Fig. (2b) we plot the correlation dimension $d_{\rm 2}$ as a function of St for different values of $\phi_{\rm m}$. Note that $d_{\rm 2}=2$ for $\mbox{St}=0,\infty$ and attains a minimum value, which corresponds to maximum clustering, around $\mbox{St}\approx 0.6$ Bec et al. (2007). We observe that for all the values of $\phi_{\rm m}$ this is indeed the case. However, the amount of clustering (smallest value of $d_{\rm 2}$) decreases with increasing $\phi_{\rm m}$. We find that for a fixed St, the maximum clustering is obtained for one-way coupled simulations where the back-reaction from the dust is ignored. Similar results have also been observed for particle-laden turbulent homogeneous shear flows Gualtieri et al. (2011, 2017b). Qualitatively, the small-scale vortices produced in presence of mass-loading will expel particles hence clustering reduces as with $\phi_{\rm m}$ increases.

3.2 Energy spectra

Next, we study the angle averaged velocity power spectrum

[TABLE]

where $\bm{u}_{\bm{m}}$ is the velocity field in the Fourier space. In the absence of dust particles ($\phi_{\rm m}=0$), steady-state two-dimensional energy spectrum [Fig. (1a)] shows inertial range scaling $E(k)\sim k^{-3.9}$ for $0.03\leq k\eta\leq 0.1$ and decays exponentially in the dissipation range ($k\eta>0.10$) Perlekar et al. (2011). We find that addition of dust particles to the gas dramatically alters the dissipation range spectrum. For $k\eta>0.1$, we observe a new power-law $E(k)\sim k^{-\xi}$ with $\xi<3.9$. Let us hold $\phi_{\rm m}=1$ fixed and increase St (Fig. (3a)): we find $\xi$ first decreases, reaches its minimum value $\xi\sim 3$ for $\mbox{St}=0.33$ [see insets to Fig. (3a)] and then increases again. For a fixed $\mbox{St}=0.67$, $\xi$ reduces monotonically as $\phi_{\rm m}$ is increased Fig. (3b).

3.3 Enstrophy budget

To understand the scaling behavior we now study the scale-by-scale enstrophy budget equation:

[TABLE]

Here

[TABLE]

is the cumulative enstrophy up to wave-number $k$,

[TABLE]

is the enstrophy flux due to non-linear terms,

[TABLE]

is cumulative dissipation rate, $-\alpha\Omega(k)$ is the contribution due Ekman friction,

[TABLE]

is the cumulative energy injected due to Kolmogorov forcing, and

[TABLE]

is the contribution because of back reaction from the dust particles to the gas.

In Fig. (4a) we plot the enstrophy budget for the gas in absence of particles ($\phi_{\rm m}=0$). Similar to earlier studies, we observe that at large scales energy injected by external forcing is primarily balanced by Ekman drag and the enstrophy flux $\Pi(k)$ decreases with increasing $k$ Boffetta and Ecke (2012).

We now show that the presence of dust particles dramatically alters the enstrophy budget in the dissipation range. In Fig. (4b) we plot the cumulative contributions of all the terms in budget for $\mbox{St}=0.67$ and $\phi_{\rm m}=1.0$. The dust particles inject enstrophy ($\mathscr{R}$) at large $k$ which is then balanced by viscous dissipation $\mathscr{D}$. We find a negligible change in shape of $\Pi$, $\mathscr{F}$ and the Ekman drag term in the inertial range. A closer look at $\mathscr{R}$ [inset Fig. (4b)] reveals that it makes a net negative contribution to budget till a wavenumber $k_{\rm c}$ after which it turns positive. Clearly, the particles extracts enstrophy from the flow at small $k$ (large scales) but injects enstrophy at large $k$ (small scales). Furthermore, for $k>k_{\rm c}$ the two dominant terms that balance each other are $\mathscr{R}$ and $\mathscr{D}$. Hence, we expect $E(k)\sim k^{-4}d[\mathscr{R}(k)]/dk$ for $k>k_{\rm c}$. In Fig. (5a,b) we show that $\mathscr{R}(k)\sim k^{\beta}$ for $k>k_{\rm c}$. Using the relation between $E(k)$ and $R(k)$, we find the scaling exponent $\xi=-5+\beta$. For $\phi_{\rm m}=1$ and $\mbox{St}=0.33,0.67,1.0$ we obtain $\xi=3.08,3.17~{}\rm{and}~{}3.36$ respectively, which is very close to the exponents that we directly determine from the energy spectra [Fig. (3)]. This proves that the dominant balance between $\mathscr{R}$ and $\mathscr{D}$ determines the new scaling exponent.

4 Conclusion

We use an Eulerian-Lagrangian formalism to study the effects of dust to gas coupling in two-dimensional turbulence. The dust are modeled as heavy inertial particles immersed in the gas. We solve gas equations on fixed Eulerian grids by incorporating the forces [Eq. (3)] due to dust . The main problem with this technique is that to have a smooth Eulerian representation of the feedback, number of particles per cell needs to be equal or greater than a certain threshold ($\approx 1$) Balachandar and Eaton (2010); Gualtieri et al. (2013). We choose $N_{\rm p}$ such that in the stationary state, i.e. after the dust have clustered, the ratio $N_{\rm p}/N^{d_{\rm 2}}\approx 1$ for almost all the St. Furthermore, we use higher order weight function for extrapolation to ensure a better approximation of back reaction in fluid grids. We obtain reasonable scaling range for nearly all the St and $\phi_{\rm m}$.

The Eulerian-Lagrangian formalism has been extensively used to understand how the interaction with dust modifies three-dimensional turbulence. Here, we shall review some of them with an emphasis on energy spectra (see Balachandar and Eaton (2010); Poelma and Ooms (2006) and references therein for more details). Refs. Squires and Eaton (1990); Boivin et al. (1998) studied the effects of dust in isotropic stationary turbulence using direct numerical simulations while similar studies in decaying turbulence were done by Refs. Ferrante and Elghobashi (2003); Sundaram and Collins (1999). The key results of these studies are: (a) particle injects energy at large $k$ and reduces it at small $k$, and (b) increasing mass loading leads to reduction of the total kinetic energy. But, the effect of St or $\phi_{\rm m}$ on the scaling of energy spectra remained unclear as these simulations were done at small or moderate resolution.

More recently, Ref. Gualtieri et al. (2015) introduced a new numerical scheme to model coupling between gas and dust wherein the disturbance in the gas due to dust particle is evaluated in a closed analytic form (by using solutions of the Stokes equation) and is incorporated into fluid grids after a certain regularization time ($\epsilon_{\rm r}$). Here $\epsilon_{r}$ is the time taken by a subgrid-scale disturbance to reach nearby fluid grid locations from an off-grid particle location. One major advantage of this method is that the number of particles need not be comparable to number of grid cells for smooth feedback. Unfortunately, the method is computationally expensive and not easily parallelizable on distributed-memory machines. By studying dust laden homogeneous shear turbulent flow using this technique, Ref. Gualtieri et al. (2017a) reported a scaling exponent of $-4$ in gas kinetic energy spectrum [for $\mbox{St}=1$ and $\phi_{\rm m}=(0.4,0.8)$]. They argued that the new scaling appears due to the balance of viscous forces with the back reaction from dust. In Gualtieri et al. (2017a), the new scaling was observed for $k\eta\geq 1$, whereas our two-dimensional study shows that it starts around $k\eta\sim 0.2$. Clearly, the crucial problem with our and similar studies is that there is, as yet, no well-established algorithm to numerically calculate the feedback in DNS.

For good reason, the most important one being difficulties in experimental realization, turbulence in flows of dust and gas has been rarely studied in two dimensions. Ref. Bec et al. (2017) using Eulerian description of dust found a scaling exponent of $-2$ in the gas energy spectra, that once again emerges due to balance of viscous dissipation against the feedback, for $\mbox{St}\ll 1$ and $\phi_{\rm m}$ between $0.1-0.4$. To numerically smoothen the caustics that invariably develops in such a computation a synthetic hyper-viscous term was added in the Eulerian description. Ekman drag coefficient was chosen such that the pure gas spectra (without dust coupling) scale with an exponent $-3.3$. Notably, the new scaling here starts at much small $k$ compared to what we find.

To summarize our main results: (a) presence of dust-gas coupling decreases clustering of dust particles, (b) a new scaling regime emerges in the kinetic energy spectrum, marked by rise in the tail, (c) scale-by-scale enstrophy budget, suggests that the new scaling is because of gas viscosity dissipating the enstrophy injected by dust at those scales. (d) dust has a net negative contribution to budget till a wavenumber $k_{\rm c}$ and injects enstrophy at higher fourier modes, (e) as the form of dust-gas coupling term varies with both $\phi_{\rm m}$ and more importantly St, the scaling exponent is non-universal and a function of both.

We conclude that once feedback from the dust to the gas is significant the spectra of gas changes, not in the inertial range but a new scaling regime emerges in the erstwhile dissipative range. The central message of this paper is that this scaling exponent is non-universal, it depends on the Stokes number of the dust and the mass loading parameter. Even in two dimensions, where we have been able to do large-scale simulations for long enough time, the appearance of the new scaling regime is not always prominent. We cannot rule out the possibility that there may be no scaling range at all. But we can and do conclude that the spectra is non-universal.

It is quite difficult to perform a DNS of similar resolution, with feedback from particles, in three dimensions. So it is unlikely that in near future we shall observer clear scaling behavior in analogous cases in three dimensions. But based on our result we speculate that the same non-universal nature of spectra will be true in three-dimensions too.

5 Acknowledgment

We thank Paolo Gualtieri for useful discussions. DM acknowledges financial support from the grant Bottlenecks for particle growth in turbulent aerosols from the Knut and Alice Wallenberg Foundation (Dnr. KAW 2014.0048) and from Swedish Research Council Grant no. 638-2013-9243 as well as 2016-05225.