Table of Contents Author Guidelines Submit a Manuscript
Mathematical Problems in Engineering
Volume 2018, Article ID 2582797, 10 pages
Research Article

GPU-Based Computation of Formation Pressure for Multistage Hydraulically Fractured Horizontal Wells in Tight Oil and Gas Reservoirs

1School of Engineering Science, University of Science and Technology of China, Hefei 230026, China
2Department of Basic Teaching and Experiment, Hefei University, Hefei 230601, China
3Engineering Technology Research Institute, Southwest Oil and Gas Field Branch, Deyang 610051, China

Correspondence should be addressed to Detang Lu; nc.ude.ctsu@ultd

Received 30 October 2017; Accepted 30 January 2018; Published 24 April 2018

Academic Editor: Suzanne M. Shontz

Copyright © 2018 Rongwang Yin et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


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.

Figure 1: Schematic diagram of the relationship between horizontal wells, fractures, and formation spaces.
Figure 2: Schematic diagram of - plane projection of MFHW.
Figure 3: Schematic diagram of - plane projection of MFHW.
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

Figure 4: Pressure distribution of the source unit.

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.

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

Table 1: Well parameters.
Figure 6: History match figure of bottom-hole pressure.
Figure 7: 90 h formation pressure distribution.
Figure 8: 360 h formation pressure distribution.
Figure 9: 980 h formation pressure distribution.
Figure 10: 2580 h formation pressure distribution.
Figure 11: 3650 h formation pressure distribution.
Figure 12: 4360 h formation pressure distribution.

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.

