A sand paving optimization design method for multi-scale complex fracture networks

Abstract

A multi-scale complex fracture network sand paving optimization design method belongs to the field of petroleum and natural gas engineering technology. It carries out three-dimensional fracture expansion numerical simulation by constructing a fine geomechanical model, and on this basis obtains the actual expansion morphology of the main fracture and branch fractures, and on this basis carries out the design of the conductivity experimental scheme, and further takes the production parameter as the optimal, obtains the conductivity ratio of the main fracture and branch fracture and the respective optimal fracture conductivity, and reversely calculates the lowest fracture sand paving concentration based on the optimal fracture conductivity, and realizes the sand paving design of the complex fracture network from the dual factors of economy and production. The present invention takes the complex fracture network of volume fracturing of unconventional oil and gas reservoirs as the research object, and optimizes the design of the conductivity of the fracture network sand paving to improve the conductivity of the complex fracture network and ensure the contribution of the complex fracture network to production. The invention has practical application value for the efficient fracturing design of unconventional reservoirs.

Description

Technical Field

[0001] The invention relates to a multi-scale complex fracture network sand paving optimization design method, belonging to the technical field of petroleum and natural gas engineering. Background Art

[0002] Volume fracturing is currently the main mode of unconventional oil and gas reservoir development. During volume fracturing, it is usually necessary to fully break up the reservoir and improve the permeability of the reformed area by forming a high-conductivity fracture network to achieve economic development. The design of volume fracturing must take into account the optimal fracturing size under existing economic and technical conditions and the appropriate fracture reformation volume in the reservoir. The fracture network formed during the volume fracturing process usually has multi-scale characteristics, that is, it is composed of main fractures with high conductivity and secondary fractures with relatively weak conductivity. Before carrying out volume fracturing, it is necessary to optimize the sanding parameters of the fracture proppant parameters to meet the long-term contribution of fracture conductivity to production.

[0003] The existing optimization design ideas for proppant laying usually use production numerical simulation to achieve maximum production or maximize economic benefits through economic evaluation. Generally speaking, the maximization of production after fracturing is used as the objective function, or multi-objective functions are used for optimization. However, these methods are mainly aimed at double-wing symmetrical fractures with relatively simple fracture morphology, which are quite different from the actual complex fracture network morphological characteristics. It is impossible to guarantee the actual effectiveness of the prediction results of these methods in complex fracture networks, thus causing field engineering problems such as sand plugging. At the same time, the designed sand laying scheme is often difficult to implement, and the proppant laying in the secondary fractures is not ideal, which significantly reduces the field effect of volume fracturing development.

[0004] The conductivity of proppant fractures is affected by a combination of factors. Generally speaking, the concentration of sand, particle size of proppant, mechanical parameters of rock and stress environment will affect the conductivity of fractures. In order to obtain the best fracture proppant conductivity that meets the formation production conditions, the API conductivity test can be completed. However, the experiment mainly focuses on the fracture conductivity of the main fracture as the main research object during the analysis and testing, ignoring the contribution and influence of branch fractures on production. However, with the continuous development of deep oil and gas, the closure stress conditions of the reservoir are relatively high, and how to improve the conductivity of branch fractures has become extremely important for the increase of the entire oil and gas production. Summary of the invention

[0005] In response to the problem that there are too many assumed parameters in the existing volume fracturing sand paving design optimization process, which is greatly different from the actual situation, the present application provides a multi-scale complex fracture network sand paving optimization design method. This method constructs a fine geomechanical model to carry out three-dimensional fracture expansion numerical simulation, and on this basis obtains the actual expansion morphology of the main fracture and branch fractures, and on this basis carries out the design of the conductivity experiment scheme. Further, taking the production parameter as the optimal, the conductivity ratio of the main fracture and branch fracture and the respective optimal fracture conductivity are obtained, and the lowest fracture sand paving concentration is obtained by inverse calculation based on the optimal fracture conductivity, so as to realize the sand paving design of the complex fracture network from the dual factors of economy and production.

[0006] The technical solution of the present invention is as follows:

[0007] A multi-scale complex crack network sand paving optimization design method includes the following steps:

