Consistent Flow Structure Evolution in Accelerating Flow Through Hexagonal Sphere Pack

Direct numerical simulation based on incompressible Navier–Stokes equations with an immersed boundary method is used to simulate accelerating porous media flow through a bed of uniform spheres arranged in hexagonal close packing order. The transient flow is realised by driving initially resting fluid by a constant pressure gradient. A wide spec-trum of Reynolds number based on the sphere diameter and volume-averaged velocity is considered, which ranges from creeping flow up to a Reynolds number of approximately 350, where turbulent flow structures are evident inside the pores. It is found that nonlinear dependence of the volume-averaged velocity with respect to the applied pressure gradient is the consequence of emergence of streamwise jets and the accompanying streamwise vortices, as previously observed for other sphere pack arrangements. Furthermore, two distinct flow modes are identified in the steady flow regime which satisfy full geometric symmetries. The flow then becomes unsteady around Reynolds number of 90 which coincides with a partial breaking of the symmetries, and pore-scale turbulence emerges once all the symmetries vanish when Reynolds number is larger than 200. For all the considered unsteady flow, independent of being turbulent or not, we observe a consistent sequence of flow structure evolution during the flow development with progressively broken symmetries albeit at widely varying instantaneous Reynolds numbers. Moreover, we show that the symmetry breaking takes place in larger pore spaces first, then propagate into smaller pores located in downstream.


Introduction
There are many real-world porous media flow problems that are transient rather than steady, for instance, flow through our respiratory system, a forest inside a developing atmospheric boundary layer, and wave-induced current inside a coral community (Lowe et al. 2008;Rogers et al. 2016). However the majority of the studies in the literature focus 1 3 on steady-state flow, with few exceptions such as Hill et al. (2001a) who numerically studied transient ordered/random sphere-pack flow being accelerated from rest, however only in the zero Reynolds number limit. Whilst the typical Reynolds number range of such flow is indeed low (i.e. viscosity-dominated), in some technical and environmental applications the governing non-dimensional number can be significantly higher than unity, meaning that the inertial influence is significant and even pore-scale turbulence may appear. For the steady-state flow, it is known that different flow regimes emerge depending on the Reynolds number, namely: steady linear, steady nonlinear, periodic/quasi-periodic nonlinear, chaotic or turbulent (Dybbs and Edwards 1984;Hill and Koch 2002). To our knowledge, how those flow regimes develop in transient (development) phase is yet to be investigated. Moreover, the precise Reynolds number limits to distinguish those flow regimes vary greatly in the literature-possibly due to the diverse porous media geometries that were considered-and they remain topics of active research. To contribute to this aspect, we consider a porous medium consisting of uniform touching spheres packed in the well-defined arrangement called hexagonal closed packing (HCP). Our motivation to choose such regular yet more complex sphere-pack geometry than some of the other popular geometries [e.g. simple cubic or face-centred-cubic (FCC)] is to address our knowledge gap between those simpler geometries and the random pack of spheres which is also frequently studied (e.g. Hill et al. 2001a;Fand et al. 1987;Ghaddar 1995;Koch and Ladd 1997;Hill et al. 2001b).
Due to the high geometrical and dynamical complexities, the porous media flow is often treated macroscopically by means of the volume-averaged quantities and models. Though there exists well-established steady-state models, the time-dependent counterparts still require significant improvement. Currently, many of the unsteady models are constructed by extending the steady-state models. However, during the development phase towards equilibrium, the porous media flow would undergo a series of qualitative transitions, which may require new modelling approaches that are significantly different from the steady-state models. Consequently one may ask: Can we simply relate the intermediate states during development at certain transient Reynolds number to the fully-developed picture of the equivalent Reynolds number? Alternatively, the transient dynamics might be governed by a significantly different set of non-dimensional parameters. To address this question directly, we investigate the possible link between the qualitative change in the bulk quantities and the corresponding pore-scale flow structures. Potentially, this microscopic approach allows us to evaluate the generality of our findings, as it can be used as an additional phenomenological indicator to associate, or differentiate, the flow through different porous media geometries.
Independent of whether we study the time-dependent flow through porous media numerically or experimentally, the key technical challenge arises from the extreme spatiotemporal resolution requirements. Experimentally capturing the fully three-dimensional flow structures, whose length-scales are typically in the order of the measurement probe size, is extremely challenging. The limited optical access imposed by the complex geometry also hinders the possibility to perform any quantitative image-based measurements. Similarly, the high computational effort required to fully resolve the fine-scale dynamics inside the complex pores had been equally challenging on the simulation front. However, the recent development in advanced numerical techniques, such as immersed boundary methods combined with the conventional Navier-Stokes solvers (Chan-Braun et al. 2011), or lattice-Boltzmann method (Hill et al. 2001b;Jin et al. 2015;Suga 2016), opened up new possibilities to simulate such complex flow phenomena with a realistic effort.
Therefore we identified our research questions to be addressed in this paper as following. At which Reynolds number nonlinearity (i.e. significant inertial footprint) emerges for the hexagonal sphere pack flow? How does the nonlinearity manifest inside the pore space as structural features? Similarly, at which Reynolds number does the chaotic or turbulent flow regime appear, and how? Are the terms "chaotic" and "turbulent" flow equivalent in the flow configuration? If not, how do they differ? How those structural features in different flow regime emerge in the developing flow? Do each of flow regime follows own unique development "path", or are those paths somehow related?
Consequently, we define the following objectives for this paper. We perform a series of fully-resolved direct numerical simulations of flow through hexagonal sphere pack at the close limit, covering different (steady-state) flow regimes by varying Reynolds number. The step between the Reynolds numbers need to be fine enough to adequately locate the limits of those flow regimes.
Due to the high computational effort required to generate such high-fidelity dataset, we limit ourselves to consider a single pressure gradient direction.We then investigate the pore-scale flow structure development in the transient flow in different flow regimes by identifying the key structural features for each flow state. These identified flow features are then compared with the equivalent structures in other more well-understood pore geometry in their steady-state.
In the subsequent section, we will detail the numerical methods and the simulation setup, then we will discuss about our simulation results in Sect. 3 including Reynolds number dependence of steady-state superficial velocity, and pore-scale flow structure analysis of linear, steady/unsteady nonlinear and turbulent flow regime. We will then conclude this paper with summary and outlook in Sect. 4.

