The Effects of Periodicity Assumptions in Porous Media Modelling

The effects of periodicity assumptions on the macroscopic properties of packed porous beds are evaluated using a cascaded Lattice-Boltzmann method model. The porous bed is modelled as cubic and staggered packings of mono-radii circular obstructions where the bed porosity is varied by altering the circle radii. The results for the macroscopic properties are validated using previously published results. For unsteady flows, it is found that one unit cell is not enough to represent all structures of the fluid flow which substantially impacts the permeability and dispersive properties of the porous bed. In the steady region, a single unit cell is shown to accurately represent the fluid flow across all cases studied


:
Ambient density ( ≈ 1) : Velocity vector u: x-direction velocity component v: y-direction velocity component x i : Spatial coordinate vector t: Time coordinate D c : Obstructing cylinder diameter L c : Unit cell side length L: Macroscopic cell side length Δp: Macroscopic cell pressure drop : Porosity p * : Dimensionless pressure drop p s : Flow driving pressure gradient 1 3 x: Element size ( = 1) t: Time-step size ( = 1) N t : Number of time steps Re p : Particle Reynolds number Re: Generic Reynolds number D L : Longitudinal dispersion D T : Transversal dispersion D * L : Dimensionless longitudinal dispersion D * T : Dimensionless transversal dispersion T:

Introduction
Transitional and turbulent flows through porous media typically take in natural processes such as flow between vegetation (Nepf 1999) and stones in rapids as well as in industrial processes including flow around electronic components (Wibron et al. 2018;Emelie et al. 2019), within metal foam heat exchangers (Jamarani et al. 2017), catalytic convertors (Jouybari et al. 2016), drying pellets (Ljung et al. 2012) and packed bed combustion (Okuyama et al. 2011). The flow regimes are also relevant when considering internal erosion in embankment dams, being one of the primary breach mechanisms (Foster et al. 2000). Therefore, a basic understanding of the flow and the forces that occur in a porous media is crucial to capture the deformation of the material which may lead to dam failure (Lundstrom et al. 2010;Hellström et al. 2010;Frishfelds et al. 2014Frishfelds et al. , 2011. Traditionally transitional flow within porous media is dealt with by measurements of global parameters such as overall pressure drop and computations of empirical correlations like the Forchheimer equation (Dixon and Cresswell 1986;Gunn 1987;Hlushkou and Tallarek 2006). Global parameters can also be derived through modelling from first principles, which often include some types of homogenisations (Kuwahara et al. 1998;Pedras and De Lemos 2000). To exemplify, the validity of the Forchheimer equation for high Reynolds number flow has been investigated, e.g. (Soulaine and Quintard 2014). Often terms added to Darcy's law are investigated so that high Reynolds number laminar flow through porous media can be predicted on a global level. Also, turbulent flow through porous media has been treated theoretically, e.g. (Skjetne and Auriault 1999). Hence, a number of relationships exist for unbounded flow through porous media.
In addition to treating the porous medium as homogeneous, simplified geometries have been studied, with analytical, numerical and experimental methods. This includes network models as, for instance, described in E (1974); Dullien and Dullien (1992). It also includes detailed studies of periodic or random arrangements of cylinders e.g. (Lundström and Gebart 1995;Hellström et al. 2010;Lundstrom et al. 2010;Fabricius et al. 2016) or spheres, e.g. (Frishfelds et al. 2014;Koch and Hill 2001;Nijemeisland and Dixon 2004). When studying periodic models, a common tactic is to define a geometric unit cell. Such an approach is likely to be valid as long as the Reynolds number is low enough, but as transitional effects appear, it will here be shown that this kind of approach may be inadequate. This observation is not new, and several unit cells are usually studied when large-scale flow structures are expected to occur in the flow Jin et al. 2015). The cases for which this is considered is often limited to studies where the pore scale of the geometry is comparable in size to the macroscopic bed scale, i.e. where it is obvious that assuming an ideal bed of infinite extent would be erroneous.
Therefore, in this study, the effect of periodic assumptions in 2D flows is investigated. Simple packings of ordered arrays of cylinders will be studied using a cascaded lattice-Boltzmann solver with the number of unit cells represented in the simulated domain being varied. As will be shown, a single unit cell is, for most packings, not enough to accurately describe the macroscopic properties of a large bed when the flow is unsteady. A 3D study is also performed comparing two similar, but dimensionally differing cases to assess how well this study translates to more physically relevant 3D flows. The study finally yields a couple of interesting results for transitional flow through porous media, in general, that are presented and discussed.

