Spheres and fibres in turbulent flows at various Reynolds numbers

The sphere motion and forces are modelled using an Eulerian-based immersed boundary method developed by Hori et al. [26]. The velocity $\mathbf{U}$ and rotation rate $\mathbf{\Omega}$ of each sphere are found by integrating the Newton-Euler equations in time

	$\displaystyle m\partial_{t}U_{i}$	$\displaystyle=\oiint_{S}\left(\rho\nu(\partial_{i}u_{j}+\partial_{j}u_{i})-p\delta_{ij}\right)n_{j}dS-F^{col}n_{i},$		(4)
	$\displaystyle I\partial_{t}\Omega_{i}$	$\displaystyle=\oiint_{S}{\varepsilon_{ijk}}\frac{c}{2}n_{j}\rho\nu(\partial_{k}u_{l}+\partial_{l}u_{k})n_{l}dS,$		(5)

where $\delta_{ij}$ is the Kronecker delta and ${\varepsilon_{ijk}}$ is the Levi-Civita symbol, $S$ is the surface of the sphere, and $\mathbf{n}$ is its normal. ${m=\tilde{\rho}\pi c^{3}/6}$ and ${I=mc^{3}/20}$ are the mass and moment of inertia of the sphere with diameter $c$ and density $\tilde{\rho}$. A soft-sphere collision force $F^{col}\mathbf{n}$ is applied in the radial direction when spheres overlap [26].

Table 2 shows our choice of parameters for the sphere-laden flows. In all cases, we use 300 spheres, giving the solid phase a volume fraction of $0.079$. The characteristic time of the spheres is $\tau_{s}=\tilde{\rho}c^{2}/(18\rho\nu)$, from which we can define the Stokes number $St\equiv\tau_{s}/\tau_{f}$ of our flows. The particle Reynolds number of the spheres is $Re_{p}=c\sqrt{\langle\mathbf{\Delta u}_{n}\cdot\mathbf{\Delta u}_{n}\rangle_{n}}/\nu$, where $\mathbf{\Delta u}_{n}$ is the difference between the velocity of the $n$th particle and the local fluid velocity, averaged in a ball of diameter $2c$ centred on the sphere. The angled brackets $\langle\rangle_{n}$ denote an average over all spheres in the simulation.

For the motion and coupling of the fibres, we also use an immersed boundary method, but this one is Lagrangian and solves the Euler-Bernoulli equation for the position $\mathbf{X}$ of the beam with coordinate $l$ along its length [27, 28, 29]

$$\frac{\pi}{4}d^{2}(\tilde{\rho}-\rho)\partial^{2}_{t}X_{i}=\partial_{l}(T\partial_{l}X_{i})+\gamma\partial^{4}_{l}X_{i}-F^{fs}_{i}+F^{col}_{i},$$

(6)

where $T$ is the tension, enforcing the inextensibility condition;

$$\partial_{t}X_{i}\partial_{t}X_{i}=1.$$

(7)

The volumetric density of the fibre is $\tilde{\rho}$, and its stiffness is $\gamma$. Note that the fluid density $\rho$ in equation 6 cancels the inertia of the fluid in the fictitious domain inside the fibre [30]. To exclude deformation effects in our study, we choose a stiffness which limits the fibre deformation below 1%. The collision force $\mathbf{F^{col}}$ is the minimal collision model by Snook et al. [31]. The fluid-solid coupling force $\mathbf{F^{fs}}$ enacts the non-slip condition at the particle surface, and an equal and opposite force ($\mathbf{f^{sf}}$ in equation 1) acts on the fluid. The spreading kernel onto the Eulerian grid has a width of three grid spaces, giving the fibre diameter $d=3\Delta x$ in units of the Eulerian grid spacing [32, equation 14]. The fibre diameter $0.4\eta<d<5\eta$ is on the order of the Kolmogorov length $\eta$ in all cases, i.e., it is smaller than the turbulent eddies in the energy-containing range of the flow. This allows us to consider the fibres as infinitesimally thin objects with a high degree of anisotropy. Table 3 shows our choice of parameters for the fibre-laden flows. In all cases, we use $10^{4}$ fibres, giving the solid phase a volume fraction of around $5.4\times 10^{-3}$. The characteristic time of the fibres is calculated using a formulation which takes their aspect ratio $\beta\equiv c/d$ into account [33],

