3.1.1 Mean velocity magnitude distribution

The compactness and confinement of the geometry creates strong acceleration and deceleration of the flow within the pore space. To visualize the flow field, contours of the normalized temporal mean velocity magnitude ( $\overline{U}_{m}$ ) in the centre $xy$ -plane are plotted in figures 2(a), 2(c), and 2(e) for each Reynolds number, where $\overline{U}_{m}$ is calculated as $\sqrt{\overline{u_{x}}^{2}+\overline{u_{y}}^{2}+\overline{u_{z}}^{2}}$ and normalized by the mean interstitial velocity $\overline{\langle u_{x}\rangle ^{\unicode[STIX]{x1D6FE}}}$ . An almost identical distribution of the normalized fields for $\overline{U}_{m}$ is observed regardless of the Reynolds number. The mean flow enters the centre pore from the left-hand side corners (regions labelled by ( $\text{i}_{1}$ ) and ( $\text{i}_{2}$ ) in figure 2 c) and accelerates to region (ii) due to the geometric constrictions; then the mean flow at the centre starts to decelerate as it encounters the particle at region (iii), similar to an impinging jet. The flow then accelerates in the spanwise directions and leaves the pore through the right-hand side corners (regions ( $\text{iv}_{1}$ ) and ( $\text{iv}_{2}$ )). The acceleration and deceleration close to the right-hand side sphere can be seen as transfer of momentum from the streamwise to the spanwise directions due to the bounding of the solid phase. The pattern of mean velocity distribution indicates how the mean flow in the pore is affected by the geometry. High Reynolds number flow features, such as wake and large scale vortex shedding in an external flow behind a spherical particle, are not observed in this low porosity bed. A weak wake-like structure (small negative mean velocity just behind the left-hand side sphere near region (ii) in figure 2 c) is present, but is confined to a small region. The closely packed particles in the low porosity FCC configuration tend to break down large scale flow structures and prevent the generation of a significant wake region. This is in contrast to the high porosity studies of square cylindrical porous structures conducted by Chu et al. (Reference Chu, Weigand and Vaikuntanathan2018), where distinct wakes and flow separation were identified.

Figures 3(a), 3(c), and 3(e) show the distributions of the mean velocity magnitude on the centre $yz$ -plane. On this slice, the mean flow is moving perpendicularly into the page and there is a distinguished region with high velocity magnitude near each corner of the centre pore. Similar to what is observed from figures 2(a), 2(c) and 2(e), the mean velocity field does not change significantly at different flow rates.

3.1.2 TKE distribution

The TKE distribution presents a strong dependency on the Reynolds number, which is different from the uniform distribution of the mean velocity field described in the previous section. The turbulent kinetic energy $k$ , defined as $(1/2)(\overline{{u_{x}^{\prime }}^{2}}+\overline{{u_{y}^{\prime }}^{2}}+\overline{{u_{z}^{\prime }}^{2}})$ , is normalized by the square of the mean interstitial velocity $(\overline{\langle u_{x}\rangle ^{\unicode[STIX]{x1D6FE}}})^{2}$ . Here the temporal fluctuation velocity $u_{x}^{\prime }$ is defined as $u_{x}-\overline{u_{x}}$ . Figures 2(b), 2(d) and 2(f) illustrate the distribution of normalized TKE on the centre $xy$ -plane. Substantial disparities in the TKE distribution can be observed at different Reynolds numbers. Firstly, the overall magnitude of the TKE increases with $Re_{H}$ . Secondly, with increasing flow rate, the TKE values in the centre region of the pore remain substantially high. It also suggests that the flow through the pore, although bounded by curved particle walls, is different from simple channel/duct flows even with complex boundaries (Orlandi, Modesti & Pirozzoli Reference Orlandi, Modesti and Pirozzoli2018) wherein the TKE reaches a peak value near the boundaries and then decreases in the centre region.

The normalized TKE distributions on the centre $yz$ -plane shown in figures 3(b), 3(d) and 3(f) indicate similar trends: the magnitude of the TKE increased with increasing $Re_{H}$ ; the distribution pattern changes noticeably as well. In this slice, some similarities with duct flows are noticed. Away from the wall, the TKE increases and reaches a peak value; and then decreases as we approach the core region of the pore. The homogeneous particle packing involving four particles aligned together seem to form a locally duct-like flow pattern in this section. However, it disappears in other sections away from the centre $yz$ -plane, owing to the three-dimensional nature of the spherical particles.

3.1.3 TKE budget

To better understand the mechanism of TKE transport at the pore-scale, the TKE budget including production, dissipation and all the transport terms in the TKE transport equation (equation (3.1), where $k$ is the turbulent kinetic energy) are explicitly computed. The balance between the gain and loss of kinetic energy in the budget has been computed locally as well as globally for the entire flowfield. The average error of the TKE budgets are within $5\,\%$ for all three cases. The transport equation for turbulent kinetic energy is given by

