# The Constructive Solution of the Fractal Composite Reservoir with Stress-Sensitivity Formation

Authors: Wei Li; Songlin Zhang; Haohan Liu; Shunchu Li
Journal: Mathematical Problems in Engineering
 (2021)

---

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


---

## Body

## 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 [17–19]; 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 Figure1.Figure 1
The general sketch of the problem under the study.The steps of this work are as follows. First, the SCM for solving the BVPs for system of modified Bessel equations is presented in Section2. 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:(1)d2y1dx2+A1xdy1dx−C1xq1−2y1=0,a≤x≤b,(2)d2y2dx2+A2xdy2dx−C2xq2−2y2=0,b≤x≤c,(3)Ey1+1+EFy1′x=a=D,(4)y1|x=b=y2|x=b,(5)y1′|x=b=δy2′|x=b,(6)Gy2+Hy2′x=c=0,where A, B, C, D, E, F, G, a, b, and c are real constants and C1,C2>0,D≠0,G2+H2≠0, and 0<a<b<c. If the above BVP described in equations (1)–(6) has a unique solution, then the solution to y1x in interval [a, b] can be presented as [19](7)y1x=D⋅1E+1/F+Φ1a⋅1F+Φ1a⋅Φx,and the solution to y2x in interval [b, c] can be presented as(8)y2x=D⋅1E+1/F+Φ1a⋅1F+Φ1a⋅φ0,11b,bΦ2bφ1,11a,b−δφ1,01a,b⋅Φ2x,where Φ2x is called the similar kernal function (SKF) of y2x in interval [b, c]:(9)Φ2x=Gφ0,02x,c+Hφ0,12x,cGφ1,02b,c+Hφ1,12b,c,b≤x≤c,and Φ1x is called the SKF of y2x in interval [a, b]:(10)Φ1x=Φ2bφ0,11x,b−δφ0,01x,bΦ2bφ1,11a,b−δφ1,01a,b,a≤x≤b,where the formula φi,jτξ,xi,j=0,1;τ=1,2 is called the generating functions (GFs) of the SKF, which can be generated by the following binary function:(11)φ0,0τξ,x=ξατIvτγτξβτξατKvτγτξβτxατIvτγτxβτxατKvτγτxβτ=−ξxατψυτ,υττξβτ,xβτ,γτ,where xατIvτγτxβτandxατKvτγτxβτ,τ=1,2, are two linearly independent solutions of equations (1) and (2). Furthermore, taking the partial derivative of φi,jτξ,x with respect to ξ and x, we have(12)φ1,0τξ,x=∂∂ξφ0,0τξ,x=−ξατ−1xατατ+βτυτψυτ,υττξβτ,xβτ,γτ−γτβτξβτψυτ+1,υττξβτ,xβτ,γτ,(13)φ0,1τξ,x=∂∂xφ0,0τξ,x=−ξατxατ−1ατ+βτυτψυτ,υττξβτ,xβτ,γτ+γτβτxβτψυτ,υτ+1τξβτ,xβτ,γτ,(14)φ1,1τξ,x=∂2∂x∂ξφ0,0τξ,x=−xξατ−1×ατ+βτυτ2ψυτ,υττξβτ,xβτ,γτ+γτβτατ+βτυτxβτψυτ,υτ+1τξβτ,xβτ,γτ−γτβτατ+βτυτξβτψυτ+1,υττξβτ,xβτ,γτ−γτβτ2xξβτψυτ+1,υτ+1τξβτ,xβτ,γτ,which are the notation(15)ψm,nξ,x,t=KmξtInxt+−1m−n+1ImξtKnxt,is linear combination of modified Bessel functions IntandKnt 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):(16)xατIvτγτxβτ,xατKvτγτxβτ,where(17)ατ=1−Aτ2,βτ=qτ2,γτ;=2Cτqτ,υτ=1−Aτqτ,τ=1,2.(2)
Secondly, assemble the GFsφi,jτξ,xi,j=0,1;τ=1,2 by equations (11)–(14).(3)
Substituteφi,jτξ,x into equations (9) and (10) and combine the boundary condition equation (5) with equation (6); we can generate the kernal functions Φ1x and Φ2x.(4)
Finally, substitute the kernal functionsΦ1x and Φ2x into equations (7) and (8), the solutions of the BVP in equations (1)–(6) are obtained.