$$\tau_{s}=\frac{\tilde{\rho}d^{2}}{18\rho\nu}\frac{\beta\ln\left(\beta+\sqrt{\beta^{2}-1}\right)}{\sqrt{\beta^{2}-1}},$$

(8)

from which we can define the Stokes number $St\equiv\tau_{s}/\tau_{f}$ of our flows. The particle Reynolds number of the fibres is $Re_{p}=d\sqrt{\langle\mathbf{\Delta u}_{n}\cdot\mathbf{\Delta u}_{n}\rangle_{n}}/\nu$, where $\mathbf{\Delta u}_{n}$ is the difference between the velocity of the midpoint of the $n$th fibre and the local fluid velocity, averaged in a ball of diameter $2c$ centred on the fibre’s midpoint. The angled brackets $\langle\rangle_{n}$ denote an average over all fibres in the simulation.

The Taylor Reynolds number is defined as $Re_{\lambda}\equiv u_{rms}\lambda/\nu$, where $u_{rms}\equiv\sqrt{\langle u_{i}u_{i}\rangle_{\mathbf{x},t}}/3$ is the root-mean-square velocity, $\lambda=u_{{rms}}\sqrt{15\nu/\epsilon}$ is the Taylor length scale, $\epsilon\equiv 2\nu\langle s_{ij}s_{ij}\rangle_{\mathbf{x},t}$ is the viscous dissipation rate, ${s_{ij}\equiv(\partial_{i}u_{j}+\partial_{j}u_{i})/2}$ is the strain-rate tensor, and angled brackets $\langle\cdot\rangle_{\mathbf{x},t}$ indicate an average over space $\mathbf{x}$ and time $t$.

The Taylor-Reynolds number is an indicator of the intensity of turbulence in the flow, and in figure 2, we show how $Re_{\lambda}$ compares for all of our cases. Figure 2a shows that increasing the solid mass fraction $M$ causes a reduction in $Re_{\lambda}$ at both high and low forcing Reynolds numbers. Both spheres and fibres produce a comparable decrease in $Re_{\lambda}$, but the reduction effect is more substantial for very heavy fibres. Looking at the trend in $Re_{\lambda}$ with the forcing Reynolds number $Re_{ABC}$ in figure 2b, we see, as expected, a monotonic increase in all three cases, i.e., single phase, sphere-laden, and fibre-laden flows.

A second way to quantify turbulence is the turbulent kinetic energy $K^{\prime}\equiv\langle u_{i}^{\prime}u_{i}^{\prime}\rangle_{\mathbf{x},t}/2$. Similarly to Oka and Goto [20], we compute $K^{\prime}$ using the fluctuating part of the fluid velocity $u_{i}^{\prime}(\mathbf{x},t)\equiv u_{i}(\mathbf{x},t)-\langle{u_{i}(\mathbf{x,t})\rangle_{t}}$. Figure 3a shows the dependence of turbulent kinetic energy on particle mass fraction. We see that adding spheres and fibres reduces $K^{\prime}$ relative to the single-phase flow, with fibres having a slightly greater attenuation effect at both high and low forcing Reynolds numbers. From figure 3b, we see that the turbulent kinetic energy shows only a weak dependence on $Re_{ABC}$, the single-phase $(M=0)$ cases remain roughly constant (around $K^{\prime}\approx 6CL$), and flows with fixed (M=1) particles show a large attenuation of $K^{\prime}$ for all $Re_{ABC}$. This agrees with the observations of Peng et al. [22], who made simulations of spherical particles in homogeneous isotropic turbulence and found that the attenuation of turbulent kinetic energy has little dependence on the Reynolds number.