Table 2: SI metric conversion factors.

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.


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


  1. D. Babu and A. S. Odeh, “Productivity of a Horizontal Well (includes associated papers 20306, 20307, 20394, 20403, 20799, 21307, 21610, 21611, 21623, 21624, 25295, 25408, 26262, 26281, 31025, and 31035),” SPE Reservoir Engineering, vol. 4, no. 04, pp. 417–421, 2013. View at Publisher · View at Google Scholar
  2. P. A. Goode and F. J. Kuchuk, “Inflow performance of horizontal wells,” SPE Reservoir Engineering, vol. 6, no. 3, pp. 319–323, 1991. View at Google Scholar · View at Scopus
  3. T. Yildiz, “Inflow Performance Relationship for Perforated Horizontal Wells,” in Proceedings of the SPE Production and Operations Symposium, Oklahoma City, Oklahoma. View at Publisher · View at Google Scholar
  4. B. Guo, A. Ghalambor, and K. Ling, “A rigorous composite-inflow-performance relationship model for multilateral wells [J],” SPE Production Operations, vol. 23, no. 02, pp. 241–248, 2008. View at Google Scholar
  5. J. Liu, P. Li, Z. Sun et al., “A new method for analysis of dual pore size distributions in shale using nitrogen adsorption measurements,” Fuel, vol. 210, pp. 446–454, 2017. View at Publisher · View at Google Scholar · View at Scopus
  6. H. A. Al-Ahmadi and R. A. Wattenbarger, “Triple-porosity models: one further step towards capturing fractured reservoirs heterogeneity,” in Proceedings of the SPE/DGS Saudi Arabia Section Technical Symposium and Exhibition, SPE-149054-MS, Al-Khobar, Saudi Arabia, May 2011. View at Publisher · View at Google Scholar
  7. S. Al Rbeawi and D. Tiab, “Productivity index and inflow performance of hydraulically fractured formations,” in Proceedings of the SPE Annual Technical Conference and Exhibition 2012: Unconventional Wisdom, ATCE 2012, pp. 2980–2998, usa, October 2012. View at Scopus
  8. W. Xie and X. Li, “Research on Fractured Horizontal Wells Productivity and Productivity Influence in Shale Gas Reservoir,” in Proceedings of the SPE/EAGE European Unconventional Resources Conference and Exhibition, Vienna, Austria. View at Publisher · View at Google Scholar
  9. V. R. Penmatcha and K. Aziz, “Comprehensive reservoir/wellbore model for horizontal wells,” in Proceedings of the 1998 SPE India Oil and Gas Conference and Exhibition, pp. 191–204, New Delhi, India, April 1998. View at Scopus
  10. D. Li, L. Zhang, J. Y. Wang, and D. Lu, “Composition-transient analysis in shale-gas reservoirs with consideration of multicomponent adsorption,” SPE Journal, vol. 21, no. 2, pp. 648–664, 2016. View at Publisher · View at Google Scholar · View at Scopus
  11. J. Guo, L. Zhang, H. Wang, and G. Feng, “Pressure transient analysis for multi-stage fractured horizontal wells in shale gas reservoirs,” Transport in Porous Media, vol. 93, no. 3, pp. 635–653, 2012. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  12. H.-T. Wang, “Performance of multiple fractured horizontal wells in shale gas reservoirs with consideration of multiple mechanisms,” Journal of Hydrology, vol. 510, pp. 299–312, 2014. View at Publisher · View at Google Scholar · View at Scopus
  13. D. Li, C. Xu, J. Y. Wang, and D. Lu, “Effect of Knudsen diffusion and Langmuir adsorption on pressure transient response in tight- and shale-gas reservoirs,” Journal of Petroleum Science and Engineering, vol. 124, pp. 146–154, 2014. View at Publisher · View at Google Scholar · View at Scopus
  14. Y. Sang, H. Chen, S. Yang et al., “A new mathematical model considering adsorption and desorption process for productivity prediction of volume fractured horizontal wells in shale gas reservoirs,” Journal of Natural Gas Science and Engineering, vol. 19, pp. 228–236, 2014. View at Publisher · View at Google Scholar · View at Scopus
  15. Y.-L. Zhao, B.-C. Shan, L.-H. Zhang, and Q.-G. Liu, “Seepage flow behaviors of multi-stage fractured horizontal wells in arbitrary shaped shale gas reservoirs,” Journal of Geophysics and Engineering, vol. 13, no. 5, pp. 674–689, 2016. View at Publisher · View at Google Scholar · View at Scopus
  16. Y. Zhao, L. Zhang, Y. Xiong, Y. Zhou, Q. Liu, and D. Chen, “Pressure response and production performance for multi-fractured horizontal wells with complex seepage mechanism in box-shaped shale gas reservoir,” Journal of Natural Gas Science and Engineering, vol. 32, pp. 66–80, 2016. View at Publisher · View at Google Scholar · View at Scopus
  17. D. Li, W. Zha, S. Liu, L. Wang, and D. Lu, “Pressure transient analysis of low permeability reservoir with pseudo threshold pressure gradient,” Journal of Petroleum Science and Engineering, vol. 147, pp. 308–316, 2016. View at Publisher · View at Google Scholar · View at Scopus
  18. X. Li, Z. Guo, Y. Wan et al., “Geological characteristics and development strategies for Cambrian Longwangmiao Formation gas reservoir in Anyue gas field, Sichuan Basin, SW China,” Shiyou Kantan Yu Kaifa/Petroleum Exploration and Development, vol. 44, no. 3, pp. 398–406, 2017. View at Publisher · View at Google Scholar · View at Scopus
  19. R. Raghavan and L. Chin, “Productivity Changes in Reservoirs With Stress-Dependent Permeability,” in Proceedings of the SPE Annual Technical Conference and Exhibition, San Antonio, Tex, USA, 2002. View at Publisher · View at Google Scholar
  20. R. N. Horne and K. O. Temeng, “Relative productivities and pressure transient modeling of horizontal wells with multiple fractures,” in Proceedings of the 9th Middle East Oil Show & Conference. Part 2 (of 2), pp. 563–574, Bahrain, March 1995. View at Scopus
  21. X. Kong, “Advanced Mechanics of Fluids in Porous Media (Second Edition)[M],” in University of Science and Technology of, China Press, 2010. View at Google Scholar