1 Introduction

Fluid flow within a porous layer confined between solid walls is of great importance in a number of natural and industrial applications such as fiber-reinforced composites manufacturing (Nordlund and Lundström 2008; Frishfelds et al. 2010; Tan and Pillai 2012), paper making (Lundström et al. 2002; Singh et al. 2015), microfluidic systems (Jeon and Shin 2009), lubrication of journal bearings (Rao et al. 2013; Almqvist et al. 2017) and different layers of fuel cells (Qin and Hassanizadeh 2015). Therefore, researchers from a wide range of disciplines have been attracted to conduct analytical, numerical and experimental studies of fluid flow within such porous media which are also known as thin porous media (TPM) (Yeghiazarian et al. 2016). According to Fabricius et al. (2016), confining walls affect the Stokesian flow within porous media if the walls are within a certain distance; however, two-dimensional arrangements of cylinders can be used to predict the flow within bounding walls outside this range (Gebart 1992; Lundström and Gebart 1995; Hellström et al. 2010a2010b; Jourak et al. 2013). It is also observed that the flow is two dimensional when the walls are very close together. Fabricius et al. (2016) sorted the TPM into three categories based on the dimensions of its pores and the macroscopic thickness of the porous layer. Flow in very thin porous media (VTPM) occurs when the thickness of the porous layer is much smaller than the pore sizes. Proportionally thin porous medium (PTPM) is attributed to a porous layer whose thickness is proportional to the pore space dimensions, and homogeneously thin porous medium (HTPM) is used when the layer thickness is much larger than the pores. Recently, Qin and Hassanizadeh (2015) presented another categorization of thin porous media: physically thin porous medium and geometrically thin porous layer. In the former one, the porous layer thickness is only one order of magnitude larger than the mean pore size of the porous medium, while in the latter the porous layer thickness is much larger than the mean pore size, but its thickness is much smaller than other dimensions.

The flow regime in most of the studies mentioned above is Stokesian which is governed by Stokes equations from which the Darcy law can be derived after averaging (Whitaker 1986; Forslund 2020). Before proceeding to the different flow regimes in porous media, it is worth noting that there are generally two main approaches to study the fluid flow within the porous media. The first one is the microscale or pore scale in which the governing equations such as continuity and momentum are directly solved within the pores of the porous media to obtain the velocity, pressure and turbulent parameters. The second approach is to replace the complex pore structures with a homogenous medium with averaged properties via volume averaging technique (VAT) where the governing equations are averaged over representative elementary volumes (REV) to obtain the volume-averaged equations. Therefore, the volume-averaged values obtained from these equations are macroscale representations of microscale phenomena. Based on this, Darcy’s law is the macroscale representation of Stokesian flow through the pores. Through visualization of flow within different types of porous media such as a bed of rods in a regular, hexagonal packing, Dybbs and Edwards (1984) presented a categorization of flow regimes in porous media based on the pore Reynolds number (Rep), which is defined using average velocity within the pores, as Darcian regime (Rep < 1), steady laminar (1 − 10 < Rep < 150), unsteady-laminar (150 < Rep < 300) and turbulent regimes (Rep > 300). The particle Reynolds number, \({\text{Re}}_{D} = \phi {\text{Re}}_{P}\), in the present study is defined as:

$${\text{Re}}_{D} = \frac{{\rho U_{D} D}}{\mu }$$
(1)

where \(U_{D} = \phi \left\langle U \right\rangle\) is the Darcy velocity, \(\phi\) is the porosity, and \(\left\langle U \right\rangle\) is the averaged velocity within the pores. The macroscale counterparts of these regimes are presented in Fig. 1 (Horton and Pokrajac 2009). The fluid velocity in most of the above-mentioned studies is too small so that the Darcy law governs the fluid flow within the thin porous medium in which the pressure drop across the medium varies linearly with the fluid velocity. According to Dybbs and Edwards (1984), the range of ReD in which Darcy law governs the fluid flow in porous media is limited to the order of 1 or smaller. The post-Darcy regime, ReD > 1 range in Table 1, is started when the deviation from the Darcy law is observed due to the inertial effects resulting in nonlinear relationships such as Forchheimer’s equation in the macroscale (Hassanizadeh and Gray 1987). The ReD demarking the upper boundary of Darcian regime is too small so that the flow regime in many industrial applications is in post-Darcy regime. Although the post-Darcy flow in porous media has been the subject of many studies in the literature due to its wide range of applications (Dehghan et al. 2014a, 2014b; Jouybari et al. 2014, 2020; Nimvari et al. 2014), the number of studies that focused on the post-Darcy flow in a thin porous media is quite scarce. Driven by this, Larsson et al. (2018) investigated the capabilities of using tomography particle image velocimetry to study transitional flow through confined porous media yielding that the microscopic flow changes considerably as ReD is increased from 45 to 950.

Fig. 1
figure 1

The geometry of the thin porous media in the present study a staggered arrangement of circular cylinders, b random arrangement of circular cylinders, c top view of the xz plane of the geometry in (a)

Table 1 Flow regime categorization in porous media along with associated ReD (Horton and Pokrajac 2009)