The dissipative anomaly is a well-studied feature in single-phase turbulent flows, first proposed by Taylor [35]. It states that, even in the limit of vanishing viscosity ($Re_{\lambda}\to\infty$), the normalised dissipation $\epsilon\mathcal{L}/u_{{rms}}^{3}$ remains finite. The integral length scale is given by,

$$\mathcal{L}=\frac{\pi}{2u_{{rms}}^{2}}\int_{0}^{\infty}\frac{E}{\kappa}\mathrm{d}\kappa,$$

(9)

where $E$ is the turbulent kinetic energy spectrum and $\kappa$ is the wavenumber. In 2005, Donzis et al. [34] parametrized the dissipative anomaly, based on the mathematical derivation of Doering and Foias [36], they fit the function

$$\frac{\epsilon\mathcal{L}}{u_{{rms}}^{3}}=A\left[1+\sqrt{1+(B/Re_{\lambda})^{2}}\right]$$

(10)

to a number of single-phase flows, obtaining $A=0.2$ and $B=92$. Figure 4 shows the dependence of the normalised dissipation for our flows. We see that as $Re_{\lambda}\to\infty$, flows with spheres and fibres at all mass fractions appear to converge to the same anomalous value of dissipation $\epsilon\mathcal{L}/u_{{rms}}^{3}{\to}0.4$. Hence, we fit equation 10 with $A=0.2$ to each mass fraction of spheres and fibres, obtaining a value for $B$ in each case, shown in the insets of figure 4. The fit to our single phase flows agrees closely with Donzis et al. [34]’s result. However, both spheres and fibres cause an increase in the value of $B$ as their mass fraction increases, with fibres producing roughly double the effect. To understand how spheres and fibres modify the dissipation in the flow, we must look at how they transport energy from large to small scales into the dissipative range.

Figure 5 shows each flow’s turbulent kinetic energy spectrum $E$. Single-phase flows exhibit the canonical Kolmogorov scaling $E\sim\kappa^{-5/3}$ for one or two decades, depending on the forcing Reynolds number. Adding spheres reduces $E$ at wavenumbers up to the sphere diameter ($\kappa<\kappa_{c}$) and increases $E$ in the dissipative range. We mention, in passing, the oscillations at the wavenumber of the sphere diameter; these are an artefact resulting from the discontinuity in the velocity gradient at the sphere boundary [19]. The addition of fibres causes a reduction in $E$ across a broad band of wavelengths ($\kappa L\lesssim 100$) and a pronounced increase in the dissipative range. As was previously seen by Olivieri et al. [10], the energy scaling $E\sim\kappa^{-\beta}$ becomes flatter as the mass fraction $M$ of fibres increases. Figure 5 also shows that both spheres and fibres cause the Taylor length scale to shift to higher wavenumbers due to the increase of the energy in the dissipative scales due to the presence of particles.

Figure 6 shows how each term in the Navier-Stokes equation interacts with the energy spectrum, as expressed by the spectral energy balance,

$$\mathcal{F}_{inj}(\kappa)+\Pi(\kappa)+\Pi_{sf}(\kappa)+{D}(\kappa)=\epsilon,$$

(11)

