Abstract

The concentric two-zone composite reservoir model is a boundary value problems (BVPs) of modified Bessel equations. In this paper, we propose a constructive method to solve the BVPs for the system of modified Bessel equations with Robin mixed outer boundary condition and apply it to solve a two-zone fractal composite reservoir seepage model with stress-sensitivity formation. By using Pedrosa variable substitution, regular perturbation technique, Laplace transform, and Stehfest numerical inversion technique, the unified expression for the solutions of the reservoir model with three outer boundary (infinite, impermeable, and constant pressure) conditions is constructed. Type curves of bottom-hole pressure and pressure derivative are drawn, and sensitivity analysis of reservoir parameters are carried out. In comparison with the traditional approach, the solutions of this model are simple and regular, with continued fraction form, the constructive method is efficient and easy to operate. The application of this method avoids the complicated and trivial derivative operation and the use of Cramer’s rule to solve the system of linear equations. It can help to better understand the relationship between the solutions of the reservoir model and the inner and outer boundary conditions. The constructive method can be applied not only to solve the fractal composite reservoir model but also to solve more general reservoir model, BVPs of fluid diffusion, heat conduction, and so on.

1. Introduction

In 1990, Chang and Yortsos [1] proposed the fractal reservoir seepage model for the first time; the fractal porosity and permeability in the form of power function were given. Aprilian et al. [2] extended Chang’s model and applied it to the interpretation of interference well testing in geothermal deposits, and the fractal dimension of mineral deposits is obtained. Domestic scholars Liu [3] and Tong and Wang [4] further expanded and analyzed the fractal dimension mathematically. According to the structural characteristics of porous media, Kong [5] established a fractal diffusion equation seepage model in different coordinate systems by introducing the capillary diameter fractal dimension and the tortuosity fractal dimension. For the stress-sensitive reservoir, the dependence of permeability on pore pressure makes the flow equation strongly nonlinear, and the quadratic gradient term in seepage equation cannot be ignored. Kikani and Pedrosa [6] defined a permeability modulus to obtain a regular perturbation sequence solution for an infinite radial seepage reservoir. Considering the influence of quadratic gradient term, Wang and Dusseault [7] developed the analytical solution of the porous-media model and obtained the solution of nonlinear diffusion equation by using Laplace transform and the law of conservation of mass. Chakrabarty et al. [8] derived the mathematical model of the nonlinear diffusion equation and analyzed the influence of the quadratic gradient term quantitatively. Bai and Roegiers [9] developed a new dual medium theory by considering the quadratic gradient term in the diffusion-convection equation. Tong and Liu [10], considering the influence of quadratic gradient term, used the Douglas-Jones predictor-corrector method to solve the model and discussed the seepage law of the fractal reservoir with dual-porous media. Dan et al. [11] established the mathematical model of one dimensional fluid-solid coupling nonlinear seepage problem, derived the differential equation of the unsteady seepage field in deformable medium, and solved the analytical solution of the nonlinear seepage equation by using the separation of variables’ method and the Fourier transform. Jiang et al. [12] established a nonlinear seepage model of the fractal dual-porosity reservoir considering pressure sensitivity and solved the model based on the finite element method. Duan et al. [13] established a multizone triple-porosity well test model for pressure transient analysis of horizontal wells extending through a carbonate reservoir; by using the point source function, Green’s function, and coupling of multiple reservoir sections, the model was solved. Fatemeh et al. [14] simulated the two-phase incompressible fluid flow through highly heterogeneous porous media by using the multiscale sinite volume (MsFV) method. Comparing the total number of six boundary conditions of two general types, the effects of the localization assumption on the accuracy of the MsFV are investigated.

From the previously performed studies, we know that most of the analytical solving methods for seepage models adopted methods such as Green’s function method, separation of variables, and integral transformation, as well as numerical algorithms such as difference, finite element, and the MsFV method. Many of the solving methods were complex and varied, and the solution formulas were also complex and not unified. Li et al. [15, 16] found that the solutions of the initial-boundary value problems for a class of second-order ODEs have similar structure of continuous fraction; then, they proposed the similar structure theory (SST) and applied it to the seepage model of the oil and gas reservoir [1719]; the similar structure solutions with continued fraction form were also obtained. After years of research, Li et al. [20, 21] proposed the similar constructive method (SCM), by which the semianalytical solution of the seepage equation can be easily obtained without complicated and trivial derivation. By using the SCM, Dong et al. [22] studied how to solve a BVP of three-region composite second-order linear homogeneous ODE; then, they [23] found a simple method for constructing the exact solutions of the nonhomogeneous mixed BVP for sets of n-interval composite second-order ODE with variable coefficients.