Porous Media Geometry
Two types of porous media consisting of ordered arrays of cylinders are considered, one is a simple cubic geometry, and the other is a staggered packing, see Fig. 1. The porosity ( ) is defined as the fraction of void space in which fluid flows through (Nield and Bejan 2017) within the domain for the arrays studied and is obtained as where D c is the diameter of the obstructing cylinders and L c is the side length of the quadratic unit cell, see Fig. 3.

Flow Regions in Porous Media
Focus here is set on three flow regions, the Darcy region, the laminar-steady region, and the unsteady region. More regions exist (Wood et al. 2020), but are not of interest in this study. The dimensionless quantity defining the ranges of these regions is the well-known Reynolds number that in this work is defined as Here U D is called the interstitial velocity being the average stream-wise velocity of the fluid in the bed and is the kinematic viscosity. Notice that multiple definitions of porous bed Fig. 1 The two types of packings investigated. A simple cubic array of mono-radii cylinders (left) and a staggered packing of mono-radii cylinders (right)

3
Re have been suggested, see Ziolkowska and Ziolkowski 1988) and Hellström et al. (2010) for comprehensive reviews. The length scale of obstructing cylinder diameter ( D C ) is chosen for this work since the packing densities studied are relatively sparse. The flow regions are defined by the dispersive properties, the relation between pressure drop and velocity as well as the flow field properties. Sometimes the unsteady region is separated into a laminar-unsteady and turbulent region; since 2D-fluid models without turbulence models are of interest to this study, these two regions are combined into the unsteady region instead.
The Darcy region is characterized by a linear relationship between the interstitial velocity and the pressure gradient across the bed, i.e. U D ∝ Δp∕L where Δp is the pressure drop and L is the length of the bed in the stream-wise direction. In the laminar-steady region, the convective pressure losses become more significant than the viscous losses and, therefore, a second-order relation is obtained giving U D ∝ √ (Δp∕L) . The unsteady region is not defined by any relation between the pressure drop and velocity, although the convective losses tend to dominate giving U D ∝ √ (Δp∕L) , i.e. the same as in the laminar-steady region, Forslund et al. (2021). This region is instead defined by in time-fluctuating velocity and pressure fields, i.e. | d dt | ≠ 0 in at least one part of the domain.

Computational Method
The fluid dynamics solver utilized in this study is an in-house lattice-Boltzmann method (LBM) solver based on the models described in Geier et al. (2006); Lycett-Brown and Luo (2014). The code is implemented in the CUDA framework (NVIDIA 2015), and the computations are carried out on a GPU architecture using 32-bit floating point precision (IEE754). The LBM method recreates the Navier-Stokes equations emergently by starting from the Boltzmann equation. The method has increased in popularity ever since it was introduced due to an increase in accuracy and stability of the numerical schemes being developed, as well as its favourable parallelizability properties (Delbosc et al. 2014;Delbosc 2015). The cascaded LBM method applied in this work was first introduced by Geier et al. (2006) in 2006, with the argument that the lack of Galilean in-variance and cross-talk between moments was the main contributor to the often seen instabilities of the multiple relaxation time (MRT) and single relaxation time (SRT) schemes. The main idea is that the moments are transformed into the rest frame before being relaxed and then being transformed back. A Galilean transformation is therefore applied to the velocity distribution prior to the calculation of the moments. For further details on this see the works of Geier et al. (2006) and Premnath et al ) who provide an in-depth derivation of the cascaded LBM model. This type of LBM model was used since it can handle large cell Reynolds number flows without diverging, a shortcoming that is present for both the MRT and BGK methods.
Since the LBM model is a type of artificial compressibility model, special considerations need to be taken to avoid compressible effects from becoming prominent. Generally, the velocities in the domain, max(|U|), should be well below the speed of sound. Since a length scale x and time scale t of unity are prescribed, the speed of sound is given by . This condition therefore equates to where |U(x i , t)| is a dimensionless measure of the magnitude of the velocity for all points in space x i and time t. To assess how significant the compressibility effects are, the range of the density can also be inspected, giving the condition For all simulations, here presented, the maximum velocity U max is kept below 0.05 and the maximum density range range is about 0.01, i.e. ≈ 1% of the density of the fluid. Hence, compressibility effects have minor to negligible impact in this work. The lattice spacing x and time step t are set to unity for all cases, when simulation inputs or values are discussed they are in lattice units and dimensions are dropped, the variable discussed gives the dimensional context.

