Estimating of Non-Darcy Flow Coefﬁcient in Artiﬁcial Porous Media

: This study conducted a radial ﬂow experiment to investigate the existence of non-Darcy ﬂow and calculate the non-Darcy “inertia” coefﬁcient; the experiment was performed on seven cylindrical perforated artiﬁcial porous media samples. Two hundred thirty-one runs were performed, and the pressure drop was reported. The non-Darcy coefﬁcient β was calculated and compared with available in the literature. The results showed that the non-Darcy coefﬁcient decreased nonlinearly and converged on a value within a speciﬁc range as the permeability increased. Nonetheless, it was found that the non-Darcy ﬂow exists even in the very low ﬂow rate deployed in this study. In addition, it has been found that the non-Darcy effect is not due to turbulence but also the inertial effect. The existence of a non-Darcy ﬂow was conﬁrmed for all the investigated samples. The Forchheimer numbers for airﬂow at varied ﬂow rates are determined using experimentally measured superﬁcial velocity, permeability, and non-Darcy coefﬁcient.


Introduction
Since it was inferred, Darcy's law has been used as the fundamental equation governing the flow of fluids through porous media. However, with the progress in research and experimentation, it became clear that Darcy's equation can predict the pressure gradient only in a limited range of flows. In the case of nonlinear behavior, the Forchheimer equation is used, which has the non-Darcy "inertia" coefficient β that can be calculated in the Laboratory. The flow of fluids through porous media plays a crucial role in understanding the interaction of fluids flow with the porous media. Since the flow in porous media differs from that of other types of flow, it was necessary to develop a different approach. Darcy's law describes the behavior of fluid flow in porous media. According to Darcy's law, the pressure gradient is linearly proportional to the velocity of the porous media fluid. The Darcy Equation (1) is an empirical equation based on experimental water flow through packed sands at low velocity, Zeng [1]: where gradP = pressure, µ = fluid viscosity, u = Darcy's velocity, the volumetric flow rate per unit flow area, and k = permeability of the medium. Many efforts have been made to derive the Darcy Equation theoretically via different approaches. Using the volumetric averaging theory, Stephen et al. [2] derived the permeability tensor for the Darcy Equation under low velocities. Following a continuum approach, Hassnizadeh et al. [3] developed a set of equations to describe the macroscopic behavior of fluid flow through porous media.
In the case of converting these equations into linear equations, a suitable equation results that work for the flow of porous media at low velocities. Through experiments carried out by many researchers and a lot of real field data, it can be concluded that Darcy's Equation works in low velocities. Excess pressure drop induced by inertial effects limits Darcy's law's applicability for modeling fluid flow through porous media at high velocities Green et al. [4]. Various terms, such as non-Darcy flow, turbulent flow, inertial flow, high-velocity flow, etc., have been used to address the linear Darcy Equation's deviation. Many efforts have been made to adjust the Darcy Equation. Forchheimer (1901) added a second order of the velocity term to represent the inertial effect and corrected the Darcy Equation into the Forchheimer Equation [5]: where ∆P is the pressure drop, µ is the fluid viscosity, β is the non-Darcy coefficient, and ρ is the fluid density. Non-Darcy behavior has shown significant influence on well performance. Researchers prefer to use the term "turbulent" or "non-Darcy flow" to describe the viscous flow at high velocities near the wellborn region, such as Kadi et al. [6]. This flow behavior is considered a non-Darcy flow rather than turbulent Belhaj et al. [7]; the gas slippage and inertial flow lead to non-Darcy behavior. If they are not taken into account, they will certainly cause measurement errors. These result from the flow of fluid particles through the throats of twisted rocks of various sizes. In a steady Darcy flow, there is an increase in pressure and no corresponding increase in the velocity of fluid flow Rushing et al. [8].
In addition, when the liquid particles pass through the throat of smaller pores, their velocity increases and slows down when they pass through large pores' throats, which leads to a dissipation of inertial energy and an increase in pressure Katz et al., [9]. Holditch et al. [10] presented a numerical model to study the non-Darcy effect on effective fracture conductivity and gas well productivity. They found that the non-Darcy effect could reduce the fracture conductivity by 20% and the gas productivity by 50%, the same conclusion reached by Guppy et al. [11] and Matins et al. [12].
In the reservoir, and especially the area near the well during injection or production, deviations from the linear Darcy's equation often occur due to the high flow velocity as the pressure differentials are large and affect the inertia force David et al. [13]. The non-Darcy coefficient in wells is usually determined by analyzing the multi-rate pressure test results, but such data are not always available. The permeability, porosity, and pore size distribution are crucial in the equation for predicting the non-Darcy coefficient developed by many researchers. In addition, the permeability heterogeneity and wettability are directly related to the effect of capillary number on viscous fingering patterns in porous media. Shiri et al. [14,15] studied the fault zone and a discontinuity; they also investigated the effects of wettability and permeability heterogeneities in the fluid front and preferential flow pathway. In the fault zone, they found that the pattern of the fault zone and the adjacent layer was different when the fault zone permeability was less or more than that of the vicinity layer, the sweep efficiency, and the fingering pattern. This phenomenon reduces the displacement efficiency by capillary trap mechanism. In addition, they concluded that the wettability difference between all the model components led to oil being cut off in wet oil regions. Faez et al. [16] studied the effect of fracture geometry on permeability; their study results concluded that increased fracture orientation would exponentially increase permeability. Namdari et al. [17] investigated the effect of the discontinuity direction on fluid flow in porous rock masses; they used a hybrid FVM-DFN and streamline simulation approach. The study results indicated that the FVM-DFN hybrid method is effective if it uses the streamlined simulation to study the fluid flow in a large model with discontinuity.
Many equations of non-Darcy coefficients based on mathematical models were presented; for example, [18][19][20][21][22][23][24][25] introduced a mathematical model to evaluate the fracture length and the non-Darcy coefficient. Using the model and the data from variable-rate tests from low-permeability hydraulically fractured wells, they were able to determine the non-Darcy coefficients.  Geertsma (1974) against the data obtained by [26][27][28][29]. They found that Equation (4) ( Table 1) was inaccurate. While suspecting tortuosity may influence data from limestone and sandstone samples measurements, Coles et al. [24] proposed two Equation (13) ( Table 1) to calculate the non-Darcy coefficient with the same porosity method, where β is expressed in 1/ft and K in MD. Comparing Equations (12) and (13) with equations developed by other investigators, the flow enters the exponents for porosity in Equations (12) and (13) Table 1 are positive instead of being negative in other equations.
Li et al. [30] used a combination of reservoir simulator and experimental procedure to investigate the non-Darcy flow in Berea sandstone cores, where nitrogen was injected at various flow rates in several different directions. Comparing differential pressures from simulations with their counterparts from experiments, they found a β for Berea sandstone.
Janicek et al. [9] suggested the non-Darcy coefficient equation for fluid flowing through sandstone, limestone, and dolomite porous media by rearranging Cornell's experimental results. Tek et al. [31] have used experimental data to generate an equation to evaluate β for any porous media system and come with a correlation, but the tortuosity was not considered (8) ( Table 1). Yuedong et al. [32] based on the dimensionless analysis method investigated samples displacement experimental data and built a new mathematical model to describe the seepage of high and highly production reservoirs and formed the following conclusions. Whether high-velocity non-Darcy flow occurs is determined by the value of the flow velocity and the combined action of the displacement medium (fluid), porous media, and external driving force. Studying the flow in porous media is understanding the flow behavior, and it is more challenging to inspect the fluid behavior in rock and simulate large-scale models. We have previously performed both single-phase and two-phase flow studies in nearwellbore regions. Single-phase flow experiments in heterogeneous core samples Shachi et al. [36] and two-phase flow modeling homogeneous core samples [37] were performed. We have also completed several studies related to foam flow in heterogeneous core samples without the presence of oil [38,39]. We have also created synthetic core samples and a perforation tunnel to conduct single-phase flow experiments and model a peroration tunnel [40][41][42][43]. This study is based on the non-Darcy flow performed radial steady pressure gradient-flow rate in large cores at a high airflow rate. The available experimental studies use relatively small samples without perforation and linear low flow rate range. This study utilizes the Radial Flow Facility "RFC." The RFC allows for performing a radial flow experiment on perforated samples, which is very rare in the literature; also in addition, the size of the samples is large enough to allow the flow to develop fully.