(3.1) $$\begin{eqnarray}\displaystyle \underbrace{\frac{\unicode[STIX]{x2202}k}{\unicode[STIX]{x2202}t}}_{\mathit{Unsteady\,term}}+\underbrace{\overline{u_{j}}\frac{\unicode[STIX]{x2202}k}{\unicode[STIX]{x2202}x_{j}}}_{\mathit{Mean\,advection}} & = & \displaystyle \underbrace{-\overline{u_{i}^{\prime }u_{j}^{\prime }}\frac{\unicode[STIX]{x2202}\overline{u_{i}}}{\unicode[STIX]{x2202}x_{j}}}_{\mathit{Production}\,(P)}-\underbrace{\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D6FE}}\overline{\frac{\unicode[STIX]{x2202}u_{i}^{\prime }}{\unicode[STIX]{x2202}x_{j}}\frac{\unicode[STIX]{x2202}u_{i}^{\prime }}{\unicode[STIX]{x2202}x_{j}}}}_{\mathit{Dissipation}\,(\unicode[STIX]{x1D716})}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}\left[\underbrace{\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D6FE}}\frac{\unicode[STIX]{x2202}k}{\unicode[STIX]{x2202}x_{j}}}_{\mathit{Viscous\,transport}}-\underbrace{\overline{\frac{1}{2}u_{i}^{\prime }u_{i}^{\prime }u_{j}^{\prime }}}_{\mathit{Turbulent\,transport}}-\underbrace{\frac{1}{\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FE}}}\overline{p^{\prime }u_{j}^{\prime }}}_{\mathit{Pressure\,transport}}\right],\quad\end{eqnarray}$$

where $P$ represents the production and $\unicode[STIX]{x1D716}$ the dissipation rate of TKE, respectively. TKE production normally serves as the source term in the TKE transport equation. The production term arising from the Reynolds stress working on the background mean shear extracts energy from the mean flow field and is typically expected to be positive. The dissipation term represents the transfer of the kinetic energy to heat under the effect of viscosity. Thus, it acts as a sink term in the TKE transport equation. If production and dissipation are in balance, the turbulence is in an equilibrium state.

Figure 4.

Distribution of the normalized net production of TKE ( $P$ ) and its normal components ( $P_{kk}$ ): the left column shows contours of normalized TKE production on the centre $xy$ -plane; and the right column shows variation of the three normal components of the production along the centre line. (a,b) $Re_{H}=300$ ; (c,d) $Re_{H}=500$ ; (e,f) $Re_{H}=1000$ .

Figures 4(a), 4(c) and 4(e) show the distribution of the net normalized TKE production ( $P$ ) in the centre $xy$ -plane for the three Reynolds numbers of 300, 500 and 1000, respectively. For the centre pore, the peak TKE production occurs just behind the left particle on either side of the centre line (regions $\text{ii}_{1}$ and $\text{ii}_{2}$ in figure 4 c). The size of the peak production region also increases with increasing Reynolds number. More importantly, it is clear from the contour plots that negative values of net production are observed near the surfaces of the particles (regions iii, $\text{iv}_{1}$ and $\text{iv}_{2}$ ) especially near the flow exit regions. In most turbulent flows, the net TKE production term mainly remains positive, indicating that the energy is extracted from mean flow and fed into turbulence. However, this term being negative for the present cases, implies the unusual condition that the energy transport process is reversed. Turbulence is providing energy to maintain the mean flowfield. It is rare to observe such phenomenon and is attributed to the complex nature of the tortuous flow paths created by the porous geometry.

To further understand such phenomenon, the axial variation of the normal components of production, $P_{11}$ , $P_{22}$ and $P_{33}$ , are plotted along the centre line for the three different Reynolds numbers in figures 4(b), 4(d) and 4(f), respectively. The normal components of production, $P_{kk}$ , computed as $-\overline{({u_{k}}^{\prime })^{2}}(\unicode[STIX]{x2202}\overline{u_{k}}/\unicode[STIX]{x2202}x_{k})$ , are dominant along the centre line. At the edge of the left sphere ( $x=0.35$ ), all components of production are zero. As we move away from the left sphere surface, the axial component $P_{11}$ is positive, whereas the spanwise components, $P_{22}$ and $P_{33}$ , are negative. Although individual components of TKE production can be either positive or negative, the net summation of these components turns out to be negative in this small wake region directly behind the left sphere. Moving further to the right, the axial component turns to be negative and the other two components become positive, all peaking at about $x=0.42$ , reaching almost zero values at the midway point $x=0.5$ , and then switch signs with a peak at about $x=0.61$ . The net production of TKE is positive in this centre region, $0.38\leqslant x\leqslant 0.61$ . Approaching the stagnation point near the right sphere, $P_{11}$ is positive but $P_{22}$ and $P_{33}$ are negative, such that the net production becomes negative. It is also noteworthy that, $P_{22}$ and $P_{33}$ are mostly identical owing to the symmetry in the flow geometry.

In the negative production regions, the spanwise components, $P_{22}$ and $P_{33}$ are negative with the streamwise component $P_{11}$ being positive. In order to reveal the underlying physics related to this phenomenon, the production term is decomposed for more detailed examination. Here $P_{ij}$ is defined as the negative product of Reynolds stress $\overline{u_{i}^{\prime }u_{j}^{\prime }}$ multiplied by the mean strain rate $\unicode[STIX]{x2202}\overline{u_{i}}/\unicode[STIX]{x2202}x_{j}$ , where the normal components of Reynolds stress $\overline{{u_{i}}^{^{\prime }2}}$ are always positive. So the positivity of the normal TKE production components is exclusively determined by the mean rate of strain, i.e. whether the mean flow is accelerating or decelerating.