To the best of the authors’ knowledge, there is no study in the literature that addresses the turbulent flow characteristics within the pores of a thin ordered porous media. However, there are several studies that have investigated the turbulent flow within the non-confined array of cylinders. Pore-scale study of turbulent flow in an array of cylinders with different cross sections has been initially carried out in the literature with an aim to calculate the macroscopic (macroscale) values of flow and turbulence parameters such as pressure drop, turbulent kinetic energy and dissipation rate (Pedras and de Lemos 2001; Nakayama and Kuwahara 1999; Yang et al. 2014; Torabi et al. 2019). However, due to the incorporation of eddy-viscosity turbulence models in these studies, they were not able to report the values of Reynolds stresses and hydrodynamic dispersion terms from a direct calculation of these parameters within the pores. In order to calculate these parameters, it is inevitable to carry out the large eddy simulations (LESs) or direct numerical simulations (DNSs). However, it is well known that the computational cost of DNS is too high, making it impossible for complex geometries such as pore structures and high Reynolds numbers. Therefore, DNS is used for some simple pore geometries of porous media at lower particle Reynolds numbers, as in Jin and Kuznetsov (2017) and Uth et al. (2016). The problem gets exacerbated in the cases investigated in the present study where the bounding walls exist in the computational domain. In this case, the number of computation meshes needed to resolve the turbulence structures in the domain highly increases which in turn increases the computational cost. In this condition, LES provides an opportunity to directly calculate the Reynolds stress and hydrodynamic dispersion components within the pores with a significantly smaller computational cost since only the large scales of turbulence are resolved in LES instead of resolving all turbulence scales, as in DNS. Therefore, LES is used in the present study to investigate the turbulence effects within the pore space of porous media confined between walls. Although few studies such as Kuwahara et al. (2006) and Kundu et al. (2016) used LES to obtain the macroscopic flow and turbulence parameters in a non-confined arrangement of cylinders, LES study of turbulent flow in a thin porous media confined between walls has not been addressed in the literature.

Motivated by the above discussion, we investigated the fluid flow through a porous medium confined between two walls in a wide range of Re starting from the laminar unsteady regime to the fully turbulent flow (see Table 1). Large eddy simulation is used to simulate the fluid flow and capture the turbulence effects within the pore space. The variation of streamlines and vortical structures with porosity, pore geometry and Re in the presence of confining walls is calculated and compared with the same variables for flow through non-confined porous media. The effects of confining walls on macroscopic quantities such as pressure drop, Reynolds stress and hydrodynamic dispersion are also compared with the non-confined cases in different flow regimes. These effects have not been addressed in the literature since most of the previous studies on thin porous media are, as previously stated, limited to the Darcy regime. Finally, a discussion is provided on the effect of variation in the channel thickness on the microscopic and macroscopic quantities in TPM.

2 Physical Model

Two different arrangements of staggered and random circular cylinders confined between two parallel walls are considered to mimic thin porous media (see Fig. 1). Regular porous media in a form of staggered or inline arrangement of cylinders have been widely used in the literature to avoid the individuality of results that occurs if a specific geometry of porous media is studied. For example, such a regular geometry of porous medium has been used to establish the universal regime limits (Khalifa et al. 2020), to close the macroscopic turbulence equations in porous media (Pedras and de Lemos 2001; Nakayama and Kuwahara 1999), to obtain the convective heat transfer coefficient in porous media (Saito and de Lemos 2006) and to obtain a correlation for the pressure drop across the porous media (Kuwahara et al. 2006). Moreover, the results obtained for a thin porous media modeled as a confined periodic array of circular cylinders have been successfully applied in different applications of thin porous media such as to study capillary filtration through the clefts between endothelial cells in continuous capillaries (Tsay and Weinbaum 1991) and to model flow through a fracture partially occupied by particles or trapped droplets or bubbles (Chauvet et al. 2012). Therefore, the staggered arrangement of circular cylinder confined between two walls is selected as a model for thin porous media in the present study. Besides, in order to investigate the effects of any possible irregularity in the pore structure of porous media, a random configuration of circular cylinders is considered here. The dimensions of the unit cell are \(2H \times H \times H\) in both cases shown in Fig. 1. The diameter of cylinders is varied to change the porosity of the porous medium, according to:

$$\phi = 1 - \left( \frac{D}{H} \right)^{2} .$$
(2)

The simulations are carried out for three porosities of the staggered arrangement shown in Fig. 1a, \(\phi = 0.5,{ 0}{\text{.65}}\) and \(0.8(s)\), and one porosity for the random distribution in Fig. 1b, \(\phi = 0.8(r)\). The (s) and (r) in \(\phi = 0.8\) denote the staggered and random arrangements, respectively. The distance between walls that confine the porous media is set to \(t \approx H\) characterizing the porous medium considered here as PTMP according to the categorization presented in Fabricius et al. (2016). However, instead of comparing the results obtained in the geometry of Fig. 1 with their two-dimensional (2D) counterparts, as done for Stokesian flow in Fabricius et al. (2016), the comparison is made with those obtained in a similar geometry with symmetry boundaries instead of walls (Semi-3 dimensional). This is because the transitional and turbulent flows are three dimensional in nature, and using a two-dimensional unit cell would deteriorate the results. Semi-3-dimensional geometries are extensively used in the porous media literature to study various properties of a porous medium such as pressure drop correlations or dispersion in the range of \({\text{Re}}_{D}\) above the Darcy limit (Kuwahara et al. 2006). Considering a thin porous material, another Re definition is of interest being based on the distance between the confining walls, t, as:

$${\text{Re}}_{t} = \frac{{\rho U_{D} t}}{\mu }.$$
(3)

Note that in this study, t is varied to study its effects on the macroscopic parameters in a thin porous material. For random geometry, the porosity cannot be obtained using Eq. (2) and is directly derived from the ratio of fluids volume to the total volume of unit cell in the geometry file.

3 Governing Equations

3.1 Filtered Equations

For the LES study, the filtered continuity, momentum and energy equations are obtained by applying a filtering operation to the Navier–Stokes equations, according to:

$$\frac{{\partial \overline{U}_{j} }}{{\partial x_{j} }} = 0$$
(4)
$$\frac{{\partial \overline{U}_{j} }}{\partial t} + \frac{{\partial \overline{U}_{i} \overline{U}_{j} }}{{\partial x_{i} }} = - \frac{1}{\rho }\frac{{\partial \overline{p}}}{{\partial x_{j} }} + \upsilon \frac{{\partial^{2} \overline{U}_{j} }}{{\partial x_{i} \partial x_{i} }} - \frac{{\partial \tau_{ij}^{r} }}{{\partial x_{i} }}$$
(5)