Formulation
Equation (2) is considered for a compressible flow, and the air is an ideal gas. The density is a function of pressure and temperature in the case of compressible flow. Assuming mass flow rate (Q m ), gas density (ρ), and volumetric flow rate (q) can be expressed as follows: where p = pressure, A = cross-sectional area of fluid flow, and v = fluid velocity. Then, it is possible to derive the following expression: M = molecular weight, z = compressibility factor, R = gas constant, and T = temperature. For gases, the equation is best expressed in terms of mass velocity Qm = ρv because the mass velocity is constant when the cross-section is constant. By substituting Equations (3)-(5) in Equation (2) for radial flow, we obtain: where r 1 is the perforation radius, r 2 sample radius, and h is the sample height, the integration of Equation (6) will result in the following equation: The last alternative to estimate is β to rearrange the Forchheimer Equation in the following form: The experimental data of ( were obtained in the laboratory and utilized linear regression, concluded that the beta factor β is constant for the range of flow rates of practical interests.

Experimental Procedure
The experimental procedure was divided into two stages. The first stage is the preparation of samples from sandstone and epoxy. The second is conducting flow experiments on the Radial Flow Cell Facility. The experimental process is briefly described as follows. Initially, the experiments started at a low flow rate, and then the flow rate increased until the non-Darcy flow occurred. The air compressor has a large flow rate and continuously

Experimental Procedure
The experimental procedure was divided into two stages. The first stage is the preparation of samples from sandstone and epoxy. The second is conducting flow experiments on the Radial Flow Cell Facility. The experimental process is briefly described as follows. Initially, the experiments started at a low flow rate, and then the flow rate increased until the non-Darcy flow occurred. The air compressor has a large flow rate and continuously works for an extended period. High-velocity flow experiments have been performed on artificial samples to simulate the airflow in the near-wellbore region. The samples are cylindrical with 15.54 cm in diameter, and the perforation has a 25.54 depth and 2.54 cm in diameter, as shown in Figure 1.

Preparation of the Samples
Seven artificial samples were created in memorial University laboratories. Both Aziz and Ahmed et al. [44] made samples in the labs to produce samples using sand and epoxy. Lately, seven synthetic sand samples were created from sand collected from Nova Scotia Canada in the Memorial University of Newfoundland laboratory. The sand was dried in the civil engineering lab using a hot-air oven at 110 C. The sieving process resulted in six different sand sizes ranging from 0.18 mm to 1.18 mm. The sand at different sizes is mixed with epoxy at different quantities. The mixture is then placed in a plastic container throughout four stages, using an electric vibrator to ensure grain distribution with the epoxy glue, and then they lift to dry and consolidate for 24 h, as shown in Figure 1. The sample dimensions are 30.48 cm high, 15.54 cm diameter, and a perforation tunnel has a 25.54 cm depth and 2.54 cm diameter Abobaker et al. [45].
There are many available ways to measure the properties of the samples. Mercury Intrusion Porosimetry (MIP) is one of the most widely used ways at the University of Memorial of Newfoundland labs to characterize and analyze the synthetic samples' pore morphology and index properties. The index properties include permeability, porosity, median pore diameter, and tortuosity, which are listed in the following Table 2.

Preparation of the Samples
Seven artificial samples were created in memorial University laboratories. Both Aziz and Ahmed et al. [44] made samples in the labs to produce samples using sand and epoxy. Lately, seven synthetic sand samples were created from sand collected from Nova Scotia Canada in the Memorial University of Newfoundland laboratory. The sand was dried in the civil engineering lab using a hot-air oven at 110 C. The sieving process resulted in six different sand sizes ranging from 0.18 mm to 1.18 mm. The sand at different sizes is mixed with epoxy at different quantities. The mixture is then placed in a plastic container throughout four stages, using an electric vibrator to ensure grain distribution with the epoxy glue, and then they lift to dry and consolidate for 24 h, as shown in Figure 1. The sample dimensions are 30.48 cm high, 15.54 cm diameter, and a perforation tunnel has a 25.54 cm depth and 2.54 cm diameter Abobaker et al. [45].
There are many available ways to measure the properties of the samples. Mercury Intrusion Porosimetry (MIP) is one of the most widely used ways at the University of Memorial of Newfoundland labs to characterize and analyze the synthetic samples' pore morphology and index properties. The index properties include permeability, porosity, median pore diameter, and tortuosity, which are listed in the following Table 2.

Performing the Flow Experiments
At this stage, seven sandstone samples were prepared; the next step is to conduct the flow experiments. RFC Figure 2 Ahamed et al. [44] has been updated to be suitable for conducting multi-phase flow experiments, as five pressure sensors have been calibrated with pressure sensors from other experiments. Two pressure sensors were placed on the inlet and outlet, and the rest was placed on the fluid mixing lines. With help from the university's technical department, an air flowmeter was repaired and successfully calibrated. In addition, another two lines were added to make the experiment ready to perform experiments on a three-phase flow. The experiment procedure begins by placing the sample in the cylinder and connecting the air compressor with the injection lines. The pressure sensors and flow meters are connected to the data acquisition to monitor the flow rate and record the pressure data during the experiment. The airflow rate was ranged from 3 to 99 L per minute. The outlet inserted into the pack measured line pressures without end effects. Note that the flow was entering the sample radially direction of gas flow to avoid spurious readings due to gas impinging on or accelerating off the probes (pitot effects).

Performing the Flow Experiments
At this stage, seven sandstone samples were prepared; the next step is to conduct the flow experiments. RFC Figure 2 Ahamed et al. [44] has been updated to be suitable for conducting multi-phase flow experiments, as five pressure sensors have been calibrated with pressure sensors from other experiments. Two pressure sensors were placed on the inlet and outlet, and the rest was placed on the fluid mixing lines. With help from the university's technical department, an air flowmeter was repaired and successfully calibrated. In addition, another two lines were added to make the experiment ready to perform experiments on a three-phase flow. The experiment procedure begins by placing the sample in the cylinder and connecting the air compressor with the injection lines. The pressure sensors and flow meters are connected to the data acquisition to monitor the flow rate and record the pressure data during the experiment. The airflow rate was ranged from 3 to 99 L per minute. The outlet inserted into the pack measured line pressures without end effects. Note that the flow was entering the sample radially direction of gas flow to avoid spurious readings due to gas impinging on or accelerating off the probes (pitot effects). The experimental run starts by installing the sample in the samples chamber (12) ( Figure 2) and then locking the lid tightly to prevent any leakage. The pressure sensors and the flowmeters are then connected to the Data Acquisition inlet, where each sensor has a specific channel. All the data for each run from the flowmeters and the pressure sensor were transferred to digital numbers and charts using Lab View. The flow starts by adjusting the air compressor at a certain pressure and putting the flow rate control valve The experimental run starts by installing the sample in the samples chamber (12) ( Figure 2) and then locking the lid tightly to prevent any leakage. The pressure sensors and the flowmeters are then connected to the Data Acquisition inlet, where each sensor has a specific channel. All the data for each run from the flowmeters and the pressure sensor were transferred to digital numbers and charts using Lab View. The flow starts by adjusting the air compressor at a certain pressure and putting the flow rate control valve at the needed flow rate. The flow enters the sample radially Figure 2; the outlet is installed at the top of the sample perforation. The flow rate range was between 3 and 99 LPM; at each run, the pressure and flow data are converted to an Excel table and then later analyzed.

Non-Darcy Flow Regime
If Darcy's law is considered assuming the air properties and k are constants and the term P 2 2 − P 2 1 plotted versus (Qm µ) would create a straight line. Figures 3 and 4 show the results of a test performed to check the nonlinearity of flow. By looking at Figures 3 and 4, the lines are not straight, and this indicates the presence of non-Darcy flow. In fact, these lines can be described by equations of the second degree.
During the flow of air through the samples, a pressure loss occurs, which results in a decrease in the volumetric flow rate, but the mass flow rate remains constant. These changes in pressure and volumetric flow lead to a change in air properties such as density and viscosity.

Non-Darcy Flow Regime
If Darcy's law is considered assuming the air properties and k are constants and the term ( − ) plotted versus ( ) would create a straight line. Figures 3 and 4 show the results of a test performed to check the nonlinearity of flow. By looking at Figures 3  and 4, the lines are not straight, and this indicates the presence of non-Darcy flow. In fact, these lines can be described by equations of the second degree. During the flow of air through the samples, a pressure loss occurs, which results in a decrease in the volumetric flow rate, but the mass flow rate remains constant. These changes in pressure and volumetric flow lead to a change in air properties such as density and viscosity.   and 4, the lines are not straight, and this indicates the presence of non-Darcy flow. In fact, these lines can be described by equations of the second degree. During the flow of air through the samples, a pressure loss occurs, which results in a decrease in the volumetric flow rate, but the mass flow rate remains constant. These changes in pressure and volumetric flow lead to a change in air properties such as density and viscosity.

Calculating Non-Darcy Coefficient β
A key to applying the Forchheimer Equation is to estimate a value for β. Methods developed to calculate β are based on experimental work, correlations, and the Forchheimer Equation.
To determine (β) experimentally, the procedure is first to measure each of the core samples' absolute permeability and then apply a series of increasing pressure differentials across each sample by flowing fluid through the core plugs at ever-increasing rates. Knowing the flow rates and pressure differentials across the sample, the inertial resistance coefficient can be directly calculated using linear regression of the Forchheimer Equation ) and Q m 2πhµ 1 r 1 − 1 r 2 were obtained in our Laboratory and utilized linear regression as in Figures 5-7. Table 3 contains the calculated values for the non-Darcy coefficient. The non-Darcy coefficient was compared with values obtained from models available in the literature. The closest result to the current study is Geertsma [27] model was the value of the difference is about 0.6 to 4.0%. It is clear from the table that there is a large discrepancy between the results obtained from this study and the results of previous studies. This discrepancy can be attributed to several reasons, the first of which is the size and nature of the samples. Studies use small samples during which the flow cannot reach full development. The other possible reason is that the current study benefits from the radial flow and the perforation, which is the outlet of the flow. As Zeng et al. [1] and others mentioned, the flow accelerates due to pressure differences when approaching the perforation, thus increasing the velocity and turbulence. The results obtained from the experiment were compared with previous studies in which the porous media used in those studies are similar to the samples of the current study.
A key to applying the Forchheimer Equation is to estimate a value for . Methods developed to calculate are based on experimental work, correlations, and the Forchheimer Equation.
To determine ( ) experimentally, the procedure is first to measure each of the core samples' absolute permeability and then apply a series of increasing pressure differentials across each sample by flowing fluid through the core plugs at ever-increasing rates. Knowing the flow rates and pressure differentials across the sample, the inertial resistance coefficient can be directly calculated using linear regression of the Forchheimer Equation ( Table 3 contains the calculated values for the non-Darcy coefficient. The non-D coefficient was compared with values obtained from models available in the litera The closest result to the current study is Geertsma [27] model was the value of the d ence is about 0.6 to 4.0%. It is clear from the table that there is a large discrepancy bet the results obtained from this study and the results of previous studies. This discrep can be attributed to several reasons, the first of which is the size and nature of the sam Studies use small samples during which the flow cannot reach full development other possible reason is that the current study benefits from the radial flow and the p ration, which is the outlet of the flow. As Zeng et al. [1] and others mentioned, the accelerates due to pressure differences when approaching the perforation, thus incre the velocity and turbulence. The results obtained from the experiment were comp with previous studies in which the porous media used in those studies are similar t samples of the current study.   Figure 8 reveals the variation of the non-Darcy according to the permeability from the modified Forchheimer plot of each sandstone sample. The plot illustrates the direct effect of permeability on the non-Darcy coefficient; with permeability increased, the non-Darcy coefficient decreases dramatically. These results show that the non-Darcy behavior is more severe in low permeability porous media. That can be explained as the pressure drop increases with the permeability decrease. Thus, the superficial velocity rises and leads to a more turbulent flow. The curvature of the aspect of decline is definite by increased permeability. Some previous studies proposed defining the equation of the non-Darcy coefficient with permeability. They attempted linear regression analysis on the experimental data set of the coefficient despite the curvature noted of transition [27].
It can be seen from Figure 9 that the non-Darcy coefficient decreases with the increase of median pore diameter. The non-Darcy coefficient decreases slowly in the bigger pore diameters, and that is because the flow cross-section area is larger in the porous media with a large pore. Figure 10 shows the non-Darcy coefficient versus the tortuosity of the samples. The graphs indicate that the non-Darcy coefficient is directly proportional to the tortuosity. The non-Darcy coefficient increases with the increase in tortuosity.
is more severe in low permeability porous media. That can be explained as the pressure drop increases with the permeability decrease. Thus, the superficial velocity rises and leads to a more turbulent flow. The curvature of the aspect of decline is definite by increased permeability. Some previous studies proposed defining the equation of the non-Darcy coefficient with permeability. They attempted linear regression analysis on the experimental data set of the coefficient despite the curvature noted of transition [27]. It can be seen from Figure 9 that the non-Darcy coefficient decreases with the increase of median pore diameter. The non-Darcy coefficient decreases slowly in the bigger pore diameters, and that is because the flow cross-section area is larger in the porous media with a large pore.  Figure 10 shows the non-Darcy coefficient versus the tortuosity of the samples. The graphs indicate that the non-Darcy coefficient is directly proportional to the tortuosity. The non-Darcy coefficient increases with the increase in tortuosity.  Tortuosity Figure 9. The effect of Median pore diameter on inertia coefficient β.

Forchheimer Number
Forchheimer number , which is used instead of the Reynolds number, which indicates when nonlinear effects occur, is given as the following: where is the Darcy velocity, the superficial velocity. for the sample is calculated for 0.0E+00 Tortuosity Figure 10. The effect of tortuosity on inertia coefficient β.

Forchheimer Number Fo
Forchheimer number Fo, which is used instead of the Reynolds number, which indicates when nonlinear effects occur, is given as the following: where v is the Darcy velocity, the superficial velocity. v for the sample is calculated for each flow rate as follows: Assuming there is no loss of air in the flow system due to leakage, reaction, or any other reasons, the mass flow rate in the air compressor is the same as that in the core sample at flow equilibrium, though the volumetric flow rates can be different due to the change in pressures. On the other hand, the mass flow rate equals the product of the density and the volumetric flow rate. Similar to the calculation of z and µ, the ρ density of air in the sample is calculated using the pressure data [1]. This behavior implies that the superficial velocity alone, as in the Reynolds number, is not a criterion for identifying the non-Darcy flow behavior; therefore, the Forchheimer number can be calculated. By analyzing the values in Table 4, it is observed that Fo increases nonlinearly with the increase in the flow rate, which is a natural result of increasing the velocity leads to an increase in the nonlinearity. The nonlinearity increases in the Fo indicate that the superficial velocity is not the only reason. However, interestingly, when comparing Fo with permeability, it is noticed that the divergence in Fo values is greater with the increase in permeability values at the same velocity. That conclusion is consistent with the fact that the non-Darcy behavior is more severe in low permeability porous media.

1.
This study conducted a radial flow experiment to investigate the existence of non-Darcy flow and calculate the non-Darcy "inertia" coefficient. Seven artificial samples were used. The flow rate of the air ranged from 3 LPM to 99 LPM, and in total, 231 run were conducted. 2.
Using the mean of pressure square difference versus Q m µ plot the non-Darcy behavior was conformed. This resulted in lines better fit to a polynomial.

3.
The non-Darcy coefficient β was calculated for each sample from the experimental results of the pressure gradient and using linear regression. The β measurement results were between 276,180.32 cm −1 and 19,589.15 cm −1 .

4.
The non-Darcy coefficient decreases with the increase in the median pore diameter and the porosity. When the median pore diameter at 25.31 µm non-Darcy coefficient β 276,180.32 cm −1 and at median pore diameter 181 µm non-Darcy coefficient β = 19,589.15 cm −1 .

5.
Forchheimer numbers for airflow at varied flow rates are determined using experimental permeability and non-Darcy coefficient data.