Fully coupled fluid-solid productivity numerical simulation of multistage fractured horizontal well in tight oil reservoirs

Fig. 1. Multi-scale spatial schematic diagram of the tight oil reservoir (modified according to reference [16]).

1.1. Deformation model of reservoir rock

Regarding the rock mass as a continuous porous medium composed of rock skeleton phases and fluid phases. Under the assumption of small deformation, the rock mass mechanics governing equation of the tight oil reservoirs is ^[25]:

(1) $\sigma _{ij,j,\eta }^{{}}\text{+}f_{i,\eta }^{{}}\text{=}\mathbf{0}$

The deformation of the rock skeleton is mainly controlled by the effective stress. According to the Terzaghi’s effective stress principle, it can be expressed as ^[26]:

(2)$\sigma _{ij,\eta }^{{}}=\sigma _{ij,\eta }^{\prime }-{{b}_{\eta }}{{p}_{\eta }}{{\delta }_{ij}}$

Considering the stress balance condition, the stress and strain in the dual-media system satisfy the following relationship ^[27]:

(3)${{\sigma }_{ij}}=\sigma _{ij,\text{I}}^{{}}=\sigma _{ij,\text{II}}^{{}}$

(4)${{\varepsilon }_{ij}}=\varepsilon _{ij,\text{I}}^{{}}+\varepsilon _{ij,\text{II}}^{{}}$

Based on the assumption of isotropic elasticity, in the state of linear elastic deformation, the expression of the reservoir elastic constitutive equation can be obtained from Eqs. (2)-(4):

(5)${{\sigma }_{ij}}=D_{ijkl,\text{I},\text{II}}^{{}}{{\varepsilon }_{kl}}-D_{ijkl,\text{I},\text{II}}^{{}}C_{klrt,\text{I}}^{{}}{{b}_{\text{I}}}{{p}_{\text{I}}}{{\delta }_{rt}}-D_{ijkl,\text{I},\text{II}}^{{}}C_{klrs,\text{II}}^{{}}{{b}_{\text{II}}}{{p}_{\text{II}}}{{\delta }_{rt}}$

According to Eqs. (1)-(5), the governing equations for rock mass deformation of the tight oil reservoir can be written as follows:

(6)$\left( D_{ijkl,\text{I},\text{II}}^{{}}{{\varepsilon }_{kl}}-D_{ijkl,\text{I},\text{II}}^{{}}C_{klrs,\text{I}}^{{}}{{b}_{\text{I}}}{{p}_{\text{I}}}{{\delta }_{rt}}- \right.{{\left. D_{ijkl,\text{I},\text{II}}^{{}}C_{klrs,\text{II}}^{{}}{{b}_{\text{II}}}{{p}_{\text{II}}}{{\delta }_{rt}} \right)}_{,j}}\text{+}{{f}_{i}}\text{=}\mathbf{0}$

Under the assumption of small deformation, the relationship between the strain component of the reservoir and the displacement is as below:

(7)${{\varepsilon }_{ij}}=\frac{1}{2}\left( {{u}_{i,j}}+{{u}_{j,i}} \right)$

1.2. Fluid flow model

Aiming at the problem of saturated oil phase flow in porous media, the continuity equation of the matrix system and the natural fracture system is the mass conservation equation ^[28]:

(8)$\left\{ \begin{align} & \frac{d{{W}_{\text{I}}}}{dt}+\nabla {{H}_{\text{I}}}-\alpha \rho \frac{{{K}_{I}}}{\mu }\left( {{p}_{\text{I}}}-{{p}_{\text{II}}} \right)={{Q}_{\text{I}}} \\ & \frac{d{{W}_{\text{II}}}}{dt}+\nabla {{H}_{\text{II}}}+\alpha \rho \frac{{{K}_{I}}}{\mu }\left( {{p}_{\text{I}}}-{{p}_{\text{II}}} \right)={{Q}_{\text{II}}} \\ \end{align} \right.$

According to the theory of elastic porous media, the change of mass in each system can be written in differential form^[25]:

(9)$\frac{d{{W}_{\eta }}}{\rho }={{b}_{\eta }}d{{\varepsilon }_{v,\eta }}+\frac{d{{p}_{\eta }}}{{{M}_{b,\eta }}}$

Due to the extremely low permeability of the matrix system of the tight oil reservoir, it has obvious characteristics of non-linear seepage flow. The equation of motion of seepage flow considering the threshold pressure gradient is as below:

(10)${{V}_{\text{I}}}=\left\{ \begin{matrix} -\frac{{{K}_{\text{I}}}}{\mu }\left( \nabla {{p}_{\text{I}}}-G \right)\begin{matrix} {} & \left| \nabla {{p}_{\text{I}}} \right|>G \\ \end{matrix} \\ \begin{matrix} \begin{matrix} {} & {} \\ \end{matrix} & 0 & {} & {} \\ \end{matrix}\begin{matrix} {} & \left| \nabla {{p}_{\text{I}}} \right|\le G \\ \end{matrix} \\ \end{matrix} \right.$

The motion equation of the fracture system is expressed by Darcy's law:

(11)${{V}_{\text{II}}}=-\frac{{{K}_{\text{II}}}}{\mu }\nabla {{p}_{\text{II}}}$

Substituting Eqs. (9)-(11) into Eq. (8) can get a complete expression of the mass conservation equation of the matrix system and the natural fracture system ^[29]:

(12)$\left\{ \begin{align} & {{b}_{\text{I}}}\rho \frac{\partial {{\varepsilon }_{v,\text{I}}}}{\partial t}+\rho \frac{1}{{{M}_{b,\text{I}}}}\frac{\partial {{p}_{\text{I}}}}{\partial t}-\nabla \left[ \frac{\rho {{K}_{\text{I}}}}{\mu }\left( \nabla {{p}_{\text{I}}}-G \right) \right]- \\ & \text{ }\alpha \rho \frac{{{K}_{\text{I}}}}{\mu }\left( {{p}_{\text{I}}}-{{p}_{\text{II}}} \right)={{Q}_{\text{I}}} \\ & {{b}_{\text{II}}}\rho \frac{\partial {{\varepsilon }_{v,\text{II}}}}{\partial t}+\rho \frac{1}{{{M}_{b,\text{II}}}}\frac{\partial {{p}_{\text{II}}}}{\partial t}-\nabla \left( \frac{\rho {{K}_{\text{II}}}}{\mu }\nabla {{p}_{\text{I}}} \right)+ \\ & \text{ }\alpha \rho \frac{{{K}_{\text{I}}}}{\mu }\left( {{p}_{\text{I}}}-{{p}_{\text{II}}} \right)={{Q}_{\text{II}}} \\ \end{align} \right.$

As the main channel of the oil flow to the wellbore, artificial fractures can improve the overall seepage capacity of tight oil reservoirs. The fractures are treated as internal flow boundaries. The fluid flow continuity equation in the fractures is^[30]:

(13)${{s}_{^{F}}}\rho \frac{\partial {{p}_{F}}}{\partial t}-\nabla \left( \rho \frac{{{K}_{^{F}}}}{\mu }\nabla {{p}_{F}} \right)={{Q}_{F}}$

According to the discrete fracture model theory, the seepage flow equation of hydraulic fractures is processed by "dimension reduction" and assembled into the seepage flow equation of dual media to form a whole-region seepage flow equation system, which can be expressed as^[8]:

(14)$\int_{{{\Omega }_{\ e}}}{Q}d{{\Omega }_{\ e}}=\sum\limits_{\eta \in \left\{ I,II \right\}}{\int_{{{\Omega }_{\eta }}}{{{Q}_{\eta }}}d{{\Omega }_{\eta }}}+\sum\limits_{z}{{{w}_{F,z}}\int_{{{\Omega }_{\ F,z}}}{Q}d{{\Omega }_{\ F,z}}}$

1.3. Fluid-solid dynamic coupling model

The process of fluid-solid coupling includes: pressure change caused by fluid flow will affect the effective stress acting on the rock skeleton, resulting in rock mass deformation; rock mass deformation will change the rock pore space, the pore throat radius, and the fracture opening, etc. It will change the permeability, porosity and fracture conductivity of the rock mass, and at the same time affect the fluid flow (Fig. 2).

Fig. 2.

Fig. 2. Schematic diagram of the multi-physical field coupling relationship of tight oil reservoirs.

The matrix porosity is defined as the ratio of the current matrix pore space volume to the current total rock volume, so the matrix porosity after the reservoir deformation is ^[29]:

(15)${{\phi }_{I}}={{\phi }_{0,I}}+{{b}_{I}}{{\varepsilon }_{v,I}}+\frac{{{b}_{^{I}}}-{{\phi }_{0,I}}}{{{M}_{s,I}}}\left( {{p}_{\text{I}}}-{{p}_{0,I}} \right)$

According to Cubic's law, the relationship between matrix permeability and porosity satisfies ^[30]:

(16)${{K}_{I}}={{K}_{0,I}}{{\left( \frac{{{\phi }_{I}}}{{{\phi }_{0,I}}} \right)}^{3}}$

Since the porosity of the natural fracture system is extremely low, this paper ignores the influence of solid deformation on its porosity, and selects the model of Jiang et al. ^[31] that considers the average stress and geomechanical parameters to calculate the permeability of the natural fracture system:

(17)${{K}_{II}}={{K}_{0,II}}\exp \left[ {{c}_{II}}\frac{1+\upsilon }{1-\upsilon }{{b}_{II}}\left( {{p}_{II}}-{{p}_{0,II}} \right) \right]$

The opening of the artificial fracture is controlled by the effective stress on the fracture surface. The effective stress on the fracture surface is equal to the normal stress acting on the surface of the fracture minus the fluid pressure in the fracture. Therefore, the expression of the fracture opening is as follows ^[30]:

(18)${{w}_{F}}={{w}_{F,0}}exp\left[ -{{c}_{F}}\left( {{\sigma }_{\beta n}}-{{p}_{F}} \right) \right]$

Calculate the permeability of artificial fracture with parallel plate fracture seepage flow equation:

(19)${{K}_{F}}=\frac{w_{F}^{2}}{12}$

Eqs. (15)-(19) constitute the fluid-solid dynamic coupling model of matrix system, natural fracture system and artificial fracture system, which fully reflects the relationship between the seepage field and the stress field.

1.4. Initial and boundary conditions

The stress-displacement condition is defined as:

(20)$u\left| {{_{\Omega }}_{_{e}}} \right.={{\bar{u}}_{s}}$

(21)$\sigma \cdot n\left| {{_{\Omega }}_{_{e}}} \right.={{\bar{F}}_{s}}$

The outer boundary condition of the flow closure is:

(22)$\frac{\partial {{p}_{\eta }}}{\partial n}=0$

The initial boundary conditions are:

(23)${{p}_{\eta }}\left| _{{{\Omega }_{\text{e}}}} \right.\left( x,y,t=0 \right)={{p}_{0,\eta }}$

(24)$u\left| _{{{\Omega }_{\text{e}}}} \right.\left( x,y,t=0 \right)={{u}_{0}}$

(25)$\sigma \left| _{{{\Omega }_{\text{e}}}} \right.\left( x,y,t=0 \right)={{\sigma }_{0}}$

2. Model solution and validation

2.1. Discretization and solution

Eqs. (6), (12) and (13) can make the description of the governing equations of rock mass deformation and fluid flow, each of which is complex partial differential equation. The multi-physic field coupling problem is highly nonlinear, so it can only solve numerical solution. In this study, the finite element method is used to discretize the control equation, and the displacement and pressure in the control equation are expressed as the interpolation functions of node variables:

(26)$\left\{ \begin{matrix} u={{N}_{u}}\bar{u} \\ {{p}_{\eta }}={{N}_{p,\eta }}{{{\bar{p}}}_{\eta }} \\ {{p}_{F}}={{N}_{p,F}}{{{\bar{p}}}_{F}} \\ \end{matrix} \right.$

Based on the finite element method and discrete fracture model, considering the influence of fluid flow in artificial fracture and sink-source term, Eq. (6) of system mechanical equilibrium control equation and Eq. (12) and Eq. (13) of mass conservation equation in system can be discretized into the following matrix form in space ^[32]:

(27)$\left[ \begin{matrix} \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & {{H}_{I}}+H_{I,\tau }^{{}} & -H_{I,\tau }^{{}} \\ \mathbf{0} & -H_{II,\tau }^{{}} & {{H}_{II}}+H_{II,\tau }^{{}} \\ \end{matrix} \right]\left[ \begin{matrix} \overline{u} \\ {{\overline{p}}_{\text{I}}} \\ {{\overline{p}}_{II}} \\ \end{matrix} \right]+\left[ \begin{matrix} {{K}_{I,II}} & -{{L}_{I}} & -{{L}_{II}} \\ {{K}_{I}} & {{S}_{I}} & 0 \\ {{K}_{II}} & 0 & {{S}_{II}} \\ \end{matrix} \right]\frac{d}{dt}\left[ \begin{matrix} \overline{u} \\ {{\overline{p}}_{\text{I}}} \\ {{\overline{p}}_{II}} \\ \end{matrix} \right]=\left[ \begin{matrix} \frac{d{{R}_{\text{u}}}}{dt} \\ {{R}_{I}} \\ {{R}_{II}} \\ \end{matrix} \right]$

where ${{K}_{I\text{,}II}}=\int_{{{\Omega }_{e}}}{{{B}^{T}}}{{D}_{I\text{,}II}}Bd{{\Omega }_{e}}$

${{L}_{\eta }}=\int_{{{\Omega }_{e}}}{{{B}^{T}}}{{D}_{I,II}}{{C}_{\eta }}{{b}_{\eta }}m{{N}_{p,\eta }}d{{\Omega }_{e}}$

${{R}_{u}}=\int_{{{\Omega }_{e}}}{N_{u}^{T}{{f}_{i}}}d{{\Omega }_{\text{e}}}+\int_{{{\Omega }_{s}}}{N_{u}^{T}{{{\bar{F}}}_{s}}}d{{\Omega }_{s}}$

${{H}_{I}}=\int_{{{\Omega }_{e}}}{\rho \frac{{{K}_{I}}}{\mu }\left[ {{\left( \nabla {{N}_{p,I}} \right)}^{T}}\nabla {{N}_{p,I}}-GN_{p,I}^{T}\nabla {{N}_{p,I}} \right]d{{\Omega }_{e}}}$

${{H}_{II}}=\int_{{{\Omega }_{e}}}{{{\left( \nabla {{N}_{p,II}} \right)}^{T}}\rho \frac{{{K}_{II}}}{\mu }}\left( \nabla {{N}_{p,II}} \right)d{{\Omega }_{e}}+$

$\int_{{{\Omega }_{F}}}{{{w}_{F}}{{\left( \nabla {{N}_{p,F}} \right)}^{T}}\rho \frac{{{K}_{F}}}{\mu }}\left( \nabla {{N}_{p,F}} \right)d{{\Omega }_{F}}$

$H_{\eta,\tau }^{{}}=\int_{{{\Omega }_{e}}}{N_{p,\eta }^{T}\alpha \rho \frac{{{K}_{I}}}{\mu }{{N}_{p,\eta }}}d{{\Omega }_{e}}$

${{K}_{\eta }}=\int_{{{\Omega }_{e}}}{N_{p,\eta }^{T}{{b}_{\eta }}\rho }{{m}^{T}}{{C}_{\eta }}{{D}_{I,II}}Bd{{\Omega }_{e}}$

${{S}_{\eta }}=\int_{{{\Omega }_{\eta }}}{N_{p,\eta }^{T}\rho \frac{1}{{{M}_{b,\eta }}}}{{N}_{p,\eta }}d{{\Omega }_{\eta }}+\int_{{{\Omega }_{F}}}{{{w}_{F}}N_{p,F}^{T}{{s}_{F}}{{N}_{p,F}}}d{{\Omega }_{F}}$

${{R}_{\text{I}}}=\int_{{{\Gamma }_{\text{m}}}}{N_{p,I}^{T}{{Q}_{\text{I}}}}d{{\Gamma }_{m}}$

${{R}_{\text{II}}}=\int_{{{\Gamma }_{\text{f}}}}{N_{p,II}^{T}{{Q}_{II}}}d{{\Gamma }_{\text{f}}}+\int_{{{\Gamma }_{F}}}{{{w}_{F}}N_{p,F}^{T}{{Q}_{F}}}d{{\Gamma }_{F}}$

The fluid-solid fully coupled finite element equilibrium equation is discretized in time domain. The three main variables of the governing equation (reservoir node displacement, matrix system pressure and natural fracture system pressure) are solved simultaneously at the same time step in an implicit method. The secondary variables related to the primary variables and their coupling relationship can be updated iteratively at each time step to obtain the values of system displacement and pressure. The specific solution flow of fluid-solid coupling is shown in Fig. 3: (1) The model initialization part includes flow initialization and geomechanical initialization, the fluid- pressure field and rock stress state at different positions in space at t=0 can be obtained. (2) Fluid-solid fully coupling. The solving process is divided into several time step increments, the time step forward method is used to solve the pressure field and displacement field, and check their convergence. If not, return to repeated iteration. (3) The displacement and pressure (main variables) obtained in Step 2 can be used to calculate porosity, permeability and fracture conductivity (secondary variables) at the same time step, which are updated iteratively at each time Step (4) Repeat steps 2-3 until the simulation time is over.

Fig. 3.

Fig. 3. Flow chart of fluid-solid coupling solution.

2.2. Model validation

In order to verify the correctness of the above fluid- solid fully coupling model, a multistage fractured horizontal well model of the tight oil reservoir has been established. The basic parameters of the model are shown in Table 1. The values of numerical simulation parameters are mainly referred to references [5], [8] and [33]. The model is set with triangular mesh, in which the matrix system and natural fracture system share a set of mesh. The dimension of fracturing fracture is reduced in the process of meshing. Due to the large pressure gradient around the fracture, grid refinement is carried out around the artificial fracture. Thus, the grid division diagram of the reservoir area in whole model is obtained (Fig. 4). The simulation results of the fully coupled rigid model are obtained by the numerical calculation program and compared with the daily oil production and cumulative oil production obtained by the finite difference model used by the commercial software Eclipse (Fig. 5). The numerical simulation results of the fully coupled rigid model in this paper is in good agreement with the simulation results of Eclipse. At the same time, in two-dimensional homogeneous, linear elastic and impermeable rocks, the opening distribution of uniform compression fractures along the fracture length satisfies the analytical solution. The fracture opening calculated by this model and Jiang et al. ^[34] model is compared with the analytical solution (Fig. 6). It can be seen that fracture opening calculated by the three models is similar. The result of our model is closer to the analytical solution. The rationality and correctness of the fluid-solid coupling model are verified from the above two aspects in this paper.

Table 1. Basic model parameters

Parameter	Value	Parameter	Value	Parameter	Value
Reservoir size	1000 m× 600 m×18 m	Matrix initial permeability	0.12×10^-3μm²	Matrix elastic modulus	18 GPa
Horizontal well length	600 m	Natural fracture initial porosity	0.08%	Natural fracture elastic modulus	12 GPa
Artificial fracture length	160 m	Natural fracture initial permeability	1.5×10^-3μm²	Poisson ratio	0.27
Artificial fracture height	18 m	Natural fracture compression coefficient	0.15 MPa^-1	Maximum horizontal principal stress	38.5 MPa
Artificial fracture opening	2.884×10^-4m	Wellbore radius	0.05 m	Minimum horizontal principal stress	35.5 MPa
Artificial fracture number	7	Initial formation pressure	30 MPa	Channeling coefficient	1
Artificial fracture compression coefficient	0.06 MPa^-1	Bottom hole pressure	20 MPa	Biot coefficient	1
Artificial fracture porosity	100%	Fluid viscosity	8 mPa•s	Fluid bulk modulus	1×10⁴MPa
Matrix initial porosity	10%	Fluid density	960 kg/m³	Rock skeleton bulk modulus	6.5×10⁴MPa

Fig. 4.

Fig. 4. Schematic diagram of the grid division of the reservoir area.

Fig. 5.

Fig. 5. Comparison of simulation results between fully coupled rigid model and Eclipse model.

Fig. 6.

Fig. 6. Comparison of fracture opening distribution calculated by different models and analytical solutions.

2.3. Model comparison

We compared the production of multi-stage fractured horizontal wells simulated and calculated by fully coupled model (Model 1), matrix + natural fracture rigid model (Model 2), artificial fracture rigid model (Model 3) and fully rigid model (Model 4) (Fig. 7). Taking the cumulative oil production of the fully coupled model as the base, we calculated the relative error between the cumulative oil production of other models and the base. The results are shown in Table 2. It can be seen that there is a large difference in the predicted production between the fully coupled model and the fully rigid model. At the early stage of development, the daily oil production calculated by the fully coupled model decreased rapidly, the cumulative production calculated by the fully coupled model was the smallest after 3000 d, and the cumulative production calculated by the fully rigid model was the highest. Ignoring the fluid-solid coupling effect between matrix and natural fracture system, the prediction error of 3000 d production of fractured horizontal wells was 0.92%, which was mainly due to the gradual deterioration of physical properties of matrix and natural fracture system during the actual production process, and the gradual weakening of its supply capacity to horizontal wellbore, which was not considered in Model 2. Therefore, the cumulative production calculated by Model 2 is larger than that of Model 1. Ignoring the fluid-solid coupling effect of the artificial fracture system, the prediction error of 3000 d production of fractured horizontal wells was 36.46%, mainly due to the strong conductivity of the artificial fracture, the gradual decrease of the pressure in the fracture and the continuous decrease of the opening in the actual production process, which was not considered in Model 3. Therefore, the cumulative production calculated by Model 3 is greater than that of Model 1. Comparing the relative error values of predicted production with Model 2-Model 4, it can be seen that the relative error (38.30%) of Model 4 is not a simple linear superposition (37.38%) of the relative errors of Model 2 and Model 3, which further shows that the effects of matrix, natural fracture and artificial fracture are coupled with each other.

Fig. 7.

Fig. 7. Production of multi-stage fractured horizontal wells calculated by different models.

Table 2. Comparison of cumulative oil production calculated by different models

Model number	Model	Cumulative oil production/m³	Relative error/%
1	Fully coupled model	31 454.76	0
2	Matrix + natural fracture rigid model	31 742.79	0.92
3	Artificial fracture rigid model	42 922.81	36.46
4	Fully rigid model (uncoupled)	43 502.75	38.30

In conclusion, the influence of fluid-solid coupling on tight oil production prediction cannot be ignored. The productivity of multi-stage fractured horizontal wells predicted by uncoupled model is higher, while the productivity predicted by fluid-solid fully coupled model is more accurate. The physical property of artificial fracture is an important factor to determine the productivity of fractured horizontal wells. It is necessary to optimize the design of artificial fracture parameters.

3. Evolution law of reservoir physical properties

Using the numerical model in the model validation part, the permeability distribution of matrix system and natural fracture system at different times of multi-stage fractured horizontal well in tight oil reservoir after put into production were calculated (Figs. 8-11). When the fluid flew into the artificial fracture and wellbore through the matrix system and natural fracture system, due to the strong conductivity and small flow resistance in the artificial fracture, there was a large fluid pressure gradient near the artificial fracture area, and the reservoir permeability changed significantly. After 2000 d of production of the horizontal well, the matrix system permeability changed from 0.120×10^-3μm² to the minimum about 0.117×10^-3μm², the permeability of natural fracture system changed from 1.500×10^-3μm² to the minimum about 1.070×10^-3μm². Because the matrix system was denser, the variation range of permeability was less than that of natural fracture system. The minimum value of reservoir permeability was distributed around the artificial fracture. With the extension of time, the range of permeability reduction gradually expanded from near-artificial fracture area to the surrounding reservoir.

Fig. 8.

Fig. 8. Permeability distribution of matrix system after 100 d of production.

Fig. 9.

Fig. 9. Matrix system permeability distribution after 2000 d of production.

Fig. 10.

Fig. 10. Permeability distribution of natural fractures after 100 d of production.

Fig. 11.

Fig. 11. Permeability distribution of natural fractures after 2000 d of production.

After 2000 d of production of multi-stage fractured horizontal wells in tight oil reservoirs, the opening of artificial fractures was minimized from 2.884×10^-4m to about 1.440×10^-4m (Figs. 11-12). Fig. 14 shows the loss extent curves of artificial fracture’s opening and conductivity. At the early stage of development, due to the rapid production of fluid in the artificial fractures, the pressure near the artificial fractures decreased rapidly. As a result, the loss extent of artificial fracture’s opening and conductivity was large. After 500 d of development, the loss extent of artificial fracture’s opening was 40.37%, and the loss extent of artificial fracture conductivity was 78.80%. Then the artificial fracture’s physical property loss extent increased slowly. After 3000 d, the artificial fracture’s opening loss extent reached 52.12%, while the conductivity loss extent reached 89.02%.

Fig. 12.

Fig. 12. Artificial fracture opening distribution after 100 d of production.

4. Optimization of fracturing parameters

4.1. Parameter sensitivity

Design Case 1-Case 9, the artificial fracture compression coefficients were 0.060, 0.062, 0.064 MPa^-1, the artificial fracture openings were 2.884×10^-4, 2.978 × 10^-4, 3.065×10^-4 m, the compression coefficients of natural fracture were 0.15, 0.25, 0.35 MPa^-1, the threshold pressure gradient was 0, 0.05, 0.10 MPa/m, respectively. The cumulative oil production of Case 1 was the base. The cumulative production of 9 cases of multi-stage fractured horizontal wells was simulated and calculated. The relative error between the cumulative oil production of the other 8 cases and the base was calculated. The results are shown in Table 3. The comparison results of the cumulative oil production after 3000 d under different parameters are shown in Fig. 13 to Fig. 18. Case 5 has the largest cumulative oil production while Case 9 has the smallest cumulative oil production. The production of horizontal well is most sensitive to threshold pressure gradient, followed by artificial fracture opening. The threshold pressure gradient slows down the propagation of pressure waves in the reservoir and inhibits the production of tight oil well. Increasing the initial conductivity of artificial fractures will help reduce the production loss. Meanwhile, the cumulative production of horizontal well decreases with the increase of compression coefficients of artificial fracture and natural fracture.

Fig. 13.

Fig. 13. Artificial fracture opening distribution after 2000 d of production.

Fig. 14.

Fig. 14. Curves of loss extent of artificial fracture opening and conductivity at different time.

Fig. 15.

Fig. 15. Cumulative oil production under different artificial fracture compressibility factors.

Fig. 16.

Fig. 16. Cumulative oil production under different artificial fracture openings.

Fig. 17.

Fig. 17. Cumulative oil yield under different natural fracture compressibility factors.

Fig. 18.

Fig. 18. Cumulative oil production under different threshold pressure gradients.

Table 3. The influence of different sensitivity parameters on production

Case number	Artificial fracture compression coefficient/MPa^-1	Artificial fracture opening/10^-4m	Natural fracture compression coefficient/MPa^-1	Threshold pressure gradient/(MPa•m^-1)	Cumulative oil production/m³	Relative error/%
1	0.060	2.884	0.15	0	31 454.76	0
2	0.062	2.884	0.15	0	30 809.01	-2.05
3	0.064	2.884	0.15	0	30 171.76	-4.08
4	0.060	2.978	0.15	0	32 368.10	2.90
5	0.060	3.065	0.15	0	33 179.36	5.48
6	0.060	2.884	0.25	0	31 629.62	-0.59
7	0.060	2.884	0.35	0	31 077.90	-1.20
8	0.060	2.884	0.15	0.05	26 887.96	-14.52
9	0.060	2.884	0.15	0.10	21 600.05	-31.33

4.2. Optimization of fracture parameters

The geometric distribution of artificial fractures is one of the important factors affecting the productivity of tight oil reservoirs. In this study, while keeping the product of the artificial fracture number and its length constant, we designed 11 cases (Case 10-Case 20). The production of multi-stage fractured horizontal wells under different fracturing parameters was calculated respectively. The cumulative oil production of Case 10 was used as the base to calculate the relative error between the cumulative oil production of other cases and the base. The case design and simulation results are shown in the Table 4. The cumulative oil production of 3000 d in different cases was compared, and the results are shown in Fig. 19. The Case 20 has the largest cumulative oil production, and Case 13 has the smallest cumulative oil production. With the increase in the number of artificial fractures, the horizontal well production continued to increase. The cumulative oil production curve of each case is close to parallel at the later stage of production. The impact of artificial fracture densification on the production of horizontal wells is mainly concentrated at the early stage of production. Fig. 20 shows the relationship between the different number of artificial fractures and the relative error of the predicted production. Take the first 5 points and the last 6 points of the output prediction relative error curve respectively. We made a fitting curve to form an intersection diagram. When the number of fractures was 8-9, the output prediction relative error curve had an inflection point. When the number of fractures was over 10, the increase in the production of horizontal wells weakened, and the influence degree of the number of fractures on the production also decreased gradually. During the optimal design of fracturing construction, it is necessary to consider the comprehensive influence of artificial fracture spacing, number and length. Only by increasing the number of fractures, we cannot achieve the expected increase in production. The optimal number of fractures is 9 under the simulation conditions in this study.

Fig. 19.

Fig. 19. Comparison of cumulative oil yield of multi- stage fractured horizontal wells in different cases.

Fig. 20.

Fig. 20. The relationship curve between the number of artificial fractures and the relative error of predicted production.

Table 4. The influence of different fracture parameters on production

Case number	Artificial fracture number	Artificial fracture length/m	Artificial fracture spacing/m	Cumulative oil production/m³	Relative error/%
10	7	160.00	100.000	31 454.76	0
11	6	186.70	120.000	29 655.21	-5.72
12	5	224.00	150.000	27 332.83	-13.10
13	4	280.00	200.000	24 330.56	-22.65
14	8	140.00	85.714	33 234.68	5.66
15	9	124.44	75.000	34 127.34	8.50
16	10	112.00	66.660	34 867.50	10.85
17	11	101.82	60.000	35 361.73	12.42
18	12	93.33	54.540	35 765.05	13.70
19	13	86.15	50.000	35 992.17	14.42
20	14	80.00	46.150	36 183.16	15.03

5. Conclusions

A fully coupled numerical model of porous medium deformation and fluid flow in multi-stage fractured horizontal wells in tight oil reservoirs has been established. The model has considered comprehensively the deformation characteristics of reservoir matrix system, the natural fractures, the artificial fractures, and the zonal fluid flow law. It can be used to simulate and predict accurately the productivity of multi-stage fractured horizontal wells in tight oil reservoirs.

The model calculation results show that the reservoir physical properties in the area near the artificial fractures become worse during the production of horizontal wells in tight oil reservoirs, especially when the loss extent of artificial fracture’s opening and conductivity reached 52.12% and 89.02%, respectively. The influence of fluid- solid coupling on productivity prediction of tight oil reservoir cannot be ignored. After 3000 d of production with horizontal wells in tight oil reservoirs, the error of predicted productivity between fully coupled model and uncoupled model was up to 38.30%. The parameter sensitivity analysis shows that the productivity of horizontal wells in tight oil reservoirs is the most sensitive to the threshold of the pressure gradient, followed by the opening of artificial fracture. By improving the initial conductivity of artificial fractures, it is helpful to improve the productivity of tight oil wells. The comprehensive influence of artificial fracture conductivity, spacing, quantity and length should be considered in the fracturing construction design of tight oil reservoir. The one-sided pursuit of increasing the number of fractures cannot achieve the expected stimulation effect.

Nomenclature

b—Biot’s correction coefficient;

B—strain matrix;

c—compression coefficient, Pa^-1;

C—flexibility matrix;

C_klrt_,I—matrix system flexibility tensor, Pa^-1;

C_klrt_,II—natural fracture system flexibility tensor, Pa^-1;

D—elastic matrix;

D_ijkl_,I,II—elastic tensor of dual media system, Pa;

f_i—gravity term vector, N/m³;

${{\bar{F}}_{s}}$—known stress at the boundary, Pa;

G—threshold pressure gradient, Pa/m;

h—artificial fracture surface height, m;

H—fluid flux, kg/(m²•s);

k—iterative step;

K—permeability, m²;

K₀—initial permeability, m²;

m—unit tensor;

M_b—Biot’s modulus, Pa;

M_s,I—matrix system rock skeleton bulk modulus, Pa;

n—boundary normal direction;

n—boundary unit normal vector;

N_p—shape function of pressure field;

N_u—shape function of displacement field;

p—pore pressure, Pa;

p₀—initial pore pressure, Pa;

$\bar{p}$—node pressure, Pa;

Q—source sink term, kg/(m³•s);

Q—flow equation;

s_F—artificial fracture storage coefficient, Pa^-1;

t—time, s;

t_s—time step;

u₀—system initial displacement, m;

u—displacement vector, m;

u_i,j—the partial derivative of u_i with respect to direction j;

$\bar{u}$—node displacement, m;

${{\bar{u}}_{s}}$—known displacement at the boundary, m;

V—fluid velocity, m/s;

w_F—artificial fracture width, m;

w_F_,z—artificial fracture width of z unit, m;

w_F,0—artificial fracture initial width, m;

W—mass of fluid per unit volume, kg/m³;

x,y—reservoir coordinate position, m;

z—fracture unit;

α—channeling coefficient, m^-2;

▽—Hamiltonian operator;

μ—fluid viscosity, Pa•s;

υ—rock Poisson's ratio;

δ—Kronecker symbol;

ε—strain tensor;

ε_ij—strain component;

ε_v—rock mss strain;

ρ—fluid density, kg/m³;

σ_βn—normal stress on fracture surface, Pa;

σ—total stress tensor;

σ₀—system initial stress, Pa;

σ_ij—stress component, Pa;

σ_ij_,_j—the partial derivative of σ_ijwith respect to direction j, Pa/m;

${{{\sigma }'}_{ij}}$—components of the effective stress tensor, Pa;

Γ_m—matrix system boundary;

Γ_f—natural fracture system boundary;

Γ_F—artificial fracture system boundary;

Ω_e—reservoir system;

Ω_F—artificial fracture system;

Ω_s—reservoir boundary;

ϕ—porosity, %;

ϕ₀—initial porosity, %;

F—artificial fracture

i,j,k,l,r,s—indicators in Voigt notation;

η—matrix system or natural fracture system;

τ—interporosity flow;

I—matrix system;

II—natural fracture system.

Reference

By original order

By published year

By cited within times

By Impact factor

[1]

SUN

Longde

, ZOU

Caineng

, JIA

Ailin

, et al.

Development characteristics and orientation of tight oil and gas in China

Petroleum Exploration and Development, 2019, 46(6): 1015-1026.

[2]

ZOU

Caineng

, ZHAI

Guangming

, ZHANG

Guangya

, et al. Formation,

distribution, potential and prediction of global conventional and unconventional hydrocarbon resources

Petroleum Exploration and Development, 2015, 42(1): 13-25.

[3]

Suyun

, ZHU

Rukai

, WU

Songtao

, et al.

Profitable exploration and development of continental tight oil in China

Petroleum Exploration and Development, 2018, 45(4): 737-748.

[4]

Guoxin

, LUO

Kai

, SHI

Deqin

. Key

technologies

engineering management and important suggestions of shale oil/gas development: Case study of a Duvernay shale project in Western Canada Sedimentary Basin

Petroleum Exploration and Development, 2020, 47(4): 739-749.

[5]

REN

, SU

, ZHAN

, et al.

Fully coupled fluid-solid numerical simulation of stimulated reservoir volume (SRV)-fractured horizontal well with multi-porosity media in tight oil reservoirs

Journal of Petroleum Science and Engineering, 2019, 174: 757-775.

DOI:10.1016/j.petrol.2018.11.080 URL [Cited within: 4]

[6]

FANG

Wenchao

, JIANG

Hanqiao

, LI

Junjian

, et al.

A numerical simulation model for multi-scale flow in tight oil reservoirs

Petroleum Exploration and Development, 2017, 44(3): 415-422.

[7]

REN

Long

, SU

Yuliang

, HAO

Yongmao

, et al.

Dynamic analysis of SRV-fractured horizontal wells in tight oil reservoirs based on stimulated patterns

Acta Petrolei Sinica, 2015, 36(10): 1272-1279.

[8]

ZHANG

, ZHANG

, TANG

, et al.

A novel fluid-solid coupling model for the oil-water flow in the natural fractured reservoirs

Physics of Fluids, 2021, 33(3): 036601.

DOI:10.1063/5.0041267 URL [Cited within: 5]

[9]

CRYER

C W

A comparison of the three-dimensional consolidation theories of Biot and Terzaghi

The Quarterly Journal of Mechanics and Applied Mathematics, 1963, 16(4): 401-412.

DOI:10.1093/qjmam/16.4.401 URL [Cited within: 1]

[10]

WANG

Ziming

, DU

Zhimin

The coupled model and numerical simulation of multiphase flow in an elastoplastic deforming oil reservoir with transformation temperature

Petroleum Exploration and Development, 2001, 28(6): 68-72.

[11]

HELMIG

, FLEMISCH

, WOLFF

, et al.

Model coupling for multiphase flow in porous media

Advances in Water Resources, 2013, 51: 52-66.

DOI:10.1016/j.advwatres.2012.07.003 URL [Cited within: 1]

[12]

TIAN

, XIONG

, WANG

, et al.

A compositional model for gas injection IOR/EOR in tight oil reservoirs under coupled nanopore confinement and geomechanics effects

SPE-193818, 2019.

[13]

LIANG

Tao

, CHANG

Yuwen

, GUO

Xiaofei

, et al.

Influence factors of single well’s productivity in the Bakken tight oil reservoir

Petroleum Exploration and Development, 2013, 40(3): 357-362.

[14]

LEI

Qun

, WENG

Dingwei

, XIONG

Shengchun

, et al.

Progress and development directions of shale oil reservoir stimulation technology of China National Petroleum Corporation

Petroleum Exploration and Development, 2021, 48(5): 1035-1042.

[15]

Guoxin

, QIN

Jianhua

, XIAN

Chenggang

, et al.

Theoretical understandings, key technologies and practices of tight conglomerate oilfield efficient development: A case study of the Mahu oilfield, Junggar Basin, NW China

Petroleum Exploration and Development, 2020, 47(6): 1185-1197.

[16]

Jianchun

Production performance analysis and numerical simulation for multistage fractured horizontal well

Qingdao: China University of Petroleum, 2017.

[Cited within: 2]

[17]

FANG

Pingliang

Study on mechanism and numerical simulation method of multiphase flow-geomechanical deformation coupling in tight oil reservoirs

Beijing: China University of Geosciences, 2017.

[18]

ZHANG

, ZHANG

, WANG

, et al.

Simulation of a multistage fractured horizontal well with finite conductivity in composite shale gas reservoir through finite- element method

Energy & Fuels, 2016, 30(11): 9036-9049.

DOI:10.1021/acs.energyfuels.6b01565 URL [Cited within: 1]

[19]

ZHANG

, ZHANG

, WANG

, et al.

Simulation of a multistage fractured horizontal well in a water-bearing tight fractured gas reservoir under non-Darcy flow

Journal of Geophysics and Engineering, 2018, 15(3): 877-894.

DOI:10.1088/1742-2140/aaa5ce URL [Cited within: 1]

[20]

LIU

, LIU

, LEUNG

, et al.

Coupled flow/geomechanics modeling of interfracture water injection to enhance oil recovery in tight reservoirs

SPE Journal, 2021, 26(1): 1-21.

DOI:10.2118/199983-PA URL

[21]

SETTARI

, WALTERS

D A

Advances in coupled geomechanical and reservoir modeling with applications to reservoir compaction

SPE Journal, 2001, 6(3): 334-342.

DOI:10.2118/74142-PA URL [Cited within: 1]

[22]

SETTARI

, MOURITS

F M

Coupling of geomechanics and reservoir simulation models

Balkema: Routledge, 1994.

[23]

ZHANG

Ruihan

Numerical simulation of multi-stage fractured horizontal well in shale gas reservoir based on multi-scale flow theory

Chengdu: Southwest Petroleum University, 2019.

[24]

ZHANG

, WU

, ZHAO

, et al.

Numerical simulation of the feasibility of supercritical CO₂ storage and enhanced shale gas recovery considering complex fracture networks

Journal of Petroleum Science and Engineering, 2021, 204: 108671.

DOI:10.1016/j.petrol.2021.108671 URL [Cited within: 1]

[25]

COUSSY

Poromechanics. New York: John Wiley & Sons Ltd., 2004.

[Cited within: 2]

[26]

TERZAGHI

Theoretical soil mechanics. New York: John Wiley & Sons Ltd., 1943.

[27]

NAIR

, ABOUSLEIMAN

, ZAMAN

A finite element porothermoelastic model for dual-porosity media

International Journal for Numerical and Analytical Methods in Geomechanics, 2004, 28(9): 875-898.

DOI:10.1002/nag.336 URL [Cited within: 1]

[28]

NAIR

R S

The poromechanics of naturally fractured rock formations: A finite element approach

Norman: The University of Oklahoma, 2003.

[29]

, LI

, ZHANG

A fully coupled thermo-hydro-mechanical, three-dimensional model for hydraulic stimulation treatments

Journal of Natural Gas Science and Engineering, 2016, 34: 64-84.

DOI:10.1016/j.jngse.2016.06.046 URL [Cited within: 2]

[30]

FAN

, LI

, SHAH

, et al.

Analysis of a fully coupled gas flow and deformation process in fractured shale gas reservoirs

Journal of Natural Gas Science and Engineering, 2015, 27(Part 2): 901-913.

DOI:10.1016/j.jngse.2015.09.040 URL [Cited within: 3]

[31]

JIANG

, YOUNIS

R M

A multimechanistic multicontinuum model for simulating shale gas reservoir with complex fractured system

Fuel, 2015, 161: 333-344.

DOI:10.1016/j.fuel.2015.08.069 URL [Cited within: 1]

[32]

PENG

, ZHANG

Engineering geology for underground rocks. Berlin, Heidelberg: Springer, 2007.

[33]

ZHAO

, LIU

, ZHANG

, et al.

Simulation of a multistage fractured horizontal well in a tight oil reservoir using an embedded discrete fracture model

Energy Science & Engineering, 2019, 7(5): 1485-1503.