Particle Tracking
Particle tracking is carried out using a corrector method to avoid wall entrapment. From testing it was concluded that this corrector step was reduced the amount of particles being trapped in wall cells. The method can be summarized in pseudocode as follows: The digit 2 in front of the difference vector is included to make sure that the particle moves further from the wall than its initial position, preventing most cases of loop entrapment. This procedure is carried out by storing the positions in two buffers that are used interchangeably, one for even iteration numbers and one for uneven numbers. All boundaries are periodic resulting in an artificially larger domain for the particles than each geometry studied.

Theory
This section shortly outlines the methods applied for the analysis of the results from the LBM simulations.

Decomposition of Velocity Fields
To examine the effect of including additional unit cells in the calculations a decomposition of the velocities is carried out for the one-by-one cell contributions, two-by-two cell contributions and so on. Fourier analysis yields that the value of some variable v can be represented in frequency form in the following manner if, for simplicity, the side lengths of the X × Y cell are set to 1 2 where N x and N y are some cut-off frequencies that define the smallest spatial length scale in respective direction which is here directly related to the mesh size, and C mn is the Fourier coefficients. Now focus is set on a X × Y domain of repeating unit cells, see Fig. 2. These X × Y lattices of unit cells will be referred to as a X × Y repetition, as in how many times the unit cell repeats in the x-and y-direction, respectively. The frequencies that can exist within a single cell ( 1 × 1 ), a 2 × 2 or an in general X × Y combination of unit cells are described by the expressions where the final expression can be seen to be equal to the original decomposition into frequencies described by Eq. (5). In the present investigation only quadratic arrays of lattices with side lengths that are a power of 2 will be considered. Spatial frequencies are sought for, and therefore the K'th periodic is defined as.
This expression can be derived iteratively starting with the v 1×1 then calculating the next v 2×2 . Note that only values of v K×K evenly divisible by X and Y can be extracted by this procedure, it is, for example, not possible to get v 3×3 from a 4 × 4 unit cell repetition or v 2×2 from a 3 × 3 unit cell repetition.

Apparent Permeability and Pressure Drop
Apparent permeability is a rescaling of the pressure based on Darcy's law (Lundstrom et al. 2010) that can be written as where R is some reference length scale, the diameter of the cylinder in the current case, and L is the cell length. This way of expressing the pressure drop is useful since it gives an indication of how the pressure deviates from Darcy's law. For the LBM simulations in the present study a momentum source/constant pressure gradient, p s = Δp x , is applied to drive the flow. This method of driving the flow is standard practice for porous media flow simulation (Lundstrom et al. 2010;Koch and Hill 2001;Koch and Ladd 1997) and gives the pressure drop across the cell as p = p s L where L is the length of the cell, hence In practice, p s may be considered a pressure difference instead of a pressure gradient since the spatial distance between the mesh nodes is x = 1 . If the pressure gradient acting across an area, p s R 2 , is written as a force gradient F s = ΔF∕ x = ΔF = F and the interstitial velocity is exchanged for the average fluid velocity averaged over the solid and fluid domain, called Darcy velocity, the same formulation as that of Koch and Ladd (1997) and Sangani et al (Sangani and Acrivos 1982) is obtained. In the simulations a constant acceleration is applied to the cell to drive the flow, and this equates to a pressure gradient since ≈ 1 everywhere. The total force acting across the fluid region is then given by.
To be able to compare to earlier results the total force F is defined as an integral over the solid phase as well.

Rest Frame Energy
The total energy present in a flow within region relative to the absolute frame is given by the integral Transforming the velocity into the rest frame of the fluid, a Galilean invariant measurement of the flow energy is obtained as Here U i is the mean cell velocity in direction i, which is defined as This is a useful expression in order to find the amount of energy being transferred into the fluid by the surrounding walls. For steady flows this relation can be directly applied, while for unsteady flows an additional averaging over the time dimension needs to be done resulting in the expression

Enstrophy
The enstrophy is defined as the area integral of the square of the vorticity as (Weiss 1991) where the area is the fluid domain. The enstrophy gives an idea of the vorticity the flow; therefore, a high value indicates higher viscous dissipation. This implies that energy is being transferred from kinetic to thermal; in porous media flow a high value for the enstrophy therefore implies a larger pressure gradient.