where $\mathcal{F}_{inj},\Pi,\Pi_{sf},$ and $D$ are the rate of energy transfer by the ABC forcing, advection, solid-fluid coupling, and dissipation, respectively. See the supplementary information of Ref. [37] for a derivation of this equation. For the single-phase flows, energy is carried by the advective term $\Pi$ from large to small scales, where it is dissipated by the viscous term $D$. When particles are added, we see that the solid-fluid coupling term $\Pi_{sf}$ acts as a “spectral shortcut” [38, 10]; it bypasses the classical turbulent cascade, removing energy from large scales and injecting it at the length scale of the particles, through their wakes. In keeping with the spectra, the power of the solid-fluid coupling (shown by the peak value of $\Pi_{sf}$) increases with mass fraction $M$ of both spheres and fibres. For spheres, the coupling $\Pi_{sf}$ dominates only at wavenumbers less than that of the sphere diameter $\kappa<\kappa_{c}$ (i.e., large length scales). Around $\kappa_{c}$, the sphere wake returns energy to the fluid, and the classical cascade resumes for $\kappa>\kappa_{c}$. For fibres instead, $\Pi_{sf}$ extends deep into the viscous range. In fact, markers on the $\Pi_{sf}$ curves show that spheres mostly inject energy around the length scale of their diameter, and fibres mostly inject energy at a length scale between their thickness and length. The images of dissipation in figure 1 support this observation; around the spheres, we see wakes comparable in size to the sphere diameter, while around the fibres, we see wakes comparable in size to the fibre thickness. We presume that the fluid-sphere coupling $\Pi_{sf}$ would extend into the viscous range if the sphere diameter was smaller. Lastly, we consider the effect of Reynolds number on energy transfers in the flow, reducing $Re_{ABC}$ limits the range of scales at which the advective flux $\Pi$ occurs for single-phase flows. However, $Re_{ABC}$ has little effect on the wavenumber range of the solid-fluid coupling term $\Pi_{sf}$, suggesting that the size of the particle wakes is governed mainly by particle geometry, with a lesser effect from fluid properties like viscosity.

Moving from wavenumber space to physical space, we show the longitudinal structure functions

$$S_{p}{(r)}\equiv\langle\left[\hat{r}_{i}u_{i}(\mathbf{x}+\mathbf{r}{,t})-\hat{r}_{i}u_{i}(\mathbf{x}{,t})\right]^{p}\rangle_{{\mathbf{x},\mathbf{\hat{r}},t}}$$

(12)

of each flow in figure 7, where $\mathbf{r}$ is a separation vector between two points in the flow, it has magnitude $r$ and direction $\mathbf{\hat{r}}$, and angled brackets show an average over space $\mathbf{x}$, direction $\mathbf{\hat{r}}$, and time t. For the second moment $(p=2)$, the single-phase flows closely follow Kolmogorov’s scaling [39] in the inertial range ($r\gg\eta$). However, when particles are added, $S_{2}$ decreases relative to the single-phase case. For spheres, the decrease occurs for $r\gtrsim c$, while for fibres, it occurs at even smaller separations $r$. Much like the $E\sim k^{-\beta}$ scalings seen in figure 5, the decreased regions in figure 7 show scalings $S_{2}\sim r^{\zeta_{2}}$, which become flatter for larger mass fractions $M$. These energy spectrum and structure-function scalings are, in fact, just observations of the same phenomenon, as the exponents are related by $\zeta_{2}=\beta-1$ for any differentiable velocity field [40, p.232]. In the higher moment structure function ($p=6$), the flattening of the $S_{6}\sim r^{\zeta_{6}}$ scaling by the particles is yet more apparent, indicating that the tails of the probability distribution of $\hat{r}_{i}u_{i}(\mathbf{x}+\mathbf{r})-\hat{r}_{i}u_{i}(\mathbf{x})$ become wider for smaller separations $r$. In other words, extreme values become more common, and the velocity field becomes more intermittent in space. In the case of single-phase homogeneous-isotropic turbulence, the intermittency of the velocity field has been linked to the non-space-filling nature of the fluid dissipation [41]. Next, we explore this link in the case of particle-laden turbulence by investigating the intermittency of dissipation $\epsilon$ in space.