where \(\tau_{ij}^{r}\) is the subgrid-scale (SGS) stress. These subgrid phenomena include the small-scale effects on the momentum and heat transfer being defined as:

$$\tau_{ij}^{r} = \overline{{U_{i} U_{j} }} - \overline{U}_{i} \overline{U}_{j}$$
(6)

The anisotropic part of the SGS stresses are modeled as:

$$\tau_{ij}^{r} = - 2\upsilon_{sgs} \overline{S}_{ij}$$
(7)

where \(\overline{S}_{ij} = \frac{1}{2}\left( {\frac{{\partial \overline{{U_{i} }} }}{{\partial x_{j} }} + \frac{{\partial \overline{{U_{j} }} }}{{\partial x_{i} }}} \right)\) and \(\upsilon_{sgs}\) are the filtered rate of the strain tensor and the SGS eddy viscosity, respectively. The SGS stress is modeled via the standard k-equation SGS model, as:

$$\frac{{\partial k_{sgs} }}{\partial t} + \frac{{\partial \overline{U}_{i} k_{sgs} }}{{\partial x_{i} }} = \frac{\partial }{{\partial x_{i} }}\left( {\left( {\upsilon + \frac{{\upsilon_{sgs} }}{{\sigma_{k} }}} \right)\frac{{\partial k_{sgs} }}{{\partial x_{i} }}} \right) + 2\upsilon_{sgs} \left| {\overline{S}_{ij} } \right|^{2} - \varepsilon_{sgs}$$
(8)
$$\upsilon_{sgs} = C_{k} \Delta \sqrt {k_{sgs} }$$
(9)

where \(k_{sgs}\) is the SGS turbulent kinetic energy, \(\varepsilon_{sgs} = C_{\varepsilon } k_{sgs}^{3/2} / \Delta\) is its dissipation rate, and the second term in the RHS of Eq. (8) is the production of SGS turbulent kinetic energy. The symbol \(\Delta\) is the filter cutoff width which corresponds to the limit below which the turbulence scales are modeled instead of being resolved. According to Versteeg and Malalasekera (2007), it is usually selected to be of the same order of the grid size. Considering \(\Delta_{x}\), \(\Delta_{y}\) and \(\Delta_{z}\) as the length, width and height of the grid cells, \(\Delta\) is calculated as the cubic root of the cell volume as:

$$\Delta = \sqrt[3]{{\Delta_{x} \Delta_{y} \Delta_{z} }}$$
(10)

Following (Pope 2001; Ansys 2012), the constants in Eqs. (8) and (9) are defined as:

$$C_{k} = 0.094, \, C_{\varepsilon } = 1.048, \, \sigma_{k} = 1.0$$
(11)

3.2 Volume-Averaged Hydrodynamic Dispersion and Reynolds Stress

The volume-averaged hydrodynamic dispersion and the Reynolds stresses that appear by volume–time averaging of the convection term in the momentum equation are presented as follows (Pedras and de Lemos 2001):

$$\begin{gathered} \frac{\partial }{{\partial x_{i} }}\left( {\phi \overline{\overline{{\left\langle {U_{i} U_{j} } \right\rangle }}} } \right) = \frac{\partial }{{\partial x_{i} }}\left( {\phi \overline{\overline{{\left\langle {\left( {\left\langle {U_{i} } \right\rangle + \tilde{U}_{i} } \right)\left( {\left\langle {U_{j} } \right\rangle + \tilde{U}_{j} } \right)} \right\rangle }}} } \right) \hfill \\ { = }\frac{\partial }{{\partial x_{i} }}\left( {\phi \overline{\overline{{\left( {\left\langle {U_{i} } \right\rangle \left\langle {U_{j} } \right\rangle + \left\langle {\tilde{U}_{i} \tilde{U}_{j} } \right\rangle } \right)}}} } \right) = \frac{\partial }{{\partial x_{i} }}\left( {\phi \overline{\overline{{\left[ {\left( {\left\langle {\overline{\overline{U}}_{i} } \right\rangle + \left\langle {U^{\prime}_{i} } \right\rangle } \right)\left( {\left\langle {\overline{\overline{U}}_{j} } \right\rangle + \left\langle {U^{\prime}_{j} } \right\rangle } \right) + \left\langle {\tilde{U}_{i} \tilde{U}_{j} } \right\rangle } \right]}}} } \right) \hfill \\ { = }\frac{\partial }{{\partial x_{i} }}\left( {\phi \left( {\left\langle {\overline{\overline{U}}_{i} } \right\rangle \left\langle {\overline{\overline{U}}_{j} } \right\rangle + \overline{\overline{{\left\langle {U^{\prime}_{i} } \right\rangle \left\langle {U^{\prime}_{j} } \right\rangle }}} + \overline{\overline{{\left\langle {\tilde{U}_{i} \tilde{U}_{j} } \right\rangle }}} } \right)} \right) \hfill \\ { = }\frac{\partial }{{\partial x_{i} }}\left( {\phi \left( {\underbrace {{\left\langle {\overline{\overline{U}}_{i} } \right\rangle \left\langle {\overline{\overline{U}}_{j} } \right\rangle }}_{I} + \underbrace {{\overline{\overline{{\left\langle {U^{\prime}_{i} } \right\rangle \left\langle {U^{\prime}_{j} } \right\rangle }}} }}_{II} + \underbrace {{\left\langle {\overline{\overline{{\tilde{U}}}}_{i} \overline{\overline{{\tilde{U}}}}_{j} } \right\rangle }}_{III} + \underbrace {{\left\langle {\overline{\overline{{\tilde{U}^{\prime}_{i} \tilde{U}^{\prime}_{j} }}} } \right\rangle }}_{IV}} \right)} \right). \hfill \\ \end{gathered}$$
(12)