Dispersion
Dispersion calculations are based on a Lagrangian approach by tracking satellite particles in the flow, as described in Zhao et al. (2010). A dimensionless dispersion measurement can be obtained by forming a dimensionless group from the interstitial velocity U D , particle group time t p and some particle cloud distribution length scale L D as A suitable variable for the particle cloud extension length scale L D may, for example, be the standard deviation (Zhao et al. 2010), and this is what will be used in this work. The particle cloud extension is separated into two parts consisting of transversal spread L D T and longitudinal spread L D L yielding two dispersion measurements as.
With this Lagrangian tracking approach a packet of particles is released at some position A and after some time has elapsed the distribution is saved at point B. The point cloud is then analyzed to obtain relevant distribution quantities.

Geometry and Boundary Conditions
The geometry of the elementary cell consists of a square domain with side lengths L c and a solid central cylinder of some diameter D c < L c , see Fig. 3.
Periodic boundary conditions are applied, and for the cubic case, the mappings are direct, while for the staggered case the outlet to inlet connection is offset in the y-direction by L c ∕2 , see Fig. 4. Note that the tiling utilized by the 2 × 2-staggered configuration can also be utilized by the 2x2-cubic configuration retaining an identical geometry. There exists freedom of choice in how to carry out these tilings while still preserving the crystal structure of the geometry. The flow is driven by a body force acting on all fluid elements. For additional details on how the cascaded Lattice-Boltzmann method with the forcing scheme is implemented see . between viscosity values. Equations (2) and (8) indicate how the variations in influence Re p and p * when the driving pressure gradient p s (in units of x∕ t 2 ) is kept at a constant value for all viscosities. Four million time steps N t = 4M calculated at each viscosity give a well-converged mean even when low-frequency flow field fluctuations are present. The variation of the pressure gradient p s between the cases studied is presented in Table 1.

Parameter Sweep
After four million time steps relevant data are saved and the viscosity is changed directly, meaning that hysteresis effects may be present. The pressure gradient is adjusted for the 1x1-periodic lattice since the resistance is significantly lower compared to the larger unit cell repetitions. This leads to very high Mach numbers with max velocities in the range of U max ≈ 0.14 resulting in to non-negligible compressibility effects if the momentum source term is not reduced.

Mesh
The mesh is a uniform quadratic lattice of computational cells that are set as either wall elements or fluid elements. The circular geometries are approximated by an edge exclusive formulation setting all cells within the circle as walls. The top right quadrant for the 64 × 64 , 128 × 128 and 256 × 256 geometries are illustrated in Fig. 5. An investigation of the mesh dependence concluded that a 128 × 128 mesh is a fair compromise between

Mesh Analysis Setup
The mesh analysis is carried out for the D c = 0.6L c packing and 1x1-representation case by varying the mesh size between 32 × 32 and 256 × 256 with doubling increments. The value of the acceleration and viscosity is varied to preserve dynamic similarity and Mach number between the cases. From Eqs. (2) and (8)

the following relation is obtained
To ensure that unsteady effects are averaged over the same amount of time, the number of time steps N t is scaled between the meshes according to the relation The simulation input variables are set up according to Table 2.

Particle Tracking Setup
At an increment of 2000 time steps a packet of L c Y number of particles is released along the line x = 0, y = 0 − L c Y , see Fig. 6 for an illustration for the D c = 0.5L c 1 × 1-repetition case. These particles are tracked for a total of 200k time steps before terminating. All particle positions are saved with an increment of 200k time steps giving a total of 20 different point clouds for each viscosity, during the simulation period of 4000k time steps. The particle positions are updated every twentieth time step; since the velocity in the computational cells has a maximum value of ≈ 0.05 in lattice units, it implies that any particle never moves more than approximately one computational cell between each update. The dispersion is calculated by a procedure similar to that presented in Zhao et al. (2010). Since the original particle distribution already has an extent in the y-direction at release, the dispersion in the transversal direction, D T , is defined as here is the standard deviation operator and (t) is the particle cloud distribution in the y-direction at time t. Equation (20) ensures that the measurement of the dispersion in the transversal direction is only relative to the initial dispersion. For the longitudinal dispersion, ( (0)) is already zero initially and does not need to be adjusted. Each particle cloud has 100 groups, and these groups contain 128 particles that are released at the same instant. Plotting the dimensionless dispersions D * T and D * L as defined by Eqs. (16) and (17) for each of the groups a discrepancy in the early particle time values is seen (t < 50000 t) , see  where * and * are arrays containing the dimensionless spread for each of the particle groups and N is the number of particle groups. The quantities D * T and D * L are now averaged over all particle clouds that were released during the N t = 4M calculation period yielding the final value for D * T and D * L . For the 2 × 2 , 3 × 3 and 4 × 4-repetition cases each released group is divided into packets of 128 particles in the analysis as indicated in Fig. 8. This gives a higher number of total particle clouds for each of the X × Y repetition cases scaling by the value Y.