## 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:(1)d2y1dx2+A1xdy1dx−C1xq1−2y1=0,a≤x≤b,(2)d2y2dx2+A2xdy2dx−C2xq2−2y2=0,b≤x≤c,(3)Ey1+1+EFy1′x=a=D,(4)y1|x=b=y2|x=b,(5)y1′|x=b=δy2′|x=b,(6)Gy2+Hy2′x=c=0,where A, B, C, D, E, F, G, a, b, and c are real constants and C1,C2>0,D≠0,G2+H2≠0, and 0<a<b<c. If the above BVP described in equations (1)–(6) has a unique solution, then the solution to y1x in interval [a, b] can be presented as [19](7)y1x=D⋅1E+1/F+Φ1a⋅1F+Φ1a⋅Φx,and the solution to y2x in interval [b, c] can be presented as(8)y2x=D⋅1E+1/F+Φ1a⋅1F+Φ1a⋅φ0,11b,bΦ2bφ1,11a,b−δφ1,01a,b⋅Φ2x,where Φ2x is called the similar kernal function (SKF) of y2x in interval [b, c]:(9)Φ2x=Gφ0,02x,c+Hφ0,12x,cGφ1,02b,c+Hφ1,12b,c,b≤x≤c,and Φ1x is called the SKF of y2x in interval [a, b]:(10)Φ1x=Φ2bφ0,11x,b−δφ0,01x,bΦ2bφ1,11a,b−δφ1,01a,b,a≤x≤b,where the formula φi,jτξ,xi,j=0,1;τ=1,2 is called the generating functions (GFs) of the SKF, which can be generated by the following binary function:(11)φ0,0τξ,x=ξατIvτγτξβτξατKvτγτξβτxατIvτγτxβτxατKvτγτxβτ=−ξxατψυτ,υττξβτ,xβτ,γτ,where xατIvτγτxβτandxατKvτγτxβτ,τ=1,2, are two linearly independent solutions of equations (1) and (2). Furthermore, taking the partial derivative of φi,jτξ,x with respect to ξ and x, we have(12)φ1,0τξ,x=∂∂ξφ0,0τξ,x=−ξατ−1xατατ+βτυτψυτ,υττξβτ,xβτ,γτ−γτβτξβτψυτ+1,υττξβτ,xβτ,γτ,(13)φ0,1τξ,x=∂∂xφ0,0τξ,x=−ξατxατ−1ατ+βτυτψυτ,υττξβτ,xβτ,γτ+γτβτxβτψυτ,υτ+1τξβτ,xβτ,γτ,(14)φ1,1τξ,x=∂2∂x∂ξφ0,0τξ,x=−xξατ−1×ατ+βτυτ2ψυτ,υττξβτ,xβτ,γτ+γτβτατ+βτυτxβτψυτ,υτ+1τξβτ,xβτ,γτ−γτβτατ+βτυτξβτψυτ+1,υττξβτ,xβτ,γτ−γτβτ2xξβτψυτ+1,υτ+1τξβτ,xβτ,γτ,which are the notation(15)ψm,nξ,x,t=KmξtInxt+−1m−n+1ImξtKnxt,is linear combination of modified Bessel functions IntandKnt 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):(16)xατIvτγτxβτ,xατKvτγτxβτ,where(17)ατ=1−Aτ2,βτ=qτ2,γτ;=2Cτqτ,υτ=1−Aτqτ,τ=1,2.(2)
Secondly, assemble the GFsφi,jτξ,xi,j=0,1;τ=1,2 by equations (11)–(14).(3)
Substituteφi,jτξ,x into equations (9) and (10) and combine the boundary condition equation (5) with equation (6); we can generate the kernal functions Φ1x and Φ2x.(4)
Finally, substitute the kernal functionsΦ1x and Φ2x 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 radiusrw≤r≤brw, and the area away from the wellbore is called outer zone if brw≤r≤R.The permeability and porosity with fractal reservoir behavior are introduced into stress-sensitive formation, and considering the wellbore storage effect and effective wellbore radiusrwe=rwe−s, we have(18)Kτr,p=Kweτeγτp−p0rrwedfτ−θτ−d,ϕτr,p=ϕweτeγτp−p0rrwedfτ−d,where γ is permeability modulus, df is fractal dimension, d is Euclidean dimension, θ is fractal index, and the subscript τ=1,2 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 p0; after well opening, it is produced with variable production rate qt. (5) The fluids and reservoir compressibility cL,cf≪γτ.To simplify the model, we need to introduce the following dimensionless variables, as shown in Table1.Table 1
Dimensionless variables (τ=1,2).
PτD=Kwe1h/1.842×10−3Bqeμ1p0−pτPressure in two zonesrDe=r/rweEffective wellbore radiusRDe=R/rweEffective radius of outer zonetD=3.6Kwe1/ϕwe1μ1Ct1rwe2tDimensionless timeTD=tD/CDDimensionless time of effective diameterCD=C/2πϕwe1Ct1hrwe2Wellbore storage coefficientqDtD=qt/qeProduction rateαD=1.842×10−3Bqeμ1rwe2/Kwe1hγτStress-sensitive coefficientω12=ϕwe1Ct1/ϕwe2Ct2Storativity ratio of the inner zone to the outer zoneM12=kwe1μ2/kwe2μ1Mobility ratio of the inner zone to the outer zoneσ=ϕwe2μ2Ct2/ϕwe1μ1Ct1kwe1/kwe2=M12/ω12Ratio ofM12 to ω12Then, 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(19)∂2p1D∂rDe2+df1−d−θ1+1rDe∂p1D∂rDe−αD∂p1D∂rDe2=rDeθ1eαDp1DCDeθ1+2S∂p1D∂TD,1≤rDe≤b,∂2p2D∂rDe2+df2−d−θ2+1rDe∂p2D∂rDe−αD∂p2D∂rDe2=rDeθ2σeαDp2DCDeθ2+2S∂p2D∂TD,b≤rDe≤RDe.The initial conditions of two zones are(20)p1DrDe,0=p2DrDe,0=0.The inner boundary conditions are(21)pwDTD=p1D1,TD,∂pwD∂TD−eαDp1D∂p1D∂rDe|rDe=1=qDTD.The interface conditions are(22)p1Db,TD=p2Db,TD∂p1D∂rDe|rDe=b=1M12∂p2D∂rDe|rDe=b.The three typical outer boundary conditions areInfinite boundary:(23)p2D∞,TD=0.Constant pressure boundary:(24)∂p2D∂rDe|rDe=RDe=0.Closed boundary:(25)p2DRDe,TD=0.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]:(26)PDrDe,tD=−1αDln1−αDηrDe,tD,where ηrDe,tD is the dimensionless dependent variable, and the regular perturbation series with respect to stress-sensitive parameter αD is defined as the series(27)η=η0+αDη1+αD2η2+⋯.If we stipulateαDη<1 and by using the regular perturbation theory, the following series is expanded as(28)11−αDη=1+αDη+αD2η2+⋯.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(29)∂2η01∂rDe2+df1−d−θ1+1rDe∂η01∂rDe=rDeθ11CDeθ2+2S∂η01∂TD,1≤rDe≤b,∂2η02∂rDe2+df2−d−θ2+1rDe∂η02∂rDe=rDeθ2σCDeθ2+2S∂η02∂TD,b≤rDe≤RDe.The initial conditions in equation (20) become(30)η01rDe,0=η02rDe,0=0.The inner boundary conditions in equation (21) become(31)η0wTD=η011,TD,∂η0w∂TD−rDedf1−d−θ1+1∂η01∂rDe|rDe=1=qDTD.The interface conditions in equation (22) become(32)η01b,TD=η02b,TD,∂η01∂rDe|rDe=b=1M12∂η02∂rDe|rDe=b.The three outer-boundary conditions in equation (23) and (24) are changed as(33)η02∞,TD=0,∂η02∂rDe|rDe=RDe=0,η02RDe,TD=0.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 TD, the following dimensionless BVPs of the modified Bessel equations are obtained:(34)d2η¯01drDe2+df1−d−θ1+1rDedη¯01drDe=rDeθ1zCDeθ2+2Sη¯01,d2η¯02drDe2+df2−d−θ2+1rDedη¯02drDe=rDeθ2zσCDeθ2+2Sη¯02,(35)η¯0wz=η¯011,z,zη¯0w−rpDedf1−d−θ1+1dη¯01drDe|rDe=1=q¯Dz,(36)η¯01b,z=η¯02b,z,dη01∂rDe|rDe=b=1M12dη02∂rDe|rDe=b,(37)η¯02∞,z=0,dη¯02∂rDe|rDe=RDe=0,η¯02RDe,z=0.


