#### Abstract

The existence of a TPG can generate a relatively high pressure gradient in the process of fluid flow in porous media in low-permeable reservoirs, and neglecting the QPGTs in the governing equations, by assuming a small pressure gradient for such a problem, can cause a significant error in predicting the formation pressure. Based on these concerns, in consideration of the QPGT, a moving boundary model of radial flow in low-permeable reservoirs with the TPG for the case of a constant flow rate at the inner boundary is constructed. Due to strong nonlinearity of the mathematical model, a numerical method is presented: the system of partial differential equations for the moving boundary problem is first transformed equivalently into a closed system of partial differential equations with fixed boundary conditions by a spatial coordinate transformation method; and then a stable, fully implicit finite difference method is used to obtain its numerical solution. Numerical result analysis shows that the mathematical models of radial flow in low-permeable reservoirs with TPG must take the QPGT into account in their governing equations, which is more important than those of Darcy’s flow; the sensitive effects of the QPGT for the radial flow model do not change with an increase of the dimensionless TPG.

#### 1. Introduction

Due to a continuously decreasing crude oil output from conventional reservoirs and a high international gasoline price in recent years, unconventional reservoirs such as low-permeable reservoirs and shale oil and gas reservoirs have become urgent development resources in the petroleum industry. Consequently, considerable attention has been paid to the relevant research on the kinematic principles for the fluid flow behavior in these unconventional reservoirs [1–5] at present. Abundant experimental and theoretical analyses [6–13] have demonstrated that the fluid flow, in low-permeable porous media, does not obey the classical Darcy’s law: the seepage velocity is not proportional to the formation pressure gradient, and there exists a threshold pressure gradient (TPG) *λ*; the fluid flow happens only if the formation pressure gradient is larger than TPG.

Much research on these relevant moving boundary models has been conducted [14–22]. The computed formation pressure distributions corresponding to these moving boundary problems of the fluid flow in the porous media with TPG show big difference from the ones based on Darcy’s law (see Figure 1): the formation pressure gradient is much steeper, it decreases until up to zero at a certain value of a dimensionless distance from a well, that is, the position of a moving boundary, and the pressure distribution curve shows a property of compact support [21]; whereas, for Darcy’s flow problem, the formation pressure drop can propagate to any infinite distance transiently according to the exact analytical solution [20], and the formation pressure gradient is much more smooth. Actually, the pressure distribution difference between Darcy’s flow and fluid flow in the porous media with TPG can be explained through the angle between the dimensionless formation pressure curve and the dimensionless distance at the place of moving boundary, as shown in Figure 1: for the case of the existence of TPG, the tangent of is equal to the derivative of pressure with respect to distance at the place of moving boundary, that is, the TPG; and, as decreases gradually and is equal to zero, it becomes a limit case, which corresponds to Darcy’s flow; that is, TPG is equal to zero. And if a constant value of TPG is given, the tangent of will not change transiently. Therefore, it can be concluded that the formation pressure gradient is relatively much higher for the fluid flow in porous media with large TPG. The main physical reason is that the presence of TPG can make the formation pressure drop propagate more slowly, which causes a larger pressure gradient in a relatively shorter pressure disturbed distance from a production well.

It is well known that the rock porosity and fluid density are dependent on the formation pressure, and their formula are commonly expressed in the nonlinear forms of exponential functions [23]. Therefore, in the mathematical modeling, the deduced governing equation, by substituting these equations of state and a kinematic equation into a mass balance equation, always contains a nonlinear quadratic pressure gradient term (QPGT), where the coefficient of this term is proportional to the fluid compressibility. In general, the nonlinear QPGT is usually neglected by assuming small fluid compressibility or a small formation pressure gradient [24]. The error by neglecting the nonlinear term may be acceptable for most routine engineering applications in the development of conventional reservoirs. However, it has been realized that the linearization is not applicable for large values of time [25, 26], and certain operations such as hydraulic fracturing, high injection or production rates of wellbore, and large pressure pulse testing can generate a high formation pressure gradient.

