The Effect of Porous Media Grain Size on non-Darcy Flow Behavior using Pore Scale Simulation

The impact of grain size of the porous media of the 3D images of a beadpack on the permeability values and non-Darcy flow behaviour is analyzed using pore-scale flow modelling utilizing the finite volume method. In this study, the sample used was a 3D image from random tight packing of spheres (beadpack) with a size of 250×250×250 voxels and a porosity of 0.36. Variation in grain size is carried out by upscaling the 3D image sample, so that this treatment affects the specific surface area value but does not change the porosity and tortuosity of the sample. Several variations of grain diameter include 100 μm, 125 μm, 150 μm, 175 μm, and 200 μm. In this study, we adopted PIMPLE algorithm, which combines SIMPLE and PISO algorithms to accomplish the Navier-Stokes equations of fluid dynamics in such porous media. The simulation results show that the permeability value exhibits an inverse relationship with the specific surface area. This correlation is in accordance with the Kozeny-Carman equation. In addition, the initiation of non-Darcy flow occurs at the reference length-based Reynolds number of 2.06 and the permeability-based Reynolds number of 0.049.


Introduction
Darcy's law plays a vital role in creeping fluid flow in porous media, which express that the pressure gradient is proportional to the flow rate [1].However, Darcy's law no longer applies as the flow rate increases.The correlation between pressure gradient and flow rate becomes non-linear due to the dominating inertial effect, and in this situation, the flow transitions into the non-Darcy regime.Forchheimer [2] added a quadratic term to the Darcy equation as a non-linear correction known as the Forchheimer equation.The supplementary term in the Forchheimer equation demonstrates a direct correlation with the β non-Darcy factor.The initiation of non-Darcy flow and determination of the β non-Darcy factor in porous media are typically established through the findings obtained from pressure variation tests.The findings have been utilized to formulate empirical relationships to determine the β non-Darcy factor [3][4][5][6].However, due to the microscopic heterogeneity in reservoir rocks, these correlations are anticipated to provide unreliable forecasts for particular samples lacking experimental data.[7,8].The structural heterogeneity of porous media is characterized by tortuosity, pore size distribution, coordination number, etc. [9,10].
Numerical simulations aim to overcome this issue and offer microscopic understanding into flow phenomena in rock pore spaces by solving fluid dynamics equations directly.However, the intricacy of pore structure and the necessity to accurately quantify flow phenomena across representative volume elements has currently restricted the applicability of this method.Several simplified media and spherical packages have been studied in various numerical investigations.Fourar et al. [11] adopted the finite element method to model and estimate the initiation of non-Darcy Newtonian fluid flow in 2D and 3D periodic spherical packages.Furthermore, Newman and Yin [12] predicted the initiation of non-Darcy flow and the determination of the dimensionless inertial resistance factor in artificial 2D media using the Boltzmann lattice method.Thauvin and Mohanty [13] used a different method, the pore network model, to describe the increasing inertial effects.They incorporated the pore size distribution and coordination number into the model, and several quantities such as permeability, β non-Darcy factor, tortuosity, and porosity were obtained.In recent years, several studies related to pore-scale simulations have been applied to actual images of porous media.Generally, the research focuses on the characteristics of permeability and non-Darcy flow based on the pore structure parameters of porous media.For example, Muljadi et al. [14] explored the influence of the heterogeneity of porous media on non-Darcy flow characteristics, namely non-Darcy onset and Forchheimer's constant.In his research, three types of 3D images from porous media were used: beadpack, bentheimer, and estaillades.Cheng et al. [15] utilized the lattice Boltzmann method to investigate the effect of pore stucture on non-Darcy flow behaviour where the pore structure involves porosity and specific surface area.Moghimi et al. [16] analyzed the phenomenon of non-Darcy flow within highly porous structures featuring various levels of geometric intricacy.They explored the influence of foam structure porosity and differences in pore number on non-Darcy flow behavior and permeability.The relationship between permeability and pore structure has been further investigated in recent research [17][18][19][20].In porous media, pore size is one of the determining factors in the characteristics of pore structure, thus clearly influencing the flow through the rock.Research related to the influence of grain size on permeability has recently been conducted using a pore-scale approach linked to the porosity of the samples [21].Furthermore, Payton dkk [22] experimentally studied the influence of grain size on nonlinear flow behavior.It is worth noting that experimental methods can affect the arrangement and distribution of grains.To date, there has been no pore-scale study analyzing how grain size specifically affects non-Darcy flow behavior.In this study, we will investigate the effects of consolidated porous media grain size on permeability values and non-Darcy flow behavior.We use a 3D image of a bead (sphere) pack porous media as our sample.The finite volume method is employed to model the fluid flow occurring within the pore space of the image.A crucial step in this method is the integration of differential equations, specifically the incompressible Navier-Stokes equations, within a three-dimensional control volume [23].Grain size variations are introduced by upscaling the 3D image sample, thereby influencing specific surface area values without altering other parameters like sample porosity and tortuosity.