Using the method described by Frisch [41, p.159], Meneveau and Sreenivasan [42] and others, we obtain the quantity $\epsilon_{l}$ by averaging the viscous dissipation within a spherical region of radius $l$ (hereafter, we refer to this averaging region as a “ball” to avoid confusion with the particles). To observe the variation of $\epsilon_{l}$ we take its square and average over many balls of radius $l$ at different locations and times, where the average is denoted using angled brackets $\langle\cdot\rangle$. Figure 8 shows how $\langle\epsilon_{l}^{2}\rangle$ depends on the radius $l$ of the balls. At large radii $(l\to\infty)$, $\epsilon_{l}$ is by definition equivalent to the bulk value $\epsilon$, and so $\langle\epsilon_{l}^{2}\rangle/\epsilon^{2}$ tends to unity. However, for $l<L$, we see that $\langle\epsilon_{l}^{2}\rangle/\epsilon^{2}>1$ in all of the cases plotted. This increase indicates that there are localised regions of high dissipation in the fluid. In other words, the dissipation fields are intermittent in space. For the single-phase flows, $\langle\epsilon_{l}^{2}\rangle$ reaches a plateau for $l<\eta$, giving an indication of the size of these flow structures. For the sphere-laden flows, we see that $\langle\epsilon_{l}^{2}\rangle$ reaches almost 100 times the bulk value and remains roughly constant for small radii $(l<0.05L)$. This suggests that the spherical particles are creating high dissipation structures in the fluid, and balls with $l<0.05L$ fit entirely inside these structures. That is, reducing the ball radius below $l=0.05L$ has minimal effect on $\langle\epsilon_{l}^{2}\rangle$ because the majority of balls are sampling either entirely inside a high dissipation structure or entirely outside. In the case of fibres, high dissipation structures also produce an increase in $\langle\epsilon_{l}^{2}\rangle$ for small $l$; here no plateau is seen, implying that even the smallest balls $(l=6\times 10^{-3}L)$ cannot fit inside the high dissipation structures created by the fibres.

To measure the fractal dimension of the high dissipation structures, we take a range of moments $-6\leq q\leq 6$ of $\epsilon_{l}$. This quantity has a scaling behaviour

$$\langle\epsilon_{l}^{q}\rangle\sim l^{\tau},$$

(13)

if the dissipation field is multifractal [41]. Indeed, for the sphere-laden flows in figure 8a, $\langle\epsilon_{l}^{2}\rangle$ has a strong dependence on $l$ in the range $0.1L\leq l\leq 0.2L$, and for the fibre-laden flows in figure 8b, $\langle\epsilon_{l}^{2}\rangle$ has a strong dependence on $l$ in the range $0.04L\leq l\leq 0.08L$. Hence, we calculate $\tau$ by making lines of best fit in these ranges, and we obtain the multifractal spectra $F(\alpha)$ with a Legendre transformation

	$\displaystyle\alpha$	$\displaystyle=\frac{d\tau}{dq}+1,$		(14)
	$\displaystyle F$	$\displaystyle=q(\alpha-1)-\tau+3.$		(15)

Figure 9 shows the multifractal spectra for all cases. Similarly to Mukherjee et al. [43], we split our analysis into separate regions: solid curves result from an ensemble average over balls not containing particles, while dashed curves result from an ensemble average over balls containing particles. We see that the multifractal spectra of the single-phase flows and the regions not containing particles have peaks at $\alpha\approx 1,F\approx 3$, showing there is a background of space-filling dissipation in these regions [41, p.163]. In addition, the presence of tails in the spectra shows that the dissipation fields are not self-similar and violate the hypotheses made by Kolmogorov [39] for a turbulent flow. This can explain the flattening of the high-order structure functions ($S_{6}$ in figure 7) relative to Kolmogorov’s predictions. We note in passing that the single-phase flows with lower Reynolds number ($Re_{ABC}=55.9$ and $Re_{ABC}=112$) show narrower tails in figure 9f, as viscous effects are present here.

The multifractal spectra of the regions containing spheres (blue dashed lines in figure 9) have peaks which are shifted to the left. This shift can be explained by considering the high dissipation in the boundary layers around the spheres. As imaged in figure 12, the boundary layers have a higher dissipation rate than the surrounding fluid, and are confined in space to relatively thin sheets. In figure 9e we sketch a ball of radius $l$ intersecting with a sphere’s boundary layer which has thickness $\delta\ll l$. The volume of the high dissipation region inside the ball is $\pi l^{2}\delta$, hence the total dissipation in the ball roughly scales as $\sim l^{2}$. When we average over the volume of the ball (${4\pi l^{3}/3}$) and take the $q$th moment, we obtain $\langle\epsilon_{l}^{q}\rangle\sim l^{2q}l^{-3q}$. Comparing with equation 13, this scales as $\tau=-q$. Applying a Legendre transformation (equations 14 and 15) this corresponds to the point $\alpha=0,F=3$, which is roughly where we find the peaks in the multifractal spectra of the regions containing spheres. This reinforces our observation that the dissipation in the regions around the spherical particles is confined to sheet-like structures.