For conventional models of fluid flow in porous media based on Darcy’s law, Odeh and Babu [27] studied the analytical solutions of the nonlinear model of slightly compressible fluid flow through porous media with consideration of QPGT by a Laplace transform and pointed out that the nonlinear solutions showed the difference in the pressure change for injection and pumping conditions. Wang and Dusseault [28] presented an analytical solution for pore pressure coupled with deformation in a porous medium by taking the quadratic gradient term into account, and the deviations from the existing solutions were identified in cases of high pressure gradients. Chakrabarty et al. [29, 30] gave the analytical dimensionless pressure solutions for radial flow systems by a Laplace transform and concluded that the linearized analysis by neglecting QPGT could lead to serious errors in some cases such as high injection rates. Braeuning et al. [31] studied the effect of QPGT on the variable-rate well-tests, and concluded that the linearization error depended on the magnitude of wellbore damage, pseudoskin, and a nonlinear flow parameter. Tong et al. [25, 26] presented the exact analytical solutions of nonlinear transient flow models and dual-porosity models including QPGTs by the generalized Weber transform and Hankel transform; it was concluded that the effect of the QPGT should be considered in large time well testing in reservoir engineering. Li et al. [32] took the effects of both QPGT and wellbore storage into account to build a mathematical flow model in a fractal multilayer reservoir and presented its analytical solution by a Laplace transform. Dewei et al. [33] obtained the analytical solution of a mathematical model of transient seepage flow including QPGT by a Laplace transform, and they demonstrated that QPGT could not be ignored for the unsteady flow model with consideration of a wellbore skin effect. Bai et al. [24] constructed a nonlinear dual-porosity model incorporating QPGT in fracture space and obtained its analytical solution using a Hankel transform, and their analysis showed that the presented nonlinear model appeared to be suitable for simulating naturally fractured reservoirs, subjected to a high injection or production rate or significant fracture compressibility. Nie et al. [34, 35] studied nonlinear flow models with a QPGT for both a double-porosity reservoir and a triple-porosity reservoir and obtained a solution through a variable substitution for linearization; it was found that the influence of QPGT was very distinct, and the parameter values of nonlinear model explanation law were much accurate than those of linear model explanation. Yao et al. [36] established the mathematical model of transient flow in a double-porosity fractal reservoir by considering QPGT and solved the model by a Laplace transform; it was demonstrated that the relative errors, caused by ignoring QPGT in the constructed model, may amount to 10%. Guo and Nie [37] discussed the origins of nonlinear flow issues in underground formations by using quadratic pressure gradients to deduce diffusion equations and concluded that the nonlinear models could describe fluid flow in underground formations realistically and accurately.

As far as we know, for the moving boundary problems of fluid flow in porous media with TPG, which are based on modified Darcy’s law [8] and have been widely involved in engineering applications for low-permeable reservoirs, the effect of QPGT has not yet been taken into account in their governing equations [14–22]. However, this nonlinear QPGT may play an important role in the temporal and spatial pressure profiles due to a relatively high formation pressure gradient, generated in fluid flow in porous media with TPG. What is more, in modern times, with the development of advanced analysis methods and improved resolutions of pressure measurement machines [29], it is necessary to study quantitatively the effect of QPGT on these moving boundary problems.

Hence, based on these concerns, in this paper, QPGT is incorporated in the governing equation for a basic moving boundary problem of radial flow in an infinite reservoir with TPG for the case of a constant flow rate at the inner boundary. Owing to the existence of the moving boundary and the nonlinearity of partial differential equations, it is really difficult to obtain its exact analytical solution. Here, we adopted a spatial coordinate transform based finite difference method to solve the nonlinear moving boundary model for the radial flow case. The numerical method has been strictly verified for solving such moving boundary problems with good accuracy by a one-dimensional flow case in the recently published literature [21]. Furthermore, the effect of this QPGT on the numerical solutions of the moving boundary problem can be discussed and analyzed quantitatively, by using the numerical results with respect to different values of a dimensionless TPG.