The net production near the exit (regions $\text{iv}_{1}$ and $\text{iv}_{2}$ in figure 4 d) is discussed first. As shown in figure 2(c), the mean flow moves in the axial direction through the centre region of the pore and is decelerated by the particle surface on the right (region iii) similar to an impinging jet, then merges with the rest and exits the pore through the crevices. During this process, the flow is accelerating in spanwise directions, and decelerating in streamwise direction. As a result, the mean strain rates in $y$ - and $z$ -directions are positive, while negative in $x$ -direction. Together with the positive normal Reynolds stresses, the sum of negative production in the spanwise directions outnumbers that in the streamwise direction, which causes the net production to be negative near the stagnation point. A similar observation was reported by Nishino et al. (Reference Nishino, Samada, Kasuya and Torii1996) who observed negative production in the stagnation zone of an impinging jet.

In the entrance region (regions $\text{i}_{1}$ and $\text{i}_{2}$ in figure 4 d), the mechanism is similar except that there is no jet impingement. The flow goes into the pore through the corners on the left side of the unit cell and merges as it approaches the centre of the pore. Thus, the flow is accelerating in the spanwise directions and a confined wake near the left sphere is created by the merging flow, which leads to the negative mean strain rate in $x$ -direction; this results in a negative net production of TKE.

In summary, the mean flow acceleration and deceleration together with normal components of Reynolds stress are mainly responsible for the negative net production. Since the normal components of Reynolds stress are positive by definition, the locations of the acceleration and deceleration of mean flow actually determines where the negative production may happen. Observing from figures 2(a), 2(c) and 2(e), the mean flow field changes little with respect to Reynolds number. Consequently, regions with $-P$ are expected to locate at the same positions in the flowfield regardless of Reynolds number.

Figure 5 shows all the terms in the TKE transport equation along the centre line for the three values of $Re_{H}$ , respectively. Each term of the TKE budget, normalized by $(\overline{\langle u_{x}\rangle ^{\unicode[STIX]{x1D6FE}}})^{3}/D_{B}$ , is plotted, except the local unsteady term, which is essentially zero for stationary flows. As explained before, close to the surfaces of the left and right spheres, the net production term ( $P$ ) becomes negative giving rise to negative production. The dissipation term ( $-\unicode[STIX]{x1D716}$ ) peaks near the particle surfaces and is negative over the entire pore along the centre line. Moreover, production is not equal to the dissipation over most of the flow domain. Thus, the majority of the turbulent flow within the pore is not at an equilibrium state. The viscous transport is mostly zero in the centre region of the pore, but increases near the particle surfaces as expected. Although the fluctuation velocity $u_{i}^{\prime }$ is zero at the wall, the fluctuating strain rate $\unicode[STIX]{x2202}u_{i}^{\prime }/\unicode[STIX]{x2202}x_{j}$ is not, resulting in the non-zero viscous transport terms.

It is also interesting to note that the contribution from the turbulent transport term to the TKE budget is significant over the entire pore. Consequently, the triple correlation, which can be taken as the fluctuation velocity components working with the Reynolds stress, is redistributing the TKE over the entire flow regime. This is different from that of a channel or a duct flow wherein the contribution of the turbulent transport term is not significant (e.g. Pope Reference Pope2000). Such observations may be attributed to high turbulence intensities over the entire pore region as opposed to the wall-bounded regions of a duct flow, where the turbulent transport term is significant only near the wall.

In the negative production regions, the pressure transport term shows dominant and positive peaks for all Reynolds numbers. It is evident that the pressure plays an important role for the redistribution of TKE near the wall. The pressure transport term moves TKE close to the wall, especially to the wall on the right-hand (downstream) side of the pore. For $Re_{H}=500$ and $1000$ , the turbulent transport remains positive, and the pressure transport terms show a negative peak. For $Re_{H}=300$ , however, near the region where the net production shows a positive peak, the turbulent transport is found to be negative whereas the pressure transport is positive. This may be attributed to the fact that the flow at $Re_{H}=300$ is not fully developed into a turbulent flow. The smaller magnitude fluctuations in the velocity and pressure fields are likely to be causing these differences.

Figure 5.

The TKE budgets along the centre line at $xy$ -plane: the left column shows the production, negative dissipation and mean advection terms; the right column the viscous transport, turbulent transport and pressure transport terms. (a) and (b) are for $Re_{H}=300$ , (c) and (d) $Re_{H}=500$ and (e) and (f) $Re_{H}=1000$ .

3.3.1 Spatial variance

Following Hill & Koch (Reference Hill and Koch2002), the spatial variances for the instantaneous velocity ( $R_{ij}$ ) and vorticity ( $O_{ij}$ ) are defined as

(3.7) $$\begin{eqnarray}\displaystyle & \displaystyle R_{ij}(t)=\langle \tilde{u} _{i}\tilde{u} _{j}\rangle ^{\unicode[STIX]{x1D6FE}}/(\langle u_{x}\rangle ^{\unicode[STIX]{x1D6FE}})^{2}, & \displaystyle\end{eqnarray}$$

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle O_{ij}(t)=\langle \tilde{\unicode[STIX]{x1D714}}_{i}\tilde{\unicode[STIX]{x1D714}}_{j}\rangle ^{\unicode[STIX]{x1D6FE}}(D_{B}/2)^{2}/(\langle u_{x}\rangle ^{\unicode[STIX]{x1D6FE}})^{2}. & \displaystyle\end{eqnarray}$$

The streamwise and spanwise components of the variances are given by the diagonal part of $R_{ij}$ and $O_{ij}$ , which are defined as

(3.9a,b ) $$\begin{eqnarray}\displaystyle R_{||}=R_{11}\quad \text{and}\quad R_{\bot }=(R_{22}+R_{33})/2, & & \displaystyle\end{eqnarray}$$