The term I represents the convection of macroscopic mean velocity. The term III in Eq. (12) denotes the dispersion arising from the space fluctuation of time mean velocity within the pores which is presented in both laminar and turbulent flows. The averaged Reynolds stress term is composed of terms II and term IV in Eq. (12), which are originated from the macroscopic velocity fluctuations and time and space fluctuations of velocity within the pores, respectively.

4 Numerical Simulation and Boundary Condition

The governing equations are solved in the open-source platform OpenFOAM with the pimpleFoam solver. This algorithm is based on a blend of the transient SIMPLE and PISO algorithms. The grid structure is generated in two steps: blockMesh (BM) and snappyHexMesh (SHM). The BM is used to generate a background mesh in the channel in which the mesh is refined close to the channel walls. Then, SHM uses the cylinders geometry file generated in OpenSCAD and the background mesh to produce the final grid structure. Advection terms are discretized using an interpolation scheme based on the upwind and central differencing methods, and the second-order backward differencing scheme is used to discretize the unsteady terms. The Courant number in the simulations is between 0 and 1 based on the time steps. Periodic boundary conditions are applied on the faces in the x- and z-directions. No-slip boundary conditions are considered at the top and bottom walls of the channel for the simulations of thin porous media, while symmetry boundary condition is used on the upper and lower boundaries in the y-direction for semi-3D cases. The macroscopic values are calculated and monitored during specific intervals of the simulations. The results are considered to be macroscopically fully developed when the difference between these values in two subsequent intervals becomes negligible.

5 Grid Independency, Verification and Validation

The grid independency is carried out on a thin porous medium with \(t = H\) and \(\phi = 0.8(s)\) for axial fluctuating velocity along the z-direction at \(Y = 0.5\) and \(X = 1\) and shown for three different cell numbers for ReD = 4000 in Fig. 2. A maximum variation of less than 5 percent is observed in the fluctuation velocity with an increase in the cell numbers from 390,000 to 750,000 in Fig. 2 suggesting the sufficiency of the grid resolution used in the present study. Since the macroscopic values of flow variables are of interest in the present study, the grid independence study is also carried out for volume-averaged turbulent kinetic energy and pressure drop. The variation of these variables is small between the two finest meshes, see Table 2. Therefore, the finest mesh number is used from now on since it also satisfies the condition of \(y^{ + } \approx 1\).

Fig. 2
figure 2

Axial fluctuating velocity along the z-direction for \(\phi = 0.8(S)\) and \(Y = 0.5\) and \(X = 1\)

Table 2 Macroscopic values for \(\phi = 0.8(S)\) in thin porous media for different cell numbers

The verification of the present solver is carried out by comparing the xy Reynolds stress component with that obtained in Moser et al. (1999) through direct numerical simulation of turbulent flow in a channel, see Fig. 3. The difference between the present results and DNS data is calculated to be 4 percent at its maximum in \({{2y} \mathord{\left/ {\vphantom {{2y} H}} \right. \kern-\nulldelimiterspace} H} = 0.1\) confirming the validity of the present solver. Further validation of the present simulations is carried out through a comparison of the friction factor calculated in a channel with the experimental and DNS data reported in the literature (Abe et al. 2001; Kuroda et al. 1989; Kim et al. 1987; Kim et al. 1990; Dean 1978), see Fig. 4. As can be seen, all numerical data follow closely the Blasius correlation including the present results for which a maximum difference of 1.5% between two results is observed at the maximum Rem investigated. In order to measure the accuracy of the simulations concerning the flow in porous media, the pressure drop calculated for \(\phi = 0.8(s)\) geometry is compared with those obtained by Kuwahara et al. (1998) and Pedras and de Lemos (2001) for turbulent flow through a 2D geometry of staggered arrangement of circular cylinders in Fig. 5. Standard and low-Re k-epsilon models are used to take into account the turbulence effects in the simulations in Kuwahara et al. (1998) and Pedras and de Lemos (2001), respectively. Apart from the trend of the present results which successfully follows that of Pedras and de Lemos (2001), a quantitative comparison is possible between the non-dimensional pressure drop interpolated in \({\text{Re}}_{H} \approx 60\) and that of Pedras and de Lemos (2001). A difference of 5 percent is obtained which further verifies the accuracy of the present results.

Fig. 3
figure 3

Comparison of the Reynolds stresses with the results from Moser et al. (1999)

Fig. 4
figure 4

Comparison of the friction factor with the experimental correlations and DNS results

Fig. 5
figure 5

Comparison of pressure drop with those of Kuwahara et al. (1998) and Pedras and de Lemos (2001) for \(\phi = 0.8(S)\)

6 Results and Discussion

To start with, an investigation of pore-scale flow and turbulence structures is carried out. Then, focus is set on the macroscopic variables such as pressure drop, hydrodynamic dispersion and Reynolds stresses. Finally, the effect of channel thickness, t, is scrutinized. In the following, the TPM and semi-3D models of porous media studied here are referred to as confined and non-confined porous media, respectively.

6.1 Flow Field and Turbulence Structures

