Bidirectional coupling numerical simulation method of fracture propagation and Lagrangian proppant migration

Abstract

The present invention discloses a bidirectional coupling numerical simulation method of crack extension and Lagrangian proppant migration, comprising the following steps: S1: defining a numerical simulation region and inputting basic parameters; S2: establishing a fracturing fluid flow control equation and a hydraulic fracture deformation control equation, and performing coupling solution to obtain flow field parameters and crack width; S3: judging whether the crack is extended; S4: placing a proppant particle package and establishing a proppant migration control equation; S5: calculating particle stress terms, drag force terms, pressure gradient terms, net buoyancy terms, and collision damping terms; S6: solving the proppant migration control equation and updating the dynamic properties of the particle package; S7: calculating fluid-particle interaction terms and updating flow field parameters; S8: repeating steps S2-S7 until a preset simulation time is reached. The present invention can realize bidirectional coupling of crack extension and proppant migration, and can more accurately realize numerical simulation of crack (deformation) extension and proppant migration at an actual engineering scale.

Description

Technical Field

[0001] The invention relates to the technical field of hydraulic fracturing in oil and gas field development engineering, and in particular to a bidirectional coupling numerical simulation method for crack extension and Lagrangian proppant migration. Background Art

[0002] Hydraulic fracturing is the core technology for developing unconventional oil and gas resources such as shale oil and gas reservoirs and tight sandstone oil and gas reservoirs. The key to hydraulic fracturing is to form propped fractures with certain fracture conductivity, thereby providing a "high-speed channel" for oil and gas to flow from the reservoir into the wellbore. In hydraulic fracturing operations, the slick water fracturing fluid system has been widely used in the transformation process of unconventional oil and gas reservoirs due to its low reservoir damage, low operating cost, and strong fracture-making ability. However, the slick water fracturing fluid has poor sand-carrying performance due to its low viscosity, which often causes insufficient support for hydraulic fractures, restricts the effect of hydraulic fracturing, and becomes a bottleneck problem in the development of unconventional oil and gas. Therefore, studying the migration behavior and laying characteristics of proppants in fractures during slick water fracturing can provide important theoretical support and technical guidance for the optimization of hydraulic fracturing construction parameters.

[0003] At present, most of the numerical simulation studies on proppant migration assume that the fracture wall is fixed (static fracture), which is quite different from the fracture extension process in actual hydraulic fracturing. Hydraulic fracturing involves the coupling between multiple fields such as fluid velocity field, fluid pressure field and solid stress field. Discrete element method, finite element method, displacement discontinuity method and variants and combinations of the above three methods are the most common numerical methods for simulating fracture (deformation) extension. Compared with discrete element method and finite element method, displacement discontinuity method has the significant advantages of less number of grids, simple data preparation and faster solution speed. More importantly, displacement discontinuity method directly solves the unknown quantity with fracture aperture, which provides convenience for solving proppant migration in fractures. In addition, multiphase particle grid method has been successfully applied to the numerical simulation of proppant migration at engineering scale (fractures of hundreds of meters in length), and has great application prospects. Therefore, the combination of displacement discontinuity method and multiphase particle grid method is an important research direction for the study of fracture extension and proppant migration.

[0004] However, most of the existing coupling models based on the displacement discontinuity method use the concentration transport model to describe the proppant migration process in the fracture. The concentration transport model is an Euler-Euler method, which will destroy the "discrete characteristics" of the particle phase. Secondly, the concentration transport model is suitable for the sand-carrying situation of high-viscosity fracturing fluid, but not for the low-viscosity fracturing fluid (slick water) working condition. Thirdly, the concentration transport model needs to use an empirical formula to fit the experimental results, and the applicability and reliability of the relevant empirical formulas remain to be verified.