(3.10a,b ) $$\begin{eqnarray}\displaystyle O_{||}=O_{11}\quad \text{and}\quad O_{\bot }=(O_{22}+O_{33})/2, & & \displaystyle\end{eqnarray}$$

where the subscript $\Vert$ represents streamwise and $\bot$ spanwise. As for the relative helicity $\cos (\unicode[STIX]{x1D6F7})$ , the variance $C$ is calculated as

(3.11) $$\begin{eqnarray}\displaystyle C(t)=\langle \cos ^{2}(\unicode[STIX]{x1D6F7})\rangle ^{\unicode[STIX]{x1D6FE}}-(\langle \cos (\unicode[STIX]{x1D6F7})\rangle ^{\unicode[STIX]{x1D6FE}})^{2}. & & \displaystyle\end{eqnarray}$$

Table 5 shows the time-averaged spatial variances for the velocities and vorticities in both streamwise and spanwise directions, as well as the relative helicity. A reasonable agreement of the variances at $Re_{H}=300$ is found between the current study and Hill & Koch (Reference Hill and Koch2002). The small differences may be attributed to (i) the time-averaging, which may have removed some of the temporal fluctuations and was not performed by Hill & Koch (Reference Hill and Koch2002), and (ii) the grid resolution, which was much coarser in Hill & Koch (Reference Hill and Koch2002) so that some very small velocities near the solid wall may not be completely resolved in their work.

Table 5.

Time-averaged spatial variances of velocity and vorticity in the streamwise and spanwise directions together with the relative helicity at three different Reynolds numbers; the overbar means time-averaging; the first line are the data reported by Hill & Koch (Reference Hill and Koch2002) at the same Reynolds number of $Re_{H}=300$ at one time instant.

Table 5 indicates that the variance values for velocity in the streamwise direction are smaller than that in the spanwise directions. This is attributed to the fact that the mean flow is strong in the streamwise direction, while almost zero in the other directions. For the vorticity components in both streamwise and spanwise directions, the variances are increasing almost linearly with $Re_{H}$ , and the values are essentially equal regardless of the direction. This is quite different for the steady inertial and weakly turbulent flows reported in Hill & Koch (Reference Hill and Koch2002). They observed that the variance of vorticity in spanwise directions is much larger than that in the streamwise direction. For steady inertial flows, the vortices are primarily initiated by shear close to the particle surfaces. Thus the cross-flow components are much stronger. For the unsteady flows, however, vortex shedding becomes important, stretching the flow structures and initiating the generation of vortices in all directions, even away from the solid surface. In the fully turbulent flows, this mechanism is greatly magnified. Finally, the variance of relative helicity roughly remains constant for the high Reynolds number turbulent flows studied here. This is also distinctly different compared to the observations of Hill & Koch (Reference Hill and Koch2002). They reported that the variance increases greatly from the Darcy flow regime to the unsteady flow regime. To understand the above trends better, the p.d.f.s of velocity, vorticity and relative helicity are computed and plotted below.

Figure 8.

The normalized p.d.f.s for velocity and vorticity for different Reynolds numbers: (a) $x$ -component of velocity, (b) $y$ -component of velocity, (c)  $x$ -component of vorticity, and (d) $y$ -component of vorticity. The dotted green line is for $Re_{H}=300$ , red $Re_{H}=500$ and blue $Re_{H}=1000$ . Here ○, ♢, and ▵ denote the data from Hill & Koch (Reference Hill and Koch2002) for $Re_{H}=300$ , $123$ , and $70$ , respectively. The black dash-dot line represents the Gaussian distribution with zero mean and unit standard deviation. The black dashed line stands for the double-exponential distribution.

3.3.2 Probability density functions

The normalized probability density functions for the velocities, vorticities and relative helicity are provided in this section. The spatial variances in table 5 are used to normalize the p.d.f.s making their standard deviations unity and allowing direct comparison with the Gaussian or double exponential distributions. For the plots in this section, the subscript ‘n’ in the titles of the abscissas, for example $\text{}_{n}u_{x}$ , means that the abscissa has been normalized by its spatial variance. Owing to the periodicity, the statistics in the $y$ - and $z$ -directions should be the same (which was verified) and hence only the p.d.f.s of velocity in the $y$ -direction are shown.

The p.d.f.s for the instantaneous velocity and vorticity components in the streamwise and spanwise directions are given in figure 8 for the three Reynolds numbers. Data from Hill & Koch (Reference Hill and Koch2002) is also included for the lower Reynolds number steady inertial or weakly turbulent flows for comparison. It can be observed that the p.d.f.s for lowest Reynolds number of $300$ studied in the present work are in reasonable agreement for the corresponding p.d.f.s from Hill & Koch (Reference Hill and Koch2002).