Numerical Methods and Simulation Set-Up
We employ our in-house fluid solver MGLET which solves the incompressible Navier-Stokes equations for the primitive variables (i.e. velocity and pressure) in the Cartesian coordinate system: where = [u, v, w] is velocity vector corresponding to the spatial coordinate = [x, y, z] , whilst p is pressure. Moreover, and are fluid density and the kinematic viscosity respectively. Those variables are stored in a staggered grid arrangement, and discretised in space using second-order central finite-volume method. An explicit three-step low-storage Runge-Kutta time integration scheme is used (Williamson 1980), which is combined with Chorin's projection method to decouple velocity and pressure computations (Chorin 1968). Consequently, a pressure Poisson equation is solved for each Runge-Kutta substep, which is done by the Strongly Implicit Procedure (SIP) iterative solver (Stone 1968).
To resolve the complex geometry of the porous media, a second-order spatially-accurate mass-conserving immersed boundary method is employed (Peller et al. 2006;Peller 2010).
For completeness, an outline of the employed immersed boundary method is given as follows. For setting the internal boundaries, grid cells that are intersected by the geometries are omitted from the computation. Instead, the velocities which belong to the interface cells are determined in two different ways: first variant is done by point-wise interpolation based Dirichlet boundary conditions and the neighbouring fluid-phase cell values; whereas second variant is also based on the boundary conditions and the adjacent fluid  cell values, but by computing the mass flux through the open face areas instead. The physical  dimension of the second variant is made consistent with the first variant by dividing them by  the open face areas, and their mass conservation property is maintained by distributing the  divergence to the neighbouring cells based on their open areas (Peller 2010). In the framework of the projection method, the first non-divergence-free variant is used to compute the advective and diffusive terms in order to determine the intermediate velocity field. Then, the interface cell values are overwritten by the divergence-free flux variant in prior to the pressure correction step. The code is MPI-parallel. The parallelisation is done based on a conventional domain decomposition, which can be combined with a local grid refinement strategy by adding grid boxes with finer resolutions in an octree-like, hierarchical and overlapping manner (Manhart 2004). The parallel implementation of the code has recently undergone intensive optimisation for the modern massively-parallel architectures, and the applied methods and improvements have been documented in Sakai et al. (2019).
In our simulations, the porous media flow domain is represented by a triply-periodic 3 ] , which is filled by spheres with uniform diameter d. The spheres are arranged in the HCP arrangement: where i, j, k = 0 … n are the sphere indices in x, y, z-direction respectively. The resulting sphere pack has the following two global geometric symmetries: reflectional symmetry along z-axis, and ∕3-rotational symmetry on xy-plane. Moreover, the sphere pack forms two distinct types of pore geometries, namely tetrahedral and octahedral pores, which are depicted in Fig. 1. On a yz-plane, where we will extensively place our attention on in the following analysis part of this paper, there exist two local perpendicular reflectional symmetry lines for each pore type (cf. Fig. 1b) in addition to the global z-axis symmetry. Those two distinct types of pores are interconnected and appear in an alternating order in the x-direction, following a split-and-merging pattern as: one octahedral pore splitting into a pair of tetrahedral pores, and merging into the subsequent octahedral pore (cf. Fig. 1c and d).
Along such pore network, those two pore geometries are separated by 0.5d in x-direction. The porosity of the sphere pack is at the close-pack limit of ≈ 0.26 (or the corresponding solid volume fraction = 1 − ≈ 0.74 ), which means all spheres are in contact to their neighbours.
To simulate the transient behaviour of our porous media flow, resting flow is accelerated by a constant pressure gradient in x-direction, namely ∇ x p = p x = C , and simulated until the amplitude of the following superficial (volume-averaged) streamwise velocity reaches to the steady-state level visually: where V is the total volume of numerical domain including the solid phase. The velocity inside the solid phase is fixed at zero, therefore the intrinsic (i.e. fluid volume averaged) (3) ⟨u⟩ s = 1 V ∭ u(x, y, z) dx dy dz . Fig. 1 Visualisations of simulation domain: a 2D projection view on a xy-plane where the pores coloured in purple and red indicate tetrahedral and octahedral pores respectively; b 3D parallel projection view on a yzplane; c streamline visualisation of the pore network; d streamwise velocity contour at u = 2 together with the same quantity plotted on a yz-plane at x∕d = 1.5 where the opacity is scaled with the magnitude velocity ⟨u⟩ i can be computed from ⟨u⟩ s , as ⟨u⟩ i = ⟨u⟩ s ∕ .For the cases with unsteady oscillations, we simulate longer than the above criterion so that we can determine the long-term behaviour. Throughout this paper, we use the Reynolds number based on the superficial velocity ⟨u⟩ s and sphere diameter d, as: where is the dynamic viscosity. Additionally, we refer to the Reynolds number based on the steady-state ⟨u⟩ s as Re steady , where ⟨u⟩ s is averaged over the statistical steady-state: To control the Reynolds number, we fix ∇ x p , d and , whereas is varied to adjust kinematic viscosity = ∕ .
Overall 16 simulations at their final grid resolutions are performed for this study, and their numerical and physical parameters are summarised in Table 1. All of the production simulations were performed on CoolMUC-2 Linux Cluster hosted by Leibniz Supercomputing Center of the Bavarian Academy of Sciences and Humanities (LRZ). The smaller cases with the d∕Δx = 160 spatial resolution (i.e. L1-6, SNL1-3) required 1 Intel Haswell node with 28 cores/node, whereas the larger cases with d∕Δx = 320 (SNL4, UNL1-2, T1-4) were performed on 32 nodes (896 MPI processes). Moreover, one additional simulation with d∕Δx = 640 was performed to evaluate the grid convergence. Due to its significantly larger problem size, this particular simulation was run on SuperMUC Phase 2 of LRZ, employing 256 Intel Haswell nodes and 7168 MPI processes. The total walltime of each simulation ranges between 1 and 9 days. Darcy permeability K of the sphere pack was calculated a-posteriori using the result of the L1 case and the following 1D version of Darcy's linear empirical relation: The resulting normalised permeability K∕d 2 = 1.731 × 10 −4 compares reasonably to the one from the well-known Kozeny-Carman equation (Kozeny 1927;Carman 1937Carman , 1956: and the FCC sphere pack result of Chapman and Higdon (1992) We employ an equidistant grid point distribution in each direction, and the grid spacing in all directions are approximately identical, i.e. Δx ≈ Δy ≈ Δz . Grid resolution of 160 finite-volume cells per diameter up to Re steady = 48 , and 320 points per diameter between 59 ≤ Re steady ≤ 347 are used. This particularly high grid resolution is required for our sphere pack simulations partially because all the spheres are in contact with their neighbours at this close-pack limit. Sufficiently resolving the flow around the contact points is one of the most computationally expensive aspects of our simulations. A visual representation of the d∕Δx = 320 grid shown in Fig. 2 highlights the required high resolution around the critical regions.
Sufficiency of the grid resolution was validated in terms of the convergence of the superficial streamwise velocity at Re = 59 (cf. Fig. 3a) and Re = 347 (cf. Fig. 3b). The grid resolutions of d∕Δx = [80, 160, 320] were considered for the former, whilst d∕Δx = [160, 320, 640] were evaluated for the higher Reynolds number case. In case of Re = 59 , the discrepancy between the finest and the subsequent resolution in terms of the peak ⟨u⟩ s was found at 0.85% , whereas 0.8% discrepancy with respect to the steadystate value was observed. From this, we concluded that d∕Δx = 160 is sufficient. On the

