Numerical Solution of a Moving Boundary Problem of One-Dimensional Flow in Semi-Infinite Long Porous Media with Threshold Pressure Gradient

A numerical method is presented for the solution of a moving boundary problem of one-dimensional flow in semi-infinite long porousmediawith threshold pressure gradient (TPG) for the case of a constant flow rate at the inner boundary. In order to overcome the difficulty in the space discretization of the transient flow region with a moving boundary in the process of numerical solution, 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.Then a stable, fully implicit finite difference method is adopted to obtain its numerical solution. Finally, numerical results of transient distance of the moving boundary, transient production pressure of wellbore, and formation pressure distribution are compared graphically with those from a published exact analytical solution under different values of dimensionless TPG as calculated from actual experimental data. Comparison analysis shows that numerical solutions are in good agreement with the exact analytical solutions, and there is a big difference of model solutions between Darcy’s flow and the fluid flow in porous media with TPG, especially for the case of a large dimensionless TPG.


Introduction
With the increase of international oil price, unconventional reservoirs such as low-permeable reservoirs, heavy oil reservoirs, and shale gas reservoirs [1] have become new development targets in the field of petroleum engineering in modern times.The relevant research on the kinematic principles of the fluid flow in unconventional reservoirs is a recently hot topic.Abundant experimental and theoretical analyses [2][3][4][5][6][7][8][9][10][11][12][13][14][15] have demonstrated that the fluid flow in low-permeable reservoirs and the Bingham non-Newtonian fluid flow in heavy oil reservoirs do not obey the conventional Darcy's law; there exists a threshold pressure gradient (TPG) or the yield stress.That is, the fluid flow happens only if the formation pressure gradient exceeds the TPG or the shear stress for the fluid is larger than the yield stress.
Many literatures have demonstrated the significance of taking into account the TPG for the fluid flow in unconventional reservoirs.For example, Zhu et al. [16] conducted experimental investigation of gas flow in water-bearing tight gas reservoirs with TPG and analytical investigation of mathematical model of low-velocity gas flow; their calculation results indicate that the peripheral reserves of the wellbore are difficult to deploy, and the reservoir energy is mainly consumed near the wellbore due to existence of TPG, which is really unlike Darcy's flow.Zhu et al. [17] presented a method for improving history-matching precision of low-permeable reservoir numerical simulation by considering TPG, and relevant applications in Units X10 and X11 of Daqing Oilfield verify the obviously improved history-matching precision of reservoir numerical simulation.Yin and Pu [18] performed study on the surfactant flooding simulation for lowpermeable oilfield in the condition of TPG; the enhanced matching degree between theoretical model and field practice was verified through a pilot test of surfactant flooding in Chao 45 Block of Daqing Oilfield.
Due to the derivative discontinuity for the modified Darcy's law [2], the mathematical model for this physics, which the oilfield-development technologies such as well testing and reservoir numerical simulation for unconventional reservoirs involve, should be built as a moving boundary problem.Some analytical or numerical investigations [19][20][21][22][23][24][25][26][27][28][29][30] have demonstrated the big difference between the mathematical models with considering moving boundary conditions and those without considering moving boundary conditions.First of all, pressure distribution curves calculated from the mathematical model with moving boundary conditions due to the TPG show compact support [29,30], which is similar to the characteristics of power-law non-Newtonian fluid flow [20,22,25] in porous media or the fluid flow in nanoporous media [1].Wu et al. [20] incorporated moving boundary conditions for theoretical study of the radial flow and displacement of a Bingham fluid in porous media with TPG by integral method and also presented a confirmed reservoir well-test-analysis method.Feng et al. [24] studied the model of unsteady radial flow in lowpermeable gas reservoirs considering the TPG by using Green function method with numerical approximation; the moving boundary can represent the single-well control radius; by using the data from Sichuan Gas Field, computation results show that the mathematical model with moving boundary conditions could correctly reflect the seepage mechanics and production performance for low-permeable gas reservoirs.Wang et al. [27,28] studied the effect of moving boundary on the transient radial flow in low-permeable reservoirs with TPG and also demonstrated the significance of taking into account the effect of the moving boundary due to the TPG for engineering applications.
The moving boundary problem of fluid flow in porous media with TPG is different from the classical Stefan problem in heat conduction theory [29]: the velocity of the moving boundary is proportional to the second derivative of the unknown pressure function with respect to the distance at this moving boundary.However, due to the strong nonlinearity of the moving boundary problem, it is very difficult to obtain its exact analytical solution, whereas, recently, in the new published paper [29], exact analytical solutions for the non-Stefan moving boundary problems of one-dimensional flow in semi-infinite long porous media with TPG are presented through a similarity transformation method.Before the occurrence of the exact analytical solutions, approximate analytical solutions [19][20][21]26] and numerical solutions [22-25, 27, 28, 30] have been the main research tools for solving the non-Stefan moving boundary problems of the fluid flow in porous media with TPG.However, these solutions generally lack a strict verification with the related exact analytical solutions.Furthermore, some problems of fluid flow in porous media with threshold pressure are much more complicated and nonlinear such as the consideration of nonlinear terms in the governing equations [31] and the coupling of stress sensitive effect in low-permeable porous media [32], and then it will be extremely hard to obtain the exact analytical solutions.Therefore, it becomes very necessary to develop a verified numerical method by the published exact analytical solution [29] in order to solve more complicated moving boundary problems of fluid flow in porous media with TPG.
The objective of this paper is to present a simple and novel method for numerical solution of the moving boundary problem of one-dimensional flow in semi-infinite long porous media with TPG for the case of a constant flow rate at the inner boundary.First, this moving boundary problem is equivalently transformed into a closed partial differential equation system with fixed boundary conditions by the spatial coordinate transformation method [33][34][35].Then a stable, fully implicit finite difference method [36] is adopted to obtain its numerical solution.Finally, the accuracy of the numerical solution is verified through graphically comparing with the published exact analytical solution.The numerical method presented here is applicable to the moving boundary problems of multidimensional flow in porous media with TPG.