[0005] In addition, although the prior art realizes the one-way coupling of proppant migration and PKN type fracture extension based on the multiphase particle grid method, the one-way coupling model does not actually strictly realize the coupling between the fluid flow field and the solid stress field, that is, the two-way coupling, and has certain limitations. In addition, the one-way coupling model is difficult to realize the numerical simulation of hydraulic fracture extension and proppant migration in reservoirs with geostress differences. Summary of the invention

[0006] In view of the above problems, the present invention aims to provide a bidirectional coupling numerical simulation method of fracture propagation and Lagrangian proppant migration, which realizes efficient bidirectional coupling of fracture propagation and proppant migration based on displacement discontinuity method and multiphase particle grid method, and can realize engineering-scale fracture propagation and Lagrangian proppant migration simulation, providing a numerical tool for optimization of hydraulic fracturing pumping program.

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

[0008] A bidirectional coupling numerical simulation method for fracture propagation and Lagrangian proppant migration includes the following steps:

[0009] S1: Define the numerical simulation area and input basic parameters;

[0010] S2: establishing a fracturing fluid flow control equation and a hydraulic fracture deformation control equation of a plane three-dimensional displacement discontinuity method, and coupling and solving the fracturing fluid flow control equation and the hydraulic fracture deformation control equation to obtain flow field parameters and fracture width;

[0011] S3: Determine whether the crack is expanding: if it is expanding, add a new displacement discontinuity unit adjacent to the existing crack tip unit; if it is not expanding, keep the crack length unchanged;

[0012] S4: placing proppant particle packages into the numerical simulation area at a predetermined time node and establishing a proppant migration control equation;

[0013] S5: calculating the particle stress term, the drag force term, the pressure gradient term, the net buoyancy term and the collision damping term on the proppant particle package;

[0014] S6: solving the proppant migration control equation and updating the dynamic properties of the proppant particle package;

[0015] S7: Calculate the fluid-particle interaction term and update the flow field parameters of step S2 according to the calculation results;

[0016] S8: Repeat steps S2-S7 until the preset simulation time is reached.

[0017] Preferably, in step S2, the fracturing fluid flow control equation includes:

[0018] The mass conservation equation:

[0019]

[0020] Where: w is the crack width, m; θ f is the fluid volume fraction, m 3 / m 3 ; t is time, s; is the divergence operator; u f is the fracturing fluid flow velocity vector, m / s; C L is the fracturing fluid loss coefficient, m / s 0.5 ; τ(x) is the time when the fracturing fluid starts to lose at a certain position x, s;

[0021] Momentum conservation equation:

[0022]

[0023] Where: μ f is the fracturing fluid viscosity, Pa·s; is the gradient operator; p is the fracturing fluid pressure, Pa; F fp is the fluid-particle interaction term, Pa / m.

[0024] Preferably, in step S2, the hydraulic fracture deformation control equation of the plane three-dimensional displacement discontinuity method is:

[0025]

[0026] Where: is the normal stress of the plane three-dimensional displacement discontinuous unit i, Pa; N is the total number of plane three-dimensional displacement discontinuous units in the target system; is the influence coefficient of the normal displacement discontinuity of the plane three-dimensional displacement discontinuity element j on the normal stress of the plane three-dimensional displacement discontinuity element i, Pa / m; is the normal displacement discontinuity of the plane three-dimensional displacement discontinuity element j, m; is the in-situ stress.

[0027] Preferably, in step S2, a Newton-Raphson iterative method is used for coupling solution.

[0028] Preferably, in step S3, when judging whether the crack is expanding, the judgment is made based on the Type I stress intensity factor of the crack tip unit; if the Type I stress intensity factor of the crack tip unit is greater than or equal to the fracture toughness of the rock, the crack expands; if the Type I stress intensity factor of the crack tip unit is less than the fracture toughness of the rock, the crack does not expand.

[0029] Preferably, the type I stress intensity factor of the seam tip unit is calculated by the following formula:

[0030]

