Abstract

A mathematical model for multistage hydraulically fractured horizontal wells (MFHWs) in tight oil and gas reservoirs was derived by considering the variations in the permeability and porosity of tight oil and gas reservoirs that depend on formation pressure and mixed fluid properties and introducing the pseudo-pressure; analytical solutions were presented using the Newman superposition principle. The CPU-GPU asynchronous computing model was designed based on the CUDA platform, and the analytic solution was decomposed into infinite summation and integral forms for parallel computation. Implementation of this algorithm on an Intel i5 4590 CPU and NVIDIA GT 730 GPU demonstrates that computation speed increased by almost 80 times, which meets the requirement for real-time calculation of the formation pressure of MFHWs.

1. Introduction

Horizontal well drilling and fracturing are key processes developed for unconventional oil and gas (such as shale gas and tight oil and gas).

The longer the horizontal wells, the more fractured the segments and the higher the corresponding monthly production and costs. Therefore, horizontal well fracturing optimization based on single well productivity prediction technology is of great significance.

The prediction of single well productivity is realized by solving the equation of porous flow and obtaining the bottom-hole pressure or output. The commonly used solution method is an analytic method based on the superposition principle. According to different permeabilities of the formation, single well production capacity can be divided into two types: capacity prediction based on the steady-state seepage equation and capacity prediction based on the transient seepage equation.

Because of the simplicity of the calculation of partial differential equations for steady-state capacity, it is widely used in the field.

Using the homogeneous horizontal well model, Babu et al. obtained the quasi-steady-state production formula for horizontal wells, which is applicable to horizontal wells in the presence of bottom water or gas-filled rectangular reservoirs [1, 2]. Yildiz et al. studied the characteristics of perforated wellhead wells in heterogeneous reservoirs and analyzed the variations in law of yield under the quasi-steady state of horizontal wells [35]. Al-Ahmadi et al. introduced a three-porosity model and established an analytical model based on the diffusion equation of the instantaneous source solution. The new technology took into account the effects of reservoir geometry, reservoir performance, fracture size, number of fractures, and spacing on production capacity [6, 7]. Xie and Li divided the horizontal well seepage area into four parts: the first and the second stage of the cylinder, the hemispherical sphere facing the heart flow, and the wellbore area flow. Based on the principle of equivalent seepage resistance, the expression of steady-state seepage capacity of fractured horizontal wells was proposed according to the formula of shale gas desorption diffusion rate [8].

The expression of steady-state capacity is simple and suitable for medium-high permeability formation prediction; however, for low-permeability oil and gas reservoirs, especially tight oil and gas reservoirs and shale gas reservoirs, the computation error is large.

With the development of computational techniques, partial differential equations, including the time factor, can be directly solved. Penmatcha and Aziz established a dynamic forecasting model of horizontal well production based on the Babu–Odeh method, which considered the influence of friction pressure drop, acceleration pressure drop, and radial flow in the horizontal well [9]. Penmatcha et al. established a model considering the basic physical properties of oil and gas reservoirs as well as some special characteristics of the shale reservoir matrix exchanging fluids with different reservoir medias and proposed a pressure transient solution for shale gas production [9, 10]. Guo et al. analyzed the effects of desorption, diffusion, viscous flow, pressure sensitivity, fracturing number of horizontal wells, and fracture spacing. Through the source function and Laplace transform combined with the numerical discretization method, the solution was obtained in the Laplace space. Furthermore, the plate curve was obtained through the Stehfest algorithm, which was applied to identify a variety of flow patterns [1114]. Through the shale gas fracturing model, Sang et al. obtained shale gas single well production and bottom flow pressure calculation formulas, by which shale gas well production was calculated [14]. Zhao et al. studied the seepage flow behaviors of multistage fractured horizontal wells in arbitrary shaped shale gas reservoirs and got the pressure response and production performance for multifractured horizontal wells [1518].

In the above study, both steady-state production and transient productivity were calculated for bottom pressure or flow, with no computation of formation pressure distribution. Since the pressure at any point during oil and gas formation can change over time, it is necessary to calculate the pressure distribution at each point of formation at different times. If the computation area is decomposed into 100,000 meshes at each time step, the computation is performed hundreds of thousands of times with existing method, and serial computation cannot meet the requirements of the oil and gas field. For this reason, this paper presents GPU-based model for computing formation pressure distribution.

2. Model

2.1. Physical Model