### 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 (τ=1,2),(38)y=η¯0τ,x=rDe,Aτ=βfτ=dfτ−d−θτ+1,C1=zCDeθ1+2S,C2=σzCDeθ2+2S,D=−q¯Dz,E=−z,F=0,δ=1M12,qτ=θτ+2,a=1,b=b,c=RDe.By using the SCM given in Section2.2, the solution of the seepage model can be easily obtained. According to equation (7), the Laplace space solution of η¯01, which is corresponding to the pressure in the inner zone, can be obtained(39)η¯01rDe,z=q¯Dz⋅1z−1/Φ11⋅Φ1rDeΦ11,1≤rDe≤b.IfrDe=1, the Laplace space solution of η¯0w, which is corresponding to bottom-hole pressure, can be obtained:(40)η¯0wz=q¯Dz⋅1z−1/Φ11.By equation (8), the Laplace space solution of η¯02, which is corresponding to the pressure in the outer zone, can be obtained:(41)η¯02rDe,z=q¯Dz⋅1z−1/Φ11⋅φ0,11b,bΦ2bφ1,111,b−1/M12φ1,011,b⋅Φ2rDeΦ11,b≤rDe≤RDe.According to equations (9) and (10), the SKF of η¯02rDe,z in the outer zone becomes(42)Φ2rDe=Gφ0,02rDe,RDe+Hφ0,12rDe,RDeGφ1,02b,RDe+Hφ1,12b,RDe,b≤rDe≤RDe,and the SKF of η¯01rDe,z in the inner zone becomes(43)Φ1rDe=Φ2bφ0,11rDe,b−1/M12φ0,01rDe,bΦ2bφ1,111,b−1/M12φ1,011,b,1≤rDe≤b.In real applications, the following special cases are more common. IfG=0, in equation (42), the outer condition is impermeable boundary; then, it is simplified as(44)Φ2rDe=φ0,12rDe,RDeφ1,12b,RDe,b≤rDe≤RDe.IfH=0, in equation (42), the outer condition is constant pressure boundary; then, it is simplified as(45)Φ2rDe=φ0,02rDe,RDeφ1,02b,RDe,b≤rDe≤RDe.If the outer boundary is infinite, the SKF is presented as follows:(46)Φ2rDe=rDe1−βfτ/2Kvτ2Cτz/θτ+2rDeθτ+2/21−βfτKvτ2Cτz/θτ+2−CτzKvτ+12Cτz/θτ+2,rDe⟶∞.Therefore, the key to solve the problem is to find the GFs,φi,jτrDe,ξi,j=0,1;τ=1,2, and construct the SKFs Φ1rDe and Φ2rDe; 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):(47)rDe1−βfτ/2Ivτ2Cτzθτ+2rDeθτ+2/2,rDe1−βfτ/2Kvτ2Cτzθτ+2rDeθτ+2/2.Compared with equation (16), we have(48)ατ=1−Aτ2,βτ=θτ+22,γτ=2Cτθτ+2,υτ=1−Aτθτ+2,τ=1,2.Step 2: secondly, assemble the GFs, φi,jτrDe,ξi,j=0,1;τ=1,2, by equation (11)–(14), as follows:(49)φ0,0τrDe,ξ=rDe1−βfτ/2Ivτ2Cτzθfτ+2rDeθfτ+2/2rDe1−βfτ/2Kvτ2Cτzθfτ+2rDeθfτ+2/2ξ1−βfτ/2Ivτ2Cτzθfτ+2ξθfτ+2/2ξ1−βfτ/2Kvτ2Cτzθfτ+2ξθfτ+2/2=−rDeξ1−βfτ/2ψυτ,υττrDeθfτ+2/2,ξθfτ+2/2,2Cτzθfτ+2.Similarly, taking the partial derivatives ofφi,jτrDe,ξi,j=0,1;τ=1,2, we have(50)φ1,0τrDe,ξ=∂∂rDeφ0,0τrDe,ξ=−rDe1−βfτ/2−1ξ1−βfτ/21−βfτψυτ,υττrDeθfτ+2/2,ξθfτ+2/2,2Cτzθfτ+2−CτzrDeβτψυτ+1,υττrDeθfτ+2/2,ξθfτ+2/2,2Cτzθfτ+2,φ0,1τrDe,ξ=∂∂ξφ0,0τrDe,ξ=−rDe1−βfτ/2ξ1−βfτ/2−11−βfτψυτ,υττrDeθfτ+2/2,ξθfτ+2/2,2Cτzθfτ+2−Cτzξθfτ+2/2ψυτ,υτ+1τrDeθfτ+2/2,ξθfτ+2/2,2Cτzθfτ+2,φ1,1τrDe,ξ=∂2∂rDe∂ξφ0,0τrDe,ξ=−rDeξ1−βfτ/2−11−βfτ2ψυτ,υττrDeθfτ+2/2,ξθfτ+2/2,2Cτzθfτ+2+Cτzβτ1−βfτξθfτ+2/2ψυτ,υτ+1τrDeθfτ+2/2,ξθfτ+2/2,2Cτzθfτ+2−Cτzβτ1−βfτrDeθfτ+2/2ψυτ+1,υττrDeθfτ+2/2,ξθfτ+2/2,2Cτzθfτ+2−fτzrDeξθfτ+2/2ψυτ+1,υτ+1τrDeθfτ+2/2,ξθfτ+2/2,2Cτzθfτ+2.The notation,(51)ψm,nη,ξ,t=KmηtInξt+−1m−n+1ImηtKnξt,is linear combination of modified Bessel functionsKnt,Int of the order n.Step 3: substitute φi,jτrDe,ξi,j=0,1;τ=1,2 into equations (42) and (43) and combine the boundary conditions in equations (36) and (37); we can generate the KFs, Φ1rDe and Φ2rDe.Step 4: finally, substitute the KFs, Φ1rDe and Φ2rDe, into equations (39) and (41), then we have the solutions of the BVPs in equations (34)–(37).