[0031] Where: K Ⅰ is the type I stress intensity factor, Pa·m 1 / 2 ; E is Young's modulus, Pa; v is Poisson's ratio; d is the crack tip unit length, m; D n is the normal displacement discontinuity of the crack tip unit, m.

[0032] Preferably, in step S4, the proppant migration control equation is:

[0033]

[0034] Where: a p is the acceleration vector of the proppant particle package, m / s 2 ;u p is the velocity vector of the proppant particle package, m / s; D p is the drag force coefficient, 1 / s; ρ p is the proppant particle density, kg / m 3 ρ f is the density of fracturing fluid, kg / m 3 ; g is the acceleration due to gravity, m / s 2 θ p is the particle volume fraction, m 3 / m 3 ; τ p is the particle stress, Pa; p > is the average velocity vector of all proppant particle packages around a certain proppant particle package, m / s; τ D is the collision relaxation time, s.

[0035] Preferably, in step S6, the dynamic properties of the proppant particle package are updated based on the forward Euler method, the leapfrog method or the prediction-correction method.

[0036] Preferably, in step S7, the fluid-particle interaction term is calculated by the following formula:

[0037]

[0038] ​Where: f is the probability distribution function of the particle phase; V p is the volume of the proppant particle package, m 3 .

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

[0040] Firstly, through the fracturing fluid flow control equation, the present invention can realize the bidirectional coupling of hydraulic fracturing crack deformation (extension) and proppant migration; secondly, the present invention retains the discrete characteristics (Lagrangian characteristics) of the particles, and can flexibly and conveniently realize more complex proppant particle injection (variable density, variable particle size proppant particle injection); finally, through the introduction of particle packing and particle stress, the present invention can realize the numerical simulation of crack (deformation) expansion and proppant migration at an actual engineering scale. BRIEF DESCRIPTION OF THE DRAWINGS

[0041] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0042] Figure 1 It is a schematic diagram of the calculation flow of the bidirectional coupling numerical simulation method of fracture extension and Lagrangian proppant migration of the present invention;

[0043] Figure 2 A schematic diagram of particle volume fraction distribution results at the end of simulation in a specific embodiment;

[0044] Figure 3 It is a schematic diagram of the particle package position distribution results at the end of the simulation in a specific embodiment. DETAILED DESCRIPTION

[0045] The present invention is further described below in conjunction with the accompanying drawings and embodiments. It should be noted that, in the absence of conflict, the embodiments in this application and the technical features in the embodiments can be combined with each other. It should be noted that, unless otherwise specified, all technical and scientific terms used in this application have the same meanings as those generally understood by those of ordinary skill in the art to which this application belongs. The words "including" or "comprising" and the like used in the disclosure of the present invention mean that the elements or objects appearing before the word cover the elements or objects listed after the word and their equivalents, without excluding other elements or objects.

[0046] like Figure 1 As shown, the present invention provides a bidirectional coupling numerical simulation method for fracture propagation and Lagrangian proppant migration, comprising the following steps:

[0047] S1: Define the numerical simulation area and input basic parameters.

[0048] In a specific embodiment, the basic parameters include displacement discontinuity unit size, fracturing fluid properties, proppant properties, rock properties and other parameters.

[0049] S2: Establishing a fracturing fluid flow control equation and a hydraulic fracture deformation control equation of a plane three-dimensional displacement discontinuity method, and coupling and solving the fracturing fluid flow control equation and the hydraulic fracture deformation control equation to obtain flow field parameters and fracture width.

[0050] In a specific embodiment, the fracturing fluid flow control equation includes:

[0051] The mass conservation equation:

[0052]

[0053] Where: w is the crack width, m; θ f is the fluid volume fraction, m 3 / m 3 ; t is time, s; is the divergence operator; u f is the fracturing fluid flow velocity vector, m / s; C L is the fracturing fluid loss coefficient, m / s 0.5 ; τ(x) is the time when the fracturing fluid starts to lose at a certain position x, s;

[0054] Momentum conservation equation:

[0055]

[0056] Where: μ f is the fracturing fluid viscosity, Pa·s; is the gradient operator; p is the fracturing fluid pressure, Pa; F fp is the fluid-particle interaction term, Pa / m.

[0057] From the above examples, it can be seen that the fluid-particle interaction term F fp The flow of the fracturing fluid is affected. The present invention improves the flow equation so that the present invention can achieve the bidirectional coupling between the hydraulic fracturing crack deformation (extension) and the proppant migration.

[0058] In a specific embodiment, the hydraulic fracture deformation control equation of the plane three-dimensional displacement discontinuity method is:

[0059]

[0060] Where: is the normal stress of the plane three-dimensional displacement discontinuous unit i, Pa; N is the total number of plane three-dimensional displacement discontinuous units in the target system; is the influence coefficient of the normal displacement discontinuity of the plane three-dimensional displacement discontinuity element j on the normal stress of the plane three-dimensional displacement discontinuity element i, Pa / m; is the normal displacement discontinuity of the plane three-dimensional displacement discontinuity element j, m; is the in-situ stress.

[0061] It should be noted that the hydraulic fracture deformation control equation of the plane three-dimensional displacement discontinuity method is the prior art, in which the normal stress influence coefficient The calculation method is an existing technology, and the specific mathematical expression is not repeated. For an open hydraulic fracture, the normal stress of the plane three-dimensional displacement discontinuous unit is equal to the sum of the fluid pressure inside the hydraulic fracture and the in-situ ground stress, and the tangential stress is equal to zero.

[0062] In a specific embodiment, the Newton-Raphson iterative method is used for coupled solution. It should be noted that the coupled solution of equations is a prior art. In addition to the coupled solution method used in this embodiment, other existing coupled solution methods can also be applied to the present invention. The flow field information obtained by the coupled solution includes flow velocity and pressure, etc.

[0063] S3: Determine whether the crack is expanding: if it is expanding, add a new displacement discontinuity unit adjacent to the existing crack tip unit; if it is not expanding, keep the crack length unchanged.

[0064] The expansion of hydraulic fractures follows linear elastic fracture mechanics. In a specific embodiment, when judging whether a fracture has expanded, the judgment is made based on the Type I stress intensity factor of the fracture tip unit; if the Type I stress intensity factor of the fracture tip unit is greater than or equal to the fracture toughness of the rock, the fracture expands; if the Type I stress intensity factor of the fracture tip unit is less than the fracture toughness of the rock, the fracture does not expand.

[0065] In a specific embodiment, the type I stress intensity factor of the seam tip unit is calculated by the following formula:

[0066]

[0067] Where: K Ⅰ is the type I stress intensity factor, Pa·m 1 / 2 ; E is Young's modulus, Pa; v is Poisson's ratio; d is the crack tip unit length, m; D n is the normal displacement discontinuity of the crack tip unit, m.

[0068] It should be noted that judging whether a crack has expanded is a prior art. In addition to the judgment method in the above embodiment, other existing methods for judging whether a crack has expanded can also be applied to the present invention. In addition, the calculation method of the type I stress intensity factor of the crack tip unit is also a prior art. In addition to the calculation method in the above embodiment, other methods for obtaining this parameter can also be applied to the present invention.

[0069] S4: placing proppant particle packages into the numerical simulation area at a predetermined time node and establishing a proppant migration control equation.

[0070] The multiphase particle grid method does not track the motion trajectory of a single physical particle, but describes the motion trajectory of a "parcel". A parcel is a collection of a certain number of physical particles with the same physical properties (size, density, acceleration, velocity and position, etc.). The present invention can achieve efficient bidirectional coupling between proppant migration and fracture propagation based on the displacement discontinuity method and the multiphase particle grid method, and can realize engineering-scale fracture propagation and Lagrangian proppant migration simulation.

[0071] In a specific embodiment, the proppant migration control equation is:

[0072]