Fig. 2
Grid visualisation on yz-plane at x∕d = 0 with d∕Δx = 320 resolution. a Overview of tetrahedral pore region; b zoom-up of the contact point in a other hand, we concluded that d∕Δx = 320 is sufficient for the Re = 347 case since the peak value discrepancy of 0.06% (Fig. 3c) and the steady-state discrepancy of 2.5% were observed between d∕Δx = 320 and d∕Δx = 640 grid resolutions. Note that the relatively large discrepancy for the steady-state value can be attributed to the short averaging period for the highest resolution case due to the high computational effort. Also note that Fig. 3c exhibits convergences that are higher than second-order, which we attribute to gradual opening of the sphere contacting regions with increasing grid resolution. Those regions are otherwise impermeable, since the current immersed boundary method blocks an entire cell when multiple boundaries intersect within the cell. We decided to use d∕Δx = 320 even from Re steady = 59 to be on the safe side for the reason which will be introduced shortly, and the same philosophy applies to the maximum CFL number which was kept under 0.16.
Note that the above grid resolutions are significantly finer than the ones used in the previous numerical simulations of sphere-pack flow in the literature. For instance, Hill and Koch (2002) Chu et al. (2018) performed direct numerical simulations of flow through a staggered array of square cylinders up to Re = 1500 using 94 cells per unit square length. More recently, He et al. (2019) simulated the flow through a close-pack of spheres in FCC arrangement up to Re = 740 based on the superficial velocity using 250 points per sphere diameter. Note that all those studies focused on the steady-state dynamics when the pore-scale flow structures are already fully-developed, whereas our developing porous media flow forms infinitely thin boundary layers next to the sphere boundaries when the flow is suddenly accelerated from rest. To resolve the thin developing boundary layers at adequate accuracy, it is essential to use such exceptionally fine grid resolutions.