The multifractal spectra of the regions containing fibres (orange dashed lines in figure 9) have peaks which are shifted even further to the left. Again, we explain this by considering the boundary layers around the particles. In figure 9f we sketch a ball of radius $l$ intersecting with a fibre’s boundary layer which has radius $\delta\ll l$. In this case, the volume of the high dissipation region is $\pi\delta^{2}l$. Thus, the average dissipation in the ball scales as $\langle\epsilon_{l}^{q}\rangle\sim l^{q}l^{-3q}$, so the resulting scaling is $\tau=-2q$, which (by equations 14 and 15) corresponds to the point $\alpha=-1,F=3$. In fact, in figure 9f, we see that the ensemble average over balls containing $M=1$ fibres produces peaks around $\alpha=-0.5$, suggesting that the wakes around the fibres are not exactly thin tubes, but more space-filling structures with a fractal dimension between one and two. This is also supported by the images in figure 12, as wakes are seen to extend behind the fibres.

The middle and upper panels of figure 9 show that (at $Re_{ABC}=224$ and $Re_{ABC}=894$) the sphere- and fibre-containing multifractal spectra tend toward the single phase multifractal spectrum as $M$ is reduced. This is expected as the particle-fluid relative velocity is reduced and the wakes become less prominent compared to the background space-filling turbulence. On the other hand, the ensembles of balls not containing particles produce multifractal spectra with low $\alpha$ tails that extend beyond those of the single-phase multifractal spectra, suggesting that the wakes of spheres and fibres extend into the bulk of the fluid and that they retain their non-space-filling nature.