## 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 radiusrw≤r≤brw, and the area away from the wellbore is called outer zone if brw≤r≤R.The permeability and porosity with fractal reservoir behavior are introduced into stress-sensitive formation, and considering the wellbore storage effect and effective wellbore radiusrwe=rwe−s, we have(18)Kτr,p=Kweτeγτp−p0rrwedfτ−θτ−d,ϕτr,p=ϕweτeγτp−p0rrwedfτ−d,where γ is permeability modulus, df is fractal dimension, d is Euclidean dimension, θ is fractal index, and the subscript τ=1,2 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 p0; after well opening, it is produced with variable production rate qt. (5) The fluids and reservoir compressibility cL,cf≪γτ.To simplify the model, we need to introduce the following dimensionless variables, as shown in Table1.Table 1
Dimensionless variables (τ=1,2).
PτD=Kwe1h/1.842×10−3Bqeμ1p0−pτPressure in two zonesrDe=r/rweEffective wellbore radiusRDe=R/rweEffective radius of outer zonetD=3.6Kwe1/ϕwe1μ1Ct1rwe2tDimensionless timeTD=tD/CDDimensionless time of effective diameterCD=C/2πϕwe1Ct1hrwe2Wellbore storage coefficientqDtD=qt/qeProduction rateαD=1.842×10−3Bqeμ1rwe2/Kwe1hγτStress-sensitive coefficientω12=ϕwe1Ct1/ϕwe2Ct2Storativity ratio of the inner zone to the outer zoneM12=kwe1μ2/kwe2μ1Mobility ratio of the inner zone to the outer zoneσ=ϕwe2μ2Ct2/ϕwe1μ1Ct1kwe1/kwe2=M12/ω12Ratio ofM12 to ω12Then, 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(19)∂2p1D∂rDe2+df1−d−θ1+1rDe∂p1D∂rDe−αD∂p1D∂rDe2=rDeθ1eαDp1DCDeθ1+2S∂p1D∂TD,1≤rDe≤b,∂2p2D∂rDe2+df2−d−θ2+1rDe∂p2D∂rDe−αD∂p2D∂rDe2=rDeθ2σeαDp2DCDeθ2+2S∂p2D∂TD,b≤rDe≤RDe.The initial conditions of two zones are(20)p1DrDe,0=p2DrDe,0=0.The inner boundary conditions are(21)pwDTD=p1D1,TD,∂pwD∂TD−eαDp1D∂p1D∂rDe|rDe=1=qDTD.The interface conditions are(22)p1Db,TD=p2Db,TD∂p1D∂rDe|rDe=b=1M12∂p2D∂rDe|rDe=b.The three typical outer boundary conditions areInfinite boundary:(23)p2D∞,TD=0.Constant pressure boundary:(24)∂p2D∂rDe|rDe=RDe=0.Closed boundary:(25)p2DRDe,TD=0.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]:(26)PDrDe,tD=−1αDln1−αDηrDe,tD,where ηrDe,tD is the dimensionless dependent variable, and the regular perturbation series with respect to stress-sensitive parameter αD is defined as the series(27)η=η0+αDη1+αD2η2+⋯.If we stipulateαDη<1 and by using the regular perturbation theory, the following series is expanded as(28)11−αDη=1+αDη+αD2η2+⋯.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(29)∂2η01∂rDe2+df1−d−θ1+1rDe∂η01∂rDe=rDeθ11CDeθ2+2S∂η01∂TD,1≤rDe≤b,∂2η02∂rDe2+df2−d−θ2+1rDe∂η02∂rDe=rDeθ2σCDeθ2+2S∂η02∂TD,b≤rDe≤RDe.The initial conditions in equation (20) become(30)η01rDe,0=η02rDe,0=0.The inner boundary conditions in equation (21) become(31)η0wTD=η011,TD,∂η0w∂TD−rDedf1−d−θ1+1∂η01∂rDe|rDe=1=qDTD.The interface conditions in equation (22) become(32)η01b,TD=η02b,TD,∂η01∂rDe|rDe=b=1M12∂η02∂rDe|rDe=b.The three outer-boundary conditions in equation (23) and (24) are changed as(33)η02∞,TD=0,∂η02∂rDe|rDe=RDe=0,η02RDe,TD=0.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 TD, the following dimensionless BVPs of the modified Bessel equations are obtained:(34)d2η¯01drDe2+df1−d−θ1+1rDedη¯01drDe=rDeθ1zCDeθ2+2Sη¯01,d2η¯02drDe2+df2−d−θ2+1rDedη¯02drDe=rDeθ2zσCDeθ2+2Sη¯02,(35)η¯0wz=η¯011,z,zη¯0w−rpDedf1−d−θ1+1dη¯01drDe|rDe=1=q¯Dz,(36)η¯01b,z=η¯02b,z,dη01∂rDe|rDe=b=1M12dη02∂rDe|rDe=b,(37)η¯02∞,z=0,dη¯02∂rDe|rDe=RDe=0,η¯02RDe,z=0.