Images and physical parameters of the sample
Fluid flow in the Darcy and Forchheimer regimes is simulated in porous media from 3D bead (sphere) pack images where the samples are based on determination of the centre coordinates of identically sized spherical particles in random tight packing [24], which are then segmented into images by Prodanovic´ and Bryant [25].The beadpack has a porosity of 0.36 with a size of 600×600×600 voxels.However, in this study, the sample used was 250×250×250 voxels, which were cropped from the original image.The image can be seen in figure 1. Variations in grain size diameter are carried out by upscaling the original image with a particular factor so that several variations are obtained, including 100 μm, 125 μm, 150 μm, 175 μm, and 200 μm.This method did not change the pore structure of the sample, such as porosity, tortuosity, and density distribution.However, this method changed the specific surface area of the sample.The characteristic length L of each sample is obtained from the grain diameter (DP), namely L = DP, while the calculation of the specific surface area of the porous media consisting of spherical grains is Sv = 6/DP [26].Information related to grain diameter, image resolution, reference length, and specific surface area are all provided in Table 1.

Numerical method
Fluid moving through porous media is assumed to satisfy incompressible Navier-Stokes equation described as [27,28]: Here, p and v represent pressure and velocity vectors, respectively.For the steady state solution the time derivative in equation ( 1) is equal to zero.Equation ( 2) represents the continuity equation, which is an implication of incompressible fluids (i.e., constant density in all infinitesimal volumes) [27] The pressure and velocity in equation ( 1) are accomplished on the pore space discretization of the beadpack, as depicted in figure 1(c), using implicit pressure with the PIMPLE algorithm (see figure 2), which combines SIMPLE and PISO algorithms [29,30].The numerical solver employed is embedded in OpenFOAM, an open-source CFD software designed for solving the Navier-Stokes equation [31].
The solution is considered convergent if e ≤ 10 -6 where . Vn is the magnitude of volumetric flowrate at the outlet of the computation domain at time level n.
In all simulations, the computational domain matches the image shape with one inlet and one outlet facing each other.The pore space is discretized with a resolution of 1:1 concerning the image resolution.We applied specific pressure values on the image's inlet and outlet boundaries, while zero velocity boundary conditions are specified on the other surfaces.The fluid used in this study is water with a viscosity of 0.001 kg/m.s and a density of 1000 kg/m 3 .

Permeability and criteria for non-Darcy flow
At extremely small Reynolds numbers (Re <<1), Darcy equation accurately describes the correlation between pressure gradient and superficial fluid velocity in porous media [1,32].
where KD represents the Darcy permeability tensor, μ represents the fluid viscosity, and v=q/A, q is the flow rate and A is the cross-sectional area of the specimen.The linear relationship shown in equation ( 3) arises from the flow of Stokes with the omission of the non-linear inertia term.As the flow rate increases, it becomes necessary to acknowledge the presence of the non-linear inertia term, Here, ρ represents the fluid density, KF represents the permeability tensor according to Forchheimer, a value that is similar to but not identical to KD, n is the unit vector in the p  direction, and β is known as the non-Darcy coefficient or Forchheimer constant.Calculation of permeability values for each case is deduced from equation (3): where p  is the pressure difference between the inlet and outlet divided by the length of the specimen, and v is velocity in the In this study, the termination condition for Darcy flow is met once the pressure gradient resulting from the linear term falls below 99% [33].To calculate the Reynolds number for the initiation of non-Darcy flow, the authors use two formulations.The initial aspect is the conventional interpretation of the Reynolds number which is established on the characteristic length L, i.e., [34] L The next expression of the Reynolds number is achieved through the establishment of L using the Brinkman screening length ( D K ), i.e., [12,35] To decide the initiation of non-Darcy flow, dimensionless apparent permeability is defined in a subsequent manner [14] * Considering the previously determined non-Darcy flow onset criteria, in this study, the initiation of the non-Darcy flow that will be employed is the juncture at which K* equals 0.99.The Kapp value is the same as the K value in equation ( 3) which is calculated for all regimes, both Darcy and non-Darcy.
Mathematically, the expression for the apparent permeability can be formulated as follows [14].
The β non-Darcy factor is derived from the gradient of the graph of equation ( 9), namely by plotting the reciprocal of the apparent permeability 1/Kapp against ρv/μ in the Forchheimer regime.