[0073] Where: a p is the acceleration vector of the proppant particle package, m / s 2 ;u p is the velocity vector of the proppant particle package, m / s; D p is the drag coefficient, 1 / s; ρ p is the proppant particle density, kg / m 3 ρ f is the density of fracturing fluid, kg / m 3 ; g is the acceleration due to gravity, m / s 2 θ p is the particle volume fraction, m 3 / m 3 ; τ p is the particle stress, Pa; p > is the average velocity vector of all proppant particle packages around a certain proppant particle package, m / s; τ D is the collision relaxation time, s.

[0074] In the above embodiment, the particle volume fraction and the fluid volume fraction in the fracturing fluid flow control equation satisfy the condition that the sum is 1.

[0075] S5: Calculate the particle stress term, drag force term, pressure gradient term, net buoyancy term and collision damping term on the proppant particle package. ​

[0076] It should be noted that the calculation methods of the particle stress term, the drag force term, the pressure gradient term, the net buoyancy term and the collision damping term are prior arts, and the specific expressions are not repeated here.

[0077] S6: Solve the proppant migration control equation and update the dynamic properties of the proppant particle package.

[0078] In a specific embodiment, the dynamic properties of the proppant particle package can be updated based on the forward Euler method, the leapfrog method or the prediction-correction method. It should be noted that the updating method of this embodiment is only a preferred updating method of the present invention, and other methods capable of updating parameters in the prior art can also be applied to the present invention.

[0079] S7: Calculate the fluid-particle interaction term and update the flow field parameters of step S2 according to the calculation results.

[0080] In a specific embodiment, the fluid-particle interaction term is calculated by the following formula:

[0081]

[0082] Where: f is the particle phase probability distribution function, which is a function of particle position, particle velocity, particle density, particle volume and time; V p is the volume of the proppant particle package, m 3 .

[0083] In the above embodiment, after the fluid-particle interaction term is calculated by formula (6), the flow field parameters of step S2 are updated according to the calculation results in the present invention, which means that when step S2 is repeated in step S8, the fluid-particle interaction term is used for the flow field update of this step.

[0084] S8: Repeat steps S2-S7 until the preset simulation time is reached.

[0085] In a specific embodiment, the bidirectional coupling numerical simulation method of fracture propagation and Lagrangian proppant migration of the present invention is used to simulate a variable-size proppant particle injection case, that is, two proppant particles of different sizes (type 1: diameter 0.3 mm; type 2: diameter 0.6 mm) are injected successively, specifically comprising the following steps:

[0086] (1) Define the numerical simulation area and input basic parameters;

[0087] In this embodiment, the basic parameters are shown in Table 1:

[0088] Table 1 Basic parameters

[0089] Fracturing fluid viscosity 1mPa·s Young's modulus of rock 35GPa Fracturing fluid density <![CDATA[1000kg / m 3 ]]> Rock Poisson's ratio 0.25 Proppant density <![CDATA[2650kg / m 3 ]]> Rock fracture toughness <![CDATA[1.5MPa·m 0.5 ]]> Proppant particle size 0.3,0.6mm Pre-liquid time 900s Proppant concentration 10％ Total simulation time 2000s Filtration coefficient <![CDATA[2.5×10 -6 m / s 0.5 ]]> Injection displacement <![CDATA[0.01m 3 / s]]>

[0090] (2) establishing the fracturing fluid flow control equations shown in equations (1)-(2) and the hydraulic fracture deformation control equations of the plane three-dimensional displacement discontinuity method shown in equation (3), and using the Newton-Raphson iterative method to couple and solve them to obtain the flow field parameters and the fracture width;

[0091] (3) Calculate the type I stress intensity factor of the crack tip unit according to formula (4) and determine whether the crack is expanding: if it is expanding, add a new displacement discontinuity unit at the adjacent position of the existing crack tip unit; if it is not expanding, keep the crack length unchanged;