There is a noticeable difference between the streamlines for confined and non-confined porous media (\(\phi = 0.8(s)\)) for \({\text{Re}}_{D} = 400,\) while the shape of streamlines is more similar for \({\text{Re}}_{D} = 4000\), see Fig. 6. The streamlines, shown in Fig. 6, are extracted based on the time mean velocity within the unit cell. As can be seen, tornado-like vortexes, based on the terminology used in Kirkil and Constantinescu (2015), are formed in the wake region of the cylinders for all ReD. While the vortex size decreases significantly with an increase in Re from 400 to 4000 for the confined case, this variation is less for the non-confined case. The presence of larger vortex downstream the confined cylinders for \({\text{Re}}_{D} = 400\) may be attributed to the larger resistance to flow where the cylinders meet the confining walls. The presence of these walls also reduces the turbulent intensity in this area, allowing the vortices downstream of cylinders to become larger in the confined case as compared to the case without bounding walls. The reduction in the size of these vortices with an increase in Re to \({\text{Re}}_{D} = 4000\) can be attributed to the enhancement of turbulence effects at higher Re which enhances the mixing resulting in a smaller vortex downstream of cylinders. However, it can still be observed that the size of the vortices for the confined porous media is larger than those in the non-confined case. Such a behavior was previously reported for a confined porous media in Khayamyan et al. (2017a2017b) through their experiments carried out on transition and turbulent flow in a random packed bed of spheres. They observed that the non-uniformity in the mean velocity profile within the pores increases at first with an increase in Re, while the mean velocity distribution within the pores becomes more uniform at higher Re. Khayamyan et al. (2017a, 2017b) also reported that such variations may also be reflected in macroscopic properties of flow such as volume-averaged hydrodynamic dispersion and pressure drop.

Fig. 6
figure 6

Comparison of streamlines, colored by mean velocity magnitude (m/s), between confined and non-confined porous media in \(\phi = 0.8(s)\) (The flow is in + x-direction)

A similar comparison between the streamlines obtained for confined and non-confined porous media is carried out in Fig. 7 for the random geometry shown in Fig. 1b for \(\phi = 0.8(r)\) where the diameter of the cylinders is half of that in the staggered arrangement. As can be seen in Fig. 7, the streamlines for the two porous media are similar regardless of ReD. Therefore, it can be concluded that the dimension of solid particles is also important in addition to Re for the development of flow structures within the thin porous media, as also reported in Fabricius et al. (2016) for Stokesian flow. In this case, the similarity between the two cases of confined and non-confined cases for all ReD of the random geometry can be attributed to the relatively smaller diameter of circular cylinder resulting in weaker three-dimensional structure formation at the base of the cylinders and, therefore, smaller vortices downstream of the cylinders.

Fig. 7
figure 7

Comparison of streamlines, colored by mean velocity magnitude (m/s), between confined and non-confined porous media in \(\phi = 0.8(r)\) (The flow is in + x-direction)

Further investigations of the effect of the confining walls for the staggered arrangement are carried out for \(\phi = 0.5\) and \(\phi = 0.65\). For the sake of conciseness, only the streamlines for \(\phi = 0.5\) are illustrated in Fig. 8 where the cylinders are made transparent in order to clearly observe the streamlines behind them. In this case, the sizes of the recirculation zones for the confined porous media are only slightly larger than those for the non-confined case. This behavior can be attributed to the smaller void space between particles for \(\phi = 0.5\) as compared to \(\phi = 0.8(s)\). In this case, the recirculation zones downstream of the cylinders in the confined porous media do not have enough space to grow, and therefore, they are slightly smaller than those observed in non-confined porous media.

Fig. 8
figure 8

Comparison of streamlines, colored by mean velocity magnitude (m/s), between confined and non-confined porous media in \(\phi = 0.5\) (The flow is in + x-direction)

In order to study the difference between turbulence structures in the confined and non-confined porous media, the iso-surfaces of non-dimensional Q-criterion are presented in Fig. 9. The Q-criterion is the second invariance of the velocity gradient tensor whose positive values denote the vortical structures within the flow field. Considering \(S\) as the symmetric part of velocity gradient tensor known as the rate of strain and \(\Omega\) as the antisymmetric part or the vorticity tensor, Q-criterion is defined as \(Q = 0.5\left( {\left\| \Omega \right\|^{2} - \left\| S \right\|^{2} } \right)\). Accordingly, the Q-criterion provides information about the balance between vorticity magnitude and shear strain rate. Positive values of Q-criterion refer to the areas in which the vorticity magnitude is greater than that of rate of strain (Zhan et al. 2019). The iso-surfaces in Fig. 9 are colored with y-velocity to show the motion of these structures. The presence of side walls in the case of confined porous media causes horseshoe-like vortices to be formed around the base of the cylinders though they are harder to observe at higher Re due to higher inertia effects causing them to break. These vortices are denoted by filled arrow in Fig. 9. Such vortical structures cannot be observed for the non-confined porous media. For \({\text{Re}}_{D} = 400\), the vortical structures associated with tornado-like vortices downstream of the cylinders can be observed for both the confined and non-confined porous media as indicated by the opened arrows. These vortices correspond to the streamlines observed downstream of the cylinders in Fig. 6 for \(\phi = 0.8(s)\) and \({\text{Re}}_{D} = 400\). Such structures would significantly enhance the heat or mass transport from the confining walls. These tornado-like vortical structures can also be observed at \({\text{Re}}_{D} = 4000\) although their size becomes smaller due to higher inertia effects and turbulence mixing. Also, due to turbulence, the legs of the horseshoe-like vortex are broken, and they are harder to observe. Notice that in order to better show the tornado-like vortices, just half of the iso-surfaces of the Q-criterion are shown in Fig. 9 for \({\text{Re}}_{D} = 4000\).

Fig. 9
figure 9

Comparison of iso-surfaces of Q-criterion streamlines, colored by y-velocity (m/s), between confined and non-confined porous media in \(\phi = 0.8(s)\) (The flow is in + x-direction)

Also, for the random geometry, horseshoe-like vortices can be observed for the confined porous medium, while such structures cannot be observed for non-confined arrangement of cylinders, see Fig. 10. The vortical structures in the wake region downstream of the cylinders are also present and elongated in the y-direction parallel to the axis of the cylinders. A detailed examination of the turbulence structures in Figs. 9 and 10 shows the presence of hairpin eddies at the interface of the high-speed flow around the cylinders and the low-speed flow in the wake region. These structures are indicated with a diamond head arrow in Fig. 10. A qualitative comparison between the vortical structures presented in Figs. 9 and 10 can be made to investigate the effect of pore geometry on the form of the turbulent structures. It is seen that the turbulence structures in confined and non-confined cases are similar for the random and staggered cases. A quantitative comparison between the cases studied will be presented in the following sections.

