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

Authors: Rongwang Yin; Qingyu Li; Peichao Li; Yang Guo; Yurong An; Detang Lu
Journal: Mathematical Problems in Engineering
 (2018)

---

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


---

## Body

## 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 [3–5]. 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 [11–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–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.


## 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 Figure1. Figures 2 and 3 show the x-y plan and x-z 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 Figure3), 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 (Figure2). Horizontal wells have N cracks, their half-lengths are xfi, and all crack spacings are changeable. The plane of the fracture is at an angle of αi to the horizontal wells of the wellbore. Cracks and horizontal wells are not symmetrical, and the length of crack is li.(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 ofx-y plane projection of MFHW.Figure 3
Schematic diagram ofx-z 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 is(1)kxiniμini∂2m∂x2+kyiniμini∂2m∂y2+kziniμini∂2m∂z2=ϕinictini∂m∂tx,y,z∈Ω,t>0.The initial conditions are(2)mx,y,z,tt=0=mini.


#### 2.2.2. External Boundary Conditions

(3) X direction:  ∂ m ∂ x x = 0 = 0 ∂ m ∂ x x = x e = 0 Y direction:  ∂ m ∂ y y = 0 = 0 ∂ m ∂ y y = y e = 0 Z direction:  ∂ m ∂ z z = 0 = 0 ∂ m ∂ z z = h = 0 .Constant flow conditions:(4)Q=∑i=1nQi,wherekx,ky indicate horizontal permeability (kx=ky=kh), (m2); kz is vertical permeability (m2); μ is fluid viscosity (pa·s) and mini is original formation pressure (pa); Ct is stratigraphic and fluid comprehensive compression coefficient (1/pa); ϕ is porosity, Qi, Q are the ith crack production and total output, respectively (m3/s); m is the normalized pressure of the tight oil and gas; the formula is as follows [19]: (5)mp=1-ϕinikini∫0pkp1-ϕpdp(the normalized pressure of the tight oil)(pa);mp=1-ϕinikini·μinizinipini∫0pkp·pμpzp1-ϕpdp(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:(6)mDxD,yD,zD,tD=2πKHhαBμmini-mx,y,z,t;tD=KHtϕμcixf2;xD=xxf,yD=yxf,xwiD=xwixf,xed=xexf,yeD=yexf,zD=zh,QiD=QiQ;xfiD=xfixf;xf=∑i=1nxfi,where (x, y, z) is the position; xfi is the ith crack half-length; (xe,ye,h) are the area length, width, and height, respectively; xwi is the location of the crack; and xwiD is the dimensionless crack location.The dimensionless equation and definite solution conditions are(7)∂2mD∂xD2+∂2mD∂yD2+KzKh∂2mD∂zD2=∂2mD∂tD;mDxD,yD,zD,tDtD=0=0(8)∂mD∂xDxD=0=∂mD∂xDxD=xeD=0∂mD∂yDyD=0=∂mD∂yDyD=yeD=0∂mD∂zDzD=0=∂mD∂xDxD=1=0∑i=1nQiD=1.Equation (8) is the mathematical model in this study. In the next section, (8) will be solved analytically.


## 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 Figure1. Figures 2 and 3 show the x-y plan and x-z 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 Figure3), 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 (Figure2). Horizontal wells have N cracks, their half-lengths are xfi, and all crack spacings are changeable. The plane of the fracture is at an angle of αi to the horizontal wells of the wellbore. Cracks and horizontal wells are not symmetrical, and the length of crack is li.(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 ofx-y plane projection of MFHW.Figure 3
Schematic diagram ofx-z 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 is(1)kxiniμini∂2m∂x2+kyiniμini∂2m∂y2+kziniμini∂2m∂z2=ϕinictini∂m∂tx,y,z∈Ω,t>0.The initial conditions are(2)mx,y,z,tt=0=mini.


### 2.2.2. External Boundary Conditions

(3) X direction:  ∂ m ∂ x x = 0 = 0 ∂ m ∂ x x = x e = 0 Y direction:  ∂ m ∂ y y = 0 = 0 ∂ m ∂ y y = y e = 0 Z direction:  ∂ m ∂ z z = 0 = 0 ∂ m ∂ z z = h = 0 .Constant flow conditions:(4)Q=∑i=1nQi,wherekx,ky indicate horizontal permeability (kx=ky=kh), (m2); kz is vertical permeability (m2); μ is fluid viscosity (pa·s) and mini is original formation pressure (pa); Ct is stratigraphic and fluid comprehensive compression coefficient (1/pa); ϕ is porosity, Qi, Q are the ith crack production and total output, respectively (m3/s); m is the normalized pressure of the tight oil and gas; the formula is as follows [19]: (5)mp=1-ϕinikini∫0pkp1-ϕpdp(the normalized pressure of the tight oil)(pa);mp=1-ϕinikini·μinizinipini∫0pkp·pμpzp1-ϕpdp(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:(6)mDxD,yD,zD,tD=2πKHhαBμmini-mx,y,z,t;tD=KHtϕμcixf2;xD=xxf,yD=yxf,xwiD=xwixf,xed=xexf,yeD=yexf,zD=zh,QiD=QiQ;xfiD=xfixf;xf=∑i=1nxfi,where (x, y, z) is the position; xfi is the ith crack half-length; (xe,ye,h) are the area length, width, and height, respectively; xwi is the location of the crack; and xwiD is the dimensionless crack location.The dimensionless equation and definite solution conditions are(7)∂2mD∂xD2+∂2mD∂yD2+KzKh∂2mD∂zD2=∂2mD∂tD;mDxD,yD,zD,tDtD=0=0(8)∂mD∂xDxD=0=∂mD∂xDxD=xeD=0∂mD∂yDyD=0=∂mD∂yDyD=yeD=0∂mD∂zDzD=0=∂mD∂xDxD=1=0∑i=1nQiD=1.Equation (8) is the mathematical model in this study. In the next section, (8) will be solved analytically.


## 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 is(1)kxiniμini∂2m∂x2+kyiniμini∂2m∂y2+kziniμini∂2m∂z2=ϕinictini∂m∂tx,y,z∈Ω,t>0.The initial conditions are(2)mx,y,z,tt=0=mini.


## 2.2.2. External Boundary Conditions

(3) X direction:  ∂ m ∂ x x = 0 = 0 ∂ m ∂ x x = x e = 0 Y direction:  ∂ m ∂ y y = 0 = 0 ∂ m ∂ y y = y e = 0 Z direction:  ∂ m ∂ z z = 0 = 0 ∂ m ∂ z z = h = 0 .Constant flow conditions:(4)Q=∑i=1nQi,wherekx,ky indicate horizontal permeability (kx=ky=kh), (m2); kz is vertical permeability (m2); μ is fluid viscosity (pa·s) and mini is original formation pressure (pa); Ct is stratigraphic and fluid comprehensive compression coefficient (1/pa); ϕ is porosity, Qi, Q are the ith crack production and total output, respectively (m3/s); m is the normalized pressure of the tight oil and gas; the formula is as follows [19]: (5)mp=1-ϕinikini∫0pkp1-ϕpdp(the normalized pressure of the tight oil)(pa);mp=1-ϕinikini·μinizinipini∫0pkp·pμpzp1-ϕpdp(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:(6)mDxD,yD,zD,tD=2πKHhαBμmini-mx,y,z,t;tD=KHtϕμcixf2;xD=xxf,yD=yxf,xwiD=xwixf,xed=xexf,yeD=yexf,zD=zh,QiD=QiQ;xfiD=xfixf;xf=∑i=1nxfi,where (x, y, z) is the position; xfi is the ith crack half-length; (xe,ye,h) are the area length, width, and height, respectively; xwi is the location of the crack; and xwiD is the dimensionless crack location.The dimensionless equation and definite solution conditions are(7)∂2mD∂xD2+∂2mD∂yD2+KzKh∂2mD∂zD2=∂2mD∂tD;mDxD,yD,zD,tDtD=0=0(8)∂mD∂xDxD=0=∂mD∂xDxD=xeD=0∂mD∂yDyD=0=∂mD∂yDyD=yeD=0∂mD∂zDzD=0=∂mD∂xDxD=1=0∑i=1nQiD=1.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.( 1 ) There are N cracks in the formation; the ith crack is divided into Mi microcells, with length dij, flow Qij, and coordinate position (xwij,ywij). These microcells are called the source (or sink) δij.( 2 ) For a given source δij, according to the principle of the line source solution, the pressure distribution on each source unit is [20](9)mDijxwDij,ywDij,zwDij,tD=∫0tDqDijτGxDxwDij,τDGyDywDij,τDGzDzwDij,τDdτ,where(10)GxDxwDij,τ=1xeD1+2∑n=1∞exp⁡-n2π2xeD2τcos2⁡nπxwDijxeDGyDywDij,τ=1yeD1+2∑n=1∞exp⁡-n2π2yeD2τcos2⁡nπywDijyeDGzDzwDij,τ=zbiD-zaiD1+4πzbiD-zaiD∑n=1∞exp⁡-n2π2τ·KzKhsin⁡nπzbiD-zaiD2cos2⁡nπzwdij.According to the assumptions of the previous physical model, the horizontal wells haveN cracks, each crack is divided into Mi 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(11)mDxDkl,yDkl,zDkl,tD=∑i=1N∑j=1MimDijxwDij,ywDij,zwDij,tD1≤k≤N,1≤l≤Mi,k≠i,l≠j.Figure 4
Pressure distribution of the source unit.( 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(12)mDwf=∑i=1N∑j=1MimDijxwDij,ywDij,zwDij,tD.There areN×Mi+1 unknowns in the above equation, respectively: N∗MiqDij and mDwf. To determine the N×Mi+1 unknowns, we must add another equation. Because the horizontal wells are produced at constant yield, there is ∑i=1N‍∑j=1MiqDij=1. The coefficient matrix is expressed as(13)ψ11ψ12⋯ψ1NMi-1ψ21ψ22⋯ψ2NMi-1⋮⋮⋮⋮⋮ψNMi1ψNMi2⋯ψNMiNMi⋮11⋯10·qD11qD12⋮qDNMimDwf=00⋮01,where(14)ψDtD=∫0tDGxDxwDij,τDGyDywDij,τDGzDzwDij,τDdτ.( 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 TDi=TDminTDmax/TDmini/Nt-1, where TDi is the dimensionless time of stepi, TDmin is the dimensionless time start, TDmax is the dimensionless end of time, and Nt is the number of equal parts of time.( 5 ) The calculation is related to the integral of the time term and three directions of the source function (x,y,z), 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. Figure5 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 calculateG 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 theN2 section.(5)
Call kernel 3 to carry out the reduction calculation on theN3 crack.(6)
Copy the calculated result in step5 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 calculatesG of each segment of the crack, which consists of three parts, namely, Gx, Gy, and Gz. The pseudocode under the CPU computation condition is as follows:Fork = 0 to N1-1Fornx = 1 to NxCalculateGx/∗ At the integration pointi, sum, calculate Gx ∗/Forny=1 to NyCalculateGy/∗ At the integration pointi, sum, calculate Gy ∗/Fornz=1 to NzCalculateGz/∗ At the integration pointi, sum, calculate Gz ∗/CalculateG/∗ Calculate the G 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 memoryunsigned int tid = blockIdx.x∗blockDim.x + threadIdx.x;unsigned int i = blockIdx.x∗(blockDim.x∗2) + threadIdx.x;double mySum = (i < n) ? g_idata[i] : 0;if (i + blockSize < n)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 synchronizationfor (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 singlethread bundle (Warp)∗/if (tid == 0) g_odata[blockIdx.x] = mySum;/∗ Write the final summation result Gx back to the global store ∗///Gy and Gz are obtained by an algorithm similar to Gx// Combine the values ofGx, Gy, and Gz to get GIn the case of each segmentG has been obtained, and the overall pressure of the underground horizontal wells of N2 segments and N3 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 isON3×N2×lgNx+lgNy+lgNz+lgN1+lgN2+lgN3, while that on the CPU is ON3×N2×N1×Nx+Ny+Nz+N2+N3. 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.


## 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. Figure5 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 calculateG 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 theN2 section.(5)
Call kernel 3 to carry out the reduction calculation on theN3 crack.(6)
Copy the calculated result in step5 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 calculatesG of each segment of the crack, which consists of three parts, namely, Gx, Gy, and Gz. The pseudocode under the CPU computation condition is as follows:Fork = 0 to N1-1Fornx = 1 to NxCalculateGx/∗ At the integration pointi, sum, calculate Gx ∗/Forny=1 to NyCalculateGy/∗ At the integration pointi, sum, calculate Gy ∗/Fornz=1 to NzCalculateGz/∗ At the integration pointi, sum, calculate Gz ∗/CalculateG/∗ Calculate the G 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 memoryunsigned int tid = blockIdx.x∗blockDim.x + threadIdx.x;unsigned int i = blockIdx.x∗(blockDim.x∗2) + threadIdx.x;double mySum = (i < n) ? g_idata[i] : 0;if (i + blockSize < n)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 synchronizationfor (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 singlethread bundle (Warp)∗/if (tid == 0) g_odata[blockIdx.x] = mySum;/∗ Write the final summation result Gx back to the global store ∗///Gy and Gz are obtained by an algorithm similar to Gx// Combine the values ofGx, Gy, and Gz to get GIn the case of each segmentG has been obtained, and the overall pressure of the underground horizontal wells of N2 segments and N3 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 isON3×N2×lgNx+lgNy+lgNz+lgN1+lgN2+lgN3, while that on the CPU is ON3×N2×N1×Nx+Ny+Nz+N2+N3. 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 Table1 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 7–12 under the time steps 90 h, 360 h, 980 h, 2580 h, 3650 h, and 4360 h, respectively.Table 1
Well parameters.
Permeability K = 2.8453 mD Well constants C = 0.9718  m 3 / M P a Wellbore epidermis S w = 1.5 Horizontal well length L = 120 0  m Boundary width X e = 1421.45 9  m Boundary height Y e = 812.26 2  m Horizontal well position X w = 91.70 9  m ,   Y w = 178.82 7  m Fracture half-length X f = 64.196 3  m Formation pressure P i = 20.0 4  M P a Formation thickness H = 3.2m Oil well radius R w = 0.1 m Porosity ϕ=0.09 Fluid viscosity μ = 4.5mPa⋅sFigure 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 Table1 (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.
M P a × 1.00 ∗ E 06 = Pa m D × 9.87 E - 14 = m2 mPa⋅s ×1.00∗ E - 03 = Pa⋅s m3/MPa ×1.00∗ E - 06 = m3/Pa C   ∗onversion factor is exact.Figures7–12 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.


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



---

*Source: 2582797-2018-04-24.xml*