The p.d.f.s for velocities seem to have similar behaviour over a range of Reynolds numbers after proper normalization. This suggests that, at least in the sense of the first order statistics, the main topology does not change much with respect to the Reynolds number, which agrees qualitatively with the observation in § 3.1.1. Consistent with observations of Hill & Koch (Reference Hill and Koch2002), more $x$ -component velocities in the upstream direction are detected for higher $Re_{H}$ . This can be seen at the right tails of p.d.f.s especially in figure 8(a). The stronger recirculation, as well as the large number of vortices generated with one of the components perpendicular to the $x$ -direction, are potentially responsible for this observation. In addition, a near-exponential decay (at a rate faster than the Gaussian distribution) is observed at the right-hand side tail of the velocity p.d.f.s, especially for lower Reynolds numbers; and the decaying rate is closer to the Gaussian distribution at higher Reynolds numbers. The p.d.f.s of instantaneous velocities shown in figure 8(b) generally follow the Gaussian distribution until about three times the standard deviation in both positive and negative abscissas, then continue to decay, but at rates noticeably larger than the Gaussian function. This indicates that the decay rates of the tails of the velocity field are Reynolds-number-dependent for the range considered in the present study. Hill & Koch (Reference Hill and Koch2002) conjectured that the exponential decay at lower Reynolds numbers is because of the distribution of the pore size and is not dependent upon the Reynolds number. However, as is shown in the present study, with increasing Reynolds numbers, the flow structures and length scales become smaller even in the main flow direction. This isotropization of the velocity field at higher $Re_{H}$ is probably the cause of the more Gaussian-distributed p.d.f.s observed.

The p.d.f.s for the vorticity components in the $x$ - and $y$ -directions for different Reynolds numbers are shown in figures 8(c) and 8(d), respectively. The p.d.f.s for steady inertial flow at $Re_{H}=123$ (♢) from Hill & Koch (Reference Hill and Koch2002) is also shown for comparison. It can be observed that the p.d.f.s for the turbulent flows studied in the present work are quite different from those for the steady inertial flow reported in Hill & Koch (Reference Hill and Koch2002). For the turbulent flows, the p.d.f.s for $\text{}_{n}\unicode[STIX]{x1D714}_{x}$ are almost decaying exponentially until about five times the standard deviation regardless of $Re_{H}$ . This trend is even more pronounced for the p.d.f.s of the vorticities in the spanwise directions, which follow the double exponential function till about seven times the standard deviation. While for Hill and Koch’s steady inertial flow, the p.d.f.s first follow the double exponential distribution until about three times the standard deviation, then start to decay in a rate closer to the Gaussian function. Such near-Gaussian decay is almost entirely absent from any of the unsteady flows studied in the current work.

Figure 9.

The p.d.f. of normalized relative helicity. The dotted green line is for $Re_{H}=300$ , red $Re_{H}=500$ and blue $Re_{H}=1000$ ; ○ illustrate the data from Hill & Koch (Reference Hill and Koch2002) for $Re_{H}=300$ .

Finally, the p.d.f.s for the relative helicity are displayed in figure 9. Note that the abscissa of figure 9 has not been scaled by the standard deviation because the p.d.f.s for relative helicity are periodic with respect to the angle ( $\unicode[STIX]{x1D6F7}$ ). At $Re_{H}=300$ , there is a distinct peak in the p.d.f. at $\cos \unicode[STIX]{x1D6F7}=0$ , indicating that the vorticity and velocity vectors are not aligned for the majority of the flow region. Comparing the p.d.f.s for $\cos \unicode[STIX]{x1D6F7}$ at $Re_{H}=300$ from the current study and Hill & Koch (Reference Hill and Koch2002), some deviations could be observed, especially at large magnitude of $\cos \unicode[STIX]{x1D6F7}$ . This difference could be due to the fact that the p.d.f. reported in Hill & Koch (Reference Hill and Koch2002) is only at one time instant, while it has been time-averaged in the present study. Unlike the velocity and vorticity, the p.d.f.s of relative helicity are strongly dependent on the Reynolds number. It can be observed from figure 9 that higher Reynolds-number flows tend to have higher probability of smaller angles between the velocity and vorticity vectors, aligning them to create helical structures. It indicates the mechanism by which the turbulent kinetic energy is transferred from large to small scale structures by helical motions and stretching.

Figure 10 shows instantaneous stream-ribbons coloured by the relative helicity for $Re_{H}=1000$ . It is observed that substantial helicity is generated in the initial part of the centre pore and propagated downstream, where the stream ribbons are coloured either by orange ( $\cos (\unicode[STIX]{x1D6F7})=1$ ) or blue ( $\cos (\unicode[STIX]{x1D6F7})=-1$ ), confirming the alignment of velocity and vorticity. With such alignment, the vorticity is advected downstream in the streamwise direction by the velocity giving rise to the helical motions.

Figure 10.

Instantaneous stream ribbons coloured by the relative helicity ( $\cos (\unicode[STIX]{x1D6F7})$ ) of flow through the centre pore for $Re_{H}=1000$ (mean flow is in the positive $x$ -direction).

3.4.1 Vortical structure visualization

To inspect how the vortical structures are influenced by inertia effects, the local swirling strength $\unicode[STIX]{x1D706}_{ci}$ is used to visualize the flows at different $Re_{H}$ as shown in figures 11(a), 11(b) and 11(c). The field variable $\unicode[STIX]{x1D706}_{ci}$ is determined by computing the imaginary component of the complex conjugate eigenvalue of the velocity gradient tensor $\unicode[STIX]{x1D608}_{ij}$ (Zhou et al. Reference Zhou, Adrian, Balachandar and Kendall1999), defined as

(3.12) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D608}_{ij}=\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{j}. & & \displaystyle\end{eqnarray}$$