Mathematical Model
The problem considered [29] involves the one-dimensional flow in a semi-infinite long porous medium 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; and the fluid and porous medium are slightly compressible.
The continuous (mass balance) equation for the onedimensional flow in the porous medium is given as follows [20,21,26,29]: where   is the dimensionless formation pressure,   is the dimensionless distance;   is the dimensionless time, and  is the dimensionless distance of the moving boundary.The initial conditions are as follows: The inner boundary condition with a constant production rate is as follows: where   is the dimensionless TPG.
The moving boundary conditions are as follows [20,21,26,29]: Equations ( 1)-( 6) form the mathematical model for the one-dimensional flow in semi-infinite long porous media with TPG for the case of a constant flow rate at the inner boundary; from (1), ( 5), and ( 6), the velocity of the moving boundary can be deduced as follows [29]:

Numerical Solution of the Problem
It is well known that numerical solution needs the space discretization for the flow region; however, due to the existence of the moving boundary in the model, the onedimensional flow region is not fixed but expands outside continuously with the time increasing according to (7).In order to overcome the difficulty in the space discretization for the transient flow region with the moving boundary, a spatial coordinate transformation method is introduced [33][34][35].Then the mathematical model with the moving boundary conditions can be transformed into a mathematical model with fixed boundary conditions.The formula of the spatial coordinate transformation is as follows: By the spatial coordinate transformation, the following two end points in the one-dimensional spatial coordinate of   can be, respectively, transformed as By (8), the dimensionless pressure   (  ,   ) at the flow interval of the moving boundary model [0, (  )] in the one-dimensional spatial coordinate of   can be equivalently transformed to the function P (  ,   ) at the fixed interval [0, 1] in the one-dimensional spatial coordinate of   , and then the differential variables in the mathematical model can be transformed, respectively, as follows: Substituting ( 10)-( 12) into (1) yields Substituting ( 10)-( 12) into ( 4)- (7), respectively, yields Substituting ( 17) into (13) yields that is, The equivalent form of ( 14) is as follows: Substituting ( 20) into (19) to cancel the variable  yields From ( 8), ( 2) can be transformed as From ( 14) and ( 15), the following equation can be deduced: Equations ( 21)-( 23) and ( 16) form a closed system of partial differential equations with the fixed boundary conditions with respect to P (  ,   ), which is equivalent to the mathematical model for the moving boundary problem of the one-dimensional flow in semi-infinite long porous media with TPG for the case of a constant flow rate at the inner boundary.From (21), it can be seen that the partial differential equations have strong nonlinearity, which indirectly reflects Mathematical Problems in Engineering the strong nonlinearity of the original moving boundary model.Owing to the strong nonlinearity of the transformed model, a stable, fully implicit finite difference method [36] is adopted to numerically solve the relevant nonlinear partial differential equations.
The one-dimensional unit interval is discretized as  spatial grid subintervals with the same length, and then the length of every subinterval Δ  is equal to 1/; the distance at the th spatial grid   is  ⋅ Δ  ; the time step is assumed to be Δ  ; replace the first derivative by first-order forward difference and the second derivative by second-order central difference, and finally the fully implicit difference equations [36] corresponding to (21) can be expressed as follows: where  denotes the index of time step.The difference equation of ( 16) is as follows: Then, by (25), the difference equation corresponding to the (−1)th spatial grid can be expressed as follows: The difference equation of ( 23) is as follows: From (22), the initial values of the unknown dimensionless formation pressure are as follows:
From ( 24), (26), and ( 27), it can be seen that there are  difference equations at the ( + 1)th time step, which contain  unknown variables P+1   ( = 0, 1, 2, . . .,  − 1); due to the strong nonlinearity of these difference equations, the Newton-Raphson iteration method [36] is adopted to obtain their numerical solutions; when P+1   ( = 0, 1, 2, . . .,  − 1) are numerically solved for, let  be equal to  + 1, and in the same manner the numerical solutions of  difference equations with respect to  unknown variables P+2   ( = 0, 1, 2, . . .,  − 1) at the ( + 2)th time step can also be numerically solved for; the rest can be deduced by analogy; then the numerical solutions of the nonlinear partial differential equations with respect to P (  ,   ) can be obtained.
The difference equation of ( 8) is The difference equation of ( 20) is Substituting ( 30) into ( 29) yields From (31), in the process of numerical solutions at every time step, the one-dimensional spatial coordinate of   can be transformed as the one of   , and then the numerical solutions of P (  ,   ) can be transformed as the numerical solutions of   (  ,   ).From (31), it can also be seen that it is just the spatial coordinate transformation, that is (8), that lets the time-dependent space discretization in the spatial coordinate of   be transformed as a time-independent space discretization in the spatial coordinate of   , which makes the numerical solutions by the finite difference method more applicable and simple.