Suppose that there exist multistage hydraulically fractured horizontal wells (MFHWs) in a fully enclosed rectangular reservoir. Horizontal well positions, crack locations at all levels, crack shapes, and reservoir boundaries are shown in Figure 1. Figures 2 and 3 show the - plan and - plan, respectively. The basic assumptions are as follows:(1)There exists stratigraphic homogeneity, but the horizontal and vertical directions of permeability are different.(2)Fluid flow meets Darcy’s law and ignores the impact of gravity.(3)The effective thickness of the formation is the same, and the horizontal wells are parallel to the top and bottom. Assuming that the upper roof is completely closed (no gas roof and bottom water, as shown in Figure 3), the height of each crack in the horizontal well is different.(4)The stratum is a rectangular closed reservoir and the horizontal wells are parallel to both sides of the rectangle (Figure 2). Horizontal wells have cracks, their half-lengths are , and all crack spacings are changeable. The plane of the fracture is at an angle of to the horizontal wells of the wellbore. Cracks and horizontal wells are not symmetrical, and the length of crack is .(5)The fluid in the reservoir is tight oil and is mined at a constant flow rate.(6)Fluid and formation are microcompressible.

2.2. Mathematical Model
2.2.1. Governing Equations

According to the description and assumptions of the physical model, in the three-dimensional Cartesian coordinate system with normalized pressure as the variable, the MFHW seepage control differential equation isThe initial conditions are

2.2.2. External Boundary Conditions

Constant flow conditions:

where indicate horizontal permeability , (m2); is vertical permeability (m2); is fluid viscosity (pa·s) and is original formation pressure (pa); is stratigraphic and fluid comprehensive compression coefficient (1/pa); is porosity, , are the th crack production and total output, respectively (m3/s); is the normalized pressure of the tight oil and gas; the formula is as follows [19]: ini is subscript, indicating the initial state.

In order to make the equation more concise and facilitate computation, the dimensionless form of the equation was used, which is defined as follows:where (, , ) is the position; is the th crack half-length; are the area length, width, and height, respectively; is the location of the crack; and is the dimensionless crack location.

The dimensionless equation and definite solution conditions are

Equation (8) is the mathematical model in this study. In the next section, (8) will be solved analytically.

3. Solution

For the abovementioned dimensionless equations, researchers have only given the expression for the pressure of the bottom of the well over time.

In this paper, the pressure distribution in the formation needs to be solved, and the existing method is not suitable for the above model. Thus, the principle of superposition combined with point source integral method was used to solve the problem.

There are cracks in the formation; the th crack is divided into microcells, with length , flow , and coordinate position . These microcells are called the source (or sink) .

For a given source , according to the principle of the line source solution, the pressure distribution on each source unit is [20]where

According to the assumptions of the previous physical model, the horizontal wells have cracks, each crack is divided into sections in addition to its own, and any one of the cracks has pressure on the other target points (Figure 4). According to the superposition principle, the pressure distribution at any point in the fracture is

Assuming that the wellbore and crack are infinite diversions, the frictional resistance is not considered between the sections within the fracture so that the pressure in the wellbore and cracks is equal everywhere. Hence, the bottom flow pressure is

There are unknowns in the above equation, respectively: and . To determine the unknowns, we must add another equation. Because the horizontal wells are produced at constant yield, there is . The coefficient matrix is expressed aswhere

According to the characteristics of the point source function solution in seepage mechanics, when the time is long, the bottom-hole pressure is proportional to the logarithmic time [21]. Therefore, the time step in this paper adopts the logarithmic time of equal step with the formula , where is the dimensionless time of step i, is the dimensionless time start, is the dimensionless end of time, and is the number of equal parts of time.

The calculation is related to the integral of the time term and three directions of the source function , and the function is characterized by polynomial summation and product. Thus, it is very suitable for parallel computing, especially for GPU-based computation.

4. Parallel Algorithm

For MFHWs, computing the pressure of all segments of each crack at each time step within an acceptable time is a major challenge. Using the parallel computing power of the GPU, a parallel strategy was designed on the CUDA platform to better solve this computationally intensive problem.

4.1. Parallel Algorithm Design

CUDA’s master-slave design pattern for carrying out parallel operation of the GPU was used. In the calculation process, in order to improve the calculation efficiency and reduce the waste of resources, the asynchronous call calculation mode was adopted. Figure 5 shows the CPU-GPU implementation flowchart.

