A Non-Linear Flow Model for Porous Media Based on Conformable Derivative Approach

: Prediction of the non-linear ﬂow in porous media is still a major scientiﬁc and engineering challenge, despite major technological advances in both theoretical and computational thermodynamics in the past two decades. Speciﬁcally, essential controls on non-linear ﬂow in porous media are not yet deﬁnitive. The principal aim of this paper is to develop a meaningful and reasonable quantitative model that manifests the most important fundamental controls on low velocity non-linear ﬂow. By coupling a new derivative with fractional order, referred to conformable derivative, Swartzendruber equation and modiﬁed Hertzian contact theory as well as fractal geometry theory, a ﬂow velocity model for porous media is proposed to improve the modeling of Non-linear ﬂow in porous media. Predictions using the proposed model agree well with available experimental data. Salient results presented here include (1) the ﬂow velocity decreases as effective stress increases; (2) rock types of “softer” mechanical properties may exhibit lower ﬂow velocity; (3) ﬂow velocity increases with the rougher pore surfaces and rock elastic modulus. In general, the proposed model illustrates mechanisms that affect non-linear ﬂow behavior in porous media.


Introduction
Ever since Henry Darcy (1865) developed his famous linear flow model (the classical Darcy's law), based on a series of sand pack experiments, the linear flow through porous media has drawn tremendous attention in various scientific and engineering field [1,2]. However, it's a common phenomenon that experiments on low velocity flow in tight porous media, deviate from the Darcy's law and the flow velocity is lower than that predicted from Darcy's law. As stated in the literature, the existence of low velocity non-Darcy flow (or low velocity non-linear flow) in tight porous media (e.g., shale gas/oil reservoirs, coalbed, or tight gas/oil reservoirs) is due to the interaction forces between the fluid and tight pores [3,4]. Many scholars have documented that, there existed threshold Reynolds number or pressure gradient, which could be used to well describe low velocity non-linear flow [5][6][7]. And they concluded that there is no flow in tight porous media when the pressure gradient is beyond the certain value (i.e., threshold pressure gradient). However, Li provided contradictory evidence for the threshold pressure gradient [8]. He suggested that the threshold pressure gradient measured in labs can be probably ascribed to the difficulty in measuring lower flow velocity, and the false phenomenon of the existence of threshold pressure gradient is strengthened by the skin effect.

Mathematical Model
In this section, the analytical low velocity non-linear flow model for porous media is detailed. The conformable derivative is used to develop the Swartzendruber model for description of low velocity non-linear flow in pores, the fractal geometry theory and modified Hertzian contact theory are used to describe the complex pore structure of porous media under stress condition.

Model Assumptions
The following assumptions are made to simplify the flow system:

1.
The porous media is composed by a bundle of capillary bundles and a single capillary with the equivalent radius r is made up of a packing of equivalent spherical grains.
The single phase flow is under isothermal and stress condition, which is fully developed and at steady state. 4.
The deformation of porous media obeys Hertzian contact theory. 5.
During the flow, the fluid has constant viscosity and density.

Conformable Derivative Approach to Swartzendruber Equation
As suggested by decades of literature, the Swartzendruber equation can well describe the non-Darcian flow in tight porous media with low permeability [13,14]. Based on the Swartzendruber equation, the following equation can be written as: where u is the flow velocity in the cross section; K is the permeability of porous media; µ is fluid viscosity; ρ is fluid density, g is the gravitational acceleration, i and I represent hydraulic gradient and threshold hydraulic gradient, respectively. According to Equation (1), the flow velocity in a single capillary with radius r can be written as: where r is the radius of the capillary. As the linear operator does not inherit all the operational behaviors from the typical first derivative, the Swartzendruber equation fails to capture the full range of non-Darcy flow behavior in porous media [13][14][15]. Fortunately, as a well-behaved and efficient method, conformable derivative approach with real order can be applied to address this problem. The conformable derivative of the flow velocity u(i): [0, ∞) → R for all i > 0 with order α ∈ (0, 1] can be defined by [10]: and the conformable derivative at 0 is given by As stated in the literature [10,14,25], the relationship between the conformable derivative and the first derivative can be written as: Equation (3b) shows that the conformable derivative coincides with the classical first derivative with a given differential order α = 1, which means the conformable derivative is a modification of classical derivative in direction and magnitude [25][26][27]. As stated in the literature [25][26][27], the physical interpretation of the conformable derivative is a modification of classical derivative indirection and magnitude. Replacing the first order derivative in Equation (2) with conformable derivative, the Swartzendruber equation can be rewritten as: By solving Equation (4) with Laplace transform and inverse Laplace transform, the flow velocity in the cross section can be determined as: is the Kummer confluent hypergeometric function [23].
Based on Equation (5), the flow rate q in the cross section can be written as: where q is the flow rate in the capillary with the radius r.

Non-Linear Flow Model
The fractal theory is used to develop the non-linear flow model for tight porous media. Based on fractal theory and modified Hertzian contact theory described in detail in [28,29], the total volumetric flow rate Q in the cross section under stress condition can be calculated by integrating the flow rate over the radius ranging from the minimum radius r min to the maximum porous radius r max [28][29][30][31][32][33]: Energies 2018, 11, 2986 4 of 11 with: In Equation (8), N is the number of pores; r 0 and r are the initial pore radius and pore radius under stress condition, respectively. f is the probability density function for pore size distribution; β is the power law index which is related to the structure of pore surface. r max0 and r min0 are the maximum and minimum pore radius of porous media at zero stress, respectively. r max and r min are the stress-dependent maximum and the stress-dependent minimum pore radius of porous media, respectively. p eff is the effective stress, E is rock elastic modulus and v is rock Poisson's ratio. D f is the fractal dimension for pore size distribution which can be determined as [32,33]: where ϕ 0 is the initial porosity of the porous media, and D f0 is the fractal dimension for pore size distribution at zero stress. Substituting Equations (6), (8) and (9) into Equation (7), the flow rate can be rewritten as: As stated in the literature [32][33][34][35], the cross sectional area of a unit cell A can be expressed as: Then, based on Equations (10) and (11), the average flow velocity u av can be written as: It is evident that the flow rate (or average flow velocity) is the function of pore structural parameters, power law index, rock elastic modulus, effective stress, and differential order α as well as hydraulic gradient and threshold hydraulic gradient.

Results and Discussion
This section aims at studying the novel analytical models in detail. In the following, we first compare our results with those from experimental data. Then, in order to analyze essential controls on non-linear flow in tight porous media, the effects of relevant parameters on average flow velocity are studied in detail.
The availability of the proposed model Equation (12) depends on its ability to adequately fit experimental data [15]. To verify our quantitative model, the measured average flow velocity versus hydraulic gradient relationship in [36] and that predicted by our proposed model are compared  Figure 1). In the experiment of Prakash K. et al. [36], the flow velocity tests were conducted on soils at an effective consolidation stress of 6.25 kPa during loading process. In our proposed model, the initial porosity of rock is 15%, the rock elastic modulus is 45 GPa, the rock Poisson's ratio is 0.23 and power law index is 3/4. Furthermore, to ensure the effective consolidation stress is 6.25 kPa, the effective stress assigned is 6.25 kPa. The values of other parameters (e.g., r max0 , r min0 , differential order α and threshold hydraulic gradient I) are listed in the Figure 1. Results displayed in Figure 1 suggest that the proposed model calculated relationship between average flow velocity and hydraulic gradient is in good agreement with that determined by experimental data [36]. Results (Figure 1) also suggest a definitive positive correlation between the average flow velocity and hydraulic gradient.
Energies 2018, 11, x FOR PEER REVIEW 5 of 11 on non-linear flow in tight porous media, the effects of relevant parameters on average flow velocity are studied in detail. The availability of the proposed model Equation (12) depends on its ability to adequately fit experimental data [15]. To verify our quantitative model, the measured average flow velocity versus hydraulic gradient relationship in [36] and that predicted by our proposed model are compared ( Figure 1). In the experiment of Prakash K. et al. [36], the flow velocity tests were conducted on soils at an effective consolidation stress of 6.25 kPa during loading process. In our proposed model, the initial porosity of rock is 15%, the rock elastic modulus is 45 GPa, the rock Poisson's ratio is 0.23 and power law index is 3/4. Furthermore, to ensure the effective consolidation stress is 6.25 kPa, the effective stress assigned is 6.25 kPa. The values of other parameters (e.g., rmax0, rmin0, differential order α and threshold hydraulic gradient I) are listed in the Figure 1. Results displayed in Figure 1 suggest that the proposed model calculated relationship between average flow velocity and hydraulic gradient is in good agreement with that determined by experimental data [36]. Results (Figure 1) also suggest a definitive positive correlation between the average flow velocity and hydraulic gradient.   [37]. In the experiment of Zhang et al. [37], the permeability tests were conducted on gap-graded sands (e.g., sand A, sand B and sand C) to determine the critical hydraulic gradient of piping in sands. Sand A, sand B and sand C were prepared by mixing a coarse sand with particle size of 3-5 mm and a fine sand with particle size of 0.25-0.50 mm at various ratios of 4:1, 5:1 and 6:1 respectively.