Reynolds Number Dependence of the Steady-State Superficial Streamwise Velocity
Although our main focus in this paper is the transient part for the porous media flow, first we discuss about the Reynolds number dependency of the steady-state superficial streamwise velocity. As it will become apparent, such information can be related to the evolution of the pore-scale flow structures that will be discussed in the following section. Steadystate is visibly achieved for all the simulations, however, we observe that the flow with Re steady ≥ 91 becomes unsteady, therefore the values from the cases with above-threshold Reynolds number need to be treated with care. To evaluate the above Re-dependence, we plot the normalised Darcy permeability using the simulated ⟨u⟩ s and other parameters that are included in Eq. (6) as a function of Re steady in Fig. 4. Figure 4a illustrates that the Darcy permeability is independent of Reynolds number between 5 × 10 −7 ≤ Re steady ≤ 1 , therefore the simulation results in this Reynolds number range are in accordance with Darcy's linear relation between the imposed pressure gradient and ⟨u⟩ s . Some previous studies suggest the existence of a so-called pre-Darcy flow regime indicated by localised increase of flow resistance (i.e. lower permeability) in comparison to the following linear Darcy flow regime (cf. Fand et al. 1987;Kececioglu and Jiang 1994;Baǧci et al. 2014;Kundu et al. 2016  supporting data, that the lower limit of Reynolds number where the Darcy's law holds is around Re ≈ 10 −5 (Fand et al. 1987). Our simulation results, however, do not support the existence of such flow regime which is associated with extremely low Reynolds number at least for this flow configuration, which agrees well with the experimental study of Farkas et al. that such flow regime cannot be observed (Farkas et al. 1999).
For Re steady > 1 , the non-dimensional pressure gradient appears to be linearly proportional to Reynolds number, which supports the quadratic Forchheimer correction (Forchheimer 1901), namely: where a can be set as − K as in Eq. (6), whereas b is referred to as Forchheimer coefficient. Notice that the gradient of the linear correlations slightly but noticeably differ in two different Reynolds number ranges, and the transition occurs somewhere between 91 ≤ Re steady ≤ 138 . The above observation indicates the existence of the Reynolds number dependence for the Forchheimer coefficient. Particularly, the range 10 ≤ Re steady ≤ 59 yields a = −1.679 × 10 4 and b = 1.028 × 10 4 , whereas the subsequent 91 ≤ Re steady ≤ 347 cases result a = −4.8435 × 10 4 and b = 7.9804 × 10 4 .
In the form of non-dimensional drag force employed in Hill et al. (2001a), Hill and Koch (2002): the current results yield: This observation agrees with the experimental finding of Fand et al. that the gradient of the linear Re-dependence of permeability changes somewhere between Re = 80 and 120, and they relate it to the onset of unsteady oscillations (Fand et al. 1987). As you will see in the following discussion, indeed the first onset of unsteady oscillation was observed for our case at Re steady = 91 , which coincides with their observation. Although for different sphere pack geometry, Hill and Koch (2002) also reported such transition around Re = 160.