Results and discussion
To analyze the non-Darcy flow behaviour, We conducted flow simulations using varying pressure gradients across the pore space of all beadpack samples.The flow velocity varied so as to cover flow regimes from Darcy to Forchheimer.The output of the simulation is the volumetric flow rate, which is then employed to calculate the permeability value of each condition.( ) where  represents the Kozeny-Carman constant and  = 5 applies to porous media with spherical particles [36].The permeability values from the simulation results and from the Kozeny-Carman equation are displayed in Table 2.The Reynolds number is obtained from the flow velocity of each sample where a pressure gradient of 0.0001 MPa/m is applied.The permeability prediction based on the Kozeny-Carman equation is similar to our results and the deviation is still below 10%.We plotted the permeability values for each sample to obtain a permeability relationship based on specific surface area, which is presented in figure 3. The figure shows a similar trend between our results and Kozeny-Carman's prediction where permeability increases with decreasing specific surface area.The dotted line in figure 3 is a curve fitting from the simulation results which states the proportional relationship between permeability and specific surface area with the order of -2.068, which is in good agreement with equation (10).The critical ReL value is in good agreement with Scheiddeger [37] and Hassanizadeh and Gray [6], but slightly smaller than the estimations of Fancher and Lewis [5] and Bear [4].Meanwhile, if we use the Reynolds number formula from Ergun, Recritical = ReL/(1- ) then we get Recritical = 3.21, where this value falls within the Ergun estimate range of 3-10 [38].This demonstrates that the predicted initiation of non-Darcy flow is in numerical alignment with experimental data acquired from similar treatments.

Conclussions
In this research, the impact of grain size of 3D image of packed particles on the initiation of non-Darcy flow has been examined through numerical modeling.The initiation of non-Darcy flow is characterized as the juncture at which the dimensionless apparent permeability, denoted as K*, equals 0.99.According to this particular criterion, our investigation indicates that the critical Reynolds numbers ReL and ReK obtained are 2.06 and 0.049 respectively.Our calculations of Darcy permeability at S1, S2, S3, S4, and S5 agree with the Kozeny-Carman equation with tolerances below

Figure 1 .
Figure 1.Appearance of the beadpack image: (a) 2-D binary image consisting of pores (black) and grains (white), (b) 3-D voxelised grain (solid) space, (c) 3-D voxelised pore space Each sample was simulated at extremely low Reynolds number (Re << 1) and the Darcy permeability was calculated based on equation(3).Darcy permeability can also be predicted from porosity (  ) and specific surface area ( v S ) which is formulated by the Kozeny-Carman equation as follows.

Figure 3 .
Figure 3.The relationship between specific surface area and permeability.

Figure 4 (
Figure 4(a) depicts the dimensionless permeability K* for samples S1, S2, S3, S4, and S5 given as a function of ReL, covering both Darcy and non-Darcy flow regimes.In Figure 4(b), the plotted data illustrates the relationship between K* and ReK.The critical ReL and ReK values, representing the transition from Darcy to Forchheimer flow are taken at K* = 0.99, as marked by the dashed red line.The trends in the ReL and ReK graphs are consistent across all samples and align with each other.Figure 4 reveals a critical ReL value of about 2.06 and a critical ReK value of about 0.049.

Figure 4 Figure 4 .
Figure 4(a) depicts the dimensionless permeability K* for samples S1, S2, S3, S4, and S5 given as a function of ReL, covering both Darcy and non-Darcy flow regimes.In Figure 4(b), the plotted data illustrates the relationship between K* and ReK.The critical ReL and ReK values, representing the transition from Darcy to Forchheimer flow are taken at K* = 0.99, as marked by the dashed red line.The trends in the ReL and ReK graphs are consistent across all samples and align with each other.Figure 4 reveals a critical ReL value of about 2.06 and a critical ReK value of about 0.049.

Table 1 .
Physical parameters of several samples used in this study after upscaling from the original image (DP = 100 μm).The grain diameter is taken as the characteristic length L

Table 2 .
Darcy permeability for the the samples calculated at low Reynolds number.