To look more closely at the flow structures produced by the particles, we compute the alignment of the unit-length eigenvectors $\mathbf{\hat{a}_{1},\hat{a}_{2},\hat{a}_{3}}$ of the strain-rate tensor ${s_{ij}\equiv(\partial_{i}u_{j}+\partial_{j}u_{i})/2}$ with the vorticity ${\boldsymbol{\omega}\equiv{\varepsilon_{ijk}}\partial_{j}u_{k}\mathbf{\hat{e}}_{i}}$ [44, 45], where $\mathbf{\hat{e}_{1},\hat{e}_{2},\hat{e}_{3}}$ are the Cartesian unit vectors, and $\varepsilon_{ijk}$ is the Levi-Civita symbol. We choose $\mathbf{\hat{a}_{1}}$ to be the eigenvector corresponding to the largest eigenvalue. For an incompressible fluid, the three eigenvalues of $s_{ij}$ sum to zero, so the largest eigenvalue is never negative. Hence, $\mathbf{\hat{a}_{1}}$ aligns with the direction of elongation in the flow. Similarly, we choose the eigenvector $\mathbf{\hat{a}_{3}}$ to correspond with the smallest eigenvalue, which is never positive, so $\mathbf{\hat{a}_{3}}$ aligns with the direction of compression in the flow. Finally, due to the symmetry of the strain-rate tensor, $\mathbf{\hat{a}_{2}}$ is orthogonal to the other two eigenvectors and, depending on the flow, there can be compression or elongation along its axis. Figure 10 shows probability-density functions $P$ of the scalar product of the unit-length vorticity $\boldsymbol{\hat{\omega}}$ with $\mathbf{\hat{a}_{1},\hat{a}_{2}}$ and $\mathbf{\hat{a}_{3}}$ for all of the flows studied. The quantity $\boldsymbol{\hat{\omega}}\cdot\mathbf{\hat{a}_{i}}$ is simply the cosine of the angle between the two vectors, and we take the modulus because $\mathbf{\hat{a}_{i}}$ and $-\mathbf{\hat{a}_{i}}$ are degenerate eigenvectors. In the single-phase flows, $P$ is maximum at $|\boldsymbol{\hat{\omega}}\cdot\mathbf{\hat{a}_{2}|}=1$, that is, we see a strong alignment of vorticity with the intermediate eigenvector $\mathbf{\hat{a}_{2}}$. This is a well-known feature of turbulent flows and has been attributed to the axial stretching of vortices [44]. Also in keeping with previous observations of single-phase turbulence, we see that the first eigenvector $\mathbf{\hat{a}_{1}}$ shows very little correlation with $\boldsymbol{\hat{\omega}}$, and the last eigenvector is mostly perpendicular to the vorticity, producing a peak in $P$ at $\boldsymbol{\hat{\omega}}\cdot\mathbf{\hat{a}_{3}}=0$. When spheres are added, vorticity aligns even more strongly with the intermediate eigenvector $\mathbf{\hat{a}_{2}}$, and becomes more perpendicular to the first and last eigenvectors $\mathbf{\hat{a}_{1}}$ and $\mathbf{\hat{a}_{3}}$. This change is consistent with the existence of a shear layer at the surface of the spheres, or in the wakes behind the spheres. To demonstrate this, we draw a laminar shearing flow and label the directions of $\mathbf{\hat{a}_{1},\hat{a}_{2},\hat{a}_{3}}$ and $\boldsymbol{\hat{\omega}}$ in the inset of figure 10. The compression $\mathbf{\hat{a}_{3}}$ and extension $\mathbf{\hat{a}_{1}}$ are in the plane of the shear flow, while the vorticity $\boldsymbol{\hat{\omega}}$ and the intermediate eigenvector $\mathbf{\hat{a}_{3}}$ are perpendicular to the shear plane; this gives rise to the $\boldsymbol{\hat{\omega}}\cdot\mathbf{\hat{a}_{1}}=0,$ $|\boldsymbol{\hat{\omega}}\cdot\mathbf{\hat{a}_{2}}|=1$ and $\boldsymbol{\hat{\omega}}\cdot\mathbf{\hat{a}_{3}}=0$ modes in the sphere-laden flows in figure 10. The larger mass fraction spheres have a greater effect, as their motion relative to the fluid is larger. Also, the spheres’ effect is more pronounced at lower $Re_{ABC}$, which can be explained by an increased thickness in the shear layer as the fluid viscosity increases. For fibre-laden flows, peaks are also seen at $\boldsymbol{\hat{\omega}}\cdot\mathbf{\hat{a}_{1}}=0,$ $|\boldsymbol{\hat{\omega}}\cdot\mathbf{\hat{a}_{2}}|=1$ and $\boldsymbol{\hat{\omega}}\cdot\mathbf{\hat{a}_{3}}=0$, but the peaks are spread over a larger range of angles, presumably due to the high curvature of the fibre surfaces, which adds variance to the local direction of vorticity and shear. At small Reynolds numbers ($Re_{ABC}=55.9$ and $112$) the peaks are widest, and $P(|\boldsymbol{\hat{\omega}}\cdot\mathbf{\hat{a}_{i}}|)$ becomes an almost uniform distribution.

The local topology of a flow can be described entirely using the three principle invariants of the velocity gradient tensor $\partial_{i}u_{j}$ [46]. The first invariant $\partial_{i}u_{i}$ is not interesting in our case because it is zero for an incompressible fluid. However, the second invariant,

$$Q=\frac{1}{4}\omega_{i}\omega_{i}-\frac{1}{2}s_{ij}s_{ij},$$

(16)

expresses the balance of vorticity and strain, while the third invariant,

$$R=\frac{1}{4}\omega_{i}s_{ij}\omega_{j}-\frac{1}{3}s_{ij}s_{jk}s_{ki},$$

(17)

is the balance of the production of vorticity with the production of strain [47]. The eigenvalues $\Lambda$ of $\partial_{i}u_{j}$ are the roots of the polynomial equation