[0008] S1: Constructing a one-dimensional geomechanical model of the reservoir: First, the corresponding rock mechanical parameters are obtained through the results of indoor core test experiments, combined with the logging interpretation results of acoustic wave time difference, gamma value, density, seismic inversion data, uniaxial and triaxial compression test results, that is, the one-dimensional geomechanical model is obtained;

[0009] The construction of a three-dimensional geomechanical model is carried out by combining the three-dimensional seismic data volume or experimental data points as the benchmark for constraints: the three-dimensional modeling of reservoir geomechanical parameters is carried out by using the sequential-Gaussian interpolation method in combination with the sedimentary phase model of the research block; the three-dimensional modeling of a series of strength parameters such as tensile strength and fracture toughness is completed by fitting formula calculation; the elastic model and Poisson's ratio are calculated as follows:

[0010]

[0011]

[0012] Where: E d is the dynamic elastic modulus of rock, GPa; μ d is the dynamic Poisson's ratio of rock; v p 、v s are the velocity of transverse and longitudinal waves, respectively, km / s; ρ is the density of rock, g / cm 3 ;

[0013] S2: Use the Karhunen-Loeve orthogonal series expansion method to invert the random three-dimensional geomechanical model, and use the simulated annealing method to optimize the inversion results of the ground stress, obtain the spatial distribution characteristics of the three-dimensional shale stress field, and check and correct the three-dimensional geomechanical model obtained in S1; preferably, the use of the simulated annealing algorithm optimization can better avoid the inversion results from falling into a local range. The specific steps are as follows:

[0014] (a) Determine the range of variation of model parameters and set an initial model m 0 , calculate the square value of the corresponding travel time residual's second norm, that is, the error equation E(m 0 );

[0015] (b) Perturb each parameter of the current model to generate a new model m 1 , calculate the corresponding error equation E(m 1 ), and calculate △E=E(m 1 )-E(m 0 );

[0016] (c) If △E < 0, the new model is accepted and assigned to the old model, that is, m 0 =m 1 ; If △E ≥ 0, then the new model m 1 Assign values ​​according to probability P = exp(-△E / T), where exp is the natural exponent and T is the temperature. A number R is randomly generated between 0 and 1. When P(△E)>R, the modification is accepted, i.e. m 0 =m 1 ;

[0017] (d) repeating steps (b) and (c) at temperature T;

[0018] (e) slowly lowering the temperature T;

[0019] (f) Repeat steps (b)-(e) until △E satisfies the specified convergence condition.

[0020] S3: Based on the modified three-dimensional geomechanical model, the commercial three-dimensional fracture extension software is used to carry out the numerical simulation of hydraulic fracturing three-dimensional fracture extension, and the Slab grid method is used to calculate the overall fracture transformation volume. This equation is mainly used to estimate the equivalent volume opening of natural fractures containing fracturing fluid. The specific calculation formula is shown in (3):

[0021]

[0022] Where: e is the crack opening, P pore is the fracture pore pressure; σ is the fracture surface normal stress; E is Young's modulus; ν is Poisson's ratio; r is the flow distance; r max is the maximum flow distance;

[0023] The actual geometric morphology of the complex fracture network after volume fracturing was further analyzed, and the relevant information of the main fracture and branch fractures, including the fracture length, width and height, was statistically obtained.

[0024] S4: Classify the fracture geometry results after processing in step S3 and import them into the production numerical simulator. Use the Oda method to set the equivalent fracture permeability. Based on the circuit series-parallel idea, the conductivity of the complex fracture network is calculated by the hydropower similarity principle. The calculation formula is as follows:

[0025]

[0026] R FN is the overall seepage resistance, R 1 —R 10 are the seepage resistance of each fracture;

[0027] The basic principle of Oda's method is as follows:

[0028]

[0029]

[0030] Where: R is the crack stiffness ratio, dimensionless; K n is the crack normal stiffness, MPa; K s is the crack shear stiffness, MPa; σ n is the normal stress on the crack surface, MPa; r is the crack size; JRC is the roughness, dimensionless; JCS is the compressive strength, MPa; is the friction angle, °;