The isosurfaces in these figures were set at 25 % of the local maximum value of $\unicode[STIX]{x1D706}_{ci}$ . With increasing $Re_{H}$ , the size of vortical structures becomes smaller. These structures seem to be stretched and aligned in the streamwise mean flow direction at lower $Re_{H}$ , however, and are broken down into much smaller structures at higher Reynolds numbers. For the cases studied in the present work, the size of these flow structures are still much smaller than the particle size. Thus, the turbulence is expected to be confined to the microscales of pores, as discussed in § 3.2. The size and connectivity of vortical structures located at distinct parts of unit cell are quite different for the different Reynolds numbers examined, suggesting a highly anisotropic flow field. To assess this further, the Reynolds stress anisotropy tensor is computed and its characteristics are analysed in the following section.

Figure 11.

The isosurface of the swirling strength at $\unicode[STIX]{x1D706}=0.25\unicode[STIX]{x1D706}_{max}$ . (a) is for $Re_{H}=300$ , (b) $Re_{H}=500$ and (c) $Re_{H}=1000$ .

3.4.2 Anisotropy characteristics

The anisotropy in the Reynolds stress distribution in the pore is important for understanding both the dynamics of the turbulent motions and the kinetic energy transport mechanism within the pore space. It is conventionally characterized by the normalized Reynolds stress anisotropy tensor, $\unicode[STIX]{x1D623}_{ij}$ , defined as

(3.13) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D623}_{ij}=\frac{\unicode[STIX]{x1D622}_{ij}}{2k}=\frac{\overline{u_{i}^{\prime }u_{j}^{\prime }}}{\overline{u_{k}^{\prime }u_{k}^{\prime }}}-\frac{1}{3}\unicode[STIX]{x1D6FF}_{ij}. & & \displaystyle\end{eqnarray}$$

Here, $\unicode[STIX]{x1D6FF}_{ij}$ is the Kronecker delta tensor and $k$ is the turbulent kinetic energy. The invariants of the tensor $\unicode[STIX]{x1D623}_{ij}$ are related to its three eigenvalues, $\unicode[STIX]{x1D706}_{i}$ . In particular, the second ( $II=3\unicode[STIX]{x1D702}^{2}$ ) and third ( $III=2\unicode[STIX]{x1D709}^{3}$ ) invariants associated with $\unicode[STIX]{x1D623}_{ij}$ can be computed from

(3.14) $$\begin{eqnarray}\displaystyle \displaystyle II=3\unicode[STIX]{x1D702}^{2} & = & \displaystyle b_{ji}\unicode[STIX]{x1D623}_{ij}/2\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D706}_{1}^{2}+\unicode[STIX]{x1D706}_{1}\unicode[STIX]{x1D706}_{2}+\unicode[STIX]{x1D706}_{2}^{2}\end{eqnarray}$$

(3.15) $$\begin{eqnarray}\displaystyle \displaystyle III=2\unicode[STIX]{x1D709}^{3} & = & \displaystyle \unicode[STIX]{x1D623}_{ij}\unicode[STIX]{x1D623}_{jn}\unicode[STIX]{x1D623}_{ni}/3\nonumber\\ \displaystyle & = & \displaystyle -\unicode[STIX]{x1D706}_{1}\unicode[STIX]{x1D706}_{2}(\unicode[STIX]{x1D706}_{1}+\unicode[STIX]{x1D706}_{2}),\end{eqnarray}$$

where $\unicode[STIX]{x1D702}$ characterizes the degree of anisotropy and $\unicode[STIX]{x1D709}$ describes its characteristic shape. These invariants are routinely used to compose the coordinate system to illustrate the anisotropy, for example, Lumley triangle (Lumley & Newman Reference Lumley and Newman1977) and turbulence triangle (Choi & Lumley Reference Choi and Lumley2001). These are both nonlinear invariant maps, but the first one is constructed based on the invariants $II$ and $III$ of $\unicode[STIX]{x1D623}_{ij}$ , and the second on the same invariants but in the transformed coordinate system $(\unicode[STIX]{x1D709},\unicode[STIX]{x1D702})$ . These maps are very useful, especially for the geometrically simple scenarios such as the turbulent channel flow. However, when dealing with a complex flow field, since the physical coordinate information is lost, they fail to identify the local region to which the anisotropic characteristics belong. In addition, they are appropriate only for visualizing small amounts of data (Emory & Iaccarino Reference Emory and Iaccarino2014).

For the current study, the barycentric componentality contour map developed by Emory & Iaccarino (Reference Emory and Iaccarino2014) is applied to diminish the drawbacks of the nonlinear invariant maps. Componentality maps are combinations of the eigenvalues that help characterize the topology of turbulence. In short, the limiting values of componentality can be thought of as being 1-component (roughly like a needle), 2-component (flat disk-like), or 3-component (roughly spherical). Using componentality, the anisotropy can be graphically visualized. This approach is briefly described as follows.

The barycentric map, first introduced by Banerjee et al. (Reference Banerjee, Krahl, Durst and Zenger2007), utilizes the fact that all the turbulence states with physical meanings should be bounded by the three limiting states of componentality (described above) in a convex combination. Thus, a Euclidian domain in an equilateral triangle, where the three limiting states with equal weights are placed at $\boldsymbol{x}_{1c}=(1,0)$ , $\boldsymbol{x}_{2c}=(0,0)$ and $\boldsymbol{x}_{3c}=({\textstyle \frac{1}{2}},\sqrt{3}/2)$ (see figure 12 a) can be established. The corresponding coordinate system $(x_{B},y_{B})$ in the map can be written as

(3.16) $$\begin{eqnarray}\displaystyle & \displaystyle x_{B}=C_{1c}x_{1c}+C_{2c}x_{2c}+C_{3c}x_{3c}, & \displaystyle\end{eqnarray}$$

