Numerical Simulation of the Non-Darcy Flow Based on Random Fractal Micronetwork Model for Low Permeability Sandstone Gas Reservoirs

Darcy’s law is not suit for describing high velocity flow in the near wellbore region of gas reservoirs. The non-Darcy coefficient β of the Forchheimer’s equation is a main parameter for the evaluation of seepage capacity in gas reservoirs. The paper presented a new method to calculate β by performing gas and con-water flow simulations with random 3D micropore network model. Firstly, a network model is established by random fractal method. Secondly, based on the network simulation method of nonDarcy flow in the literature of Thauvin and Mohanty, a modified model is developed to describe gas non-Darcy flow with irreducible water in the porous medium. The model was verified by our experimental measurements. Then, we investigated the influence of different factors on the non-Darcy coefficient, including micropore structure (pore radius and fractal dimension), irreducible water saturation (Swi), tortuosity, and other reservoir characteristics. The simulation results showed that the value of the non-Darcy coefficient decreases with the increase in all: the average pore radius, fractal dimension, irreducible water saturation, and tortuosity. The non-Darcy coefficients obtained by the fractal method of microparameters are estimated more precisely than the conventional methods. The method provides theoretical support for the productivity prediction of non-Darcy flow in gas reservoirs.


Introduction
The Darcy law describes the flow of the subsurface fluid and shows the linear relationship between the pressure gradient and the volume flow (Darcy's velocity). However, Darcy's law only applies to laminar and viscous flows. At high velocities, such as high-velocity flows near the wellbore, fluid flow behaves as an inertial flow, and Darcy's law fails to follow this high-velocity seepage law. Forchheimer proposed a classical Forchheimer equation with a non-Darcy coefficient β to estimate non-Darcy effects in Equation (1) and accounted for the impact of both viscous and inertial effects.
Since then, the non-Darcy coefficient β has attracted a lot of research. Coefficient β can be obtained from analytical methods, physical experiments, and numerical simulations. Geertsma [1] conducted an experiment to establish the quantitative relationship between non-Darcy coefficient and porosity and permeability as well as influence of the immobile fluids on non-Darcy coefficients. Frederick and Graves [2] also studied experimentally the seepage process of gases in porous media considering the influence of both mobile and immobile water on nonlinear seepage of gases. Reid et al. [3] researched high-speed non-Darcy flow process of gas-water fluids in porous media and determined that the presence of irreducible and mobile water affects the gas flow in porous media, both of which reduce the gas permeability and increase the non-Darcy coefficient. However, at the same time, the effect of movable water on gas permeability and non-Darcy coefficient was found to be more significant than that of irreducible water.
Thauvin and Mohanty [4] reported the nonlinear seepage phenomena of gas by establishing the pore structure model considering the influence of pore characteristics on the non-Darcy coefficient. Firoozabadi [5], Jones [6], and others obtained the non-Darcy coefficient β for different rocks by physical methods. A 2D random pore network developed by Neeman et al. [7] was found that the porous medium structure has a strong effect on the flow properties. Later, Chukwudozie et al. [8] used LBM to estimate the permeability, tortuosity, and β factor of porous rock.
Dou [9] used analytic methods and determined that non-Darcy β coefficient is a function of pore size, permeability, friction coefficient, and drag coefficient. Li [10] expanded his analytical research and determined that the β is also dependent on the fluid density and fluid viscosity. Wang et al. [11] obtained the relationship between β, permeability, and porosity in high-pressure gas reservoirs by statistical analysis of data from multiple sets of rock samples under different overburden tests. Ren et al. [12] established the nonlinear seepage law affected by irreducible water in low permeability reservoirs based on the laboratory seepage experiment results.
Thus, just from the short overview of the previous literature mentioned above, one can see that a lot of theoretical and experimental work was conducted to obtain different forms of non-Darcy coefficient expression, of which the main form is the quantitative relationship between permeability, porosity, and the non-Darcy coefficient. However, due to the different methods applied as well as various porous media and fluids, the non-Darcy coefficient expressions are quite different [13][14][15]. There are two orders of magnitude difference in the results calculated by different expressions with the same parameters. The fundamental cause of this difference is the intrinsic complexity of micropore structure.
The goal of our work is to obtain an inner relationship of the non-Darcy flow in porous media by a random fractal pore-scale network model. We carried out some flow experiments at different pressure drops and showed the validity of the Forchheimer's equation. We vary certain pore-scale parameters, study the way it affects various flow parameters, and identify the correlation between different parameters. The mathematical characterization was carried out to confirm the essence of the nonlinear flow. This study is limited to single-phase flow in isotropic porous media.

Methodology
The section mainly presents the dynamic simulation method of non-Darcy flow using random fractal network model.

Construction of Network Model by Random Fractal
Theory. Since the first pore network by Fatt [16], a pore numerical model is a versatile studying platform for a variety of subjects [17]. Piri and Blunt [18] developed a quasistatic random network model for researching numerous fluid configurations for two-and three-phase flow. Even since then, the complex geometry of the pore space may be reconstructed [19,20]. Fang et al. [21] established a fundamental network model, and Mohammadi et al. [22] modified the model to calculate Scri. Consequently, the construction of pore networks is extremely effective for the research of various type of flowing process in porous medium [23][24][25].
The pore sizes and the throats in the pore network model are also the key to be determined for researchers [26,27]. It is difficult to describe the irregular pores in traditional geometry. Various studies have shown that the intrinsic properties of the pore meet the statistical significance of fractal distribution [28][29][30][31]. We have recently established a constructing digital core method based on fractal theory [32]. The fractal characteristics of pore size are calculated by capillary pressure curve measured from the mercury injection experiment [33]. We then developed a procedure for constructing a 3D pore-throat structure.
Our last paper [27,32] gives the detailed derivation procedure for the fractal relationship between the pore radius (r), fractal dimension (D), minimum pore radius (r min ), and random number (ζ).
Porous media are considered to be consisted of large amounts of pore throats which interconnect large and small pores and shrinkage intervals. These pores and pore throats can be described as "points" and "edges" in a 3D network model. The pore nodes are randomly distributed, and the pore throats are connected with a pore as the pore's coordination number which is generally between 1 and 6. Based on the node position coordinates, the pore's coordination number, and the principle that the closest pores should be priorly connected, one can first establish the main structural framework of a 3D PNM by connecting the "points" and "edges" and then obtain the geometrical parameters of pores and throats including its radius and length, shape factor, and volume to form a PNM with the topology structure and geometric characteristics of actual core pores.

Numerical Simulation of Gas Phase Non-Darcy Flow.
The dynamic model is adopted using the Matlab program. The pore-throat structure shown in Figure 1 was used as a unit Right pore Left pore Half pore length 2Rt 1 t x Figure 1: Dimension of the throat structure in the pores.

Geofluids
to analyze the pressure drop caused by the fluid flowing into the pores [34]. The fluid flow corresponds to the Hagen-Poiseuille equation during the laminar flow in a porous medium. Its pressure drop is a primary mechanism to overcome the viscous drag. In our model, in addition to considering the viscous drag, the tortuosity and the pressure drop of inertial resistance caused by the gradual pore-throat radius gradient were taken into account. Based on the pore level law presented by Thauvin and Mohanty [4], the viscous pressure drop between the centers of two adjacent pores can be described by the following equation in single-phase Darcy (slow or Stokes) flow, In the near wellbore region, the velocity is higher, and Reynolds number increases to the order of Re~10 in typical gas wells. In porous media, the inertial terms become so important that it affect the additional pressure loss. It was also taken into account: The pressure drop resulting from the narrowing of the tunnel can be expressed as: The pressure drop caused by widening of the tunnel is expressed by the following equation: Based on the above equations, the total pressure drop across the structural unit can be obtained as: The pressure and velocity relationship in pore level is combined in the network model to calculate macroscopic flow parameters. Equation (7) is shown below: wherein a and b are shown, respectively: Mass balance is applied in one of the principal directions of the pore. The relation between velocity, u in a throat, and the pressure drop between the two adjacent body centers is given by Equation (8). To solve this system of equations, we used the Newton-Raphson method. The calculations are repeated at several different pressure drops to verify the Forchheimer equation. The calculations are also repeated for pressure drops applied in each of the pore directions.
In this work, numerical simulations in pore scale are run at specified pressure gradients; meanwhile, the boundary conditions in the network consist of fixed inlet and outlet pressures.

Model Validation
We designed six groups of non-Darcy experiments for low permeability sandstone samples. To ensure that the experiments are carried out in porous media, CT scanning was used to scan the samples to avoid the existence of microfractures in the stage of preparation and the end. Finally, two real cores with #1 and #2 were selected.
3.1. Establishment of Pore Network Model. The original data and calculated fractal feature parameters of the rock samples in the certain gas reservoir are given in Table 1 by the capillary pressure curve Figure 2. The first step gets the information of the pore and throat size distribution using Equation (2). Through the relationship between porosity and the pore number, the porosity of rock samples #1 and #2 is 15.2% and 17.2%, respectively, so the corresponding pore numbers are 1053 and 1268. The second step is to generate a pore network model that best fits the original shape of the rock sample, the best model that accurately represents the true capillary pressure characteristics of the core sample.
The size of the random fractal micropore network established in this paper is 3:00 mm × 3:00 mm × 3:00 mm, as shown in Figure 3, in which the pore numbers are 1053 and 1268; the throat numbers are 2476 and 2982, respectively; the average coordination numbers are 4.362 and 4.614; and the tortuosities are 1.46 and 1.62. In this model, we assume 3 Geofluids that the shape of pore and throat is triangle, and the value ranges from 0 to 0.0418. The shape factor of pore and throat was given by the Weibull distribution.

3.2.
Simulation of the Non-Darcy Experiment. Nitrogen and distilled water corresponded to the gas and liquid phases, respectively. The experiments were performed at room temperature. The core was placed in a holder and fixed with a pressure ring. The outlet pressure was fixed at a certain value. Inlet pressure was regulated to change the pressure gradient as needed. Gas flow was monitored and adjusted by flow meters. The experimental results are shown in Table 2.
We proved the validity of the new method by comparing it with the experiments result. The porous network in our model was filled with water and gas as well as the experiments. Water in the pores was in the film form; thus, it did not participate in the flow. The irreducible water in the network model was established according to the principle of small channel occupancy. The gas state was chosen to be constant because of the very small pressure changes under the simulated conditions. The network model had four closed sides and two open ones (the entrance and exit sides). There was no fluid flow through the four closed sides. The flow process was in quasisteady state, and the flow parameters only changed with space but not over time [35]. The model ambient conditions were room temperature. The nitrogen viscosity is 0.011 mPa·s, and the density is 150 kg/m 3 .
Based on the two sets of network models established earlier in the text, we simulated the single-phase nonlinear gas flow and obtained the average gas flow rate under different pressure gradients. The relationship between gas flow rate and pressure gradient was obtained. The experimental results were compared with the simulated ones (see Figure 4).
As can be seen from Figure 4, the flow rate gradually increases with the pressure gradient increase, showing a significant nonlinear relationship. The experimental results are relatively close to the simulated ones.
According to Darcy's law, the corresponding gas permeability at each flow rate can be obtained using the following equation: As the gas velocity increases, the gas permeability calculated from the Darcy's formula decreases mainly because the gas shows nonlinear seepage flow type, and the flow resistance increases together with the flow rate. The calculated permeability was less than the intrinsic absolute permeability of the model.
According to the Forchheimer equation, the seepage velocity and pressure gradient of a fluid can be expressed as a binomial equation, in which the coefficient of the binomial first power phase is the ratio of the fluid viscosity μ to its permeability k. The quadratic phase coefficient is the product of non-Darcy coefficient and fluid density. The quantitative relationship between the pressure gradient and the flow rate can be determined using the above equations. The pressure gradient and the flow rate are obtained from the flow simulation.
where k g is the gas permeability for different gas velocities, k is the absolute permeability of the network model, β is a non-Darcy coefficient, ρ is the gas density, μ is the gas viscosity, and v is the gas flow rate. The slope of the linear relationship between 1/k g and ρv/μ in Cartesian coordinates is the non-Darcy coefficient β. It was compared with the calculated result of Geertsma formula, which is expressed as: In the formula, S wi represents the irreducible water saturation, and ∅ represents the porosity.
The non-Darcy coefficients calculated by the network model and the experiments are shown in Table 3. However, the results of the two samples calculated by the Geertsma formula are close due to the similar values. The difference between two samples in the internal microscopic structure causes the different coefficients. The non-Darcy coefficients obtained by the fractal method of microparameters are estimated more precisely than the conventional methods. It shows that the non-Darcy flow network model is effective.

Factors Affecting the Non-Darcy
Coefficient by the Fractal Method Core sample #1 Core sample #2 Figure 2: Capillary pressure curves of the two groups of rock samples. 4 Geofluids

The Impact of Irreducible Water Saturation.
To study the influence of irreducible water saturation, we simulated the nonlinear seepage flow under five irreducible water saturation levels (S wi = 0, 0.1, 0.2, 0.35, and 0.5) with different fractal dimensions and obtained the non-Darcy coefficients under different conditions. Several data points from the literature [36] were added to the graph as well (see Figure 5). Non-Darcy coefficient increased with irreducible water saturation [37][38][39] was mainly because the irreducible water filled the corners of the pores, only partially leaving the rest of the space for the gas. The indirectly reduced throat radius decreases the flow passage and aggravates the nonlinear flow of gas. Gas flow radius has a quantitative relationship with the gas saturation S g . Since there are only gas and water phases in the pores, S g = 1 − S w . We used the power function to fit the relationship between non-Darcy coefficient and irreducible water saturation: where β is the non-Darcy coefficient in m -1 ; a and b are parameters (b is~-4.5 according to the fitting).  Figure 6. The fractal dimension reflects the complexity of pore radius distribution. The larger the fractal dimension, the greater the heterogeneity of pore radii and the larger the non-Darcy coefficients are. The relationship between the non-Darcy coefficient and the fractal dimension can be fitted using exponential function: where β is the non-Darcy coefficient in m -1 ; a and b are the parameters (the value of b obtained from the fitting was 4.5).
4.3. The Impact of Average Pore Size. To analyze the effect of pore size on non-Darcy coefficients, the value of r min was changed during the reconstruction of the rock sample data by using the random fractal method. Five threedimensional network models with average pore radii of r = 1, 3, 5, 7, and 9 μm were established. The simulation of nonlinear seepage of gas under different irreducible water conditions was carried. The corresponding non-Darcy coefficients are shown in Figure 7.

Geofluids
The relationships of non-Darcy coefficient and the average pore radius are consistent under different irreducible water saturations (see Figure 7): the non-Darcy coefficient decreases with the increase of the average pore radius, which agrees with the literature data [36]. This relationship can be fitted using the following power function: where β is non-Darcy coefficient in m -1 ; a and b are parameters (the value of b obtained from the fitting was~-4.5).

The Impact of Tortuosity.
To analyze the effect of tortuosity on non-Darcy coefficients, the value of tortuosity was changed, while the other parameters remained unchanged. Five sets of three-dimensional network models with tortuosity of τ = 1, τ = 1:3, τ = 1:5, τ = 2:0, and τ = 2:5 were established and obtained non-Darcy coefficients are shown in Figure 8.
The relationship between the non-Darcy coefficient and the tortuosity behaves similarly, and the non-Darcy coefficient increased with the increase of the tortuosity. We fitted this relationship using the following power function: where β is non-Darcy coefficient in m -1 ; a and b are parameters (the value of b obtained from the fitting was~1.52).

Fractal Characterization of Non-Darcy Coefficients.
According to the analysis in the previous section, the non-Darcy coefficient can be characterized as a function of the average pore radius, fractal dimension, tortuosity, and irreducible water saturation. The mathematical characterization form is shown below. β = αD 4:5 τ 1:52 where α is the proportionality constant, and its value depends on the fluid properties (viscosity, density) as well as other parameters. Non-Darcy coefficients and their related parameters were inserted into Equation (17) to obtain α for different conditions. The value of α was calculated to be equal to 3 × 10 8 through trial and error calculations. Therefore, the   Figure 5: Curves of non-Darcy coefficient and irreducible water saturation at different fractal dimension. 6 Geofluids mathematical expression of non-Darcy coefficient of this model is: β = 3 × 10 8 × D 4:5 τ 1:52 The non-Darcy coefficient calculated by Equation (18) and the corresponding experimental results are shown in Table 3. We can see that the order of magnitude of non-Darcy coefficient calculated by the experimental data was similar to Equation (18). Equation (18) we conducted by the fractal method of microparameters reflects the influence of rock internal structure on seepage flow. It is more precise than the conventional methods that are obtained by the physical parameters of porosity and permeability. The mathematical characterization was carried out to confirm the essence of the nonlinear flow. This study is limited to single-phase flow in isotropic porous media.

Conclusion
(1) Fractal theory was used to analyze the microscopic pore size distribution obtained by two sets of real core mercury injection methods to determine the fractal dimensions and pore sizes. The constructed pore network models are in good agreement with real cores (2) The gas single-phase nonlinear seepage model was established. Besides the viscous resistance, both the tortuosity and the gradual change of the pore-throat radius can produce extraresistance. Based on the three-dimensional network model, the relation curve of flow rate and pressure gradient was obtained by the simulation of the single-phase seepage with irreducible water, and these results were compared with the non-Darcy tests.  Permeability tensor (Darcy) L: Path length from one body center to another (cm) L t : Pore-throat length (cm) Δp: Pressure drop (atm) Δp b : Pressure drop due to bending in the path (atm) Δp c : Pressure drop due to contraction (atm) Δp e : Pressure drop due to expansion (atm) Δp s : Pressure drop for Stokes' flow (atm) P: Pressure (atm) q: Flow rate (cm 3 /s) r b : Body radius (cm), r b1 is the radius of pore 1, r b2 is the radius of pore 2 r t : Throat radius (cm) r min : Minimum value of the radius (cm) r max : Maximum value of the radius (cm) r: The pore radius (cm) D: Fractal dimension ζ: Random number u: Interstitial velocity (cm/s) υ: Superficial velocity (cm/s) β: Non-Darcy coefficient (cm -1 ) μ: Viscosity (g·cm -1 ·s -1 ) φ: Porosity (%) ρ: Density (g·cm -3 ) τ: Tortuousity.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.