Steady-State Pore-Scale Flow Structures
In the following, we investigate the flow structures which appear inside the pores for various Reynolds numbers at which the flow develops to a steady-state. In other words, the considered Reynolds number range covers from the very low (Darcy) limit up to the steady nonlinear regime, which we identify to be up to Re steady = 59 . We select four representative Reynolds numbers, namely Re steady = 10 −5 and 10 −2 (linear), 10 and 59 (steady nonlinear) and discuss their flow structures in detail. Due to the complex nature of the pore geometry in the hexagonal sphere pack, we study the two-dimensional cross-sectional view of the flow structures first. Then, when it is appropriate, we will supplement those 2D findings with the corresponding three-dimensional view.
At Re ≪ 1 , the flow behaves according to Stokes equation where viscous forces dominate over inertial forces. Based on the Re-dependence of the streamwise superficial velocity shown in Fig. 4, we can expect that at Reynolds numbers larger than approximately ten nonlinearity emerges and the steady-state flow structure undergoes a series of transitions, which is in a good agreement with the previous observations as following.
For instance, using the measurements of flow inside complex rod porous media structure, Dybbs and Edwards (1984) distinguished between the linear and the steady nonlinear flow regime by development of boundary layers from the porous media surfaces and the subsequent appearance of "inertial cores", or streamwise jets existing outside of the boundary layers. The transition to the steady nonlinear flow regime takes place somewhere between Re = 1 and 10, and persists until Re = 150 where their flow starts to oscillate.
Hill et al. also observed such inertial core in their simple cubic and FCC sphere pack flow. Moreover, they showed that the critical Reynolds number for the onset of unsteady oscillations is sensitive to the direction in which the pressure gradient is applied. They conjectured the directional dependence is related to whether the pressure gradient axis is in alignment to the geometric symmetries of their pore geometry (Hill et al. 2001a). In the subsequent work, Hill and Koch showed the limit of steady flow is around Re ≈ 60 for their FCC sphere pack when the pressure gradient is applied in the symmetry line. Moreover, they showed that, on the verge of onset of unsteady oscillations, there exists a steady-state eight-vortex mode with ∕2 rotational symmetry, which consists of four counter-rotating vortex pairs (cf. their Fig. 9) (Hill and Koch 2002). Figure 5 shows the steady-state flow fields of the four representative Reynolds numbers in terms of streamwise velocity and vorticity on a cross-sectional plane at x∕d = 0 . The cross-sectional plane features two larger quadrilateral pores and four smaller triangular pores, which are the intersections of the cross-sectional plane with the octahedral and tetrahedral pores in 3D, which were introduced in Sect. 2. To avoid confusion, we will continue referring those two-dimensional pores by the corresponding three-dimensional terminologies. As we detailed in Sect. 2, those cross-sectional pores satisfy the global reflectional symmetry about the horizontal axis along z∕d = √ 6∕3 ≈ 0.816 , as well as two additional reflectional symmetry lines local to the pore cross-sections that are perpendicular to each other (cf. also Fig. 1b). Lastly, the cross-sectional plane passes through only the spheres' centroid and the contact points between two spheres, therefore there is no inclination in x-direction. Consequently in the Stokes flow limit, there is no cross-flow velocity components imposed by curvature of the geometry.
The 2D flow structures of the L3 ( Re steady = 10 −5 ) and the L4 ( Re steady = 10 −2 ) cases are shown in Fig. 5a and b, and it can be seen that the flow feature is dominated by individual high-velocity jets that are filling the both pores rather uniformly. Although the boundary layers are already visibly developed at these Reynolds numbers, it is reasonable to argue that the flat velocity distribution is in accordance with the conjecture of Dybbs and Edwards (1984), that the inertial cores do not exist in this low-Re flow regime. In Fig. 5b, those high-velocity jets are surrounded by streamwise vortices which are first observed at this Reynolds number, whereas at lower Reynolds numbers only the streamwise jets are detected (cf. Fig. 5a). The octahedral-pore jets are surrounded by a four-vortex structure in weak intensity, whilst the tetrahedral-pore jets are accompanied by pairs of counter-rotating vortices with significantly higher intensity. Notice that thin vorticity layers exist between those vortex pairs and the sphere surface. The existence of the finite streamwise vorticity indicates the development of cross-streamwise velocity components already at this Reynolds number. Finally, it is important to note that all the flow features satisfy both global and local symmetries that were discussed earlier.
By increasing the Reynolds number to Re steady = 10 , the four-vortex mode in the octahedral pores has intensified and additional thin near-wall vorticity layers have established also inside the octahedral pores. Consequently, the profile of the octahedral jets becomes noticeably rounder, indicating that the stream can no longer follow the curvature of the pore geometries due to the elevated inertial influence. Recall that for Re steady ≥ 10 , the intensity of streamwise superficial velocity starts to deviate from the linear Darcy's law, indicating an emergence of nonlinearity. Figure 6 illustrates the corresponding three-dimensional view of the flow structure topology, where iso-surfaces of arbitrary positive streamwise velocity, positive and negative streamwise vorticity are depicted. The flow direction is from the bottom-left corner to the top-right corner, and only a part of the numerical domain is shown for better visibility. The entrance to this sub-domain is an octahedral pore, which is followed by a pair of tetrahedral pores in the narrow space between the touching spheres. From this view, it is visible that tetrahedral-pore shear layers-which exist as early as Re steady = 10 −2 -evolve into counter-rotating four vortices in the subsequent octahedral pores. This observation on the intense vortex formation coinciding with the emergence of nonlinearity is in a reasonable agreement with the explanation of Nield and Bejan (2006) that the nonlinear transition is associated with the occurrence of first (intense) eddies. The intensity of the vorticity in the tetrahedral pores also increases, however its structural features are maintained. Finally, the symmetry property remains unchanged.
By further increasing Reynolds number to Re steady > 10 , the original four octahedralpore vortices move apart from each other, consequently forming four additional vortices around the pore centre, resulting a new eight-vortex system (cf. Fig. 5c for the example of the SNL4 case with Re steady = 59 ). The first appearance of the eight-vortex mode on this particular plane is observed at Re steady = 36 , however, the precise limit between the fourand eight-vortex modes is still to be determined. Emergence of this eight-vortex system is accompanied by the drastically transformed streamwise jets, which are now split into a pair of smaller shell-like structures residing between the original outer and the newlyformed internal vortices. The separated jets can be described in terms of the alternating arrangement of the two pore geometries detailed in Sect. 2, where the flow in an octahedral pore is split into two in the following tetrahedral pore section before it merges again in the subsequent octahedral pore. When the Reynolds number is sufficiently low, the dominant viscous effect is strong enough to unify the two jets from a tetrahedral pore pair inside the octahedral pore downstream via diffusion. On the contrary, with increasing Reynolds number-therefore less viscous influence available-the diffusion process is no longer sufficient to unify the split jets, therefore they remain being separated throughout. Consequently, low-velocity cores are formed inside the shell-like streamwise jets. The corresponding cross-sectional pressure distribution shows us the low-velocity cores are coupled with higher pressure regions, which may trigger cross-flow from the cores towards the sphere surfaces. Meanwhile, the vortical structure inside the triangular pores evolves into a four-vortex mode. The corresponding three-dimensional view is shown in Fig. 7, and reveals that the appearance of the eight-vortex mode is the consequence of wall shear layers in octahedral pores evolving into additional smaller vortices, extending until the subsequent octahedral pores (cf. Fig. 7b). In fact, the emergence of the four-vortex system in tetrahedral pores is also a footprint of these additional vortices, whose filament-like tails extend until the following tetrahedral pores. These intensified octahedral-pore shear-layers eventually forming the eight-vortex can be described as a consequence of the streamwise jets being closer to the sphere walls due to the aforementioned splitting.
Lastly, first appearance of reversed flow inside the octahedral pores can be observed in this flow mode (cf. Fig. 7a). Through a closer inspection, it becomes visible that there are two regions where such reverse flow can be observed: One in the close vicinity of the sphere surfaces on the downstream side due to flow separation, and the other in between the split jets on the upstream side. Once again, all those transformed flow structures still follow the aforementioned global and local reflectional symmetries.

