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

Abstract

A numerical method is presented for the solution of a moving boundary problem of one-dimensional flow in semi-infinite long porous media with 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.

1. 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–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 low-permeable 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–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 low-permeable 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–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–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.

2. 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 one-dimensional flow in the porous medium is given as follows [20, 21, 26, 29]: 𝜕²𝑃_𝐷𝜕𝑥²_𝐷=𝜕𝑃_𝐷𝜕𝑡_𝐷,0≤𝑥_𝐷≤𝛿𝑡_𝐷, (1) 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: 𝑃_{𝐷𝑡𝐷=0}𝛿(0)=0,=0. (2)(3)

The inner boundary condition with a constant production rate is as follows: 𝜕𝑃_𝐷𝜕𝑥_{𝐷𝑥𝐷=0}=−1+𝜆_𝐷, (4) where 𝜆_𝐷 is the dimensionless TPG.

The moving boundary conditions are as follows [20, 21, 26, 29]: 𝜕𝑃_𝐷𝜕𝑥_{𝐷𝑥𝐷=𝛿(𝑡𝐷)}=−𝜆_𝐷,𝑃_{𝐷𝑥𝐷=𝛿(𝑡𝐷)}=0. (5)(6)

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]: 𝜕𝛿𝜕𝑡_𝐷=1𝜆_𝐷⋅𝜕²𝑃_𝐷𝜕𝑥²_{𝐷𝑥𝐷=𝛿(𝑡𝐷)}. (7)

3. 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 one-dimensional 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–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: 𝑦_𝐷𝑥_𝐷,𝑡_𝐷=𝑥_𝐷𝛿𝑡_𝐷,0≤𝑥_𝐷≤𝛿𝑡_𝐷. (8)

By the spatial coordinate transformation, the following two end points in the one-dimensional spatial coordinate of 𝑥_𝐷 can be, respectively, transformed as 𝑦_𝐷0,𝑡_𝐷=0𝛿𝑡_𝐷𝑦=0,_𝐷𝛿𝑡_𝐷,𝑡_𝐷=𝛿𝑡_𝐷𝛿𝑡_𝐷=1. (9)

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 𝑃_𝐷(𝑦_𝐷,𝑡_𝐷) 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: 𝜕𝑃_𝐷𝜕𝑥_𝐷=𝜕𝑃_𝐷𝜕𝑦_𝐷⋅1𝛿,𝜕²𝑃_𝐷𝜕𝑥²_𝐷=𝜕²𝑃_𝐷𝜕𝑦²_𝐷⋅1𝛿²,𝜕𝑃_𝐷𝜕𝑡_𝐷=𝜕𝑃_𝐷𝜕𝑡_𝐷+𝜕𝑃_𝐷𝜕𝑦_𝐷⋅−𝑥_𝐷𝛿²⋅𝜕𝛿𝜕𝑡_𝐷=𝜕𝑃_𝐷𝜕𝑡_𝐷−𝜕𝑃_𝐷𝜕𝑦_𝐷⋅𝜕𝛿𝜕𝑡_𝐷⋅𝑦_𝐷𝛿. (10)(11)(12)

Substituting (10)–(12) into (1) yields 𝜕²𝑃_𝐷𝜕𝑦²_𝐷⋅1𝛿²=𝜕𝑃_𝐷𝜕𝑡_𝐷−𝜕𝑃_𝐷𝜕𝑦_𝐷⋅𝜕𝛿𝜕𝑡_𝐷⋅𝑦_𝐷𝛿,0≤𝑦_𝐷≤1. (13)

Substituting (10)–(12) into (4)–(7), respectively, yields 𝜕𝑃_𝐷𝜕𝑦_𝐷⋅1𝛿_𝑦𝐷=0=−1+𝜆_𝐷,𝜕𝑃_𝐷𝜕𝑦_𝐷⋅1𝛿_𝑦𝐷=1=−𝜆_𝐷,𝑃_{𝐷𝑦𝐷=1}=0,𝜕𝛿𝜕𝑡_𝐷=1𝜆_𝐷⋅𝜕²𝑃_𝐷𝜕𝑦²_{𝐷𝑦𝐷=1}⋅1𝛿². (14)(15)(16)(17)

