Steady-State Transitions in Ordered Porous Media

Previously performed experiments on flow through an ordered porous media cell with tomographic particle image velocimetry reveal a complex three-dimensional steady-state flow pattern. This flow pattern emerge in the region where inertial structures have been previously reported for a wide range of packings. The onset of these steady-state inertial flow structures is here scrutinized for three different types of packing using a finite difference method. It is concluded that the onset of the flow structure coincides with a symmetry break in the flow field and discontinuities in the pressure drop, volume averaged body forces and heat transfer. A quantity for identifying the transition is proposed, namely the pressure integral across the solid surfaces. It is also shown that the transition can both increase and decrease the heat transfer dependent on the actual geometry of the porous medium.


Introduction
Porous media flow takes place in many natural and industrial processes including flow through soils, aquifers, oil and gas reservoirs (Bear 1972), textiles (filters, for instance), biological tissues and plants. Other examples are flow through fuel cells (Farzaneh et al. 2021), cooling of batteries Moosavi et al. (2021), flow during composites manufacturing (Lu et al. 2017;Tan and Pillai 2012a, b, c),industrial wicks (Zarandi et al. 2022) flow during paper-making, drying of iron ore pellets (Burstrom et al. 2018), internal erosion in embankment dams (Frishfelds et al. 2011), heat exchangers (Odabaee and Hooman 2012) and in-tissue drug delivery. This long list of applications (although it is not by far complete) points out that it is of utmost importance to increase the understanding of different flow structures present in various kinds of porous media. Increased knowledge will enable an optimization of the properties of the porous media as well as the parameters governing the flow. As to Reynolds number the flow in porous media is usually subdivided into four regions, these are called the Darcy region, laminar steady region, laminar unsteady region and turbulent region. The regions are defined by the following flow conditions being observed (Seguin et al. 1998a, b): • Darcy -p, i ∝ u i , • Laminar steady - Where p, i denotes p x i . In addition other categories for flow through porous media can be identified including Newtonian or Non-Newtonian behaviour, single and multi-phase flow, wicking and saturated flow and free or confined flow, e.g. (Xu and Pillai 2017;Larsson et al. 2018;Jouybari and Lundstrom 2021). In this study Newtonian, single-phase, incompressible, low Knudsen number, saturated confined flow is considered.
The different regions as to Reynolds number have been considered in a number of studies with focus on different types of porous media. To exemplify, for turbulent flow through ordered arrays the particle shape and the type of array are crucial for how important different turbulent quantities are (Yang et al. 2014) and the heat flux due to thermal dispersion and turbulence is small as compared to the convective contribution (Jouybari et al. 2020). Plots from the same study also show that the heat flux caused by turbulence is limited to single pores confirming that the size of turbulent structures in porous media is generally limited to the pore scale, leading to the pore-scale prevalence hypothesis (Jin et al. 2015;Uth et al. 2016). Only for a limited range of porosities > 0.93 and Darcy numbers, the turbulence structures can extend to nearby pores (Rao and Jin 2022). Randomly packed beds of spheres is another type of porous medium that has been scrutinized. With focus on turbulence, results in (Patil and Liburdy 2013) show that the characteristics of the turbulence is nearly the same between pores. In (Nguyen et al. 2018) comparisons between velocities and the Reynolds stresses yield an increase in turbulent intensity and the mixing within the pores when the Reynolds number is increased. In (Khayamyan et al. 2017a(Khayamyan et al. , 2017b also flow through randomly packed beds is studied but with focus on the laminar unsteady region. From Stereoscopic Particle Image Velocimetry it is, for instance, shown that the spatial variation in the time averaged absolute velocity increases with particle Reynolds number up to about 410 and then the variation decreases . This is explained with a transition from the laminar unsteady to the turbulent region. A similar observation is reported in Johns et al. (2000) who analyzed inertial transitions in packed beds by using magnetic resonance imaging. The spatial variance of velocity decreases considerably at a local (pore) Reynolds number of about 30. The mechanisms is, however, likely to be different to that in  since the decrease with Reynolds number in Johns et al. (2000) takes place in the laminar steady regime. Scrutinizing the results in , in detail (Figures 9, 10), the spatial variation in velocity actually decreases when the Reynolds number increases between 20 and 40 (before the variation increases up to Reynolds number 410). This is in the laminar steady flow region for this geometry and therefore in-line with the results in Johns et al. (2000) indicating that the epsilon shaped effect is a generic feature.
In the above discussion, a few values of the Reynolds number are mentioned. It should however be stated that the transitional Reynolds numbers are crucially dependent on the geometry studied. This will also be the case for the onset of the flow structures examined in this study. For porous media flow, the different types of Reynolds numbers used in the literature can complicate comparisons to actual Reynolds numbers even more.
This study focuses on a type of inertial flow structure observed in the laminar steady flow region for saturated Newtonian single phase flow through structured porous media. The epsilon-shaped effect studied inhere, was first observed in a previous work by the authors of this article when studying the flow through a porous medium consisting of pillars placed between two parallel horizontal and two parallel vertical plates see Fig. 1  ). The effect is named the epsilon effect due to the epsilon-shaped inertial cores that are formed for this particular geometry. Based on data available it is concluded that the structure is caused by wall interactions and that it remains steady as long as the Reynolds number is kept constant. Numerical simulations later confirmed that a similar, albeit distinct, flow structure could arise independent of wall interactions. In this study, a finite difference method (FDM) solver is utilized and the onset of the effect is scrutinized to yield a greater insight into the underlying mechanisms. The results show that for high-tortuosity porous media, an inertial-viscous transition occurs that is characterized by a breakdown of the high pressure region on the impinging surfaces. The transition manifests as a symmetry break in the flow field with inertial cores that are directed around the obstructing geometry. It is also shown that the transition occurs in packings of spheres.

Theory
The full Navier-Stokes equations on tensor form reads where u i is the velocity vector, is the density, is the kinematic viscosity and F i is some external body force term like gravity on volume density form. For steady-state flows the transient flow field variation is zero, i.e. (1) The epsilon effect as observed experimentally by tomographic PIV measurements, the green planes indicates the position of the slices in the porous cell. Inlet is to the left and outlet to the right with the flow direction along the positive x-axis. The pillars seen are placed between two parallel horizontal and two parallel vertical walls. Further details are available in  From left to right the above terms are referred to as the inertial term, static pressure term, viscous drag term and the source term. Integrating this expression over the fluid domain Ω gives the average contribution of each term These quantities are referred to as spatially averaged terms and is the average force balance in the fluid domain. In this study periodic boundary conditions are applied, this implies that the entire fluid domain consists of closed steady-state streamlines. Therefore the spatially averaged value of the convective term ∫ Ω u j u i,j is zero, and Eq. (3) is reduced to If the pressure gradient term is significantly larger than the viscous term, the flow is said to be inertially limited. Another useful quantity is the volume averaged force magnitudes defined by These quantities give an idea of how much the average local force balance varies between the computational elements.

Reynolds Number
The Reynolds number characterizes the flow-dynamics of the hydrodynamic phase, here defined as Where U int is the interstitial velocity defined as the average streamwise velocity in the fluid, D p is the obstruction diameter and is the kinematic viscosity as introduced earlier. This definition of the Reynolds number will be used in the analysis. For comparison to earlier works the Darcy-velocity-based Reynolds number is used, defined as The value of U Darcy is related to the interstitial velocity via U int = U Darcy ∕ , where is the porosity.

Prandtl Number
The Prandtl number is defined as the ratio of the viscous diffusion to the thermal diffusion

Dimensionless Pressure Drop
The dimensionless pressure drop, p * , can be interpreted as a dimensionless measure of the deviation from Darcy's law for non-Stokesian flow. If p * is constant with respect to the Reynolds number then the pressure drop increases linearly with velocity. This expression can be written as Lundstrom et al. (2010) Here R is a reference length scale, which in this work is set to the radius of the obstructions, |F i | is the norm of the force density, i.e. the body force acting on each unit volume, and is the dynamic viscosity.

Nusselt Number
The Nusselt number Nu is a dimensionless measure of the heat transfer between the solid phase and the fluid phase defined as Foudhil et al. (2012) In this expression T is the temperature, n is the normal pointing out from the wall at the fluid solid boundary Ω , L is some length scale, here set to the obstruction diameter L = 2R . T fw is the wall temperature and T fm is the bulk temperature defined by Where u x is the streamwise velocity. For all cases discussed in this article T fw ≡ 0 which simplifies the expression for the Nusselt number to The average Nusselt number taken across the entire boundary Ω can then be written as

Numerical Method
An in-house finite difference, artificial compressibility solver with explicit time-stepping is utilized for the numerical analysis. The code is available under a free license at https://gitlab. com/c8383/gpu-fdm and further technical details regarding the implementation of the model is available in the repository. In addition to this, the discretization is described in detail in Appendix A. There are alternative models known to the authors such as LBM (Abbaszadeh et al. 2017;) and FVM methods (Lundstrom et al. 2010); the model in this work was chosen based on the simplicity of the discretization scheme which provides full insight into the flow field update procedure, which is essential to this study. Since the aim is to solve the equations for a fluid, Navier-Stokes equations are applied. The momentum balance hence becomes The standard continuity condition ( u i,i = 0 ) is replaced by an artificial compressibility condition (Kwak and Kiris 2011) Where is the compressibility factor. A known problem for collocated grids are checkerboarding (Ferziger et al. 2002), in this solver the checkerboarding is reduced by a pressure field smoothing operation which increases neighboring cell interactions, the expression for the smoothing operation is where S p is some smoothing factor. The thermal distribution is governed by the convectiondiffusion equation where T is the temperature, is the thermal diffusivity, and S is a source term. The density is set to unit = 1ML −3 for all calculations.

Temporal and Spatial Discretization
The spatial discretization is second-order accurate and utilizes central derivatives for all relevant quantities, i.e. if is some field variable then the second-and first-order derivatives are given by The computational kernel consists of seven neighboring interacting elements as visualized in Fig. 2. The temporal discretization is first-order accurate and utilizes a simple Eulerforward scheme, i.e. the value of the field variable at time t + 1 is only dependent on field variables at time t Here f is a general function which includes any field values at time t within the computational kernel. To ensure that the compressibility effects are minimal the Mach number, defined as cell velocity over sound velocity U∕U c , should be kept at a low value (Kwak and Kiris 2011).

Computational Grid
The nodes are spaced equidistantly in all directions ( Δx = Δy = Δz = 1 ) across a lattice of size n x × n y × n z . The nodes are then set as either wall elements or fluid elements, i.e. the geometry is approximated using the stepwise approximation as described in Ferziger et al. (2002).

Verification of Numerical Method
To ensure that the FDM model is accurate it is verified against two cases. For the thermal part the model is compared against an analytical solution for a heated channel case derived by Rybiński and Mikielewicz (2014). For the porous media part the dimensionless drag is compared to the simulations of Koch and Ladd (1997).

Thermal Plane Channel Flow
The thermal channel verification case, see Fig. 3 for the geometry, consists of a plane channel of dimensions a, b with constant temperature walls ( T w ), and a driving force ( F x ). The temperature gradient in the streamwise direction is specified as some value dT dx = T G . For this case the following analytical solution for the velocity was derived by Rybiński and Mikielewicz (2014) Where C n = (2n − 1) and K = a∕b . The analytical expression for the temperature is given by The calculation input values are summarized in Table 1, N t is the amount of time-steps before data is saved. The FDM model agrees perfectly with the analytical expression for the cases of velocity and temperature, see Fig. 4. Fig. 3 Geometry of the plane channel used for thermal verification case 1 3

Porous Media Flow
The porous media verification case consists of an ordered cubic array of cylinders with the REV being a single quadratic elementary cell of side length L and cylinder diameter D p , see Fig. 5. The flow is driven by a force gradient, F , which results in an average streamwise Darcy velocity, U Darcy , that can be used to calculate the dimensionless drag,  Table 2, N t is the amount of time-steps before data is saved for each viscosity, N n is the amount of viscosities, (n) is the viscosity for sweep number n starting from n = 0 and ending at n = N n − 1 . The FDM model agrees excellently with the computations carried out by Koch and Ladd (1997) throughout the inertial region and into the unsteady region ( Re Darcy > 150 ) as well, see Fig. 6.

Results
Three different types of porous media are investigated, see Fig. 7 for the geometries and the meshes. These are a staggered array of cylinders, in which the epsilon-shaped, inertial flow structure of interest here was first discovered, a staggered array of rods with a quadratic cross section, and a body-centered cubic (BCC) packing of spherical obstructions. The boundary conditions for the cells are periodic in all directions x, y and z where the wall elements are defined using a stepwise approximation, marked to the right in Fig. 7. The obstruction diameter D p is the width of the obstruction and is the resulting porosity from the stepwise approximation. The flow is driven by a body force in the x-direction, F x , and the viscosity is varied to investigate how the properties of the bed changes with Reynolds number. For all cases, the temperature of the wall elements is set to zero and a constant temperature source, S, is applied to all fluid elements. This method of obtaining the Nusselt number is described by Chu et al. (2019) who use a varying temperature source to enforce a constant heat flux, for a steady-state temperature field the source term S can be specified as a constant. An N t amount of time steps are calculated for each viscosity to allow all flow field variables to reach a steady-state. See Table 3 for a summary of the simulation parameters for each case. A grid sensitivity analysis is presented in Appendix B. The variation of the dimensionless pressure drop p * with Re is presented in Fig. 8, it can be seen how the value of p * tends to a constant value for low Re Darcian flow.

Flow Field Changes
For all geometries investigated a discontinuity is observed in a multitude of variables at a Reynolds number specific to each packing type. This discontinuity is accompanied by a dramatic change in the flow field as exemplified with the velocity fields in Fig. 9. For the staggered cylinder and staggered rod cases the flow field exhibits a symmetry break with the emergence of complex 3D flow structures. The value of Re where the structure appears coincides with the value in the experimental work presented in , appearing at an Re of ≈ 80 . For the BCC packing there are two separate transitions, one occurring at Re = 150 and a second one occurring at Re ≈ 230 − 250 . Since the transitions are triggered by an increase in Re it is assumed that they are caused by inertial forces starting to dominate over viscous forces.

Fig. 4 Comparison between analytical (Rybiński and Mikielewicz 2014) and FDM solution for thermal channel flow
As mentioned in the introduction (Johns et al. 2000) analyzed inertial transitions in packed beds by using magnetic resonance imaging and found an inertial transition characterized by a decrease in the spatial standard deviation of velocity with Reynolds number. Indications of such a behaviour can also be seen in . Hence, for comparison to these observations plots of the standard deviation of velocity taken over the fluid phase for the current geometries are presented in Fig. 10. For the staggered cylinder and staggered rod cases a clear reduction of the spatial variation of velocity can be seen at the transition as indicated by the vertical lines. For the BCC  packing a clear and dramatic reduction is observed for the first transition (red line) and there is also a minor change in the shape of the curve at the second transition (blue line). This indicates that the inertial transition calculated for these geometries is of the same type that has been observed earlier, characterized by a decrease in the spatial standard deviation.

Pressure Drop and Force Balance
In Fig. 11 the dimensionless pressure drop, as defined by Eq. (9), is plotted against Re for the staggered cylinder, staggered quadratic rod and BCC cases. At Re ≈ 85, 65, 250 , respectively, the pressure drop variation exhibits a sharp but small discontinuity accompanied by the aforementioned changes in the flow field that can be seen in Fig. 9. The transitions correlate with changes in the signed volume averaged forces as can be seen in Fig. 12 for the staggered cubic rods and BCC cases. However, the near non-existent change for the staggered cylinder case indicates that the streamwise volume-averaged ratio of pressure gradient to viscous forces is not a good general indicator of transition. The transitions also correlate with a significant change in the volume averaged norm of the convective and pressure gradient forces as defined by Eq. 5, see Fig. 13. For this case the BCC packing does not follow the same trend as the other cases, indicating that this is not a general feature of the transition. Fig. 6 Comparison between the computational results in Koch and Ladd (1997) and results from the FDM code for porous media flow (top). There is an excellent agreement. Visualization of flow regions (bottom), from left to right, Stokes region, inertial steady and inertial unsteady

Quantifying the Transition Onset
It is expected that the deviation of the flow around the impinging solid surfaces reduces the pressure build-up at the front of the obstruction, therefore it can be hypothesized that a Fig. 7 The two columns to the left show the side 3D-views and top views of the geometries investigated for the transition. From top to bottom a staggered packing of cylinders (a), a staggered packing of quadratic rods (b) and a BCC packing of spheres (c). The meshes for respective case can be seen in the two columns to the right (side 3D-views and top views) useful quantity for recognizing the onset is to watch for a transition in the pressure integral across the solid surfaces. This integral can be written where p is the pressure and Ω is the solid-liquid boundary. The value of p solid changes significantly at the transition, see Fig. 14. This discontinuity in pressure surrounding the solid surfaces indicate that the inertial cores no longer impinge the solid structures to the same extent as previously. For the BCC packing, the first transition reduces the value of p solid while the second transition increases it. This indicates that the first transition for the BCC packing may be of a fundamentally different type compared to the other transitions.

Impact on Heat Transfer Properties
Since the transition coincides with an increase of the resistance and small-scale structures of the flow, it can be expected that this would yield better heat transfer properties, similar to what is seen for the transition to turbulence in porous media (Abraham et al. 2011). From observing Fig. 15 it is clear that the heat transfer is significantly impaired . 10 The spatial standard deviation of absolute velocity in the critical Re range scaled by interstitial velocity U int for the staggered cylinder case (a), the staggered quadratic rod case (b) and the BCC packing (c) for both the staggered cylinder and BCC case but very slightly improves for the staggered rod case. It is therefore not possible to draw any general conclusions about how the transition will impact the heat transfer properties as it is geometry dependent. It can, however, be concluded that the transition can, significantly, affect the heat transfer.

Conclusions and Future Work
The flow fields of steady-state inertial transitions which occur across a wide range of high-tortuosity porous media has been scrutinized and a quantity which characterizes the transition has been proposed namely the pressure integral across the solid surfaces. The transition can, significantly, affect the heat transfer but also has some impact on all variables investigated which indicates a significant change in general flow characteristics. The BCC geometry case indicates that the transition can occur in systems similar to packed beds, which should be investigated explicitly in future studies. In addition to this a deepened study should be carried out on how the tortuosity of the porous bed impacts the possibility of the transition occurring. Fig. 11 The dimensionless pressure drop variation in the critical Re range for the staggered cylinder case (a), the staggered quadratic rod case (b) and the BCC packing (c) Fig. 12 The fraction of pressure gradient forces to viscous forces in the x-direction of the critical Re range for the staggered cylinder case (a), the staggered quadratic rod case (b) and the BCC packing (c) Where the last term p ,i depends on the index i and is equal to x,y,z+1 − p t x,y,z−1 .

Fig. 13
The volume averaged norm of the force for the cylinder case (a), quadratic rod case (b) and BCC structure (c)

Pressure Correction Discretization
The expression to obtain the updated pressure correction p t+1 x,y,z is based on a second-order accurate central derivative discretization of (15) with an advection term added, the expression becomes

Advection Diffusion Equation Discretization
The expression to obtain the updated species concentration or thermal value depending on interpretation of the variable T is based on a second-order accurate central derivative discretization of (17), the expression becomes

Pressure Smoothing Discretization
The pressure smoothing as described by (16) increases numerical stability and reduces checkerboarding. For this operation a second-order accurate central derivative discretization is used for the spatial differences and a first-order accurate Euler-forward temporal discretization Where the variable W x,y,z is the wall coordinate which is 1 if the element is a wall and 0 otherwise.
x,y,z . Fig. 16 The meshes, from top to bottom the staggered cylinder case, the staggered quadratic case and the BCC packing of spheres

Order of Operations and Wall Treatment
The order of operations is as follows, for the main kernel the momentum update, pressure correction update and advection diffusion update are applied to the state of the field variables at time t + 1 using the variables at time t. Following this n sweeps of the pressure smoothing operation is applied, for the cases presented in this article n = 1 . Wall elements are treated as constant value elements for all field variables except for the pressure correction equation.

Appendix B-Impact of Grid Resolution
In this section a grid sensitivity analysis is presented, four grid sizes are investigated for each geometry denoted by the short-side side length element count L, see Fig. 16. The numerical values used for the calculations are presented in Tables 4, 5 and 6. Due to the step-wise approximation of the geometry the porosity varies between the cases. The results converge well across all cases presented in Fig. 17 which are the dimensionless pressure drop p * , the surface integral of the pressure across the boundaries p solid and the Nusselt number Nu avg .

Data availability
The datasets generated and analysed during the current study are available from the corresponding author on reasonable request. The source code for the solver used in this study is publically available via the link https://gitlab.com/c8383/gpu-fdm.

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.
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://creativecommons.org/licenses/by/4.0/.