#### 2. Mathematical Modeling by Considering QPGT

The problem considered involves radial flow in an infinite reservoir with TPG for the case of a constant flow rate at the inner boundary; the porous medium is homogeneous, isotropic, and isothermal; the single-phase horizontal flow does not have any gravity effect; for a basic problem, skin effect and wellbore storage are not considered here; finally, the fluid and porous medium are slightly compressible.

The fluid density is as follows [20]:where is the fluid density, kg·m^{−3}; is the initial fluid density, kg·m^{−3}; is the initial pressure, MPa; is the pressure, MPa; and is the compression coefficient of the fluid, MPa^{−1}.

The porosity of the porous medium is as follows [20]:where is the porosity of the porous medium, fraction; is the initial porosity, fraction; and is the compression coefficient of the porosity, MPa^{−1}.

The modified Darcy’s law for the fluid flow in the porous medium with TPG is as follows [8]:where is the permeability of the porous medium, 10^{−3}· m^{2}; is the fluid viscosity, mPa·s; is the radial distance, ; is the seepage velocity, m·s^{−1}; and is the TPG, MPa·m^{−1}.

The continuous equation for the radial flow in the porous medium is as follows [15]:where is time, h and is the moving boundary, m.

Substituting (1)–(3) into (4), the governing (mass balance) equation for the radial flow in the low-permeable reservoir, considering the nonlinear QPGT, can be deduced as follows (see Appendix):where is the total compression coefficient, MPa^{−1}.

The initial condition is as follows:

The inner boundary conditions with constant flow rate are where is the constant flow rate, m^{3}·d^{−1}; is the reservoir thickness, m; is the wellbore radius, m; and is the volume factor, dimensionless.

The moving boundary conditions are:

The moving boundary conditions (8) physically mean that the seepage flow only happens at the area near the wellbore inside the moving boundary, where the formation pressure gradient is larger than TPG; however, outside the moving boundary, the formation pressure gradient is smaller than TPG (equal to zero), so there is no seepage flow behavior, and the formation pressure also keeps the initial pressure; in addition, on the moving boundary, the pressure gradient is just equal to TPG; with the time increasing, the moving boundary also moves outside gradually. The existence of moving boundary is the main difference between the models of fluid flow in porous media with TPG and the classical Darcy’s flow models.

Equations (5)–(8) together form the moving boundary problem of radial flow in the infinite reservoir with TPG for the case of a constant flow rate at the inner boundary, in consideration of QPGT.

Define the following dimensionless variables: where is the dimensionless radial distance, is the dimensionless time, is the dimensionless pressure, is the dimensionless compressibility, is the dimensionless TPG, and is the dimensionless moving boundary. Consider

Equations (10)–(15) form the dimensionless mathematical model, considering QTPG, for the radial flow in the infinite reservoir with TPG for the case of a constant flow rate at the inner boundary.

From (15), we have

Differentiating two sides of (16), with respect to , we have

Substituting (14) into (17) yields

Letting on both sides of (10) and then substituting (14) yield

By (18) and (19), the velocity of the moving boundary can be written as follows:

From (20), it can be seen that the existence of the QPGT (the dimensionless compressibility is not equal to zero) can slow down the velocity of the moving boundary of the radial flow in low-permeable reservoir with TPG.

#### 3. Numerical Method

Due to the existence of the moving boundary in the dimensionless mathematical model, the flow region is not fixed and expands outward continuously with time increasing. In order to overcome this difficulty in the space discretization for the transient flow region with the moving boundary, a spatial coordinate transformation method is used as follows [21, 38]:

Through (21), the dynamic flow region for the moving boundary problem can be transformed to a fixed region , and the dimensionless pressure can be transformed as a new unknown function of two variables and , that is, equivalently; correspondingly, the differential variables can be transformed, respectively, as follows:

Substituting (22)–(24) into (10) yields

Substituting (22)–(24) into (13)–(15) and (20), respectively, yields

From (26), we have