Fig. 10
figure 10

Comparison of iso-surfaces of Q-criterion streamlines, colored by y-velocity (m/s), between confined and non-confined porous media in \(\phi = 0.8(s)\) (The flow is in + x-direction)

6.2 Pressure Drop

As discussed above, Fabricius et al. (2016) showed that the confining walls affect the pressure drop for PTMP when it is compared to the pressure drop obtained for HTPM in the Darcy regime. For the current case, i.e., the non-Darcian regimes, the difference between non-dimensional pressure drops calculated in the confined and non-confined porous media for \({\text{Re}}_{D} = 40\) is larger than other ReD studied and this difference decreases with an increase in ReD in all cases studied, see Fig. 11. For \({\text{Re}}_{D} = 40\) in Fig. 11, it is found that the difference between pressure drops calculated for the confined and non-confined cases increases with an increase in porosity for staggered geometry, i.e., 19%, 27%, 42% and 18.5% for \(\phi = 0.5\), \(\phi = 0.65\), \(\phi = 0.8(s)\) and \(\phi = 0.8(r)\), respectively. This behavior can be due to a higher ratio of viscous drag to the total drag across the porous media at lower ReD. As a result, the surface friction over the confining wall plays a more important role at lower ReD, while at higher Re the form drag or Forchheimer term in Darcy–Forchheimer correlation is more significant. Another phenomenon causing the variation of pressure drop with ReD in Fig. 11 is the difference between flow structures in the confined and non-confined cases, as discussed earlier. The formation of the tornado-like structures in the wake region of the cylinders in the confined porous media with a size larger than the vortex structures observed in the case of non-confined porous medium as well as the presence of horseshoe vortex in the confined cases leads to an increase in pressure drop at lower ReD and higher porosities. Similarly, the smaller difference in the case of \(\phi = 0.8(r)\) in comparison with \(\phi = 0.8(s)\) for \({\text{Re}}_{\rm D} = 40\) can be explained. A smaller cylinder diameter for \(\phi = 0.8(r)\) compared to that in \(\phi = 0.8(s)\) results in a smaller difference between the wake regions downstream of cylinders leading to a smaller difference between confined and non-confined porous media. With an increase in ReD and the enhancement of mixing due to turbulence effects, the difference between the wake regions deteriorates and the pressure drops for confined and non-confined cases become closer to each other.

Fig. 11
figure 11

Non-dimensional pressure drop as a function of porosity and ReD for confined (thin) and non-confined porous media a \(\phi = 0.5,\) b \(\phi = 0.65,\) c \(\phi = 0.8(s)\) and \(\phi = 0.8(r)\)

6.3 Hydrodynamic Dispersion

Figure 12 shows the variation of the normal components of the hydrodynamic dispersion (term III in Eq. (12)) non-dimensionalized as follows:

$$\frac{{\phi \left\langle {\overline{\overline{{\tilde{U}}}}_{i} \overline{\overline{{\tilde{U}}}}_{j} } \right\rangle }}{{U_{\rm D} U_{\rm D} }}$$
(13)
Fig. 12
figure 12

Components of non-dimensional volume-averaged hydrodynamic dispersion

for different porosities and ReD.

It can be observed that the normal component in the x-direction (the main flow direction), UU, has the largest values independent of the setup. The WW component of the hydrodynamic dispersion is also in the same order of magnitude as UU though the values are lower. The relatively large values of UU and WW as compared to VV are because of the tortuosity of flow path in x- and z-directions due to the presence of cylinders. However, there is no obstacle when flow moves along the y-direction in which the tortuosity is close to 1. In spite of the small values of the VV component, some interesting observations can be made on the difference between them in confined and non-confined cases. The higher values of the VV dispersion in the confined cases can be due to the horseshoe vortices at the bases of cylinders as well as larger tornado-like vortices in the wake of cylinders in these types of porous media, as discussed in Sect. 6.1. In the case of non-confined porous media, these structures are observed to be absent or weaker. The shear components of non-dimensional hydrodynamic dispersion in Eq. (13) are also calculated, though they are not shown here, and their values are found to be negligible. Looking at Fig. 12a, c, a similar trend can be observed for the variation of hydrodynamic dispersion with ReD in \({\text{Re}}_{\rm D} > 400\) for all porosities, while the trend observed for \({\text{Re}}_{\rm D} < 400\) is different in different porosities. This can be attributed to the transition from laminar to turbulent regime in the pores of porous media in the range of \({\text{Re}}_{\rm D} < 400\) which agrees well with the \({\text{Re}}_{\rm D} = 300\) reported in Dybbs and Edwards (1984) for transition to fully turbulent regime in porous media.

In order to investigate the variation of the normal components with porosity and ReD, their values are gathered in a non-dimensional dispersive kinetic energy (\(\kappa = 0.5(\langle {\overline{\overline{{\tilde{U}}}} \overline{\overline{{\tilde{U}}}} } \rangle + \langle {\overline{\overline{{\tilde{V}}}} _{i} \overline{\overline{{\tilde{V}}}} _{j} } \rangle + \langle {\overline{\overline{{\tilde{W}}}} \overline{\overline{{\tilde{W}}}} } \rangle /U_{D}^{2} )\)). As shown in Fig. 13, its values slightly increase with an increase in ReD followed by a decrease with a further increase in Reynolds number. The variation in \(\kappa\) can be explained by the turbulence effects within the porous media, as reported in the experimental observation of Khayamyan et al. (2017a2017b). The increase in dispersive kinetic energy with an increase in \({\text{Re}}_{\rm D} = 40\) is due to higher velocity (inertia) and larger vortices at higher ReD which causes the fluid particles to be more dispersed within the pores. However, with a further increase in \({\text{Re}}_{\rm D}\), the turbulence structures start to grow which enhances the mixing within the pores causing the dispersion components to be slightly damped. It is also obvious from Fig. 13 that the dispersive kinetic energy is slightly larger for the thin porous media as compared to the unconfined one which can be attributed to the presence of confining walls.