Geometry
The 3D-comparison geometry consists of a unit cell with side lengths L c and a central cylinder of radius D c with a wall at the top and bottom, see Fig. 9.

Setup
The Darcy region is not interesting for the 3D-comparison study since initial simulations yielded no deviations between the unit cell repetitions in this region. Therefore the viscosity range is limited compared to the other sweeps. Twenty five sweeps are still performed, but the total amount of time steps ( N t ) for each viscosity is lower. The simulation input variables are set up according to Table 3.

Code Runtime and Hardware
The in-house CUDA program was used to perform the sweep, and the calculations took a total of ≈ 30 hours on an NVIDIA Quadro RTX4000 GPU.

Impact of mesh resolution
The mesh study is carried out with a porosity of ≈ 0.7 ( D c = 0.6L c ) for one unit cell and the mesh resolutions {32 × 32, 64 × 64, 128 × 128, 256 × 256} . The dimensionless pressure drop p * is presented for a cubic cylinder packing, see Fig. 10. It is clear that the difference between the two finest grids is small with an average difference in predicted pressure drop of ≈ 1.8% , for the 128 × 128 compared to the 256 × 256 . The 64 × 64 grid overestimates the pressure drop, and the 32 × 32 underestimates it. Therefore the 128 × 128 grid unit cell is a fair compromise between performance and accuracy.

General Flow Behaviour
In this section, the general flow behaviour is discussed and a comparison to earlier work on pressure drop in cubic arrays of cylinders is made. Some obvious flow patterns spatially ranging over more than one unit cell are also shown. The sweep data are analyzed more quantitatively in Sects. 5.3-5.6.

Flow Regions and Comparison to Earlier Work
The three flow regions for the 1 × 1 cubic array (Darcy, laminar-steady and unsteady) for the D c = 0.3L c , D c = 0.5L c and D c = 0.7L c cases are presented in Fig. 11 by streamlines of the time-averaged velocity field. For less dense packings, D c ≈ 0.1L c − 0.5L c , the unsteadiness starts in combination with the recirculation zone development. This leads to a relatively short laminar-steady region that continues until the flow becomes unsteady behind the cylinders, i.e. vortex shedding begins. For the denser packing, D c = 0.7L c , only smaller flow structures can exist implying that a higher Re p is necessary for these structures to have an impact on the macroscopic pressure drop. This also means that large recirculation zones may exist in the wakes of the cylinders without vortex shedding. The Darcian-flow region is characterized by a linear relation between the flow through the bed and the pressure drop across it (Lundstrom et al. 2010); in other words, the following relation between U D and Δp s holds true (23) U D ∝ Δp s . → p * = const. Since Re p ∝ U D this region is characterized by a straight vertical line between Re p and p * as marked in green in Fig. 12 where simulations were carried out on a 4 × 4 cubic array. In the laminar steady regime inertia effects make a relatively small but noticeable impact on the pressure drop; this regime is marked with blue. The transition between the Darcian and laminar regime in Fig. 12 is marked by visual inspection. As already discussed in introduction the onset of the laminar region is dependent on the geometry. As the flow becomes unsteady, there is a sharp increase of the pressure losses as marked in red. Unsteady onset is defined as regions where the average velocity fields ̄ (x i ) do not agree with the instantaneous velocity field (t, x i ) , i.e. the condition (t, x i ) ≠̄ (x i ) for any t. From Fig. 12 it is also obvious that the packing density has a large impact on when unsteady effects are of importance for the pressure losses. Using the definitions of Sangani et al. (1982) it is possible to compare the results of these LBM simulations to their analytic expression for a dense packing that may be written as and for the non-dense packing case the following expression is used instead here c max is the maximum possible value of c before the cylinders touch, i.e. 4 . As mentioned earlier the expression p * can be written as F U if the numerator R 2 p s is interpreted as a force. Sangani (1982) and Koch and Ladd (1997) chose to define this force as an integral of the acceleration across both solid and fluid phase. In Table 4  For the case of D c = 0.7L c Sangani expression (24) over-estimates the Darcy resistance for border-line dense packings; therefore, the difference of ≈ 5% is expected. Koch and Ladd (1997) fitted an expression for the dimensionless pressure drop across dense packings of cubic arrays of cylinders around the onset of the Forchheimer region yielding . where is the void space of the bed. The fit described by (26) is compared to the simulations in Fig. 13. The expression closely agrees with the densest packing in the transition region between the Forchheimer and Darcy region, see Fig. 13. The expression then quickly diverges since the relationship p * ∝ Re 2 only holds true over a short range. The poorer agreement for the sparser packings is anticipated since (26) is derived for dense packings.