Figure 2.
A comparison between the experimental data [37], and results of the proposed model.   Figure 1. A comparison between the experimental data [36], and results of the proposed model. Figure 2 provides a comparison of the relationship between average flow velocity and hydraulic gradient predicted by the proposed model with experimental data of [37]. In the experiment of Zhang et al. [37], the permeability tests were conducted on gap-graded sands (e.g., sand A, sand B and sand C) to determine the critical hydraulic gradient of piping in sands. Sand A, sand B and sand C were prepared by mixing a coarse sand with particle size of 3-5 mm and a fine sand with particle size of 0.25-0.50 mm at various ratios of 4:1, 5:1 and 6:1 respectively. on non-linear flow in tight porous media, the effects of relevant parameters on average flow velocity are studied in detail.
The availability of the proposed model Equation (12) depends on its ability to adequately fit experimental data [15]. To verify our quantitative model, the measured average flow velocity versus hydraulic gradient relationship in [36] and that predicted by our proposed model are compared (Figure 1). In the experiment of Prakash K. et al. [36], the flow velocity tests were conducted on soils at an effective consolidation stress of 6.25 kPa during loading process. In our proposed model, the initial porosity of rock is 15%, the rock elastic modulus is 45 GPa, the rock Poisson's ratio is 0.23 and power law index is 3/4. Furthermore, to ensure the effective consolidation stress is 6.25 kPa, the effective stress assigned is 6.25 kPa. The values of other parameters (e.g., rmax0, rmin0, differential order α and threshold hydraulic gradient I) are listed in the Figure 1. Results displayed in Figure 1 suggest that the proposed model calculated relationship between average flow velocity and hydraulic gradient is in good agreement with that determined by experimental data [36]. Results (Figure 1) also suggest a definitive positive correlation between the average flow velocity and hydraulic gradient.   [37]. In the experiment of Zhang et al. [37], the permeability tests were conducted on gap-graded sands (e.g., sand A, sand B and sand C) to determine the critical hydraulic gradient of piping in sands. Sand A, sand B and sand C were prepared by mixing a coarse sand with particle size of 3-5 mm and a fine sand with particle size of 0.25-0.50 mm at various ratios of 4:1, 5:1 and 6:1 respectively.