Transition to Unsteady Flow and Local Symmetry Breaking
Let us now proceed into the unsteady nonlinear flow regime, where the applied pressure gradient and the superficial velocity relation continues to be nonlinear whilst the flow becomes unsteady in addition.
Early experimental study of Jolls and Hanratty (1966) detected the onset of velocity fluctuations in flow through a random pack of spheres at around Re = 110 , using a masstransfer measurement method combined with flow visualisation via dye injection. Similarly with their measurement of flow through complex rod structure, Dybbs and Edwards (1984)  In contrast, Hill and Koch (2002) observed first appearance of periodic oscillations, only in the streamwise direction, at Re ≈ 60 in their FCC sphere pack simulations. They showed that the periodic oscillations are the consequence of their eight vortex mode rotating around x-axis, alternating between clockwise and anticlockwise rotations, which still maintains ∕2 rotational symmetry but looses reflectional symmetries along y-and z-direction at this point. Using a sequence of cross-sectional velocity field visualisation, they described the alternation of rotational direction in terms of periodic merging and coalescing processes of the base-state eight-vortex mode, which is triggered by initial perturbation to the steady-state base mode. Subsequently at Re ≈ 90 , their flow field looses even the rotational symmetry, giving a rise to quasi-periodic oscillations in all spatial directions via bifurcation of frequencies. It is important to clarify the fact that Hill et al. started their simulations from fully-developed unsteady flow fields to avoid the long development time required especially for lower near-critical Reynolds number cases.
We observe the first appearance of subtle flow oscillations at Re steady = 91 (the UNL1 case, cf. Fig. 8b), which coincides with the first break-down of the aforementioned eightvortex mode as we will discuss in the following. At Re steady = 138 , more noticeable flow oscillations develop (cf. Fig. 8c).
At both Reynolds numbers, the fully symmetric eight-vortex mode with the shelllike streamwise jets appears almost immediately during the initial flow development (cf. Fig. 9a). The cross-sectional profile of the eight-vortex mode is remarkably similar to the one we observe at Re steady = 59 although those instantaneous Reynolds numbers are significantly different ( Re(t) = 59 and 141 for Figs. 5c and 9a), thus this  Fig. 5c. Colour scheme for the iso-surface is as of Fig. 6, in addition to negative streamwise velocity coloured in purple in a. a, c To improve the visibility, only one side of the shell-like positive streamwise velocity jet is shown non-dimensional number can be ruled out as a governing parameter of this process. Similarly, the time-scale required for the eight-vortex mode to appear cannot be scaled by the viscous time-scale we apply in this paper.
Subsequently, the four inner vortices become unstable and merge into two vortices giving a rise to the overall six-vortex mode (cf. Fig. 9b). This new flow structure has a smaller degree of symmetry than the previously-observed four-and eight-vortex modes: While the former has two reflectional symmetry lines, the six-vortex mode is only -rotationally symmetric around the centre of the octahedral pore. Notice, whilst the local symmetry in the octahedral pores is broken as above, the two reflectional symmetries in the triangular pores are still maintained. The vortex merging process continues inside the octahedral pores and leads to the two inner vortices merging into one inner vortex, forming a new five-vortex mode, which again has the above local -rotational symmetry (cf. Fig. 9c). Notice the similarity of this vortex merging process to the one observed by Hill and Koch (2002). Meanwhile in the tetrahedral pores, a similar symmetry breaking finally takes place at this point: starting from two symmetric vortices and ending in three asymmetric vortices. This subsequent, not simultaneous, local symmetry breaking indicates the emergence of instability in the more spacious octahedral pores propagating towards much narrower therefore more restrictive tetrahedral pores downstream. This instability propagation process agrees well with the aforementioned observation by Dybbs and Edwards (1984). Transformation from the reflectional symmetry to the rotational symmetry can also be identified in the corresponding 3D view shown in Fig. 10. Conversely, the global symmetry around z∕d = √ 6∕3 is still maintained (cf. Fig. 8). Therefore, we characterise this unsteady nonlinear regime with the reduced local symmetry together with the preserved global symmetry.
Finally, in addition to the streamwise velocity and vorticity contours, we also plot vortical structures visualised by 2 -criterion (Jeong and Hussain 1995) in Fig. 10d. Those extracted vortices are mainly quasi-streamwise vortices.
For completeness, development of superficial streamwise velocity of this flow regime is depicted in Fig. 11a. It is visible that the local flow oscillations which appeared with the local symmetry breaking are reflected in the bulk flow oscillations. Moreover, although the global symmetry is still maintained, the superficial streamwise velocity signal for Re steady = 138 already exhibits chaotic fluctuations, which indicates that it is not sufficient to measure such bulk quantities to detect the transition to turbulence.