The specific implementation steps are as follows:(1)On the host side, initialize the cracks, fractures, and other related data.(2)Copy the processed data into GPU memory.(3)Call kernel 1 to calculate of each segment of the crack. If the calculation ends, go to the next step, or continue to call kernel 1 to calculate.(4)Call kernel 2 to carry out the reduction calculation on the 2 section.(5)Call kernel 3 to carry out the reduction calculation on the 3 crack.(6)Copy the calculated result in step to the host-side memory and release the GPU storage space.

The details of the parallel strategy implementation are discussed in the next section.

4.2. Kernel 1 Algorithm Implementation

Kernel 1 mainly calculates of each segment of the crack, which consists of three parts, namely, , , and . The pseudocode under the CPU computation condition is as follows:For = 0 to  For = 1 to    Calculate     At the integration point i, sum, calculate  For to    Calculate     At the integration point i, sum, calculate  For to    Calculate     At the integration point i, sum, calculate    Calculate     Calculate the of a section of the crack

Through the GPU reduction algorithm, the parallel algorithm of kernel 1 is obtained. The core pseudocode is as follows:__global__ void cal_oneG(double g_idata, double g_odata, unsigned int n)  extern __shared__ double sdata   // Allocate shared memory   unsigned int tid = blockIdx.x blockDim.x + threadIdx.x;   unsigned int i = blockIdx.x(blockDim.x2) + threadIdx.x;   double mySum = (i < n) ? g_idata[i] : 0;   if ( + blockSize < )     mySum += g_idata[i+blockSize];   sdata[tid] = mySum;/ Copy the data from the global storage area to thelocal storage area, and conduct a reduction /   __syncthreads(); // Thread synchronization   for (unsigned int s=blockDim.x/2; s>32; s=1)      if (tid < s)         sdata[tid] = mySum = mySum + sdata[tid + s];     __syncthreads();   / The sum calculation is continued within a single thread bundle (Warp) / if (tid == 0) godata[blockIdx.x] = mySum;/ Write the final summation result back to the global store / // and are obtained by an algorithm similar to  // Combine the values of , , and to get

In the case of each segment has been obtained, and the overall pressure of the underground horizontal wells of 2 segments and 3 cracks can be obtained by kernel 2 and kernel 3, respectively. Kernel 2 and kernel 3 algorithm design is similar to that of kernel 1.

4.3. Algorithm Analysis

By comparison, the time complexity on the GPU is , while that on the CPU is . GPU can improve the efficiency of the calculation exponentially. In addition, in order to improve the accuracy, double precision calculation was used to reduce calculation errors in the GPU kernel function.

5. Case Study

The bottom flow pressure was analyzed and calculated using the measured data of a MFHW in western China. When the parameters of Table 1 were selected, the bottom flow pressure history match chart shown in Figure 6 was obtained. The abovementioned data were also used to calculate the formation pressure distribution, and the pressure distribution is shown in Figures 712 under the time steps 90 h, 360 h, 980 h, 2580 h, 3650 h, and 4360 h, respectively.

The parameters shown in Table 1 (all unit conversions are shown in Table 2) were substituted into the equation, and the comparison of calculated results and measured data is shown in Figure 6. The figure shows that the measured and calculated pressures show good agreement with an average error of 2.5%, which is within an acceptable range.

Figures 712 show that the pressure around the crack is low and the descent speed is fast, indicating that the effect of fracturing is good. In Figures 7 and 8, the flow is concentrated in the vicinity of the crack, which is the linear flow stage at this time. From Figures 911, the flow gradually expands outward, at which point the flow behaves as a linear flow to the radial flow transition phase. In Figure 12, the flow is about to reach the boundary, where the flow is a combination of radial flow and boundary control flow. The flow patterns are in good agreement with seepage mechanics laws from Figures 712. The results of the pressure distribution are in accordance with expectations, and the correctness of the algorithm and software was also verified.

6. Conclusion

In this paper, the seepage equation of the MFHW of tight oil and gas was proposed by considering seepage and pressure-sensitive effects, and a GPU-based parallel computing method was designed and implemented. The test results show that computation speed and accuracy greatly improved. The computation speed for formation pressure distribution on the GPU is nearly 80 times that of the single-core CPU, and the calculation error is approximately 2.5%.

Conflicts of Interest

There are not any conflicts of interest related to this paper.

Acknowledgments

This work is supported by CAS Strategic Priority Research Program (Grant no. XDB10030402), CNPC-CAS Strategic Cooperation Research Program (Grant no. 2015A-4812), National Natural Science Foundation of China (Grant no. 41672114), National Science and Technology Major Project (Grant no. 2017ZX05009005-002), and Foundation of Anhui Educational Committee (Grant no. KJ2015B1105917).