[0031] The calculation formula of fracture surface permeability under different openings is as follows:

[0032]

[0033] Where: K Ω is the permeability of the polygon with area Ω in the fracture sheet, 10 -3 μm 2 ; K 0 is the permeability of the crack root, 10 -3 μm 2 ;w 0 is the opening of the crack root, m; w(x) is the opening of the crack at point x, m; x is the internal node coordinate of the crack, m; Ω is the crack area, m 2 ;

[0034] Based on the fracture network expansion model, numerical simulation of production was carried out considering the difference in conductivity between main fractures and branch fractures, and the daily production and cumulative production of the fractured wells under different fracture morphology conditions were obtained;

[0035] S5: Taking the optimal production as the constraint objective function, the genetic algorithm is used to find the optimal fracture conductivity under different reservoir conditions, and the conductivity of the differentiated reservoir with the optimal production as the target is obtained;

[0036] S6: Use API standard conductivity test equipment to test the fracture conductivity, and analyze the change rules of the fracture conductivity of the main fracture and branch fracture under different particle sizes, different sand concentrations, and different closing stress conditions. The permeability of the propped fracture can be calculated according to the following formula:

[0037]

[0038] The fracture conductivity of the proppant filling layer is expressed as:

[0039]

[0040] Where: b f is the width of the filled crack, i.e. the simulated crack width, cm; K g is the gas permeability of the proppant filling layer, μm 2 ;Q 0 is the flow rate, kPa; p 0 is the atmospheric pressure under standard conditions, kPa; μ is the test gas viscosity, mPa·s; P 1 is the upstream pressure, kPa; P 2 is the downstream pressure, kPa;

[0041] S7: Using the test rules of S6, according to the optimal conductivity ratio K of the main seam and branch seam opt The fracturing sand-paving scheme under mixed particle size conditions is designed to further obtain the minimum sand-paving concentration and minimum particle size combination that meets the optimal conductivity conditions, thereby completing the sand-paving design work of an efficient fracture network with optimal economy.

[0042] The beneficial effects of the present invention are:

[0043] The present invention utilizes the Karhunen-Loeve orthogonal series method to perform three-dimensional inversion of geostress, which saves costs to a certain extent and improves the calculation efficiency and accuracy of three-dimensional geostress.

[0044] The present invention can set up a conductivity test plan for the entire fracture network according to the actual fracture morphology through three-dimensional fracture extension numerical simulation, derive the conductivity equation of the complex fracture network, and effectively combine numerical simulation and experimental methods to achieve the optimization of the conductivity of the fracture network system.

[0045] The present invention takes into account the differences in the flow conductivity of the main and branch seams, and on this basis, conducts numerical simulation of production, forms differentiated optimization of the main and branch seams, and designs the minimum sand concentration, achieving dual optimization of production and economy. Taking the complex seam network of volume fracturing in unconventional oil and gas reservoirs as the research object, the optimization design of the seam network sand flow conductivity is carried out to improve the flow conductivity of the complex seam network and ensure the contribution of the complex seam network to production. The invention has practical application value for the efficient fracturing design of unconventional reservoirs. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] Figure 1 A technical flow chart of a sand paving optimization design method for a multi-scale complex fracture network;

[0047] Figure 2 It is a schematic diagram of the vertical geostress field of the reservoir three-dimensional geomechanical calculation results obtained by numerical simulation;

[0048] Figure 3 It is a schematic diagram of the maximum horizontal geostress field of the reservoir three-dimensional geomechanical calculation results obtained by numerical simulation;

[0049] Figure 4 This is a schematic diagram of the minimum horizontal geostress field of the reservoir three-dimensional geomechanical calculation results obtained by numerical simulation;

[0050] Figure 5 Schematic diagram of the simulation results of complex fracture network crack extension under different complexity conditions; (a) is a simple fracture, (b) is a complex multiple fracture, and (c) is a complex fracture network;

[0051] Figure 6 This is the main seam conductivity test analysis result diagram;

[0052] Figure 7 This is the result of the secondary fracture conductivity test analysis;