Global Symmetry Breaking and Emergence of Turbulence
In this final part of the discussion, we analyse the flow cases in which the symmetries of the pore-scale flow structures, which are imposed by the pore geometry itself, are now entirely broken. Additionally, the bulk quantities, such as superficial velocity, shows fully chaotic signals. By this definition, we are tempted to call this type of flow "turbulence", which we distinguish it from the term "chaotic flow". Those two terminologies are often used  Fig. 7. d Vortex visualisation using 2 -criterion (Jeong and Hussain 1995). The surface of the iso-contours are coloured by the sign of x (red for positive and blue for negative). a, c To improve the visibility, only one side of the shell-like positive streamwise velocity jet is shown interchangeably, mainly due to the difficulty involved in measuring the highly dynamic pore-scale flow structures precisely. However, as it was demonstrated in the previous section, chaotic oscillations of spatially-averaged flow quantities do not necessarily mean the pore-scale flow field is turbulent. On the other hand, the total breakdown of symmetries is also not the solely-defining indicator, as Hill and Koch (2002) showed for their FCC close pack simulations that the bulk flow quantities can exhibit quasi-periodic oscillations even after the geometrical symmetries vanish. Only after Reynolds number exceeding approximately 100, they reported chaotic fluctuations of the flow quantities. Conversely, Dybbs and Edwards (1984) observed first appearance of chaotic oscillations for Re > 300.
By increasing the Reynolds number of our sphere-pack flow, we first observe that the flow structures continue to experience the identical sequence of topology evolution (cf. Fig. 12a and b for an example from the Re steady = 347 case), from eight-vortex mode to six-vortex mode, before proceeding into more complex and less symmetric flow structures (cf. Fig. 12c). Eventually the global symmetry around z∕d = √ 6∕3 breaks between Re steady = 138 and 209, indicating an emergence of "turbulence" (cf. Fig. 8c and d). Notice that the critical point is significantly lower than the value from Dybbs and Edwards, whereas noticeably higher than the value from Hill and Koch. Figure 13 visualises an instantaneous flow field of the T4 ( Re steady = 347 ) case. Conversely to the previous flow regime, the pore-scale vortices in this flow are no longer exclusively quasi-streamwise, but aligned in more diverse directions especially in the octahedral pores (cf. Fig. 13d). Notice that the dimensions of those vortices become significantly finer in comparison to the previous visualisation (Fig. 10d), as it would be expected due to the higher Reynolds number. As the result of the turbulent fluctuations induced by those vortices, the other flow structures are highly distorted and displaying the existence of a wide range of spatio-temporal scales inside the pore spaces. Turbulent energy spectra from two turbulent cases (T1 and T4) in Fig. 14 quantify such range of scales existing inside the pores. In addition, the figure also shows the effect of increasing Reynolds number that expands the scale range towards finer scale, whilst more noticeable inertial range, which is characterised by the −5∕3-exponent, starts to appear. To establish more clear scale separation between the pore-scale and the much finer dissipative scales, further increase in Reynolds number is required which is out of the scope of this paper.
The primary flow separations in the downstream of the sphere contact points are prominent at this Reynolds number (cf. Fig. 15a), and the d-periodic translational symmetry  . a-c Colour scheme for the iso-surface is as of Fig. 7. d Vortex visualisation using 2 -criterion (Jeong and Hussain 1995). The surface of the iso-contours are coloured by the sign of x (red for positive and blue for negative). a, c To improve the visibility, only one side of the shell-like positive streamwise velocity jet is shown 1 3 in x-direction is clearly broken. Note that observation of such symmetry breaking in the streamwise direction was only possible due to the current numerical domain being twice as long as the minimum repeating unit in x-direction. The above separation leads to a clear distinction between the high-speed, highly ordered flow in the tetrahedral pores and the highly turbulent flow in the octahedral pores. It is also interesting to notice that the secondary separations can be found only in the top part of Fig. 15a, where the red-coloured jets fail to trace the curvature of the spheres. Figure 15b depicts the time-averaged velocity field. Although the averaged field has not yet converged completely, it is readily visible that the aforementioned one-sided secondary separation is an instantaneous phenomenon, since the top-bottom reflectional symmetry about the horizontal line through the contact points has been recovered in this view. Similarly, the translational symmetry in x-direction is also present. In Fig. 15c, the corresponding mean turbulent kinetic energy (TKE) field is shown. This field coincides well with the earlier observation from the instantaneous fields that the flow in the octahedral pores is highly turbulent, whereas in the tetrahedral pores the flow is rather quiescent. The maximum energy is located approximately at the end of the recirculation zones.
Finally, a cross-sectional view of the mean streamwise velocity and vorticity fields are plotted in Fig. 16, in which the pore-scale flow structures from the acceleration phase reappears by the averaging operation (compare with, for example, Fig. 12a). The distribution of the streamwise velocity inside the octahedral pores is, however, visibly more diffused than the counterpart from the acceleration phase (cf. Fig. 16a). This difference can be accounted for the emergence of turbulence inside the octahedral pores that scatters around the momentum, whereas inside the tetrahedral pores where the flow is significantly more quiescent, the velocity profile remains almost unchanged. This turbulent diffusion process consequently reduces the streamwise velocity gradient inside the octahedral pores, resulting a reduction in the streamwise vorticity with respect to the eight-vortex state during the acceleration phase.