Fig. 13
figure 13

Dispersion kinetic energy as a function of porosity and ReD

6.4 Reynolds Stress

The volume-averaged Reynolds stress in porous media is expressed as the sum of terms II and IV in Eq. (12), as:

$$\left\langle {\overline{\overline{{U^{\prime}_{i} U^{\prime}_{j} }}} } \right\rangle = \underbrace {{\overline{\overline{{\left\langle {U^{\prime}_{i} } \right\rangle \left\langle {U^{\prime}_{j} } \right\rangle }}} }}_{II} + \underbrace {{\left\langle {\overline{\overline{{\tilde{U}^{\prime}_{i} \tilde{U}^{\prime}}}}_{j} } \right\rangle }}_{IV}$$
(14)

where term II refers to Reynolds stress due to the fluctuation of macroscopic velocity and term IV originates from dispersion due to both time and space fluctuations of velocity within the pores. It should be noted that term II on the RHS of Eq. (14) is found to be negligible based on the volume-averaged values of fluctuating velocities over the unit cell, \(\left\langle {U^{\prime}_{i} } \right\rangle \approx 0\) for \(i = 1,2,3\), and therefore, only term IV remains. Anyway, the term on the LHS of Eq. (14) is calculated from the pore-scale simulations of the present study and reported in the following non-dimensional form in Fig. 14:

$$\frac{{\phi \left\langle {\overline{\overline{{U^{\prime}_{i} U^{\prime}_{j} }}} } \right\rangle }}{{U_{\rm D} U_{\rm D} }}$$
(15)
Fig. 14
figure 14

Components of non-dimensional volume-averaged Reynolds stress

It can be observed that unlike the turbulence in the regular arrangement of spheres studied in Jouybari and Lundström (2020) and Jin and Kuznetsov (2017), the turbulence in this geometry is not isotropic and the normal component of Reynolds stress in the y-direction (along the length of the cylinders) is significantly lower than those in x- and z-directions. Also, it is interesting to observe that the level of normal Reynolds stresses in the non-confined porous media is higher than that in confined porous medium, see Fig. 14. This can be attributed to the damping effect of confining walls which seems to overwhelm the turbulence production around the cylinders for the confined cases. Also, the Reynolds stress trends presented in Fig. 14 support the previous conclusion on the transition to the fully turbulent regime in \({\text{Re}}_{\rm D} < 400\). Looking at Fig. 14, it can be observed that the level of normal components of Reynolds stress remains nearly constant for \({\text{Re}}_{\rm D} > 400,\) while it decreases almost in all porosities with a decrease in ReD for \({\text{Re}}_{\rm D} < 400\) showing a variation to fully turbulent regime in this range of ReD. A comparison of vortex structures between \({\text{Re}}_{\rm D} = 100\), see Fig. 20 in Appendix A, with those obtained for \({\text{Re}}_{\rm D} = 400\) in Fig. 9 confirms the variation in flow regime.

Similar to that found for dispersion, the level of shear components of Reynolds stress averaged over the unit cells is found to be too small. Figure 15 illustrates the non-dimensional volume-averaged turbulent kinetic energy (TKE), \(k^{*} = 0.5( {\langle {\overline{\overline{{U^{\prime}U^{\prime}}}} } \rangle + \langle {\overline{\overline{{V^{\prime}V^{\prime}}}} } \rangle + \langle {\overline{\overline{{W^{\prime}W^{\prime}}}} } \rangle } )/U_{D}^{2}\). It can be observed that the dissipation of turbulence caused by confining walls is more significant at lower ReD and higher porosities. The difference between \(k^{*}\) in confined and non-confined cases decreases with a decrease in porosity and an increase in ReD. Also, the TKEs calculated in the random geometry are closer to each other almost in all ReD showing the lower effects of confining walls in this geometry.

Fig. 15
figure 15

Turbulent kinetic energy as a function of porosity and ReD

6.5 Effect of Channel Thickness, t

In this section, a case with a half-channel thickness of that considered in previous simulations is considered for the case of \(\phi = 0.8(s)\). Considering a fixed \({\text{Re}}_{\rm D}\), the \({\text{Re}}_{t,0.5H}\) in the following cases are half of those in previous ones, \({\text{Re}}_{t,H}\). As can be observed in Fig. 16, the trend observed for streamlines for \({\text{Re}}_{t,0.5H}\) and \(\phi = 0.8(s)\) is similar to that observed in Fig. 6 for \({\text{Re}}_{t,H}\). However, it seems that the difference between the size of recirculating regions behind the cylinders is larger for the thinner case. Also, a comparison between the maximum values of the velocity for \({\text{Re}}_{\rm D} = 400\) in Figs. 6 and 16 between the confined and non-confined cases reveals that the maximum value of the velocity increases more in the confined porous media for the thinner case. This can be due to the larger wake region behind the cylinders in the confined case of Fig. 16 due to which the flow is more squeezed when it passes through the pore space. Figure 17 illustrates the iso-surfaces of the Q-criterion extracted for the \(\phi = 0.8(s)\) and \({\text{Re}}_{t,0.5H}\). The horseshoe vortex and the tornado-like vortices for \({\text{Re}}_{t,0.5H}\) are denoted by the filled and open arrows. Analogous to that observed for the shape of streamlines between confined and non-confined cases for \({\text{Re}}_{\rm D} = 400\) (see Fig. 16), a significant change in the form of vortical structures can be seen for \({\text{Re}}_{t,0.5H}\) and \({\text{Re}}_{\rm D} = 400\) in Fig. 17. While the Q-criterion shows the presence of elongated vortical structures in the confined case, the vortical structures in the non-confined case are broken into smaller ones for \({\text{Re}}_{\rm D} = 400\).

