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

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.


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 [3][4][5].Al-Ahmadi et al. introduced a threeporosity 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 [11][12][13][14].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 [15][16][17][18].
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.

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.Cartesian coordinate system with normalized pressure as the variable, the MFHW seepage control differential equation is

Mathematical Model
The initial conditions are

External Boundary Conditions
direction: Constant flow conditions: where   ,   indicate horizontal permeability (  =   =  ℎ ), (m 2 );   is vertical permeability (m 2 );  is fluid viscosity (pa⋅s) and  ini 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 (m 3 /s);  is the normalized pressure of the tight oil and gas; the formula is as follows [19]: (the normalized pressure of the tight oil) (pa); (the normalized pressure of the tight gas) (pa); . ( 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: Equation ( 8) is the mathematical model in this study.In the next section, (8) will be solved analytically.

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.
(1) 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)   .
(2) 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 (3) 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  ×   + 1 unknowns in the above equation, respectively:  *     and   .To determine the  ×   + 1 unknowns, we must add another equation.Because the horizontal wells are produced at constant yield, there is The coefficient matrix is expressed as where (4) 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   =   min (  max /  min ) /(  −1) , where   is the dimensionless time of step i,   min is the dimensionless time start,   max 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.

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.

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 (5) to the host-side memory and release the GPU storage space.
The details of the parallel strategy implementation are discussed in the next section.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.

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  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.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 9-11, 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 7-12.The results of the pressure distribution are in accordance with expectations, and the correctness of the algorithm and software was also verified.

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

Figure 1 :
Figure 1: Schematic diagram of the relationship between horizontal wells, fractures, and formation spaces.

Figure 6 :
Figure 6: History match figure of bottom-hole pressure.

Table 2 :
show that the pressure around the crack is low and the descent speed is fast, indicating that the effect of fracturing is good.In Figures7 and 8, the flow is concentrated in the vicinity of the crack, which is the linear flow stage at this time.From Figures 9-11, the flow gradually expands outward, at which point the flow behaves as a linear flow to the radial flow transition phase.In Figure12, 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 Figures7-12.The results of the pressure distribution are in accordance with expectations, and the correctness of the algorithm and software was also verified.
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.