Relations Between Spatially Averaged Flow Variables for the 1 × 1 Cell
In the unsteady region the mean flow variables of the 1 × 1 cubic cell tend to oscillate significantly, see   . 13 The dimensionless pressure drop using Sangani et al. (1982) formulation F U and the particle Reynolds number Re p plotted for the numerical simulations with regions marked according to the specifications in Fig. 12 the simulation ( = 4M ). The behaviour is similar to that of transitional pipe flow, implying that the friction factor is significantly lower for the steady region as compared to the unsteady region. This causes a build-up of velocity resulting in a higher Reynolds number which causes a higher friction factor which dissipates the flow energy and brings the system back to the lower Re p steady region. The transfer of flow momentum to the porous bed can be observed in Fig. 14 as rest-frame energy increasing when the interstitial velocity and total energy decrease. The enstrophy embodies the vorticity of the flow and is therefore a measure of how much energy the vorticity is dissipating. Huang et al. (2013) investigated a similar system experimentally showing that the friction factor exhibits a sharp discontinuity at the transition Re p . Similar results are also obtained by Khayamyan et al. (2014) when studying transitional flow through a pore-doublet model.

Higher-Order Periodic Effects
The 1 × 1-unit cell repetition geometry puts limitations on the kinds of spatial structures that can exist in the flow. Spatial structures in the velocity field larger than a single cell have been observed experimentally by Forslund et al. (2021in the shape of an alternating lateral velocity. This effect also occurs in these 2D-simulations at the packing D c = 0.7L c as presented in Fig. 15. The effect presented in Fig. 16 is an example of another large-scale average flow structure that cannot be represented by a single unit cell. The middle channels have different core pressures and velocities as compared to the channel in the bottom and top of the figure.
As is apparent from these examples the packing densities each have varying largescale flow features which cannot be contained in a 1 × 1-representation. Since the onset of unsteady effects is limited by the largest fluid scales that can exist in the simulated domain, it is expected that the unsteady onset should occur at lower Re p for the 4 × 4 , 3 × 3 and 2 × 2 representations as compared to the 1 × 1 representations. This is also seen in Fig. 17, where the unsteady onset is defined as in Sect. 5.2.1. For the packing D c = 0.3L c the unsteady onset differs by almost a whole order of magnitude between the 1 × 1 and other representations. As the packing becomes denser, the difference becomes progressively smaller.

Pressure Drop
Using (2) and (8) the dependence of the dimensionless pressure drop across the cells to Re p can be calculated. See Fig. 18 where this relationship is plotted for different geometry representations and porosities.

Sparse Packings
From Fig. 18 it is clear that the cubic 1 × 1-representation underestimates the pressure drop significantly in the unsteady region for cases with D c ≤ 0.5L c . This is also true for the staggered packing when 0.3L c ≤ D c ≤ 0.5L c . In these domains and in reality flow structures larger than a single unit cell are free to exist since the packing density puts little constraint on the scales possible. From the unsteady onset analysis presented in Fig. 17 it is apparent that independent of how these structures form physically they are Fig. 15 Streamlines of the 4 × 4 -average flow field for D c = 0.7L c and Re p ≈ 1000 , the streamlines are coloured by the average transversal velocity Fig. 16 Streamlines of the 2 × 2 -average flow field for D c = 0.4L c and Re p ≈ 1000 , the streamlines are coloured by the a scaled average density, which is linearly correlated with the pressure larger in size than a single unit cell. This has been a topic of some controversy but as demonstrated by Jin et al. (2015) macroscopic turbulence is possible for sparse 3D packings, although, as is pointed out, these scales generally do not extend far beyond the pore scale.
Experimental work on imaging the flow structures of interacting cylinders has been reviewed by Zhou et al. (2016) and experimentally investigated by Ziada (2006) providing clear visualizations of anti-symmetric vortex interactions in the shape of pairs of counter-rotating vortices. This anti-symmetric vortex interaction cannot be represented in the 1 × 1 cell; therefore, the fluid dynamics of the vortex interactions will differ substantially between the case of a single unit cell as compared to two or more repetitions. As discovered by Zhou et al, the mechanism by which the drag coefficient is varied in systems of two interacting cylinders is by the interactions of the vortices being shedded by said cylinders. In other words, depending on the flow conditions and cylinder arrangement different interactions between the shedded vortices will arise and cause differing forces to act upon the cylinders. Zhou et al. listed several possible modes of interaction between two cylinders in a system similar to that presented in this work. It is therefore likely that the lack of such vortex interactions in the 1 cell case is the reason for the large deviation in dimensionless pressure drop displayed in Fig. 18  for both packings. It is also expected that absence of the interacting vortices will delay the unsteady onset as disclosed in Fig. 17.

Dense Packings
For the comparatively denser packings D c = 0.6L c to D c = 0.8L c the differences between the 1 × 1-representation and higher unit cell repetition geometries are less in comparison with the sparse packings. Large-scale flow structures are though still present and clearly seen in the flow presented in Fig. 15. In this region the 1 × 1-representation switches from underestimating the pressure drop to overestimating it for the highest Re p region and for both packings. This switch is easiest seen in Fig. 18 for D c = 0.7L c and is also visible for D c = 0.6L c . For the cubic packing and the D c = 0.7L c case a pressure discontinuity is observed at Re p ≈ 200 for all representations except for the 1 × 1-case. The discontinuity is due to the same alternating lateral velocity effect presented in Fig. 15.

Dispersion Analysis
The analysis of the dispersion yields that in the Darcy and laminar steady regions there is no difference between the unit cell repetitions of the simulated domain, while the 1 × 1 case deviates for the unsteady cases, see Figs. 19 and 20 which contains the dispersion information for all packings and Reynolds numbers. Hence, the behaviour of both the longitudinal and transversal dispersion mirrors the other results presented in this article, cf. Fig.18, for instance. The results in Figs.19 and 20 will now be scrutinized in more detail.

Darcy and Laminar Steady Region
In the region where Re p is small the values of D * T go to zero for both the staggered and cubic cases for all porosities/packing densities. Since the Peclet number is infinite, the separated flow regions do not intermingle and all cases form a single central channel which the particles reside in. For the cubic packing the longitudinal dispersion D * L tends to increase with increasing packing density. An analytical expression for the lateral dispersion around a dilute random array of cylinders was derived by Eames and Bush as (1999) where is the porosity and C is a constant dependent on the packing type, ≈ 0.74 for a random packing case. With the scaling employed earlier (Eq. (17)) this expression becomes Here the averaged over time is excluded. It is seen that the dispersion should scale with the volume of the solid material and the diameter of the cylinder yielding a D * L ∝ D 3 c relation. As the lattice becomes denser the dispersion starts scaling by D * L ∝ D c as discovered Brady & Brenner (Brady and Brenner 1989). The longitudinal dispersion in the Darcy region is presented for all the cells in Fig. 21. The figure shows that in this respect all the packings studied can be considered dense since the D * L ∝ D c relation approximately holds except for a slight power law behaviour for D c = 0.1L c − 0.3L c .

Unsteady Region
The results for the dispersion in the unsteady region agree with the other results from the unsteady onset analysis (Fig. 17), and the pressure drop (Fig. 18). The differences in dispersion between the 1 × 1-repetition and higher-order cases tend to be comparatively larger for the cubic and sparse ( D c = 0.1L c − 0.5L c ) packings. Although the agreement is better between the repetitions for the staggered cases, it is still inaccurate. It should be noted that the 2 × 2 , 3 × 3 and 4 × 4-repetitions tend to deviate from each other for both the staggered and cubic cases. Even more-so than it did for the case of pressure, this indicates that accurate dispersion calculations in the unsteady region may demand geometric representations in excess of 4 × 4-repetitions for a physically correct solution.

Decomposition Analysis
The velocities can be decomposed according to the procedure described in theory section. Following this decomposition the energy contained in each spatial scale is calculated according to the expression Here u X×Y is the x-direction velocity in the cell decomposition X × Y , v X×Y is the corresponding y-direction velocity, and is the fluid domain. Since incompressibility is assumed, ≈ 1 , and the density can be omitted from the numerical integration. From the pressure and dispersion analysis it is reasonable to expect that where the 1 × 1-representation differs from the higher-order ones large-scale structures should exist. These structures should also contain flow energy. This is visible in the energy decomposition since the 1 × 1 -mode does not contain all flow energy. Quantities at each level are defined by the fraction of flow energy contained. The results from this procedure are presented in Fig. 22.
Comparing the energy percentage in the 2 × 2 and 4 × 4 mode it is clear that the fraction for the D c = 0.6L c and D c = 0.8L c is an order of magnitude lower compared to the  other packing densities at Re p ≈ 1000 , with Dc = 0.7Lc as an exception. Looking at the pressure drop for these packings in Fig. 18 the 1 × 1-representation deviates the least from the higher-order representations. A corresponding pattern of approximate agreement can be seen in the longitudinal and transversal dispersion in Fig. 19. This indicates that large-scale flow features do not impact the pressure drop or dispersion significantly for these packings, at least in the flow region investigated.
Another clear pattern in the 4 × 4 energy fraction is that the lower density packings 0.1Lc − 0.5L c all show a local peak before going into the steady region. The flow field at D c = 0.5L c at this peak is presented in Fig. 23. The flow field makes it clear why this effect can only be represented by a 4 × 4 cell, in this case the large-scale structures take on the shape of varying channel velocities.

Comparison to a 3D case
The 3D comparison is carried out on a unit cell of size 64 3 elements. From the mesh study it is known that this will impact the accuracy of the solution, but the results should still be indicative of the behaviour in a 3D cell as compared to a 2D cell. It is known from the comprehensive review of Boffeta et al. (2012) that 2-dimensional turbulence/unsteady flow Fig. 22 Energy fraction of the 1 × 1 , 2 × 2 and 4 × 4 modes for the cubic packing. Note the different y-axis scaling for the 1 × 1 as compared to the other cases is fundamentally different from 3-dimensional turbulence/unsteadiness. To see how well the insights gathered from this 2D study translates to more physically relevant 3D flows a study on one porosity case has been carried out on a MRT-D3Q27 model based on the work of Suga et al. (2015).
Re p is varied in the same manner by changing the viscosity of the fluid. The main result is that the dimensionless pressure drop p * against Re p has a similar behaviour between the 2D and 3D cases, see Fig. 24.
In the Forchheimer region both the 2D and 3D cases follow a similar trend characterized by overlapping pressure drops for all geometric representations. This indicates that for the steady region a 1 × 1-representation is sufficient to describe the flow in 3D-porous materials. The flow becomes unsteady at Re p ≈ 200 , but while the 1 × 1-representation for the 2D case immediately deviates from the other representations, this does not happen for the 3D case. At Re p ≈ 500 the 3D 1 × 1-representation starts to significantly deviate from the other representations, which is due to the same alternating lateral velocity effect that is observed in Fig. 15 taking place in the 3D case.

Conclusions and Future Work
It is clear from the results for the pressure drop, dispersion and unsteady onset that largescale flow features not representable within a 1 × 1 cell has a significant impact on the macroscopic properties of the porous bed. Therefore pore-scale simulations of unsteady flows require a domain larger than a single unit cell to accurately model the porous media. It is also shown that when the flow is steady only a 1 × 1 cell is necessary to represent all the flow structures, at least for the cases of staggered and cubic arrays of cylinders.
Future work includes a more detailed investigation of the dispersive disagreements between the 2 × 2 , 3 × 3 and 4 × 4-repetitions for some packings and flow conditions. A clear example of where the higher-order representations deviate are the 3 × 3 and 4 × 4 /2 × 2 pressure drops for D c = 0.7L c in Figs. 18 and 24. These results indicate that the periodicity of the lattice may impact the periodicity of the structures that can appear in it. This raises the question of what amount of unit cell repetitions that should be used depending on Fig. 23 Streamlines of the 4 × 4 -average flow field for D c = 0.5L c and Re p ≈ 100 , the streamlines are coloured by the absolute velocity scaled by the interstitial velocity U D the packing that is being modelled, or, put more succinctly, how many unit cells should be modelled to obtain an accurate solution for a given packing and for a given variable.
Funding Open access funding provided by Lulea University of Technology. This work was made possible by funding from the Swedish Research Council, grant 2017-04390.
Data Availability Relevant data is available from the corresponding author upon reasonable request.
Code Availability Relevant source code is available from the corresponding author upon reasonable request.

Conflicts of interest The authors report no conflict of interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.

Fig. 24
Dimensionless pressure variation comparison between the 2D case (top) and 3D case (bottom)