Substituting (17) into (13) yields 𝜕²𝑃_𝐷𝜕𝑦²_𝐷⋅1𝛿²=𝜕𝑃_𝐷𝜕𝑡_𝐷−𝜕𝑃_𝐷𝜕𝑦_𝐷⋅1𝜆_𝐷⋅𝜕²𝑃_𝐷𝜕𝑦²_{𝐷𝑦𝐷=1}⋅𝑦_𝐷𝛿³; (18) that is, 𝜕²𝑃_𝐷𝜕𝑦²_𝐷𝜕𝑃⋅𝛿=^𝐷𝜕𝑡_𝐷⋅𝛿³−𝜕𝑃_𝐷𝜕𝑦_𝐷⋅𝑦_𝐷𝜆_𝐷⋅𝜕²𝑃_𝐷𝜕𝑦²_{𝐷𝑦𝐷=1}. (19)

The equivalent form of (14) is as follows: 𝜕𝑃𝛿=−^𝐷𝜕𝑦_{𝐷𝑦𝐷=0}⋅11+𝜆_𝐷. (20)

Substituting (20) into (19) to cancel the variable 𝛿 yields 𝜕²𝑃_𝐷𝜕𝑦²_𝐷⋅𝜕𝑃_𝐷𝜕𝑦_{𝐷𝑦𝐷=0}−11+𝜆_𝐷²𝜕𝑃_𝐷𝜕𝑦_{𝐷𝑦𝐷=0}³⋅𝜕𝑃_𝐷𝜕𝑡_𝐷−1+𝜆_𝐷𝜆_𝐷⋅𝑦_𝐷⋅𝜕𝑃_𝐷𝜕𝑦_𝐷⋅𝜕²𝑃_𝐷𝜕𝑦²_{𝐷𝑦𝐷=1}=0. (21)

From (8), (2) can be transformed as 𝑃_𝐷𝑦_𝐷,𝑡_{𝐷𝑡𝐷=0}=0. (22)

From (14) and (15), the following equation can be deduced: 𝜕𝑃_𝐷𝜕𝑦_{𝐷𝑦𝐷=0}=1+𝜆_𝐷𝜆_𝐷𝜕𝑃_𝐷𝜕𝑦_{𝐷𝑦𝐷=1}. (23)

Equations (21)–(23) and (16) form a closed system of partial differential equations with the fixed boundary conditions with respect to 𝑃_𝐷(𝑦_𝐷,𝑡_𝐷), 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 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: 𝑃_𝐷^𝑗+1_𝑖+1−2𝑃_𝐷^𝑗+1_𝑖+𝑃_𝐷^𝑗+1_𝑖−1Δ𝑦_𝐷²⋅𝑃_𝐷^𝑗+1₁−𝑃_𝐷^𝑗+1₀Δ𝑦_𝐷−11+𝜆_𝐷²𝑃_𝐷^𝑗+1₁−𝑃_𝐷^𝑗+1₀Δ𝑦_𝐷³⋅𝑃_𝐷^𝑗+1_𝑖−𝑃^𝑗_𝐷𝑖Δ𝑡_𝐷−1+𝜆_𝐷𝜆_𝐷⋅𝑖⋅Δ𝑦_𝐷⋅𝑃_𝐷^𝑗+1_𝑖+1−𝑃_𝐷^𝑗+1_𝑖Δ𝑦_𝐷⋅𝑃_𝐷^𝑗+1_𝑁−2−2𝑃_𝐷^𝑗+1_𝑁−1Δ𝑦_𝐷²(),=0𝑖=1,2,…,𝑁−2 (24) where 𝑗 denotes the index of time step.

The difference equation of (16) is as follows: 𝑃_𝐷^𝑗+1_𝑁=0. (25)

Then, by (25), the difference equation corresponding to the (𝑁−1)th spatial grid can be expressed as follows: −2𝑃_𝐷^𝑗+1_𝑁−1+𝑃_𝐷^𝑗+1_𝑁−2Δ𝑦_𝐷²⋅𝑃_𝐷^𝑗+1₁−𝑃_𝐷^𝑗+1₀Δ𝑦_𝐷−11+𝜆_𝐷²𝑃_𝐷^𝑗+1₁−𝑃_𝐷^𝑗+1₀Δ𝑦_𝐷³⋅𝑃_𝐷^𝑗+1_𝑁−1−𝑃^𝑗_𝐷𝑁−1Δ𝑡_𝐷+1+𝜆_𝐷𝜆_𝐷⋅(𝑁−1)⋅Δ𝑦_𝐷⋅𝑃_𝐷^𝑗+1_𝑁−1Δ𝑦_𝐷⋅𝑃_𝐷^𝑗+1_𝑁−2−2𝑃_𝐷^𝑗+1_𝑁−1Δ𝑦_𝐷²=0. (26)

