1 Introduction
The collision rate of heavy particles suspended in a turbulent flow is of interest in many research areas, such as rain formation in clouds (Shaw et al.
Reference Shaw, Reade, Collins and Verlinde1998) or dust grain dynamics in astrophysical environments. One contribution to the collision rate is the radial distribution function (Sundaram & Collins Reference Sundaram and Collins1997; Wang, Wexler & Zhou Reference Wang, Wexler and Zhou2000), which is a measure of the non-homogeneity or local concentration of the particle field. Maxey (Reference Maxey1987) showed that heavy particles tend to cluster in regions of low vorticity and high strain as a result of their inertia, a mechanism later called preferential concentration. Preferential concentration is most pronounced when the particles have a Stokes number
$St$
around unity, which represents the non-dimensional response of the particle to the flow field. The Stokes number is defined as
$St={\it\tau}_{p}/{\it\tau}_{{\it\eta}}$
, where
${\it\tau}_{{\it\eta}}=({\it\nu}/\langle {\it\epsilon}\rangle )^{1/2}$
is the Kolmogorov time scale of the flow and
${\it\tau}_{p}=2{\it\rho}_{p}r^{2}/(9{\it\rho}_{f}{\it\nu})$
the particle relaxation time;
${\it\nu}$
represents the viscosity of the carrier fluid,
$r$
the particle radius and
$\langle {\it\epsilon}\rangle$
is the mean dissipation of the flow;
${\it\rho}_{p}$
is the density of the particles and
${\it\rho}_{f}$
of the carrier fluid.
Squires & Eaton (Reference Squires and Eaton1991) used the flow classification of Hunt, Wray & Moin (Reference Hunt, Wray and Moin1988) to divide the flow into eddies, streams and convergence zones, and showed that particles tend to move through the flow via streaming zones and tend to avoid eddies. Streaming zones typically end in convergence zones, where an increase in particle concentration can be found (Squires & Eaton Reference Squires and Eaton1991; Eaton & Fessler Reference Eaton and Fessler1994).
Chong, Perry & Cantwell (Reference Chong, Perry and Cantwell1990) introduced a different classification of the local flow based on the
$PQR$
invariants, i.e. the first, second and third invariants of the velocity gradient tensor of the flow. In incompressible flows, mass conservation implies that
$P$
equals zero.
$Q$
is a measure of the relative intensity between strain and enstrophy and
$R$
is helpful in differentiating between regions of stretching and compression. Several studies have investigated the role of local flow topology in particle behaviour (Rouson & Eaton Reference Rouson and Eaton2001; Bijlard et al.
Reference Bijlard, Oliemans, Portela and Ooms2010) in wall bounded flows and show that particles with Stokes numbers of the order of unity preferentially concentrate into streamwise flowing low-speed streaks.
By investigating the average flow pattern in local coordinate systems, Elsinga & Marusic (Reference Elsinga and Marusic2010) identify the average flow patterns responsible for some apparently universal aspects of small-scale turbulent motion, such as a preferential alignment of the vorticity vector with the eigenvector corresponding to the intermediate eigenvalue of the strain rate tensor as well as the characteristic teardrop shape of the joint probability density function of
$Q$
and
$R$
. Those aspects have been observed in several turbulent flows, such as turbulent boundary layers, turbulent channel flows and homogeneous and isotropic turbulence.
Collisions in turbulent flows occur in regions of the flow where the dissipation rate is higher than flow averaged. Perrin & Jonker (Reference Perrin and Jonker2014) constructed the following physical picture of a collision in cloud-like turbulence using direct numerical simulation (DNS): enstrophy makes particles preferentially concentrate in quiescent flow regions, thereby increasing the particle velocity coherence, i.e. decreasing relative velocities between particles, needed to approach each other. Strongly clustered particles thus have a low collision probability, until a dissipative event accelerates the particles towards each other. The novelty in their paper was to identify the need for dissipation for particles to collide. Dissipation, however, is a scalar quantity and does not provide information about the shape of the local flow field. In this paper we will therefore investigate the effect of the shape of the local flow structure on the collision chance.
The local flow topology can be categorized based on the eigenvalues
${\it\lambda}_{1}$
,
${\it\lambda}_{2}$
and
${\it\lambda}_{3}$
of the local velocity gradient tensor of the velocity field
$\boldsymbol{u}$
, i.e.
$\mathscr{J}=\partial \boldsymbol{u}/\partial \boldsymbol{x}$
(see figure 1);
$\boldsymbol{x}$
represents the position vector. Physically, complex eigenvalues can be associated with vortices (Zhou, Adrian & Balachandar Reference Zhou, Adrian and Balachandar1996; Zhou et al.
Reference Zhou, Adrian, Balachandar and Kendall1999), whereas real eigenvalues correspond to convergence/divergence zones. Here we choose the following classification for incompressible flows. When the eigenvalues of the velocity gradient tensor are complex (two out of the three eigenvalues form a conjugate pair) rotation is present in the local flow and the following holds:
$$\begin{eqnarray}\displaystyle & \text{Re}({\it\lambda}_{1})=\text{Re}({\it\lambda}_{2}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \text{Im}({\it\lambda}_{1})=-\text{Im}({\it\lambda}_{2}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\it\lambda}_{3}=-2\text{Re}({\it\lambda}_{1}). & \displaystyle\end{eqnarray}$$
Following the nomenclature of Chong et al. (Reference Chong, Perry and Cantwell1990), complex eigenvalues represent a focus. The real part of the eigenvalues
${\it\lambda}_{1}$
and
${\it\lambda}_{2}$
in such a case determines whether the focus is inward spiralling (a) or outward spiralling (b). By looking only at the eigenvectors
$e_{1}$
and
$e_{2}$
, we could also classify the focus as stable (a) or unstable (b).
Graphical representation of the four atoms of flow. (a) Stable focus, (b) unstable focus, (c) stable saddle node, (d) unstable saddle node.
In the case of only real eigenvalues of the velocity gradient tensor, rotation is absent in the local flow structure and the eigenvalues can always be rearranged such that
${\it\lambda}_{1}<{\it\lambda}_{2}<{\it\lambda}_{3}$
. The local flow structure forms a saddle node and the sign of the eigenvalue
${\it\lambda}_{2}$
determines if the node is stable and converging along the eigenvectors
$e_{1}$
and
$e_{2}$
(c) or unstable and diverging along the eigenvectors
$e_{2}$
and
$e_{3}$
(d). Blackburn, Mansour & Cantwell (Reference Blackburn, Mansour and Cantwell1996) show a graphical representation in the
$QR$
plane of the four categories (a–d) as presented in figure 1.
Defining the local flow categories is based on a point-wise evaluation of the velocity gradient tensor; we interpret these four categories as fundamental building blocks of the flow, and refer to them as atoms of flow since they are considered elementary and they are the smallest constituent unit of the flow.
This categorization is closely related to the commonly used categorization of Chong et al. (Reference Chong, Perry and Cantwell1990) which uses the
$\mathit{PQR}$
invariants. It is a matter of taste; the use of eigenvalues can be somewhat more intuitive since an eigenvalue directly represents the local flow direction along its eigenvector. The eigenvalues
${\it\lambda}_{1}$
,
${\it\lambda}_{2}$
and
${\it\lambda}_{3}$
for incompressible flows can directly be translated into the
$\mathit{PQR}$
invariants with the following identities:
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}P={\it\lambda}_{1}+{\it\lambda}_{2}+{\it\lambda}_{3}=0,\\ Q={\it\lambda}_{1}{\it\lambda}_{2}+{\it\lambda}_{1}{\it\lambda}_{3}+{\it\lambda}_{2}{\it\lambda}_{3},\\ R={\it\lambda}_{1}{\it\lambda}_{2}{\it\lambda}_{3}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$
An outline of the paper is as follows. Using DNS of homogeneous and isotropic turbulence, we will first categorize the flow into the four atoms of flow and investigate the effect of the local flow topology on the presence of particles and the occurrence of collision (§ 2). We will use this information to find in the DNS the collision kernel associated with each of these atoms of flow. In § 3 we develop a simplified framework based on the atoms of flow to gain more conceptual insight into the effect of an individual flow structure on the occurrence of collisions.
2 DNS framework
The dynamics of the turbulent flow field are governed by the Navier–Stokes equations. If the density of the particles
${\it\rho}_{p}$
is high compared to the density of the carrier fluid
${\it\rho}_{f}$
and if the particle radius
$r$
is small compared to the Kolmogorov scale
${\it\eta}=({\it\nu}^{3}/\langle {\it\epsilon}\rangle )^{1/4}$
of the flow, the full equations of motion of particles in turbulence (Gatignol Reference Gatignol1983; Maxey & Riley Reference Maxey and Riley1983) can be reduced to
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\boldsymbol{v}(t)}{\text{d}t}=\frac{\boldsymbol{u}[\boldsymbol{x}(t),t]-\boldsymbol{v}(t)}{{\it\tau}_{p}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\boldsymbol{x}(t)}{\text{d}t}=\boldsymbol{v}(t), & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{v}(t)$
is the particle velocity vector,
$\boldsymbol{x}(t)$
the particle position vector and
$\boldsymbol{u}[\boldsymbol{x}(t),t]$
the flow velocity field at the particle position. Gravity is omitted in this paper since it adds complexity to the problem in a delicate way. The combined effect of turbulence and gravity is not merely an addition of separate phenomena (Woittiez, Jonker & Portela Reference Woittiez, Jonker and Portela2009; Bec, Homann & Ray Reference Bec, Homann and Ray2014; Gustavsson, Vajedi & Mehlig Reference Gustavsson, Vajedi and Mehlig2014).
2.1 Numerical details of the DNS
In this paper, an in-house-developed DNS code has been used, which solves the Navier–Stokes equations using pseudo-spectral methods. We use a triple-periodic computational domain, in which time stepping is restricted by the Courant–Friedrich–Lewy criterion using a Courant number
$C$
of 0.1. Since a turbulent system is inherently dissipative, energy is injected at the lowest wavenumber. To this end, we employ a forcing scheme similar to that used by Woittiez et al. (Reference Woittiez, Jonker and Portela2009) to add kinetic energy to the largest scales. This energy has been set to
$(0.25{\it\epsilon}_{t}L)^{2/3}$
;
${\it\epsilon}_{t}$
denotes the target mean dissipation rate of the simulation, which is the mean dissipation rate we aim for (it is not necessarily exactly equal to the mean dissipation rate of the actual simulation), and
$L$
denotes the physical size of the computational domain. Time integration is performed using a third-order Adams–Bashforth scheme. Both advection and diffusion are treated explicitly, and the
$3/2$
-rule is used to fully deal with aliasing errors (see e.g. Canuto, Hussaini & Quarteroni Reference Canuto, Hussaini and Quarteroni2007).
The equations of motion (2.1) and (2.2) are updated using a second-order Runge–Kutta scheme. The velocity of the flow field at the particle position is computed using trilinear interpolation.
A collision routine checks for collisions using the algorithm of Chen, Kontomaris & McLaughlin (Reference Chen, Kontomaris and McLaughlin1998), which uses cell indexing and linked lists to check only particle pairs that could collide within one time step. The cost of this algorithm is
$O((27N_{p}^{2}/2N_{x}N_{y}N_{z}))$
, where
$N_{x}$
,
$N_{y}$
and
$N_{z}$
define the size of the computation domain in the
$x$
,
$y$
, and
$z$
direction, respectively and
$N_{p}$
is the number of particles present in the computational domain. To ensure that all collisions are detected, the maximum travel distance of the particles is restricted to half a grid distance by using a dynamically adaptive time step for the particles. Collisions are only detected but not enacted (i.e. ghost particle approach). Since particles do not interact with each other, the statistics are independent of the number of particles used. See Perrin & Jonker (Reference Perrin and Jonker2014) for more information on the numerical aspects.
2.2 DNS results
DNS allows us to calculate the eigenvalues of the local velocity gradient tensor at each grid point and classify the local flow structure into one of the four categories. The local topology of the flow categorized in this way can then be conditionally sampled on the positions of the particles and the occurrence of collisions. The results shown in figure 2 are obtained from simulation O1, see table 1 for the numerical details. The solid line represents the volume fraction
$V_{n}$
defined by
$$\begin{eqnarray}V_{n}=\frac{1}{N_{x}N_{y}N_{z}}\mathop{\sum }_{\boldsymbol{x}}I_{n}\{\mathscr{J}(\boldsymbol{x})\},\end{eqnarray}$$
where
$I_{n}$
is an indicator operator that returns 1 if
$\mathscr{J}(\boldsymbol{x})$
is categorized as flow configuration (atom)
$n$
, and zero otherwise. As such,
$V_{n}$
represents the fraction of grid boxes that are associated with atom
$n\in \{1,2,3,4\}$
. In the same vein we can base a fraction on particle locations
$\boldsymbol{x}_{p}$
$$\begin{eqnarray}P_{n}(St)=\frac{1}{N_{p}}\mathop{\sum }_{\boldsymbol{x}_{p}}I_{n}\{\mathscr{J}(\boldsymbol{x}_{p})\}\end{eqnarray}$$
and collision locations
$\boldsymbol{x}_{c}$
$$\begin{eqnarray}C_{n}(St)=\frac{1}{N_{c}}\mathop{\sum }_{\boldsymbol{x}_{c}}I_{n}\{\mathscr{J}(\boldsymbol{x}_{c})\},\end{eqnarray}$$
where
$N_{p}$
and
$N_{c}$
are the number of particles and collisions, respectively.
Occurrence of particles and collisions for different values of the Stokes number inside each of the four atoms. The black line represents the relative volume
$V_{n}$
of each atom. The red line represents the relative occurrence
$P_{n}(St)$
of the particles and the blue line represents the relative occurrence
$C_{n}(St)$
of collisions. The green line shows the normalized collision kernel
${\it\Gamma}_{n}(St)$
from (2.6). The results are obtained from DNS simulation O1, see table 1.
Overview of the simulations. The simulation denoted O1 investigates the effect of the eigenvalues. Each simulation is shown together with the dimensions of the domain
$L$
, the number of grid points
$N_{x}$
, the Taylor-based Reynolds number
$Re_{{\it\lambda}}$
, the mean dissipation rate of the flow
$\langle {\it\epsilon}\rangle$
, the Stokes number
$St$
, the number of particles
$N_{p}$
, and
$N_{cat}$
, the number of different particle categories. In all simulations only monodispersed collisions are considered.
Figure 2 shows that the most often occurring flow structures are the stable focus and the unstable saddle node, which is in agreement with previous results (Blackburn et al. Reference Blackburn, Mansour and Cantwell1996). The figure shows a clear preference (relative to the flow) for particles to preferentially concentrate in regions with no complex eigenvalues (i.e. low vorticity, enstrophy), which is in good qualitative agreement with Squires & Eaton (Reference Squires and Eaton1991), Eaton & Fessler (Reference Eaton and Fessler1994) and Rouson & Eaton (Reference Rouson and Eaton2001). Collisions occur in regions with even less rotation and higher levels of convergence. The DNS clearly shows that converging regions tend to favour both the presence of particles and the occurrence of collisions. Since converging regions are often associated with dissipation, our results are in qualitative agreement with the results of Perrin & Jonker (Reference Perrin and Jonker2014), which showed that collisions tend to occur in regions of higher dissipation rates than average.
Based on the previous results, we can compute a normalized collision kernel
${\it\Gamma}_{n}(St)$
as function of the Stokes number for each category:
$$\begin{eqnarray}{\it\Gamma}_{n}(St)=\frac{C_{n}(St)}{P_{n}^{2}(St)V_{n}}.\end{eqnarray}$$
This normalized collision kernel is derived from the definition of the collision kernel and is a measure for the number of collisions given a certain volume and number of particles. Figure 2 (right-hand axis) shows the normalized collision kernel (green line) for the different flow categories as a function of the Stokes number. The stable saddle node, which is the least often occurring flow configuration, shows a huge collision kernel, indicating that in a relative sense a large number of collisions occur. By contrast, the most often occurring atom, the stable focus, has the smallest collision kernel. The unstable focus and the unstable saddle node show a similar collision kernel. In the next section we will try to understand better the effect of individual flow structures on the collision kernel. To this end we will develop a simplified analytical model to compute the collision kernel for a large range of different flow configurations.
Graphical impression of a simulation done with the atoms model for two different flow configurations. By releasing a large number of particle pairs, the collision probability of a flow structure can be computed.
3 Conceptual framework
To understand better the effect of the local flow structures or atoms of flow on the collision kernel, we construct an idealized framework, from now on referred to as the atoms model. By releasing a huge number of particle pairs near an idealized atom of flow (see figure 3), we can compute the collision kernel associated with that flow configuration. This model differs from turbulence since the atoms of flow are assumed to be static. Since we are investigating individual atoms, this model is not suited for particles with
$St\gg 1$
, where the contribution of particle collisions from caustics is dominant. Caustics are singularities in the particle flow field, which result from the detachment between the flow and the particles (see e.g. Falkovich, Fouxon & Stepanov Reference Falkovich, Fouxon and Stepanov2002; Wilkinson & Mehlig Reference Wilkinson and Mehlig2005). The higher the Stokes number, the more the particle motion becomes detached from the local flow field, and the smaller the correlation between flow structures and collision rate (Voßkuhle et al.
Reference Voßkuhle, Pumir, Lévêque and Wilkinson2014).
First the theoretical framework of the atoms model is explained in § 3.1 and how to investigate the collision kernel of a flow structure, followed by the numerical details, implementation and validation in § 3.2. Results obtained with the atoms model will be shown in § 3.3. We will investigate the effect of the velocity gradient tensor on collisions, followed by the effect of dissipation and enstrophy on collisions in § 3.4.
3.1 Motion of heavy particles and collisions statistics
Making the equation of motion for heavy particles (2.1) and (2.2) non-dimensional using the Kolmogorov scales and linearizing the instantaneous fluid velocity
$\boldsymbol{u}$
around a point
$\boldsymbol{x}$
(e.g. Rouson & Eaton Reference Rouson and Eaton2001 and Reeks, Fabbro & Soldati Reference Reeks, Fabbro and Soldati2006) yields the non-dimensional equation of motion for heavy particles as influenced by non-dimensional velocity gradient tensor
$\unicode[STIX]{x1D645}$
:
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\boldsymbol{X}(T)}{\text{d}T}=\boldsymbol{V}(T), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\boldsymbol{V}(T)}{\text{d}T}=-\frac{1}{St}[\boldsymbol{V}(T)-\unicode[STIX]{x1D645}\boldsymbol{X}(T)], & \displaystyle\end{eqnarray}$$
where
$T$
,
$\boldsymbol{V}(T)$
,
$\boldsymbol{X}(T)$
are the rescaled time, the rescaled particle velocity vector and the rescaled particle position vector, respectively. The Jacobian
$\unicode[STIX]{x1D645}$
of the flow is the rescaled velocity gradient tensor of the flow,
$\unicode[STIX]{x1D645}=\partial \boldsymbol{U}/\partial \boldsymbol{X}$
, with
$\boldsymbol{U}$
the scaled flow velocity vector. The flow field in the vicinity of such an atom is now determined by the characteristics of
$\unicode[STIX]{x1D645}$
.
In general
$\unicode[STIX]{x1D645}$
possesses nine degrees of freedom. For incompressible flows, the continuity equation,
$\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{U}=0$
, removes one degree of freedom. By neglecting gravity, we can assume that the problem is rotational invariant, reducing the degrees of freedom by three. Finally, two additional degrees of freedom can be discarded by only considering orthogonal eigenvectors. This assumption is based on the preferential alignment of the vorticity vector and the second eigenvector of
$\unicode[STIX]{x1D645}$
(Ashurst et al.
Reference Ashurst, Kerstein, Kerr and Gibson1987; Blackburn et al.
Reference Blackburn, Mansour and Cantwell1996; Elsinga & Marusic Reference Elsinga and Marusic2010). With above considerations, we can confine ourselves to the following class of Jacobians:
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D645}=\left[\begin{array}{@{}ccc@{}}a & b & 0\\ -b & d & 0\\ 0 & 0 & -(a+d)\end{array}\right]. & & \displaystyle\end{eqnarray}$$
The Jacobian of the local flow structure can be related to the different atom classifications (figure 1) defined by its eigenvalues
${\it\lambda}_{1}$
,
${\it\lambda}_{2}$
and
${\it\lambda}_{3}$
, which are related to the parameters
$a$
,
$b$
and
$d$
by
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}{\it\lambda}_{1}=\displaystyle \frac{a+d}{2}+\sqrt{\left(\displaystyle \frac{a-d}{2}\right)^{2}-b^{2}},\\ {\it\lambda}_{2}=\displaystyle \frac{a+d}{2}-\sqrt{\left(\displaystyle \frac{a-d}{2}\right)^{2}-b^{2}},\\ {\it\lambda}_{3}=-(a+d),\end{array}\right\} & & \displaystyle\end{eqnarray}$$
where
${\it\lambda}_{1}$
and
${\it\lambda}_{2}$
may be complex, in which case they are each other’s complex conjugate. The Jacobian elements
$a$
and
$d$
are related to the convergence and divergence of the flow field, whereas
$b$
is related to the amount of rotation. We now have a framework to track particles inside such an atom defined by
$a$
,
$b$
and
$d$
, using (3.1) and (3.2).
The flow structure of such an atom can be related to the dimensionless dissipation rate and enstrophy by the following relation:
$$\begin{eqnarray}\displaystyle & {\it\epsilon}=2(\unicode[STIX]{x1D645}^{S}\boldsymbol{ : }\unicode[STIX]{x1D645}^{S}),\quad \text{where }\unicode[STIX]{x1D645}^{S}={\textstyle \frac{1}{2}}(\unicode[STIX]{x1D645}+\unicode[STIX]{x1D645}^{\text{T}}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\it\Omega}=2(\unicode[STIX]{x1D645}^{A}\boldsymbol{ : }\unicode[STIX]{x1D645}^{A}),\quad \text{where }\unicode[STIX]{x1D645}^{A}={\textstyle \frac{1}{2}}(\unicode[STIX]{x1D645}-\unicode[STIX]{x1D645}^{\text{T}}), & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D645}^{S}$
and
$\unicode[STIX]{x1D645}^{A}$
are the symmetric and anti-symmetric part of the velocity gradient tensor, respectively;
$\boldsymbol{ : }$
represents the Frobenius inner product, which is the component-wise inner product of two tensors. By substituting
$\unicode[STIX]{x1D645}$
from (3.3), the above expressions can be reduced to
$$\begin{eqnarray}\displaystyle & {\it\epsilon}_{\mathit{atom}}=4(a^{2}+d^{2}+ad), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\it\Omega}_{\mathit{atom}}=4b^{2}. & \displaystyle\end{eqnarray}$$
It can be sometimes more convenient to enforce a certain value for the dissipation rate and enstrophy, and find the corresponding values for
$a$
,
$b$
and
$d$
. The conversion between enstrophy and
$b$
is given by
$b=\sqrt{{\it\Omega}_{\mathit{atom }}}/2$
.
Obtaining a combination of
$a$
and
$d$
from the dissipation rate cannot be done uniquely. Given a certain value of
${\it\epsilon}_{\mathit{atom}}$
, we look for a combination of the values of
$a$
and
$d$
satisfying (3.7). Since no parametrization exists to describe the relation between
$a$
and
$d$
we assume all combinations of
$a$
and
$d$
have an equal probability of occurring. A model relating
$a$
and
$d$
could improve the accuracy of our results. This can be achieved using the following transformation:
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle a=\frac{{\it\zeta}_{1}}{\sqrt{3}}-{\it\zeta}_{2},\\ \displaystyle d=\frac{{\it\zeta}_{1}}{\sqrt{3}}+{\it\zeta}_{2},\end{array}\right\} & & \displaystyle\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}{\it\zeta}_{1}={\textstyle \frac{1}{2}}\sqrt{{\it\epsilon}_{\mathit{atom }}}\text{cos}({\it\phi}),\\ {\it\zeta}_{2}={\textstyle \frac{1}{2}}\sqrt{{\it\epsilon}_{\mathit{atom }}}\text{sin}({\it\phi}),\end{array}\right\} & & \displaystyle\end{eqnarray}$$
and where we assume that
${\it\phi}$
is a uniformly distributed between 0 and
$2{\rm\pi}$
.
The equations of motion (3.1) and (3.2) for the particles can be solved analytically, by rewriting them as a single second-order differential equation in the eigenspace of
$\unicode[STIX]{x1D645}$
, i.e. taking
$\boldsymbol{X}={\it\xi}(t)\boldsymbol{e}$
one arrives at
$$\begin{eqnarray}St\,\ddot{{\it\xi}}+\dot{{\it\xi}}-{\bf\lambda}{\it\xi}=0.\end{eqnarray}$$
The general solution is then given by
$$\begin{eqnarray}\boldsymbol{X}(T)=\mathop{\sum }_{i=1}^{3}[A_{i}\text{e}^{{\it\alpha}_{i}^{+}T}+B_{i}\text{e}^{{\it\alpha}_{i}^{-}T}]e_{i},\end{eqnarray}$$
where
$e_{i}$
are the eigenvectors of
$\unicode[STIX]{x1D645}$
and
${\it\alpha}^{\pm }$
are given by
$$\begin{eqnarray}{\it\alpha}_{i}^{\pm }=\frac{-1\pm \sqrt{1+4St{\it\lambda}_{i}}}{2St}.\end{eqnarray}$$
The constants
$A_{i}$
and
$B_{i}$
follow from the initial conditions
$\boldsymbol{X}_{0}$
and
$\boldsymbol{V}_{0}$
.
Consider two particle trajectories
$\boldsymbol{X}_{1}(T)$
and
$\boldsymbol{X}_{2}(T)$
, both with (non-dimensional) radius
$R$
; a collision occurs when
$$\begin{eqnarray}|\boldsymbol{X}_{1}(T)-\boldsymbol{X}_{2}(T)|<2R.\end{eqnarray}$$
In principle, given the initial conditions of two particles, the occurrence of a collision can be computed with (3.14). From a computational point of view however, it is more efficient (due to the exponent in (3.12)) to numerically solve the equations of motion (3.1) and (3.2) directly, which we will do in the next section. Nevertheless, the structure of the solution and in particular the relation between
${\it\alpha}^{\pm }$
and the eigenvalues of
$\unicode[STIX]{x1D645}$
(3.13) will prove to be useful for understanding the collision behaviour.
3.2 Model implementation and validation
The collision probability of particles in the vicinity of such an atom of flow is computed by releasing consecutively a large number of particle pairs, tracking their trajectories using (3.1) and (3.2) and detecting collisions.
The code has been implemented in a multi-GPU (graphics processing unit) framework using a combination of the C++ programming language, CUDA and the MPI-communication protocol. The large number of threads present on a GPU makes it very well suited to compute the behaviour of independent particle pairs in parallel. Equations (3.1) and (3.2) are integrated using a second-order Adams–Bashforth scheme with a user-defined time step to update the particle positions. For the simulations the time step has been set to
$1/100$
th.
A basic validation of the code has been performed in the case of identical particles with infinite mass and non-dimensional radius
$R$
. For the test, one particle is placed in the centre of the atom with zero velocity. The other particle is released at a distance
$\boldsymbol{X}_{0}$
, where
$|\boldsymbol{X}_{0}|\gg R$
and with a random velocity
$\boldsymbol{V}_{0}$
. For
$N_{p}$
particles released, the collision chance is given by
$N_{c}/N_{p}=(R/|\boldsymbol{X}_{0}|)^{2}$
. It has been found that the code statistically converges well to the predicted theoretical value.
For all simulations, particle pairs are released inside a sphere of radius one (i.e. one Kolmogorov scale), with a small random initial velocity, and allowed to travel for 15 Kolmogorov time scales. When the distance between the two particles is smaller than the sum of their radii, a collision is detected and the simulation stops. A graphical impression of such a simulation is shown in figure 3. Particle pairs that have not collided within 15 time scales are most likely swung far away and it is very unlikely they will ever collide.
3.3 Results for the atoms model
The particle trajectories in the presence of the flow structures are updated using (3.1) and (3.2). The particles dynamics are therefore governed by the Stokes number
$St$
and the velocity gradient tensor
$\unicode[STIX]{x1D645}$
of the flow. In this paper we focus on particles with a Stokes number of one. Those particles exhibit a high level of clustering, but still a high correlation exists between flow structure and collision kernel (see figure 2). The second important parameter is the velocity gradient tensor
$\unicode[STIX]{x1D645}$
. We perform a phase-space analysis of the Jacobian to investigate its effect on the particle dynamics. We use a phase-space in which we draw a random set of
$a$
,
$b$
and
$d$
, all between the values of
$-20$
and 20.
Figure 4(a) shows the relation between the eigenvalues and the dissipation rate. The larger the real part of the eigenvalues (in absolute sense), the higher the dissipation rate. Note that each quadrant corresponds to one flow structure category from figure 1. For example the upper left quadrant for which the real part of the complex eigenvalues is negative represents a stable focus (see figure 1
a). Figure 4(b) shows the relation between the eigenvalues and the enstrophy. A correlation can be observed between the imaginary part of the eigenvalues and the enstrophy. Figure 5(a) shows the number of collisions as a function of the eigenvalues. Since we release two particles inside each atom, all with identical volume, the number of collisions is proportional to the collision kernel. For computing the collision kernel, we draw
$10^{8}$
different configurations, and for each configuration we release
$10^{5}$
particle pairs and compute whether a collision occurs. Note that different combinations of
$a$
and
$d$
can yield the same dissipation rate.
(a) Effect of the eigenvalues on the collision probability. The dashed line indicates the critical amount of rotation above which particle are swung outwards due to their momentum, even though the flow spirals inwards. (b) Effect of the dissipation rate and enstrophy on the collision probability.
Let us consider each of the four quadrants separately (which represent all four atoms) and keep in mind conservation of mass (
${\it\lambda}_{1}+{\it\lambda}_{2}+{\it\lambda}_{3}=0$
). For the stable saddle node (see figure 1
c), the larger the convergence, the more collisions occur, especially when the difference between
${\it\lambda}_{1}$
and
${\it\lambda}_{2}$
is small. In such cases, the flow in the atom is symmetrically converging, which appears to be very favourable for collisions. The white voids are a consequence of rearranging the three eigenvalues such that
${\it\lambda}_{1}<{\it\lambda}_{2}<{\it\lambda}_{3}$
. Next consider the stable focus in the top left quadrant (see figure 1
a). Large negative real parts (strong converging motion) also favour collisions, as long as a critical amount of rotation is not exceeded. When this critical amount is exceeded, particles are swung outwards due to their momentum, even though the flow structure is spiralling inwards. Since the third flow direction, dictated by
${\it\lambda}_{3}$
, is directed outwards due to continuity, particles are expelled outwards in all three flow directions. Using the analytical solution of the particles paths in (3.12) and (3.13), the critical amount of rotation can be calculated above which particles are swung outwards (i.e.
$\text{Re}(a^{+})>0$
), although the flow is flowing inwards (i.e.
$\text{Re}({\it\lambda})<0$
). This leads to the following inequality:
$$\begin{eqnarray}\text{Re}(\sqrt{1+4St{\it\lambda}})>1,\end{eqnarray}$$
which can be simplified to
$$\begin{eqnarray}-\!\text{Re}({\it\lambda})<St\,\text{Im}({\it\lambda})^{2}.\end{eqnarray}$$
This relation is shown as a dashed line in figure 5(a).
The unstable focus in the top right quadrant and the unstable saddle node in the bottom right quadrant behave very similarly: both atoms have a small but non-zero chance of producing collisions. When
${\it\lambda}_{1}$
and
${\it\lambda}_{2}$
have a positive real part, particles are swung outwards in two out of three directions, reducing the chances of a collision. The amount of rotation in such a case is of small influence. Collisions can only be induced by the third and negative eigenvalue
${\it\lambda}_{3}$
which makes particles converge towards each other.
The above results are qualitatively in good agreement with the results obtained with the DNS in figure 2. The stable saddle node shows the highest collision kernel. Both the atoms model and the DNS show that the unstable saddle node and the unstable focus have very similar collision kernels. The stable focus possesses the lowest collision kernel in the DNS. Considering that in turbulence flow structures of this category possess more rotation than convergence (e.g. Bijlard et al. Reference Bijlard, Oliemans, Portela and Ooms2010; Elsinga & Marusic Reference Elsinga and Marusic2010), the results are in qualitative agreement with the atoms model, which predicts zero collisions for such a situation.
3.4 Effect of dissipation and enstrophy
Perrin & Jonker (Reference Perrin and Jonker2014) sketched the conceptual picture that particles cluster in regions of low vorticity, but this process increases their velocity coherence, which decreases their collision chances. A dissipative event is needed to decorrelate the particle motion and make them collide. Collisions therefore correlate on average with regions of higher values of dissipation rate. In this section we investigate the effect of dissipation and enstrophy on the collision probability within the atoms of the flow model. To achieve this, we select a value for
${\it\epsilon}_{\mathit{atom}}$
and
${\it\Omega}_{\mathit{atom}}$
from a uniform distribution between 0 and 25 and compute values for
$a$
,
$b$
and
$d$
as explained in § 3.1;
$3\times 10^{6}$
configurations are computed, each with
$10^{5}$
pairs. Figure 5(b) shows the number of collisions obtained from this computation.
Results in figure 5 from the conceptual atoms model show that dissipation is strongly correlated with the number of collisions and that enstrophy is strongly anti-correlated with the number of collisions. When the enstrophy is higher than the dissipation, almost no collisions are occurring, similar to the behaviour observed in figure 5. The sensitive dependence of collisions on dissipation is in agreement with results found by Perrin & Jonker (Reference Perrin and Jonker2014) in turbulent flows. Results presented in the current paper imply that it is the converging aspect of dissipation which favours collisions.
4 Concluding remarks
Using the eigenvalues of the velocity gradient tensor, we divided the flow into four main categories, see figure 1. Physically, complex eigenvalues can be associated with vortices, whereas real eigenvalues correspond to convergence zones. This categorization is very close to the categorization of Chong et al. (Reference Chong, Perry and Cantwell1990) using the PQR invariants of the velocity gradient tensor.
Applying this classification to DNS, the results show a clear preference for particles to preferentially concentrate in regions with no complex eigenvalues and low enstrophy. This behaviour has already been extensively observed in previous studies. The DNS also shows that converging regions not only favour the presence of particles but also have a higher collision kernel. As a consequence, most of the collisions occur in regions with only real eigenvalues.
We have developed a conceptual framework to understand the effect of an individual flow structure on the collision efficiency. Linearization of the equations of motion of heavy particles and scaling with the Kolmogorov scales yields the equations motion of particles in the vicinity of a local flow structure. Since we are investigating individual atoms, this model is not suited for particles with
$St\gg 1$
, where the contribution of particle collisions from caustics is dominant. By releasing particle pairs, we investigated the effect of the local flow configuration on the collision probability of particles. In regions where two out of three eigenvalues have a negative real part (converging motion), most of the particles tend to move to the centre of the flow structure and have a high chance of colliding. Above a critical amount of rotation collisions become extremely rare, because particles move outwards due to their momentum, while the flow spirals inward. This critical amount of rotation can be computed analytically. When two of the eigenvalues have a positive real part (independent of the presence of an imaginary part), particles tend to move away from the centre of the flow structure. In such a case, the third flow direction has a negative eigenvalue and allows some particles to converge to the centre of the flow structure, and collide.
Next we have associated the velocity gradient tensor with a dissipation rate and enstrophy to investigate their effect on the collision chances. The results we obtained were in agreement with previous studies. Dissipation (convergence) favours collisions, whereas vortices (rotation) makes particles less likely to collide.
Acknowledgement
This work is part of the research programme of the Foundation for Fundamental Research on Matter (FOM), which is part of the Netherlands Organisation for Scientific Research (NWO). This work was carried out on the Dutch national e-infrastructure with the support of SURF Foundation.