Substituting (29) and (30) into (25) to cancel the variables and yields

From (21), (11) can be transformed into

From (26) and (27), the following equation can be obtained [21]:

Equations (31)–(33) and (28) together form a closed system of partial differential equations with the fixed boundary conditions with respect to , which is equivalent to the dimensionless mathematical model, that is, (10)–(15). The new partial differential equation (31) shows strong nonlinearity, indirectly indicating the strong nonlinearity of the original, untransformed moving boundary problem incorporating the QPGT. Here, a stable, fully implicit finite difference method [21, 39] is used to obtain its numerical solutions. The first derivative is replaced by a first-order forward difference; the second derivative is replaced by the second-order central difference [21], and then the difference equation for (31) can be written as follows:where denotes the total number of spatial grid subintervals with the same length; is the length of a grid subinterval, which is equal to 1/; denotes the index of the spatial grid from the well; denotes the index of a time step; and denotes the time step size [21].

From (28), we have

Then, from (34) and (35), the difference equation corresponding to the th spatial grid can be expressed as follows:

The difference equation of (33) is as follows:

From (32), the initial conditions are obtained as follows:

Equations (34), (36), and (37) together form difference equations at the th time step and also contain unknown variables . The Newton-Raphson iterative method [21, 39] is used to numerically solve these nonlinear difference equations; when are numerically solved for, let be equal to , and, in the same manner, the numerical solutions of difference equations with respect to unknown variables at the th time step can also be numerically solved for; the rest can be deduced by analogy [21, 39]. Eventually, numerical solutions for the transformed nonlinear partial differential equation system with respect to can be obtained.

The difference equation of (21) is

The difference equation of (30) is

Substituting (40) into (39) yields

By (41), numerical solutions of can be transformed as the ones of in the process of numerical solutions at every time step [21]. From (41), it can be further concluded that the presented spatial coordinate transformation lets the time-dependent space discretization in the spatial coordinate of be transformed as a time-independent space discretization in the new spatial coordinate of and then makes the numerical solutions by the finite difference method more applicable and simpler [21].

#### 4. Results and Discussions

##### 4.1. Effect of QPGT under Different Values of TPG

Figures 2–4 show the effect of the QPGT on numerical solutions of the model, with respect to the dimensionless formation pressure distribution, dimensionless transient wellbore pressure, and dimensionless transient distance of moving boundary, respectively, under different values of the dimensionless TPG. From Figures 2–4, it can be clearly seen that with an increase of the dimensionless TPG, the effect of the QPGT on the mathematical model solutions become more and more obvious; the reason can be explained as follows: the bigger the TPG, the steeper the formation pressure gradient, which generates the solution deviations, resulting from neglecting the QPGT in the governing equation, more seriously. Moreover, from Figures 2–4, it can also be indicated that the dimensionless transient wellbore pressure for the case of neglecting the QPGT is larger than that for the case of considering the QPGT; and the dimensionless transient distance of moving boundary for the case of neglecting the QPGT is also larger than that for the case of considering the QPGT. The QPGT can slow down the velocity of the moving boundary in low-permeable reservoir with TPG.

Tables 1 and 2 show the selected data of calculated relative errors from the already computed dimensionless formation pressure distribution and dimensionless transient wellbore pressure between the two cases considering and not considering the QPGTs, under different values of dimensionless TPG, which correspond to Figures 2 and 3, respectively. The relative error is equal to the absolute error divided by the magnitude of the exact value, that is, the solutions of the models considering the QPGT. In order to clearly figure out the change behavior of relative errors, the curves regarding the relative errors by the data from Tables 1 and 2 are also plotted in Figures 5 and 6, respectively.

From Figure 5, it can be clearly indicated that, for any case with the same value of the dimensionless distance, the larger the dimensionless TPG, the larger the relative error of the dimensionless formation pressure; furthermore, the larger dimensionless TPG can make the relative errors grow more quickly with an increase of dimensionless distance. For the cases with the values of the dimensionless TPG 0.034, 0.067, and 0.1, the relative errors have largely exceeded 5% in the whole region.