$$\Lambda^{3}+Q\Lambda+R=0.$$

(18)

Figure 11 shows the joint probability distributions $p(R,Q)$ for all of our flows. We also plot a line where the discriminant of equation 18 is zero: $27R^{2}/4+Q^{3}=0$. Above this line, $\partial_{i}u_{j}$ has one real and two complex eigenvalues and vortices dominate the flow. Below this line, all three eigenvalues are real, and the flow is dominated by strain. The single-phase distributions have tails in the top-left and bottom-right quadrants; these correspond to stretching vortices and regions where the flow compresses along one axis, respectively [46]. When spheres are added, both $Q$ and $R$ are reduced, as the non-slip boundary condition at their surfaces dampens the flow structures. Incidentally, a similar $Q$ and $R$ reduction has been seen with droplets [48]. However, we see that the addition of spheres causes $p(R,Q)$ to increase in the shaded regions in the left panels of figure 11. These regions are defined by

$$\frac{27R^{2}}{4}+Q^{3}<-\frac{4\epsilon^{3}}{\nu^{3}}.$$

(19)

These shaded regions correspond to a strain-dominated flow. In fact, for pure strain, we have $Q=-s_{ij}s_{ij}/2$ and $R=0$. Hence we can substitute for $Q$ and $R$ in equation 19 to estimate that the dissipation is

$${2\nu s_{ij}s_{ij}=-4\nu Q<4^{4/3}\epsilon\approx 6\epsilon}$$

(20)

in the straining flow. That is, the dissipation in regions around the spheres is roughly six times greater than the mean dissipation. Fibres also create straining regions in the flow, seen by the increased probability density functions in the lower parts ($27R^{2}/4+Q^{3}<0$) of the right-hand panels in figure 11. As the fibre mass fraction is increased, the variance of $R$ is reduced, showing the production of strain and vorticity is suppressed by the fibres. At high fibre mass fractions the distributions are approximately symmetrical in $R$, in other words, $Q$ and $R$ are uncorrelated (much like the decorrelation of vorticity and strain observed for fibre-laden flows in figure 10). This suggests that the small fibre diameter produces many small scale structures which disrupt the typical turbulent flow. For both particle types, increasing the Reynolds number $Re_{ABC}$ reduces the effect of the particles.

In figure 12 we visualise isosurfaces where the local fluid dissipation is six times its bulk value, i.e., where $2\nu s_{ij}s_{ij}=6\epsilon$. As predicted by equation 20, the dissipation is greater than $6\epsilon$ in the near-particle regions. And in keeping with our measurements of the fractal dimension of the dissipation field (figure 9), the high dissipation region near the sphere is approximately sheet-like and two-dimensional, while the high dissipation regions near the fibres are approximately tube-like and one-dimensional. In both cases, wakes extend behind the particles, which contribute to the fractal dimension of the dissipation fields. We see the single-phase dissipation structures take various shapes and sizes, giving rise to the multifractal behaviour discussed above. In figure 12 we colour the dissipation isosurfaces according to the local value of the flow topology parameter [49],

$$\Sigma\equiv\frac{2s_{ij}s_{ij}-\omega_{i}\omega_{i}}{2s_{ij}s_{ij}+\omega_{i}\omega_{i}},$$

(21)

which is essentially the second invariant of the velocity gradient tensor $Q$ normalised so that it is bounded between -1 and 1. Regions where $\Sigma=-1$ (coloured green) correspond to pure rotational flow, regions where $\Sigma=0$ (coloured white) correspond to pure shear, and regions where $\Sigma=1$ (coloured purple) correspond to pure straining flow. From figure 12 we see fibres produce straining flow on the upstream side, and shearing flow in wakes on the downstream side, this manifests in figure 11 as an increase in the spread of $Q$ for the fibre-laden flows. At the surface of the sphere we see a mostly shearing flow, which supports our observations of strong shear-vorticity alignment $(|\boldsymbol{\hat{\omega}}\cdot\mathbf{\hat{a}_{2}}|=1)$ in figure 10.