With the help of SCM proposed by Prof. Li, we directly give the solutions and steps of solving the BVPs for system of modified Bessel equations in this paper and then apply the solving method to a two-zone fractal composite reservoir model with stress-sensitivity formation. As long as some approximate simplifications and corresponding variable substitutions are made, the reservoir model is identical to the BVPs for system of modified Bessel equations; therefore, the solution of the reservoir model can be constructed directly according to the steps given above, which can avoid the complicated and trivial derivative operation and the solving of systems of linear equations in the traditional solution method. Moreover, the solution of the model is a continued fraction form which is simple, regular, and easy to be constructed; it can help to better understand the relationship between the solutions of the reservoir model and the inner and outer boundary conditions. The general sketch of the problem under our study is shown in Figure 1.

The steps of this work are as follows. First, the SCM for solving the BVPs for system of modified Bessel equations is presented in Section 2. Next, a two-zone fractal composite reservoir model with stress-sensitivity formation is established and simplified in Section 3.1; in Section 3.2, the variable substitution is made and the solutions of the model are constructed following the steps of the SCM in Section 2. After that, the influences of sensitive parameters on bottom-hole pressure and pressure derivative are analyzed. Finally, the conclusions are drawn.

2. The Proposal of the Constructive Method

2.1. The BVPs for System of Modified Bessel Equations

In order to fully illustrate the high efficiency of the SCM in solving the stress-sensitive fractal composite reservoir model, it is necessary to introduce the following BVPs for system of modified Bessel equations:where A, B, C, D, E, F, G, a, b, and c are real constants and , and . If the above BVP described in equations (1)–(6) has a unique solution, then the solution to in interval [a, b] can be presented as [19]and the solution to in interval [b, c] can be presented aswhere is called the similar kernal function (SKF) of in interval [b, c]:and is called the SKF of in interval [a, b]:where the formula is called the generating functions (GFs) of the SKF, which can be generated by the following binary function:where , are two linearly independent solutions of equations (1) and (2). Furthermore, taking the partial derivative of with respect to and x, we havewhich are the notationis linear combination of modified Bessel functions of order n.

2.2. The Constructive Method and Procedure

The steps to construct the solution of BVPs in equations (1)–(6) are as follows:(1)Firstly, find the two linearly independent solutions of the differential equations in equations (1) and (2):where(2)Secondly, assemble the GFs by equations (11)–(14).(3)Substitute into equations (9) and (10) and combine the boundary condition equation (5) with equation (6); we can generate the kernal functions and .(4)Finally, substitute the kernal functions and into equations (7) and (8), the solutions of the BVP in equations (1)–(6) are obtained.

3. The Application of SCM in Fractal Composite Reservoir with Stress-Sensitivity Formation

3.1. The Mathematical Model

In the process of drilling and completion, due to the inflow of drilling fluid or mud pollution, two flow areas with different seepage properties will be formed around the wellbore, which is called the concentric circle composite reservoir. The flow area near the wellbore is called inner zone if the wellbore radius , and the area away from the wellbore is called outer zone if .

The permeability and porosity with fractal reservoir behavior are introduced into stress-sensitive formation, and considering the wellbore storage effect and effective wellbore radius , we havewhere is permeability modulus, is fractal dimension, is Euclidean dimension, is fractal index, and the subscript represent inner and outer zone, respectively. Suppose we assume the following. (1) The reservoir is horizontal, the pore medium is homogeneous and isotropic, the thickness of the production layer is uniform, and the single-phase fluid flows into a well radially. (2) The viscosity of slightly compressible fluid is constant and obeys isothermal Darcy’s law. (3) The influence of gravity and capillary force is ignored. (4) Before well opening, the pressure in each part of the reservoir is equal to the original formation pressure ; after well opening, it is produced with variable production rate . (5) The fluids and reservoir compressibility .

To simplify the model, we need to introduce the following dimensionless variables, as shown in Table 1.

Then, we can establish the following dimensionless fractal composite reservoir radial flow mathematical model with stress-sensitivity formation.