Figure 2.
A comparison between the experimental data [37], and results of the proposed model.   respectively. To ensure the porous media simulated in the proposed model exhibits the same physical properties as the porous media tested in the experiments of [37], the initial void ratio e 0 and initial porosity ϕ 0 (i.e., ϕ 0 = e 0 /(1 + e 0 )) applied in the proposed model are the same with that of the sands in the experiments and the effective stress assigned is 0 MPa. The values of other parameters (e.g., r max0 , r min0 , differential order α and threshold hydraulic gradient I) are listed in the Figure 2. It can be seen from Figure 2 that our predicted values agree well with the corresponding experimental data [37].
As the average flow velocity is related to the threshold hydraulic gradient, differential order α, effective stress, rock elastic modulus, power law index, pore structural parameters, and the hydraulic gradient. We will then study the effects of these relevant parameters (e.g., threshold hydraulic gradient, differential order α, effective stress, rock elastic modulus, power law index and initial porosity) on average flow velocity in detail to analyze essential controls on non-linear flow in tight porous media. We plotted average flow velocity vs. hydraulic gradient with different threshold hydraulic gradient I and different differential order α (Figure 3a,b). For sand A, the hydraulic conductivity and void ratio are 2.6 × 10 −3 m/s and 0.56 respectively. The hydraulic conductivity and void ratio of sand B are 3.6 × 10 −3 m/s and 0.60 respectively. In addition, for sand C, the hydraulic conductivity and void ratio are 4.0 × 10 −3 m/s and 0.63 respectively. To ensure the porous media simulated in the proposed model exhibits the same physical properties as the porous media tested in the experiments of [37], the initial void ratio e0 and initial porosity φ0 (i.e., φ0 = e0/(1 + e0)) applied in the proposed model are the same with that of the sands in the experiments and the effective stress assigned is 0 MPa. The values of other parameters (e.g., rmax0, rmin0, differential order α and threshold hydraulic gradient I) are listed in the Figure 2. It can be seen from Figure 2 that our predicted values agree well with the corresponding experimental data [37].
As the average flow velocity is related to the threshold hydraulic gradient, differential order α, effective stress, rock elastic modulus, power law index, pore structural parameters, and the hydraulic gradient. We will then study the effects of these relevant parameters (e.g., threshold hydraulic gradient, differential order α, effective stress, rock elastic modulus, power law index and initial porosity) on average flow velocity in detail to analyze essential controls on non-linear flow in tight porous media. We plotted average flow velocity vs. hydraulic gradient with different threshold hydraulic gradient I and different differential order α (Figure 3a,b). The parameters assigned in the proposed model were listed in the Figure 3. As suggested by the results (Figure 3a), average flow velocity decreases as threshold hydraulic gradient increases. In addition, Figure 3b shows that larger different differential order leads to smaller average flow velocity with smaller hydraulic gradient, however, on the contrary, when the hydraulic gradient increases to a certain value, average flow velocity increases as differential order increases.  The parameters assigned in the proposed model were listed in the Figure 3. As suggested by the results (Figure 3a), average flow velocity decreases as threshold hydraulic gradient increases.
In addition, Figure 3b shows that larger different differential order leads to smaller average flow velocity with smaller hydraulic gradient, however, on the contrary, when the hydraulic gradient increases to a certain value, average flow velocity increases as differential order increases.  (Figure 4), average flow velocity decreases as effective stress increases. This may be attributed to the decrease of pore radius which is resulted from the pore compaction. Therefore, fluid flow behavior in tight porous media is strongly stress-dependent.
Energies 2018, 11, x FOR PEER REVIEW 7 of 11 Figure 4 illustrates the average flow velocity vs. hydraulic gradient with different effective stress. The parameters assigned in the proposed model were listed in the Figure 4. As suggested by the results (Figure 4), average flow velocity decreases as effective stress increases. This may be attributed to the decrease of pore radius which is resulted from the pore compaction. Therefore, fluid flow behavior in tight porous media is strongly stress-dependent. Plotting average flow velocity vs. hydraulic gradient with different rock elastic modulus, power law index and initial porosity ( Figure 5) was also useful. For the necessary calculations, the parameters assigned in the proposed model were listed in the Figure 5. As suggested by the results (Figure 5a), the average flow velocity increases as rock elastic modulus increases. We interpret this result to indicate that the larger rock elastic modulus decreases the contact surface radius of a given particle, leading to reduced pore volume compressibility and larger hydraulic conductivity. Therefore, rock types of "soft" lithology can yield lower flow velocity. As suggested by the results (Figure 5b), the average flow velocity increases as power law index increases. Correspondingly, we suggest that the larger power law index β implies rougher pore surfaces, leading to only a limited number of pores being compressed and the larger hydraulic conductivity. Figure 5c shows the average flow velocity decreases as rock initial porosity decreases, which is expected. Correspondingly, we suggest that the smaller initial porosity implies narrower pore radius, leading to the smaller hydraulic conductivity.  Plotting average flow velocity vs. hydraulic gradient with different rock elastic modulus, power law index and initial porosity ( Figure 5) was also useful. For the necessary calculations, the parameters assigned in the proposed model were listed in the Figure 5. As suggested by the results (Figure 5a), the average flow velocity increases as rock elastic modulus increases. We interpret this result to indicate that the larger rock elastic modulus decreases the contact surface radius of a given particle, leading to reduced pore volume compressibility and larger hydraulic conductivity. Therefore, rock types of "soft" lithology can yield lower flow velocity. As suggested by the results (Figure 5b), the average flow velocity increases as power law index increases. Correspondingly, we suggest that the larger power law index β implies rougher pore surfaces, leading to only a limited number of pores being compressed and the larger hydraulic conductivity. Figure 5c shows the average flow velocity decreases as rock initial porosity decreases, which is expected. Correspondingly, we suggest that the smaller initial porosity implies narrower pore radius, leading to the smaller hydraulic conductivity.
Energies 2018, 11, x FOR PEER REVIEW 7 of 11 Figure 4 illustrates the average flow velocity vs. hydraulic gradient with different effective stress. The parameters assigned in the proposed model were listed in the Figure 4. As suggested by the results (Figure 4), average flow velocity decreases as effective stress increases. This may be attributed to the decrease of pore radius which is resulted from the pore compaction. Therefore, fluid flow behavior in tight porous media is strongly stress-dependent. Plotting average flow velocity vs. hydraulic gradient with different rock elastic modulus, power law index and initial porosity ( Figure 5) was also useful. For the necessary calculations, the parameters assigned in the proposed model were listed in the Figure 5. As suggested by the results (Figure 5a), the average flow velocity increases as rock elastic modulus increases. We interpret this result to indicate that the larger rock elastic modulus decreases the contact surface radius of a given particle, leading to reduced pore volume compressibility and larger hydraulic conductivity. Therefore, rock types of "soft" lithology can yield lower flow velocity. As suggested by the results (Figure 5b), the average flow velocity increases as power law index increases. Correspondingly, we suggest that the larger power law index β implies rougher pore surfaces, leading to only a limited number of pores being compressed and the larger hydraulic conductivity. Figure 5c shows the average flow velocity decreases as rock initial porosity decreases, which is expected. Correspondingly, we suggest that the smaller initial porosity implies narrower pore radius, leading to the smaller hydraulic conductivity.