[0092] (4) placing a proppant particle package into the numerical simulation area at a predetermined time node, and establishing the proppant migration control equation shown in equation (5);

[0093] (5) calculating the particle stress term, drag force term, pressure gradient term, net buoyancy term and collision damping term of the proppant particle package;

[0094] (6) solving the proppant migration control equation and updating the dynamic properties (including acceleration, velocity and position) of the proppant particle package using the forward Euler method;

[0095] In this embodiment, the mathematical expression of the forward Euler method is:

[0096]

[0097] Where: u n+1 is the proppant particle velocity at time step n+1, m / s; u n is the proppant particle package velocity at time step n, m / s; a n is the acceleration of the proppant particle package at time step n, m / s 2 ; Δt is the time step, s; x n+1 is the position of the proppant particle package at time step n+1, m; x n is the position of the proppant particle package at time step n, m.

[0098] (7) Calculate the fluid-particle interaction term according to formula (6), and update the flow field parameters of step S2 according to the calculation results;

[0099] (8) Repeat steps (2)-(7) until the preset simulation time is reached.

[0100] In this embodiment, the particle volume fraction distribution at the end of the numerical simulation is as follows: Figure 2 As shown in Figure 2, the hydraulic fracture width distribution is as follows: Figure 3 shown.

[0101] In summary, the present invention provides a bidirectional coupling numerical simulation method of fracture extension and Lagrangian proppant migration, which solves the fracture (deformation) extension and proppant migration processes based on the planar three-dimensional displacement discontinuity method and the multiphase particle grid method respectively; adopts the improved fracturing fluid flow control equation to describe the fracturing fluid flow in the fracture, and realizes the bidirectional coupling of fracture (deformation) extension and proppant migration; retains the discrete properties of proppant particles, and can flexibly and conveniently realize more complex proppant particle injection (variable density, variable particle size proppant particle injection); due to the introduction of particle packs and particle stress, it can realize the optimization of hydraulic fracturing construction pumping program on an engineering scale. Compared with the prior art, the present invention has significant progress.

[0102] The above description is only a preferred embodiment of the present invention and does not limit the present invention in any form. Although the present invention has been disclosed as a preferred embodiment as above, it is not used to limit the present invention. Any technician familiar with this profession can make some changes or modify the technical contents disclosed above into equivalent embodiments without departing from the scope of the technical solution of the present invention. However, any simple modification, equivalent change and modification made to the above embodiments according to the technical essence of the present invention without departing from the content of the technical solution of the present invention still fall within the scope of the technical solution of the present invention.

Claims

1. A bidirectional coupling numerical simulation method for fracture propagation and Lagrangian proppant migration, characterized in that: The following steps are involved: S1: Define the numerical simulation area and input basic parameters; S2: establishing a fracturing fluid flow control equation and a hydraulic fracture deformation control equation of a plane three-dimensional displacement discontinuity method, and coupling and solving the fracturing fluid flow control equation and the hydraulic fracture deformation control equation to obtain flow field parameters and fracture width; The fracturing fluid flow control equation includes: The mass conservation equation: Where: w is the crack width, m; θ f is the fluid volume fraction, m 3 / m 3 ; t is time, s; is the divergence operator; u f is the fracturing fluid flow velocity vector, m / s; C L is the fracturing fluid loss coefficient, m / s 0.5 ; τ(x) is the time when the fracturing fluid starts to lose at a certain position x, s; Momentum conservation equation: Where: μ f is the fracturing fluid viscosity, Pa·s; is the gradient operator; p is the fracturing fluid pressure, Pa; F fp is the fluid-particle interaction term, Pa / m; S3: Determine whether the crack is expanding: if it is expanding, add a new displacement discontinuity unit adjacent to the existing crack tip unit; if it is not expanding, keep the crack length unchanged; S4: placing proppant particle packages into the numerical simulation area at a predetermined time node and establishing a proppant migration control equation; S5: calculating the particle stress term, the drag force term, the pressure gradient term, the net buoyancy term and the collision damping term on the proppant particle package; S6: solving the proppant migration control equation and updating the dynamic properties of the proppant particle package; S7: Calculate the fluid-particle interaction term and update the flow field parameters of step S2 according to the calculation results; S8: Repeat steps S2-S7 until the preset simulation time is reached.