[0053] Figure 8 The calculation result of the critical value of the optimal main seam conductivity is shown in Figure 2.

[0054] Fig. 9 Schematic diagram of production comparison results under different main fracture and secondary fracture conductivity ratio conditions. DETAILED DESCRIPTION

[0055] The present invention will be further described below by way of embodiments in conjunction with the accompanying drawings, but is not limited thereto.

[0056] Embodiment 1:

[0057] A multi-scale complex crack network sand paving optimization design method includes the following steps:

[0058] S1: Constructing a one-dimensional geomechanical model of the reservoir: First, the corresponding rock mechanical parameters are obtained through the results of indoor core test experiments, combined with the logging interpretation results of acoustic wave time difference, gamma value, density, seismic inversion data, uniaxial and triaxial compression test results, that is, the one-dimensional geomechanical model is obtained;

[0059] The construction of a three-dimensional geomechanical model is carried out by combining the three-dimensional seismic data volume or experimental data points as the benchmark for constraints: the three-dimensional modeling of reservoir geomechanical parameters is carried out by using the sequential-Gaussian interpolation method in combination with the sedimentary phase model of the research block; the three-dimensional modeling of a series of strength parameters such as tensile strength and fracture toughness is completed by fitting formula calculation; the elastic model and Poisson's ratio are calculated as follows:

[0060]

[0061]

[0062] Where: E d is the dynamic elastic modulus of rock, GPa; μ d is the dynamic Poisson's ratio of rock; v p 、v s are the velocity of transverse and longitudinal waves, respectively, km / s; ρ is the density of rock, g / cm 3 ;

[0063] S2: Use the Karhunen-Loeve orthogonal series expansion method to invert the random three-dimensional geomechanical model, and use the simulated annealing method to optimize the inversion results of the ground stress, obtain the spatial distribution characteristics of the three-dimensional stress field, and check and correct the three-dimensional geomechanical model obtained in S1; the use of simulated annealing algorithm optimization can better avoid the inversion results falling into a local range. The specific steps are as follows:

[0064] (a) Determine the range of variation of model parameters and set an initial model m 0 , calculate the square value of the corresponding travel time residual's second norm, that is, the error equation E(m 0 );

[0065] (b) Perturb each parameter of the current model to generate a new model m 1 , calculate the corresponding error equation E(m 1 ), and calculate △E=E(m 1 )-E(m 0 );

[0066] (c) If △E < 0, the new model is accepted and assigned to the old model, that is, m 0 =m 1 ; If △E ≥ 0, then the new model m 1Assign values ​​according to probability P = exp(-△E / T), where exp is the natural exponent and T is the temperature. A number R is randomly generated between 0 and 1. When P(△E)>R, the modification is accepted, i.e. m 0 =m 1 ;

[0067] (d) repeating steps (b) and (c) at temperature T;

[0068] (e) slowly lowering the temperature T;

[0069] (f) Repeat steps (b)-(e) until △E satisfies the specified convergence condition.

[0070] S3: Based on the modified three-dimensional geomechanical model, the commercial three-dimensional fracture extension software (FracMan) is used to carry out the numerical simulation of three-dimensional fracture extension of hydraulic fracturing. The Slab grid method is used to calculate the overall transformation volume of the fracture. This equation is mainly used to estimate the equivalent volume opening of the natural fracture containing fracturing fluid. The specific calculation formula is shown in (3):

[0071]

[0072] Where: e is the crack opening, P pore is the fracture pore pressure; σ is the fracture surface normal stress; E is Young's modulus; ν is Poisson's ratio; r is the flow distance; r max is the maximum flow distance;

[0073] The actual geometric morphology of the complex fracture network after volume fracturing was further analyzed, and the relevant information of the main fracture and branch fractures, including the fracture length, width and height, was statistically obtained.

[0074] S4: Classify the fracture geometry results after processing in step S3 and import them into the production numerical simulator (EDFM or GOHFER). Use the Oda method to set the equivalent fracture permeability. Based on the circuit series-parallel idea, the conductivity of the complex fracture network is calculated by the hydropower similarity principle. The calculation formula is as follows:

[0075]