The basic flow equations which obey the law of conservation of mass are

The initial conditions of two zones are

The inner boundary conditions are

The interface conditions are

The three typical outer boundary conditions areInfinite boundary:Constant pressure boundary:Closed boundary:

The mathematical model composed of equations (19)–(25) has strong nonlinearity due to the existence of the quadratic gradient term of the derivative of pressure. To facilitate the solution, we make the following Pedrosa’s substitution [24]:where is the dimensionless dependent variable, and the regular perturbation series with respect to stress-sensitive parameter is defined as the series

If we stipulate and by using the regular perturbation theory, the following series is expanded as

Then, equations (19)–(25) are reduced to the initial-boundary value problem of zero-order approximation equations satisfying the two zone composite fractal reservoir model with stress-sensitivity formation.

The control equations in equation (19) become

The initial conditions in equation (20) become

The inner boundary conditions in equation (21) become

The interface conditions in equation (22) become

The three outer-boundary conditions in equation (23) and (24) are changed as

The equations form equations (29) to (33) is a system of linear partial differential equations; if we apply the Laplace transform on the above equations with respect to time , the following dimensionless BVPs of the modified Bessel equations are obtained:

3.2. Solving of the Constructive Method

It is evident that the above seepage model in equations (34)–(37) is identical to the BVPs in equations (1)–(6), advanced in Section 2.1; if we make the following variable substitution (),

By using the SCM given in Section 2.2, the solution of the seepage model can be easily obtained. According to equation (7), the Laplace space solution of , which is corresponding to the pressure in the inner zone, can be obtained

If , the Laplace space solution of , which is corresponding to bottom-hole pressure, can be obtained:

By equation (8), the Laplace space solution of , which is corresponding to the pressure in the outer zone, can be obtained:

According to equations (9) and (10), the SKF of in the outer zone becomesand the SKF of in the inner zone becomes

In real applications, the following special cases are more common. If , in equation (42), the outer condition is impermeable boundary; then, it is simplified as

If , in equation (42), the outer condition is constant pressure boundary; then, it is simplified as

If the outer boundary is infinite, the SKF is presented as follows:

Therefore, the key to solve the problem is to find the GFs, , and construct the SKFs and ; we only need to follow the steps given in Section 2.2.Step 1: find the two linearly independent solutions of each equation in equation (34):Compared with equation (16), we haveStep 2: secondly, assemble the GFs, , by equation (11)–(14), as follows:Similarly, taking the partial derivatives of , we haveThe notation,is linear combination of modified Bessel functions of the order .Step 3: substitute into equations (42) and (43) and combine the boundary conditions in equations (36) and (37); we can generate the KFs, and .Step 4: finally, substitute the KFs, and , into equations (39) and (41), then we have the solutions of the BVPs in equations (34)–(37).

4. Type-Curve Analysis

Using the following formula and Stehfest numerical inversion technique, the real space solutions of the corresponding reservoir pressure and bottom-hole pressure distribution can be obtained:

With the aid of computer programs, type curves of pressure and pressure derivative of two-zone fractal composite reservoir with stress-sensitivity formation are drawn.

Figure 2 shows the curves of bottom-hole pressure and pressure derivative of the stress-sensitive fractal composite reservoir with infinite outer boundary affected by parameter . It can be seen from the figure that when the value of the parameter is large, the type curves of pressure and pressure derivative converge on a straight line with slope 1 in pure wellbore flow period, and if the fractal dimension , the pressure derivative curves does not coincide with the horizontal line of 0.5 as in the homogeneous reservoir, and both the pressure and pressure derivative curves show an increasing trend with the increase of parameter .

Figure 3 shows the type curves of bottom-hole pressure and pressure derivative of the fractal composite reservoir with the infinite outer boundary affected by stress-sensitivity coefficient . As can be seen from the figure, the pressure does not depend on the parameter during the pure wellbore flow period. With the increase of time, the larger the value of becomes, the more obvious the curves diverges and the more significant the rise of the curves of pressure derivative is. Especially when approaches 0.1, the growth rate of the curves of both the pressure and pressure derivative increases significantly and shows an upward trend.