Verification of Numerical Solutions
4.1.Exact Analytical Solutions.The exact analytical solution of the mathematical model for the one-dimensional flow in semi-infinite long porous media with TPG for the case of a constant flow rate at the inner boundary is as follows [29]: where  can be determined by the following equation with respect to the TPG [29]: It has been proven [29] that, for a given value of   , there is one and only one positive root of (33).
The distance of the moving boundary can be expressed as follows [29]: When   = 0, the exact analytical solution of the mathematical model for the one-dimensional Darcy's flow in semiinfinite long porous media for the case of a constant flow rate at the inner boundary is as follows [29]: 4.2.Assignment of the Values of Dimensionless TPG   .The formula for the dimensionless TPG   is as follows [29]: where  is the permeability of porous media,  is the TPG,  is the fluid viscosity, and ]  is the average seepage velocity.
Here, selected Prada and Civan' s actual experimental data of eight samples [2] for measuring the one-dimensional flow in three types of porous media (brown sandstone, sandpack, and shaly sandstone) with TPG is used to calculate the values of   ; the experimental data and specific calculation process are shown in Table 1.
From Table 1, it can be concluded that the range of the values of dimensionless TPG   is between 0.242 and 0.852 for the eight samples of three types of porous media.Without loss of generality, the values of dimensionless TPG   are set 0.242, 0.553, and 0.852 for the following numerical tests, respectively.