[0076] R FN is the overall seepage resistance, R 1 —R 10 are the seepage resistance of each fracture;

[0077] The basic principle of Oda's method is as follows:

[0078]

[0079]

[0080] Where: R is the crack stiffness ratio, dimensionless; K n is the crack normal stiffness, MPa; K s is the crack shear stiffness, MPa; σ n is the normal stress on the crack surface, MPa; r is the crack size; JRC is the roughness, dimensionless; JCS is the compressive strength, MPa; is the friction angle, °;

[0081] The calculation formula of fracture surface permeability under different openings is as follows:

[0082]

[0083] Where: K Ω is the permeability of the polygon with area Ω in the fracture sheet, 10 -3 μm 2 ; K 0 is the permeability of the crack root, 10 -3 μm 2 ;w 0 is the opening of the crack root, m; w(x) is the opening of the crack at point x, m; x is the internal node coordinate of the crack, m; Ω is the crack area, m 2 ;

[0084] Based on the fracture network expansion model, numerical simulation of production was carried out considering the difference in conductivity between main fractures and branch fractures, and the daily production and cumulative production of the fractured wells under different fracture morphology conditions were obtained;

[0085] S5: Taking the optimal production as the constraint objective function, the genetic algorithm is used to find the optimal fracture conductivity under different reservoir conditions, and the conductivity of the differentiated reservoir with the optimal production as the target is obtained;

[0086] S6: Use API standard conductivity test equipment to test the fracture conductivity, and analyze the change rules of the fracture conductivity of the main fracture and branch fracture under different particle sizes, different sand concentrations, and different closing stress conditions. The permeability of the propped fracture can be calculated according to the following formula:

[0087]

[0088] The fracture conductivity of the proppant filling layer is expressed as:

[0089]

[0090] Where: b f is the width of the filled crack, i.e. the simulated crack width, cm; K g is the gas permeability of the proppant filling layer, μm2 ;Q 0 is the flow rate, kPa; p 0 is the atmospheric pressure under standard conditions, kPa; μ is the test gas viscosity, mPa·s; P 1 is the upstream pressure, kPa; P 2 is the downstream pressure, kPa;

[0091] S7: Using the test rules of S6, according to the optimal conductivity ratio K of the main seam and branch seam opt The fracturing sand-paving scheme under mixed particle size conditions is designed to further obtain the minimum sand-paving concentration and minimum particle size combination that meets the optimal conductivity conditions, thereby completing the sand-paving design work of an efficient fracture network with optimal economy.

Claims