Fig. 16
figure 16

Comparison of streamlines, colored by mean velocity magnitude (m/s), between confined and non-confined porous media in \(\phi = 0.8(s)\) and for \({\text{Re}}_{t,0.5H}\) (The flow is in + x-direction)

Fig. 17
figure 17

Comparison of Q-criterion iso-surfaces, colored by y-velocity (m/s) in \(\phi = 0.8(s)\) and for \({\text{Re}}_{t,0.5H}\) (The flow is in + x-direction)

In Fig. 18, the effect of variation in the channel thickness on the macroscopic values of pressure drop in confined and non-confined porous media is presented. It is found that the non-dimensional pressure drop increases with a decrease in channel thickness in confined porous media, while it remains nearly constant for the non-confined case. This result agrees well with the observation made in Fig. 17 which shows that the differences between vortical structures are more significant in smaller channel thickness. Meanwhile, the results of the pressure drop in non-confined cases suggest that the macroscopic pressure drop in this type of porous media is not dependent upon the unit cell thickness despite the presence of three-dimensional effects in the flow field and a smaller unit cell can be incorporated to perform the numerical simulations. Similar to the trend observed in Fig. 17, the difference between the confined and non-confined macroscopic non-dimensional pressure drop decreases with an increase in ReD due to the dominance of turbulent effects at higher Re numbers causing the flow pattern in two confined and non-confined cases to become more similar to each other, as can be observed in Fig. 17.

Fig. 18
figure 18

Effect of channel thickness on the pressure drop for confined (thin) and non-confined porous media

In order to scrutinize the effect of variation in the channel thickness on the hydrodynamic dispersion and Reynolds stress within confined and non-confined porous media, the dispersive and turbulent kinetic energies are plotted in Fig. 19. As can be seen, the turbulent kinetic energy increases for all ReD with a decrease in channel thickness in the non-confined porous media, while, in confined cases, the major increase is observed for \({\text{Re}}_{\rm D} > 400\). With a reduction in channel thickness, the size of turbulent structures along the cylinder length is limited by the symmetry boundary conditions which may interrupt the energy cascade processes between large coherent structures to the small structures. As a result, the values of TKE in the thinner channel are higher than those in the thicker channel also for non-confined porous media. It seems that such an effect is observable for \({\text{Re}}_{\rm D} > 400\) in the confined cases and for Re numbers smaller than this value, the turbulence effects are not significant within the flow field which can be due to the dissipation caused by the presence of confining walls. Meanwhile, it can be seen in Fig. 19a that as the turbulence effects become important, the turbulent kinetic energy values become closer to each other in confined and non-confined cases in both \({\text{Re}}_{t,H}\) and \({\text{Re}}_{t,0.5H}\). The effect of turbulence on the dispersion can be observed in Fig. 19b by computing the dispersive kinetic energy values. It can be observed that the dispersive kinetic energy values in \({\text{Re}}_{t,0.5H}\) for non-confined case increase for \({\text{Re}}_{\rm D} > 400\) and become closer to that of confined porous media which can be considered as a result of enhancing turbulence effects, as observed in Fig. 19a. Also, the higher turbulence effects in non-confined cases (see Fig. 19a) causes the dispersive kinetic energy in non-confined cases to be smaller than those in confined porous media (see Fig. 19b). As discussed above, the turbulence effects tend to increase the uniformity of flow within the pores through enhancing the mixing which results in lower dispersion effects.

Fig. 19
figure 19

Comparison of a turbulent kinetic energy, b dispersive kinetic energy for the confined and non-confined porous media in \({\text{Re}}_{t,H}\) and \({\text{Re}}_{t,0.5H}\)

7 Conclusion

Numerical simulations of fluid flow within porous media are carried out in a wide range of flow regimes from nonlinear laminar flow to the turbulent regime for different porous media confined with solid walls or not. It is observed that the streamline patterns and vortical structures, specifically at the bases of cylinders, in the confined porous media are different from that in non-confined porous media and this difference becomes smaller as the ReD is increased. This can be attributed to higher turbulence effects resulting in an enhancement in mixing within the pore at higher ReD. A larger pressure drop is observed for the confined cases as compared to the unconfined ones due to the presence of confining walls and a different pore-scale flow. This difference decreases with an increase in ReD because of a decrease in the contribution from the viscous drag to the total drag over the confining walls and the higher level of similarity between pore-scale flow at higher ReD. An investigation of dispersive kinetic energy revealed its higher values for the confined cases compared to those in non-confined porous media. It is also observed that the averaged values of dispersive kinetic energy initially increase up to \({\text{Re}}_{\rm D} = 100\). With a further increase in ReD, the level of averaged dispersive kinetic energy decreases which can be attributed to higher turbulence mixing at higher ReD. This result is in agreement with the experimental observation of Khayamyan et al. (2017a2017b) for packed beds of spheres. Also, the results show a higher level of turbulent kinetic energy in non-confined porous media which may be attributed to the presence of the confining walls in confined porous media and their dissipation effects on turbulence. As a result of the higher turbulent kinetic energy in the non-confined cases, the flow parameters become more uniform resulting in lower dispersive kinetic energy in these cases. Finally, a discussion is provided on the effect of channel thickness, t, on the flow and turbulence parameters by carrying out the numerical simulations in a channel with a thickness of half of that considered in the previous simulations. It is observed that the difference between the streamlines and vortical structures in two confined and non-confined media increases which results in a higher difference between the pressure drop calculated in a thinner channel.