Figure 4 shows the type curves of bottom-hole pressure and pressure derivative with the infinite outer boundary of the two-zone fractal composite reservoir affected by the formation stress-sensitivity parameter . The figure shows that the type curves of both the bottom-hole pressure and pressure derivative in the first three flow periods (pure wellbore flow period, inner radial flow transition period, and inner radial flow period) are not affected by parameter . With the increase of time, the pressure derivative curves begin to diverge from the 4th flow period; the larger is, the more obvious the pressure derivative curves rise; especially, when it is greater than 0.01, the curves of pressure and pressure derivative increase significantly.

Figure 5 shows the type curves of bottom-hole pressure and pressure derivative with the infinite outer boundary of the two-zone fractal composite reservoir affected by . It can be seen from the figure that the parameter affects almost all the flow periods except the first pure wellbore flow period. The larger the value of is, the lower the position of pressure derivative curves are in the 2nd to 5th flow periods and the later the 6th flow period appears.

Figure 6 shows the type curves of bottom-hole pressure and pressure derivative with the infinite outer boundary of the two-zone fractal composite reservoir under the joint influence of and formation stress-sensitivity coefficient . As can be seen that and have little effect on the first three flow periods, starting from the 4th flow period, the larger the value of is, the lower the pressure derivative curves are at the 4–6th flow periods. If we take as a fixed value 1 or 10, then three groups of values of stress-sensitivity coefficient are taken as 0.013, 0.01, and 0.005, respectively. The results show that the sensitive coefficient has a greater influence on the curves of pressure dynamic when is small and has a smaller influence on the curves of pressure dynamic when is large.

5. Summary and Conclusions

(i)Permeability modulus is a measure of stress-sensitivity coefficient. Permeability and porosity are not only functions of spatial variables, but also functions of pressure.(ii)The introduction of the regularized perturbation method and the perturbation series zero-order approximation method makes the original nonlinear seepage equation with quadratic gradient term can be approximated to a linear seepage equation. After linearization and Laplace transformation, the seepage model in Section 3 and the BVPs in Section 2 are the same problem in mathematics, and the differences are only in the signs. Under the premise of satisfying the engineering requirements, the SCM can be widely applied to the nonlinear reservoir seepage model with the pressure gradient term.(iii)Compared with the BVPs in Section 2, the solution of the new mathematical model can be constructed by substituting the corresponding variables. From the constructive solutions of the model, it can be seen that the inner boundary determines the similar structure, the outer boundary determines the kernel functions, and the three outer boundary solutions can be unified as Robin-type boundary. The solution under different boundary conditions is only a special case.(iv)This paper proves the feasibility and practicability of the SCM, which provides a new way for the solution of the seepage model and the well test analysis method.

Abbreviations

:Reservoir pressure,
:Inital pressure,
:Wellbore pressure,
:Radial distance,
:Wellbore radius,
:Crude oil volume factor
:Skin effect
:Variable production rate,
:Reference production,
:Permeability,
:Total compressibility,
:Wellbore storage coefficient,
:Reservoir thickness,
:Production time, h
:Euclid dimension
:Fractal dimension
:Interface radius, m
:Radius of region, m
:Laplace transform parameter
Greek Symbols
:Fractal parameter
:Permeability modulus
:Dependent variable
:Storativity ratio
:Fluid viscosity,
:Porosity, %
:Mobility ratio
:Fractal diffusion exponent
:Ratio
Subscripts
:Effective
:Dimensionless
:Dimensionless effective variable
:Total
:Effective wellbore
0:Zero-order term
:Number
:Inner zone and outer zone
Superscripts
:Laplace transform
:Inverse Laplace transform
′:Derivative.

Data Availability

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

Disclosure

This manuscript has been submitted as preprint to Elsevier, “The Constructive Solution of the Fractal Composite Reservoir with Stress-Sensitivity, Mathematical Problems in Engineering, 2021.”

Conflicts of Interest

The authors declare that they have no conflicts of interest in connection with the work submitted.

Acknowledgments

This work was supported by the Research Project of Sichuan Education Department (Grant nos. 15ZB0450 and 18ZB0394), Talent Introduction Project of Xihua University 2020 (Grant no. Z201076), Vocational Education Reform and Innovation Project of “Science, Innovation and Education” of the Ministry of Education (Grant no. HBKC217128), and Team and Project Funds of Yibin Vocational and Technical College (Grant nos. ybzysc20bk05, ybzy21cxtd-06, and ZRKY21ZDXM-03).