## 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 (τ=1,2),(38)y=η¯0τ,x=rDe,Aτ=βfτ=dfτ−d−θτ+1,C1=zCDeθ1+2S,C2=σzCDeθ2+2S,D=−q¯Dz,E=−z,F=0,δ=1M12,qτ=θτ+2,a=1,b=b,c=RDe.By using the SCM given in Section2.2, the solution of the seepage model can be easily obtained. According to equation (7), the Laplace space solution of η¯01, which is corresponding to the pressure in the inner zone, can be obtained(39)η¯01rDe,z=q¯Dz⋅1z−1/Φ11⋅Φ1rDeΦ11,1≤rDe≤b.IfrDe=1, the Laplace space solution of η¯0w, which is corresponding to bottom-hole pressure, can be obtained:(40)η¯0wz=q¯Dz⋅1z−1/Φ11.By equation (8), the Laplace space solution of η¯02, which is corresponding to the pressure in the outer zone, can be obtained:(41)η¯02rDe,z=q¯Dz⋅1z−1/Φ11⋅φ0,11b,bΦ2bφ1,111,b−1/M12φ1,011,b⋅Φ2rDeΦ11,b≤rDe≤RDe.According to equations (9) and (10), the SKF of η¯02rDe,z in the outer zone becomes(42)Φ2rDe=Gφ0,02rDe,RDe+Hφ0,12rDe,RDeGφ1,02b,RDe+Hφ1,12b,RDe,b≤rDe≤RDe,and the SKF of η¯01rDe,z in the inner zone becomes(43)Φ1rDe=Φ2bφ0,11rDe,b−1/M12φ0,01rDe,bΦ2bφ1,111,b−1/M12φ1,011,b,1≤rDe≤b.In real applications, the following special cases are more common. IfG=0, in equation (42), the outer condition is impermeable boundary; then, it is simplified as(44)Φ2rDe=φ0,12rDe,RDeφ1,12b,RDe,b≤rDe≤RDe.IfH=0, in equation (42), the outer condition is constant pressure boundary; then, it is simplified as(45)Φ2rDe=φ0,02rDe,RDeφ1,02b,RDe,b≤rDe≤RDe.If the outer boundary is infinite, the SKF is presented as follows:(46)Φ2rDe=rDe1−βfτ/2Kvτ2Cτz/θτ+2rDeθτ+2/21−βfτKvτ2Cτz/θτ+2−CτzKvτ+12Cτz/θτ+2,rDe⟶∞.Therefore, the key to solve the problem is to find the GFs,φi,jτrDe,ξi,j=0,1;τ=1,2, and construct the SKFs Φ1rDe and Φ2rDe; 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):(47)rDe1−βfτ/2Ivτ2Cτzθτ+2rDeθτ+2/2,rDe1−βfτ/2Kvτ2Cτzθτ+2rDeθτ+2/2.Compared with equation (16), we have(48)ατ=1−Aτ2,βτ=θτ+22,γτ=2Cτθτ+2,υτ=1−Aτθτ+2,τ=1,2.Step 2: secondly, assemble the GFs, φi,jτrDe,ξi,j=0,1;τ=1,2, by equation (11)–(14), as follows:(49)φ0,0τrDe,ξ=rDe1−βfτ/2Ivτ2Cτzθfτ+2rDeθfτ+2/2rDe1−βfτ/2Kvτ2Cτzθfτ+2rDeθfτ+2/2ξ1−βfτ/2Ivτ2Cτzθfτ+2ξθfτ+2/2ξ1−βfτ/2Kvτ2Cτzθfτ+2ξθfτ+2/2=−rDeξ1−βfτ/2ψυτ,υττrDeθfτ+2/2,ξθfτ+2/2,2Cτzθfτ+2.Similarly, taking the partial derivatives ofφi,jτrDe,ξi,j=0,1;τ=1,2, we have(50)φ1,0τrDe,ξ=∂∂rDeφ0,0τrDe,ξ=−rDe1−βfτ/2−1ξ1−βfτ/21−βfτψυτ,υττrDeθfτ+2/2,ξθfτ+2/2,2Cτzθfτ+2−CτzrDeβτψυτ+1,υττrDeθfτ+2/2,ξθfτ+2/2,2Cτzθfτ+2,φ0,1τrDe,ξ=∂∂ξφ0,0τrDe,ξ=−rDe1−βfτ/2ξ1−βfτ/2−11−βfτψυτ,υττrDeθfτ+2/2,ξθfτ+2/2,2Cτzθfτ+2−Cτzξθfτ+2/2ψυτ,υτ+1τrDeθfτ+2/2,ξθfτ+2/2,2Cτzθfτ+2,φ1,1τrDe,ξ=∂2∂rDe∂ξφ0,0τrDe,ξ=−rDeξ1−βfτ/2−11−βfτ2ψυτ,υττrDeθfτ+2/2,ξθfτ+2/2,2Cτzθfτ+2+Cτzβτ1−βfτξθfτ+2/2ψυτ,υτ+1τrDeθfτ+2/2,ξθfτ+2/2,2Cτzθfτ+2−Cτzβτ1−βfτrDeθfτ+2/2ψυτ+1,υττrDeθfτ+2/2,ξθfτ+2/2,2Cτzθfτ+2−fτzrDeξθfτ+2/2ψυτ+1,υτ+1τrDeθfτ+2/2,ξθfτ+2/2,2Cτzθfτ+2.The notation,(51)ψm,nη,ξ,t=KmηtInξt+−1m−n+1ImηtKnξt,is linear combination of modified Bessel functionsKnt,Int of the order n.Step 3: substitute φi,jτrDe,ξi,j=0,1;τ=1,2 into equations (42) and (43) and combine the boundary conditions in equations (36) and (37); we can generate the KFs, Φ1rDe and Φ2rDe.Step 4: finally, substitute the KFs, Φ1rDe and Φ2rDe, 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:(52)PDrD,tD=−1αDln1−αDL−1η¯0rD,tD,PwDtD=−1αDln1−αDL−1η¯wz.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.Figure2 shows the curves of bottom-hole pressure and pressure derivative of the stress-sensitive fractal composite reservoir with infinite outer boundary affected by parameter CDe2s. It can be seen from the figure that when the value of the parameter CDe2s 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 df<2.0, 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 CDe2s.Figure 2
Type curves of pressure dynamic in the fractal composite reservoir with stress-sensitivity formation affected by parameterCDe2S.Figure3 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 αD. As can be seen from the figure, the pressure does not depend on the parameter αD during the pure wellbore flow period. With the increase of time, the larger the value of αD becomes, the more obvious the curves diverges and the more significant the rise of the curves of pressure derivative is. Especially when αD approaches 0.1, the growth rate of the curves of both the pressure and pressure derivative increases significantly and shows an upward trend.Figure 3
Type curves of pressure dynamic with the infinite outer boundary in the fractal composite reservoir affected by stress-sensitivity coefficientαD.Figure4 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 αD. 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 αD. With the increase of time, the pressure derivative curves begin to diverge from the 4th flow period; the larger αD 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 4
Type curves of pressure dynamic with the infinite outer boundary of the two-zone fractal composite reservoir affected by stress-sensitivity coefficientαD.Figure5 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 5
Type curves of pressure dynamic with the infinite outer boundary of the two-zone fractal composite reservoir affected by parameterσ.Figure6 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 αD. As can be seen that δ and αD 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 αD are taken as 0.013, 0.01, and 0.005, respectively. The results show that the sensitive coefficient αD 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.Figure 6
Type curves of pressure dynamic in the two-zone fractal composite reservoir with stress-sensitivity formation affected by parametersδ=1/M12 and αD.


## 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 Section3 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 Section2, 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.



---

*Source: 4566248-2021-09-10.xml*