The difference equation of (23) is as follows: 𝑃_𝐷^𝑗+1₁−𝑃_𝐷^𝑗+1₀Δ𝑦_𝐷=−1+𝜆_𝐷𝜆_𝐷𝑃_𝐷^𝑗+1_𝑁−1Δ𝑦_𝐷. (27)

From (22), the initial values of the unknown dimensionless formation pressure are as follows: 𝑃⁰_𝐷𝑖=0,𝑖=0,1,2,…,𝑁−1. (28)

From (24), (26), and (27), it can be seen that there are 𝑁 difference equations at the (𝑗+1)th time step, which contain 𝑁 unknown variables 𝑃_𝐷^𝑗+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 𝑃_𝐷^𝑗+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 𝑃_𝐷^𝑗+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 𝑃_𝐷(𝑦_𝐷,𝑡_𝐷) can be obtained.

The difference equation of (8) is 𝑥^𝑗+1_𝐷𝑖=𝑦_𝐷𝑖⋅𝛿^𝑗+1. (29)

The difference equation of (20) is 𝛿^𝑗+1𝑃=−^𝐷𝑗+11−𝑃_𝐷^𝑗+1₀Δ𝑦_𝐷⋅11+𝜆_𝐷. (30)

Substituting (30) into (29) yields 𝑥^𝑗+1_𝐷𝑖𝑃=−𝑖⋅^𝐷𝑗+11−𝑃_𝐷^𝑗+1₀1+𝜆_𝐷. (31)

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 𝑃_𝐷(𝑦_𝐷,𝑡_𝐷) 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.

4. 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]: 𝑃_𝐷𝑥_𝐷,𝑡_𝐷=21+𝜆_𝐷×𝜃𝑡_𝐷^1/2𝑒^{−𝑥2𝐷/4/𝑡𝐷}+𝑥_𝐷/2𝜋^1/2erf𝑥_𝐷/2/𝑡_𝐷𝑒^−𝜃2+𝜋^1/2(𝜃)−𝑥𝜃erf2^𝐷,𝑥_𝐷∈[0,𝛿], (32) where 𝜃 can be determined by the following equation with respect to the TPG [29]: 𝑒^−𝜃2𝑒^−𝜃2+𝜋^1/2(𝜃)=𝜆𝜃erf1+^𝐷𝜆_𝐷. (33)

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]: 𝛿=2⋅𝜃⋅𝑡_𝐷. (34) When 𝜆_𝐷 = 0, the exact analytical solution of the mathematical model for the one-dimensional Darcy’s flow in semi-infinite long porous media for the case of a constant flow rate at the inner boundary is as follows [29]: 𝑃_𝐷𝑥_𝐷,𝑡_𝐷=−𝑥_𝐷+2𝑡_𝐷𝜋exp−𝑥²_𝐷4𝑡_𝐷+𝑥_𝐷𝑥erf^𝐷2𝑡_𝐷,𝑥_𝐷∈[0,∞]. (35)

4.2. Assignment of the Values of Dimensionless TPG 𝜆_𝐷

The formula for the dimensionless TPG 𝜆_𝐷 is as follows [29]: 𝜆_𝐷=𝑘𝜇𝜆𝜈_𝑤, (36) 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.

Table 1 

Calculation of dimensionless TPG 𝜆_𝐷 through experimental data [2] of eight samples.

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.

4.3. 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.

Figure 1 

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

Figure 2 

Comparison of transient production pressure of wellbore under different values of dimensionless TPG.

Figure 3 

Comparison of formation pressure distribution when 𝑡𝐷 = 5,000 under different values of dimensionless TPG.

From Figures 1–3, it can be seen that the numerical solutions are in good agreement with the exact analytical 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).

5. 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.

Conflict of Interests

Jun Yao, Wenchao Liu, and Zhangxin Chen 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 project (Grant no. 11102237) sponsored by the Natural Science Foundation of China (NSFC) and the Program for Changjiang Scholars and Innovative Research Team in University (Grant no. IRT1294). In particular, Wenchao Liu also would like to express his great gratitude to the China Scholarship Council (CSC) for the generous financial support of the research.