(3.17) $$\begin{eqnarray}\displaystyle & \displaystyle y_{B}=C_{1c}y_{1c}+C_{2c}y_{2c}+C_{3c}y_{3c}, & \displaystyle\end{eqnarray}$$

where the weights $C_{ic}$ are defined as

(3.18) $$\begin{eqnarray}\displaystyle & \displaystyle C_{1c}=\unicode[STIX]{x1D706}_{1}-\unicode[STIX]{x1D706}_{2}, & \displaystyle\end{eqnarray}$$

(3.19) $$\begin{eqnarray}\displaystyle & \displaystyle C_{2c}=2(\unicode[STIX]{x1D706}_{2}-\unicode[STIX]{x1D706}_{3}), & \displaystyle\end{eqnarray}$$

(3.20) $$\begin{eqnarray}\displaystyle & \displaystyle C_{3c}=3\unicode[STIX]{x1D706}_{3}+1. & \displaystyle\end{eqnarray}$$

Note that any $C_{ic}$ resides in the range $[0,1]$ . This setting is coincident with the red-green-blue (RGB) colour system, where the $[1,0,0]$ is coded for red, $[0,1,0]$ green, and $[0,0,1]$ blue. Every colour could be expressed by a linear combination of these three triplets. As a result, the barycentric componentality coefficients $C_{ic}$ could as well be used as the coefficients for the RGB colour function in the form of

(3.21) $$\begin{eqnarray}\displaystyle \left[\begin{array}{@{}c@{}}R\\ G\\ B\end{array}\right]=C_{1c}\left[\begin{array}{@{}c@{}}1\\ 0\\ 0\end{array}\right]+C_{2c}\left[\begin{array}{@{}c@{}}0\\ 1\\ 0\end{array}\right]+C_{3c}\left[\begin{array}{@{}c@{}}0\\ 0\\ 1\end{array}\right]\!. & & \displaystyle\end{eqnarray}$$

Therefore, the different shapes of turbulence can be mapped onto the equilateral triangle (see figure 12 a) with the help of the RGB colour system. In particular, one-dimensional turbulence is identified by red, two-component turbulence by green, and the isotropic turbulence by the colour blue.

Figure 12.

(a) Reynolds stress anisotropy componentality contour in the barycentric colour map. (b) Lumley triangle showing limits of the invariants $\unicode[STIX]{x1D709}$ and $\unicode[STIX]{x1D702}$ and characteristic spheroids representing the shape of the turbulence.

Figure 13.

The column on the left-hand side shows the contours for anisotropy at the centre slice using barycentric componentality map; and the right-hand side illustrates the anisotropy by the characteristic ellipsoids with a zoom-in view at the centre pore (the dashed square in (a)): (a,b) $Re_{H}=300$ , (c,d) $Re_{H}=500$ , and third (e,f) $Re_{H}=1000$ .

Distribution of the anisotropy stress tensor on the centre $xy$ -slice of the porous unit cell described by the barycentric componentality contour map is shown in figures 13(a), 13(c) and 13(e) for the three Reynolds numbers, respectively. In addition, a new technique, first introduced by Hamilton & Cal (Reference Hamilton and Cal2015), is applied as well to visualize the turbulence structures more intuitively. The algorithm involves using the eigenvalue decomposition of the normalized Reynolds stress anisotropy tensor to extract the corresponding eigenvalues and eigenvectors. The characteristic triaxial ellipsoids (figure 12 b), with the shape determined by the eigenvalues and orientation by the eigenvectors, can then be drawn to directly represent the turbulence structure as described by Lumley triangle. Figures 13(b), 13(d) and 13(f) provide zoom-in views on the centre pore region (shown as the dashed square in figure 13 a) depicting the turbulence with such ellipsoids. Note, in order to avoid the clustering, only one in six of the data points in each direction are displayed.

Figures 13(a), 13(c) and 13(e) show that some characteristics of the anisotropy stress tensor are common to all three Reynolds-number flows. The characteristic ellipsoids in figures 13(b), 13(d) and 13(f) help to visualize the types of structures present in different regions. Firstly, in the regions close to the contact points of the spheres on the left-hand side, the near one-dimensional turbulence (red-ish) dominates as the flow is squeezed and accelerated by the particles to form a jet-like structure. At higher Reynolds numbers, the extent of this one-dimensional turbulence is reduced. As the accelerated flow extends into the pore space, transition to 2-D turbulence (green) is observed. A weak wake-like region is present behind the left sphere. In the regions of peak TKE production (yellow-dashed circles in figures 13 a, 13 c and 13 e), the boundary between red (one-dimensional (1-D) turbulence) and green (2-D turbulence) becomes less discernible with increasing $Re$ . The dominant colour transits from light red-ish (figure 13 a, $Re_{H}=300$ ), to a mixture of green and red (figure 13 c, $Re_{H}=500$ ) and finally turns to purple-ish (figure 13 e, $Re_{H}=1000$ ). Moreover, the ellipsoids in the red-dashed circles in figures 13(b), 13(d) and 13(f) indicate that the shape of the turbulence becomes more like a 2-D disk when the Reynolds number increases from 300 to 500; and then develops into the cigar shape at $Re_{H}=1000$ . The flow squeezed through the boundaries of the particles is stretched mainly by the Reynolds stress in the main flow direction giving rise to cigar-shaped structures, whereas near the centre of the pore and on the left side, the structures are more stretched by other components of the Reynolds stress leading to the 2-D turbulence. The anisotropy distribution depends on the distance from the wall, as expected. Similar to duct flows, near-wall shear leads to higher anisotropy and elongated structures. However, in the present configuration, this phenomenon is not as pronounced owing to the three-dimensional geometry of the pore and the tortuous path of the flow through the FCC packing.