Conclusions
In this study, we developed a novel non-linear flow model for tight porous media. The model is based on Swartzendruber equation and conformable derivative approach and as well as the modified Hertzian contact theory and fractal geometry, and allowed us to analyze essential controls on nonlinear flow in tight porous media. An advantage of this model is that it lacks empirical constants, and, more importantly, every parameter in the model has specific physical significance. Predictions from the proposed analytical model exhibit similar variation trends as experimental data, suggesting validity of the model to predict the average flow velocity. Moreover, we analyzed resulting model predictions in detail, to confirm the model's robustness in this context. Results of the new model show the following salient conclusions: 1. The proposed models indicate that average flow velocity in tight porous media is a function of microstructural parameters of the pore space, rock lithology and differential order α as well as hydraulic gradient and threshold hydraulic gradient. 2. The parametric study reveals that average flow velocity increases with the rougher pore surfaces and rock elastic modulus, and decreases with increasing effective stress. "Softer" rock lithology may yield lower average flow velocity. 3. This non-linear model presented here considers microstructural parameters of pore space and rock lithology; we have shown that its forecasted values are robust, at least compared to experimental data, and thus may be useful for performance predictions of non-linear flow behavior in tight porous media. Results also reveal more information about the details of specific Average flow velocity (m/s)

Conclusions
In this study, we developed a novel non-linear flow model for tight porous media. The model is based on Swartzendruber equation and conformable derivative approach and as well as the modified Hertzian contact theory and fractal geometry, and allowed us to analyze essential controls on non-linear flow in tight porous media. An advantage of this model is that it lacks empirical constants, and, more importantly, every parameter in the model has specific physical significance. Predictions from the proposed analytical model exhibit similar variation trends as experimental data, suggesting validity of the model to predict the average flow velocity. Moreover, we analyzed resulting model predictions in detail, to confirm the model's robustness in this context. Results of the new model show the following salient conclusions: 1.
The proposed models indicate that average flow velocity in tight porous media is a function of microstructural parameters of the pore space, rock lithology and differential order α as well as hydraulic gradient and threshold hydraulic gradient.

2.
The parametric study reveals that average flow velocity increases with the rougher pore surfaces and rock elastic modulus, and decreases with increasing effective stress. "Softer" rock lithology may yield lower average flow velocity.

3.
This non-linear model presented here considers microstructural parameters of pore space and rock lithology; we have shown that its forecasted values are robust, at least compared to experimental data, and thus may be useful for performance predictions of non-linear flow behavior in tight porous media. Results also reveal more information about the details of specific parameters (and therefore mechanisms) that affect non-linear flow behavior in porous media. The new model presented in this work can be used to depict the non-linear flow in tight porous media, and may provide meaningful applications for design and development of tight reservoirs. In addition, as the model takes effective stress into account, it is also useful for performance predictions of the coupled flow deformation behavior (stress sensitivity) in tight porous media.

Conflicts of Interest:
The authors declare no conflicts of interest. Initial maximum values at zero stress min minimum values min0

Nomenclature
Initial minimum values at zero stress