1. A multi-scale complex crack network sand paving optimization design method, characterized in that: The following steps are involved: S1: Constructing a one-dimensional geomechanical model of the reservoir: First, the corresponding rock mechanical parameters are obtained through the results of indoor core test experiments, combined with the logging interpretation results of acoustic wave time difference, gamma value, density, seismic inversion data, uniaxial and triaxial compression test results, that is, the one-dimensional geomechanical model is obtained; The construction of a three-dimensional geomechanical model is carried out by combining the three-dimensional seismic data volume or experimental data points as the benchmark for constraint: the three-dimensional modeling of reservoir geomechanical parameters is carried out by using the sequential-Gaussian interpolation method in combination with the sedimentary phase model of the research block; the three-dimensional modeling of strength parameters is completed by using the fitting formula calculation; the elastic model and Poisson's ratio are calculated as follows: Where: E d is the dynamic elastic modulus of rock, GPa; μ d is the dynamic Poisson's ratio of rock; v p 、v s are the velocity of transverse and longitudinal waves, respectively, km / s; ρ is the density of rock, g / cm 3 ; S2: The random three-dimensional geomechanical model is inverted using the Karhunen-Loeve orthogonal series expansion method, and the inversion results of the ground stress are optimized using the simulated annealing method to obtain the spatial distribution characteristics of the three-dimensional shale stress field, and the three-dimensional geomechanical model obtained in S1 is checked and corrected; S3: Based on the modified 3D geomechanical model, the commercial 3D fracture extension software is used to carry out the numerical simulation of hydraulic fracturing 3D fracture extension. The Slab grid method is used to calculate the overall fracture transformation volume. This equation is used to estimate the equivalent volume opening of the natural fracture containing fracturing fluid. The specific calculation formula is shown in (3): Where: e is the crack opening, P pore is the fracture pore pressure; σ is the fracture surface normal stress; E is Young's modulus; ν is Poisson's ratio; r is the flow distance; r max is the maximum flow distance; The actual geometric morphology of the complex fracture network after volume fracturing was further analyzed, and the relevant information of the main fracture and branch fractures, including the fracture length, width and height, was statistically obtained. S4: Classify the fracture geometry results after processing in step S3 and import them into the production numerical simulator. Use the Oda method to set the equivalent fracture permeability. Based on the circuit series-parallel concept, calculate the conductivity of the complex fracture network. The calculation formula is shown as follows: R FN is the overall seepage resistance, R1—R 10 are the seepage resistance of each fracture; The principle of Oda's method is as follows: Where: R is the crack stiffness ratio, dimensionless; K n is the crack normal stiffness, MPa; K s is the crack shear stiffness, MPa; σ n is the normal stress on the crack surface, MPa; JRC is the roughness, dimensionless; JCS is the compressive strength, MPa; is the friction angle, °; The calculation formula of fracture surface permeability under different openings is as follows: Where: K Ω is the permeability of the polygon with area Ω in the fracture sheet, 10 -3 μm 2 ; K0 is the permeability of the crack root, 10 -3 μm 2 ; w0 is the opening of the crack root, m; w(x) is the opening of the crack at point x, m; x is the internal node coordinate of the crack, m; Ω is the crack area, m 2 ; Based on the fracture network expansion model, numerical simulation of production was carried out considering the difference in conductivity between main fractures and branch fractures, and the daily production and cumulative production of the fractured wells under different fracture morphology conditions were obtained; S5: Taking the optimal production as the constraint objective function, the genetic algorithm is used to find the optimal fracture conductivity under different reservoir conditions, and the conductivity of the differentiated reservoir with the optimal production as the target is obtained; S6: Use API standard conductivity test equipment to test the fracture conductivity, and analyze the change rules of the fracture conductivity of the main fracture and branch fracture under different particle sizes, different sand concentrations, and different closing stress conditions; S7: Using the test rules of S6, according to the optimal conductivity ratio K of the main seam and branch seam opt The fracturing sand-paving scheme under mixed particle size conditions is designed to further obtain the minimum sand-paving concentration and minimum particle size combination that meets the optimal conductivity conditions, thereby completing the sand-paving design work of the efficient fracture network.

2. The multi-scale complex fracture network sand paving optimization design method according to claim 1 is characterized in that: In step S2, the specific steps of optimizing using the simulated annealing algorithm are as follows: (a) Determine the range of variation of model parameters, set an initial model m0, and calculate the square value of the corresponding second norm of the travel time residual, that is, the error equation E(m0); (b) Perturb each parameter of the current model to generate a new model m1, calculate the corresponding error equation E(m1), and calculate △E=E(m1)-E(m0); (c) If △E<0, the new model is accepted and assigned to the old model, that is, m0=m1; if △E≥0, the new model m1 is assigned with probability P=exp(-△E / T), where exp is the natural exponent and T is the temperature; a number R is randomly generated between 0 and 1, and when P(△E)>R, the modification is accepted, that is, m0=m1; (d) repeating steps (b) and (c) at temperature T; (e) slowly lowering the temperature T; (f) Repeat steps (b)-(e) until △E satisfies the specified convergence condition.

3. The multi-scale complex fracture network sand paving optimization design method according to claim 1 is characterized in that: In step S6, the permeability of the propped fracture is calculated according to the following formula: The fracture conductivity of the proppant filling layer is expressed as: Where: b f is the width of the filled crack, i.e. the simulated crack width, cm; K g is the permeability of the supporting fracture, μm 2 ; Q0 is the flow rate, kPa; p0 is the atmospheric pressure under standard conditions, kPa; μ is the test gas viscosity, mPa·s; P1 is the upstream pressure, kPa; P2 is the downstream pressure, kPa.

Images