Numerical Tests.
By the numerical method presented in the paper, the mathematical model for the one-dimensional flow in semi-infinite long porous media with TPG for the case of a constant flow rate at the inner boundary is numerically solved, where  = 160, Δ  = 10, and   = 0.242, 0.553, and 0.852, respectively.By (33), the values of  corresponding to the three values of   can be computed by Newton-Raphson iteration method [36] as 1.005, 0.7732, and 0.6599 for the exact analytical solutions.In order to make clear the difference of model solutions between Darcy's flow, when   = 0, and the fluid flow in porous media with the TPG, exact analytical solutions for the Darcy's flow, that is, (35), are also added for the purpose of comparison analysis.Under different values of   , the numerical results are graphically compared with those from the exact analytical solution, that is, (32)- (34).Figures 1-3 show the comparison curves between the numerical solutions and the exact analytical solutions with respect to the transient distance of the moving boundary, the transient production pressure of wellbore, and the formation pressure distribution when   = 5,000, respectively.
From Figures 1-3, it can be seen that the numerical solutions are in good agreement with the exact analytical Mathematical Problems in Engineering solutions for the moving boundary problem.Furthermore, through lots of numerical tests, it is found that the accuracy of the numerical solutions can be further improved by increasing the value of .It demonstrates the correctness and validity of the presented numerical method here for solving the moving boundary problem of the one-dimensional flow in semi-infinite long porous media with TPG.Besides, from Figures 2 and 3, it can also be concluded that the TPG has a big influence on the model solutions; the larger the value of the dimensionless TPG, the larger the difference of model solutions between Darcy's flow and fluid flow in porous media with TPG; in Figure 3, the formation pressure distribution curves corresponding to non-zero dimensionless TPG indicate an instructive characteristics of model solutions: they have compact supports [1], which are different from that of Darcy's flow.Therefore, it is very necessary to take into account the effect of TPG for the fluid flow in porous media with TPG, for mathematical modeling and engineering applications, especially for the case of a large dimensionless TPG that is, (36).

Conclusions
The utility of the presented numerical method can be attributed to its simple approach of spatial coordinate transformation in numerically solving the moving boundary problem effectively by the stable, fully implicit finite difference method.The accuracy of the numerical solutions is verified by graphically comparing with the published exact analytical solutions.Besides, the accuracy can be further improved by increasing the value of , that is, decreasing the length of the spatial grid subintervals.The numerical results also show the big difference of model solutions between Darcy's flow and the fluid flow in porous media with TPG, especially for the case of a large dimensionless TPG.In comparison with previous numerical methods for solving such moving boundary problems, the verified numerical method presented here is more reliable and convenient to implement by the finite difference method in developing well testing software and reservoir simulator.The numerical method will be applicable to the moving boundary problems of multi-dimensional flow in porous media with TPG, which will be our future research topic.The presented research supports theoretical foundations for technologies of well testing and numerical simulation in developing low-permeable reservoirs and heavy oil reservoirs.

Figure 1 :
Figure 1: Comparison of transient distance of moving boundary under different values of dimensionless TPG.

Figure 2 :Figure 3 :
Figure 2: Comparison of transient production pressure of wellbore under different values of dimensionless TPG.