From Figure 6, it can be clearly indicated that, for any case with the same value of dimensionless time, the larger the dimensionless TPG, the larger the relative error of the dimensionless transient wellbore pressure; and the larger dimensionless TPG can make the relative errors grow more quickly with an increase of dimensionless time. For the cases with the values of the dimensionless TPG 0.067 and 0.1, the relative errors have largely exceeded 5% in the whole time.

From Figures 5 and 6, it can also be concluded that the relative errors corresponding to the case with the value of the dimensionless TPG , which is relatively close to Darcy’s flow, stay a lower level for the change of relative errors of the dimensionless transient wellbore pressure with the increase of dimensionless time, whereas the change of relative errors of the dimensionless formation pressure distribution with the increase of dimensionless distance is still very considerable.

In conclusion, the QPGT has a great effect on the mathematical model solutions of radial flow in low-permeable reservoirs with TPG, especially for a large value of the dimensionless TPG; the larger the dimensionless distance and dimensionless time, the bigger the effect of the QPGT on the dimensionless formation pressure and the dimensionless transient wellbore pressure. In summary, the mathematical models of radial flow in low-permeable reservoirs with TPG should take the QPGT into account in their governing equations, which is more important than those of Darcy’s flow.

##### 4.2. Sensitive Effect of QPGT: The Parameter of Dimensionless Compressibility

Figures 7, 8, and 9 show the sensitive effect of the QPGT on the formation pressure distribution, transient wellbore pressure, and transient distance of moving boundary, respectively. From these figures, it can be seen that the larger the value of the dimensionless compressibility, the smaller the values of the dimensionless formation pressure, the dimensionless transient wellbore pressure, and the dimensionless transient distance of the moving boundary. Moreover, the sensitive effects of the QPGT on the formation pressure distribution, transient wellbore pressure, and transient distance of the moving boundary for the radial flow model do not change with an increase of the dimensionless TPG.

#### 5. Conclusions

The existence of a TPG can lead to relatively steep formation pressure gradients, and, thus, it is not appropriate to neglect the QPGT for the fluid flow in porous media with the TPG. In consideration of the QPGT in the governing equation, a numerical method is presented for investigation of the effect of QPGT on a moving boundary problem of radial flow in low-permeable reservoirs with TPG for the case of a constant flow rate at the inner boundary. Numerical result analysis shows that the QPGT plays an important role in the model solutions with respect to the formation pressure distribution, transient wellbore pressure, and transient distance of the moving boundary. The mathematical models of radial flow in low-permeable reservoirs with TPG should take the QPGT into account in their governing equations, which is more important than those of Darcy’s flow. Besides, the sensitive effects of the QPGT on the numerical solutions for the radial flow model do not change with an increase of the dimensionless TPG. The presented research supports theoretical foundations for further improving technologies of well testing and numerical simulation in developing low-permeable reservoirs.

#### Appendix

Equation (1) can be rewritten as follows:

Differentiating two sides of (A.1) with respect to , we have

Equation (A.2) can be rewritten as follows:

In the same manner as above, from (1) and (2), the following equations can also be deduced as

The left-hand side of (4) can be expanded as follows:

Because and , the small term in the right-hand side of (A.5) can be neglected, and then (A.5) can be rewritten as follows:

In (A.6), the QPGT is retained for the deduction of the governing equation.

Expanding the right-hand side of (4) yields

Substituting (A.4) into the right-hand side of (A.7) yields

Substituting (A.6) and (A.8) into (4), the governing equation in consideration of the QPGT can be obtained as follows:

Equation (A.9) can be equivalently simplified, by canceling the variable in its both sides, as follows:

#### Conflict of Interests

Wenchao Liu and Jun Yao declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

The authors would like to acknowledge the funding by the Projects Grants nos. 11102237 and 51404232, both sponsored by the Natural Science Foundation of China (NSFC), and the Program for Changjiang Scholars and Innovative Research Team inUniversity (Grant no. IRT1294).