2. The bidirectional coupling numerical simulation method of fracture propagation and Lagrangian proppant migration according to claim 1, characterized in that: In step S2, the hydraulic fracture deformation control equation of the plane three-dimensional displacement discontinuity method is: Where: is the normal stress of the plane three-dimensional displacement discontinuous unit i, Pa; N is the total number of plane three-dimensional displacement discontinuous units in the target system; is the influence coefficient of the normal displacement discontinuity of the plane three-dimensional displacement discontinuity element j on the normal stress of the plane three-dimensional displacement discontinuity element i, Pa / m; is the normal displacement discontinuity of the plane three-dimensional displacement discontinuity element j, m; is the in-situ stress.

3. The bidirectional coupling numerical simulation method of fracture propagation and Lagrangian proppant migration according to claim 1, characterized in that: In step S2, a Newton-Raphson iterative method is used to perform a coupled solution.

4. The bidirectional coupling numerical simulation method of fracture propagation and Lagrangian proppant migration according to claim 1, characterized in that: In step S3, when judging whether the crack is expanding, the judgment is made according to the type I stress intensity factor of the crack tip unit; if the type I stress intensity factor of the crack tip unit is greater than or equal to the rock fracture toughness, the crack expands; if the type I stress intensity factor of the crack tip unit is less than the rock fracture toughness, the crack does not expand.

5. The bidirectional coupling numerical simulation method of fracture propagation and Lagrangian proppant migration according to claim 4, characterized in that: The type I stress intensity factor of the seam tip unit is calculated by the following formula: Where: K Ⅰ is the type I stress intensity factor, Pa·m 1 / 2 ; E is Young's modulus, Pa; v is Poisson's ratio; d is the crack tip unit length, m; D n is the normal displacement discontinuity of the crack tip unit, m.

6. The bidirectional coupling numerical simulation method of fracture propagation and Lagrangian proppant migration according to claim 1, characterized in that: In step S4, the proppant migration control equation is: Where: a p is the acceleration vector of the proppant particle package, m / s 2 ;u p is the velocity vector of the proppant particle package, m / s; t is the time, s; D p is the drag force coefficient, 1 / s; u f is the fracturing fluid flow velocity vector, m / s; ρ p is the proppant particle density, kg / m 3 ; ρ f is the density of fracturing fluid, kg / m 3 ; g is the acceleration due to gravity, m / s 2 ; is the gradient operator; p is the fracturing fluid pressure, Pa; θ p is the particle volume fraction, m 3 / m 3 ; τ p is the particle stress, Pa; p > is the average velocity vector of all proppant particle packages around a certain proppant particle package, m / s; τ D is the collision relaxation time, s.​ 7. The bidirectional coupling numerical simulation method of fracture propagation and Lagrangian proppant migration according to claim 1, characterized in that: In step S6, the dynamic properties of the proppant particle package are updated based on the forward Euler method, the leapfrog method or the prediction-correction method.

8. The bidirectional coupling numerical simulation method of fracture propagation and Lagrangian proppant migration according to claim 1, characterized in that: In step S7, the fluid-particle interaction term is calculated by the following formula: Where: F fp is the fluid-particle interaction term, Pa / m; f is the particle phase probability distribution function; V p is the volume of the proppant particle package, m 3 ; ρ p is the proppant particle density, kg / m 3 ; D p is the drag force coefficient, 1 / s; u f is the fracturing fluid flow velocity vector, m / s; u p is the proppant particle package velocity vector, m / s; is the gradient operator; p is the fracturing fluid pressure, Pa.

Images