Summary
We performed a series of fully-resolved direct numerical simulations of accelerating porous media flow through hexagonal close pack spheres. The considered Reynolds numbers range between 5.6 × 10 −7 and 347, which covers steady linear, steady nonlinear, unsteady nonlinear and turbulent flow regimes. Our special focus was placed on the pore-scale flow structures inside the porous medium, such as vortices and jets, and their evolution process during the accelerating phase especially.
We confirmed that the emergence of the nonlinearity takes place when the superficial-velocity-based Reynolds number is between 1 and 10. We found the rise of the nonlinearity is related to the formation of pore-scale vortical structures, whilst two distinct types of flow structures in the steady nonlinear regime, both of which preserve global and local reflectional symmetries on the cross-stream plane, were observed. With increasing Reynolds number, the symmetries local to the individual pores start to break  (Cabral and Leedom 1993) of the velocity vector parallel to a local symmetry plane defined by − √ 3 3 y + √ 6 3 z = − 1 2 . c time-averaged turbulent kinetic energy on the same plane. Data is from T4 ( Re steady = 347 ) case and flow starts to oscillate, which marks a rise of unsteady nonlinear flow regime. Based on a sequence of flow structure evolution, we showed that the symmetry breaking takes place in larger pore spaces before propagating into the downstream narrower spaces. At this point, the volume-averaged flow quantities exhibit chaotic fluctuations despite the fact that the global symmetry remains unbroken. Eventually, all the symmetries, both local and global, are broken and the emergence of pore-scale turbulence is achieved.
During the initial development phase, it was found that steady/unsteady nonlinear and chaotic flow exhibit a consistent sequence of flow structure development, although the applicable Reynolds numbers vary significantly (from Re steady = 36 and 347). Moreover, we found that the appearance of such consistent flow structure sequence cannot be explained by the viscous time-scale that was used for normalisation throughout this paper. Determining the governing parameters of the flow structure development needs to be carried out in our future studies.
The key limitation of the current study is the single sphere-pack geometry which is exposed to the pressure gradient in one direction only. This is a consequence of the high computational efforts necessary to simulate such extensive range of Reynolds number covering all four distinct flow regimes that are known to us today. Although certain level of generality can be assumed for the current findings based on both qualitative and quantitative agreements with the literature, it would still be highly interesting to extend the current study to different sphere geometries and/or different pressure gradient directions as in Hill et al. (2001b), Hill and Koch (2002). By doing so, our understanding in the flow structure development in the transient porous media flow will be more complete, and might provide us, for instance, an important basis of effective flow control.
Acknowledgements Open Access funding provided by Projekt DEAL.