As the flow from the two corners of the left sphere meets near the centre of the pore, one transitions from 2-D (green-ish) to isotropic turbulence (blue). The extent of this isotropic region increases with increasing Reynolds number. This is also evident by small spherical shapes in the centre pore region towards the right-hand side. As the flow then impinges on the right sphere boundary, a region of disk-like structure is observed representing the isotropic 2-D turbulence (green region) very close to the particle surface. The anisotropy characteristics close to the stagnation zone is consistent with the results reported by Nishino et al. (Reference Nishino, Samada, Kasuya and Torii1996). As the flow turns and exits through the right corners, it is again squeezed to give the 1-D (red region) turbulence.

3.4.3 $Q-R$ invariants space

The velocity gradient tensor has been thoroughly studied because it contains enriched geometric and statistical information of the flowfield, such as shapes of fluid volume, rate of deformation and intermittency (Meneveau Reference Meneveau2011). Particularly, much attention has been paid to the second and third invariants (Cantwell Reference Cantwell1992; Martin et al. Reference Martin, Ooi, Chong and Soria1998), defined based on this second order, traceless tensor as

(3.22) $$\begin{eqnarray}\displaystyle & \displaystyle Q=-{\textstyle \frac{1}{2}}\unicode[STIX]{x1D608}_{im}\unicode[STIX]{x1D608}_{mi}, & \displaystyle\end{eqnarray}$$

(3.23) $$\begin{eqnarray}\displaystyle & \displaystyle R=-{\textstyle \frac{1}{3}}\unicode[STIX]{x1D608}_{im}\unicode[STIX]{x1D608}_{mn}\unicode[STIX]{x1D608}_{ni}. & \displaystyle\end{eqnarray}$$

It is convenient to study the velocity gradient tensor through the investigation of two such invariants because, firstly, this methodology greatly reduces the dimension of the problem (from eight to two); and secondly, it can provide a representation of the local topology of fluid motion. One common observation among different types of turbulence is the sheared ‘tear drop’ shape distribution of the $(R,Q)$ -space (see Chong, Perry & Cantwell (Reference Chong, Perry and Cantwell1990), Cantwell (Reference Cantwell1993)). Also, there exists a distinguished line, the Vieillefosse tail defined as $Q=-(3/2^{2/3})R^{2/3}$ , along which the higher probability of data points are clustered especially at the $R>0$ and $Q<0$ quadrant. The Vieillefosse tail separates the $(R,Q)$ -space into two parts. The one above the curve represents the focal or vortical fluid motions and the other, below the curve, denotes the non-focal, dissipative motions. In this section, the $(R,Q)$ -space is examined for different regions in the pore based on the anisotropy distribution illustrated in § 3.4.2.

Figure 14.

Scatter plots of the $(R,Q)$ -space at $Re_{H}=1000$ corresponding to data taken from different regions within the pore indicated in figure 13(f): (a) isotropic region corresponding to the green triangle area, (b) 2-D turbulence within the blue rectangle region, and (c) cigar-shaped turbulence structures corresponding to peak TKE production region represented by the two red circles. The dashed curves are the Vieillefosse tail.

The $(R,Q)$ -space scatter plots are shown in figure 14 for the highest Reynolds number $Re_{H}=1000$ . Data points from three different regions in figure 13(f), corresponding to the isotropic (green-dashed triangle), 2-D turbulence (blue-dashed rectangle) and cigar-shaped turbulence structures (red-dashed circles), which is also a region of high production of turbulent kinetic energy in figure 4(e), are used to draw figures 14(a), 14(b) and 14(c), respectively. The sheared ‘tear drop’ shape distribution of the $(R,Q)$ -space is generally recognized for all three regions and the Vieillefosse tail is followed in the $R>0$ and $Q<0$ quadrant. This is consistent with extensive observations of different types of turbulence observed from both numerical (Nomura & Post Reference Nomura and Post1998; Ooi et al. Reference Ooi, Martin, Soria and Chong1999; Lüthi, Tsinober & Kinzelbach Reference Lüthi, Tsinober and Kinzelbach2005) and experimental (Gultitski et al. Reference Gultitski, Kholmyansky, Kinzelbach, Luthi, Tsinober and Yorish2007; Elsinga & Marusic Reference Elsinga and Marusic2010) work. However, there are discrepancies in the scatter plots for the three regions depicted. The most noticeable difference is that considerably more data points above the Vieillefosse tail at the $R<0$ and $Q>0$ quadrant are observed for 2-D and cigar-shape turbulence, meaning that more vortical motions are present in these two regions. The fluid motions represented by the cigar-shape turbulence are mostly helical motions, which aggressively promote the production of TKE (see figure 4). Whereas for the isotropic turbulence region, the data points in the $(R,Q)$ -space are mostly clustered near the origin and follow the Vieillefosse tail in the $R>0$ and $Q<0$ quadrant, as expected. On the other hand, only small differences in data under the Vieillefosse tail are observed. As discussed in § 3.1.3, the TKE production peaks in particular regions within the pore; however, the dissipation rate of TKE is mostly non-zero over the entire pore space. This explains why there are little deviations of data points in the $(R,Q)$ -space below the Vieillefosse tail in different regions.