A calculation method for long-term permeability damage of artificially filled fractures after hydraulic fracturing

The permeability damage of artificial filling fractures after hydraulic fracturing is simulated by deep filter grid model and finite difference method, which solves the simulation difficulties in the existing technology, and achieves efficient field application and production capacity prediction.

CN115952746BActive Publication Date: 2025-07-29QINGDAO INST OF MARINE GEOLOGY
View PDF 0 Cites 0 Cited by

Patent Information

Application Number
CN202211539916.5
Authority / Receiving Office
CN · China
Patent Type
Patents(China)
Current Assignee / Owner
Filing Date
2022-12-02
Publication Date
2025-07-29
Estimated Expiration
2042-12-02

AI Technical Summary

Technical Problem

The prior art is difficult to efficiently simulate the long-term permeability damage of artificially filled cracks after hydraulic fracturing at engineering scales, and cannot meet the needs of on-site fracturing design and numerical simulation of reservoirs.

Method used

The deep filter grid model is used to determine the capture coefficient by calculating the macroscopic fracture flow velocity distribution and particle trajectory, and the finite difference method is used to simulate particle blockage and calculate the fracture permeability damage.

Benefits of technology

It effectively simulates crack blocking behavior under long-term mining conditions under engineering scale, and is suitable for on-site fracturing design and production capacity prediction.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN115952746B_ABST
    Figure CN115952746B_ABST
Patent Text Reader

Abstract

The present invention discloses a calculation method for long-term permeability damage of artificial filled fractures after hydraulic fracturing. First, the flow velocity distribution law of macroscopic fractures is determined, then the capture coefficient is determined based on the particle trajectory model, and finally the calculation of fracture permeability damage caused by particle invasion is realized. This solution is solved based on the finite difference method and can efficiently simulate the fracture plugging behavior at the engineering scale. Compared with other technical solutions based on the discrete element method, this solution determines the particle capture efficiency of the unit based on the particle trajectory in the fracture pore space, avoiding the generation and trajectory simulation of solid particles. Therefore, it can simulate the fracture damage behavior under long-term mining conditions with less computing resources, is more suitable for on-site application, and has a wider practical application value.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] The present invention belongs to the technical fields of hydraulic fracturing in oil and gas field development and reservoir numerical simulation, and particularly relates to a calculation method for long-term permeability damage of an artificial fracture after hydraulic fracturing. Background Art

[0002] Since its birth in the 1930s of the last century, the hydraulic fracturing technology has developed rapidly. At present, it has become a key technology for realizing the efficient development of complex oil and gas resources such as low-permeability and tight reservoirs, and has been widely used in the "sand control - production increase integration" well completion design of reservoirs such as extra-high water cut and unconsolidated sandstone reservoirs. In the process of hydraulic fracturing design, the conductivity of the artificial fracture is a key parameter affecting the success or failure of hydraulic fracturing construction and is also the main factor determining the fracturing transformation effect. The designed conductivity of the fracture consists of the effective fracture width and the fracture permeability. Under the complex flow and stress of the reservoir after fracturing, both of them will deviate from the designed values, resulting in the actual fracture conductivity falling short of expectations. Among them, the damage of the fracture width is mainly caused by the crushing and embedding of proppants, and this part can be determined through long-term conductivity experiments in the laboratory; the damage of the fracture permeability is caused by the invasion and blockage of reservoir fine particles, which has the characteristics of large scale, long period, and many influencing factors, and it is difficult to study by means of laboratory experiments.

[0003] Existing research on the permeability damage induced by particle invasion in fractures is carried out based on the discrete element method. This method requires modeling each proppant and reservoir invasion particle, with large computational resource requirements, and the simulated geometric and time scales are quite different from the real fractures. It can only carry out qualitative mechanism research and cannot provide guidance for the field. Although the prediction of long-term fracture permeability damage is of great significance for the evaluation of reservoir transformation effect, the prediction of post-fracture productivity, the selection of wells and layers for refracturing, etc., there is currently no systematic calculation method for long-term fracture permeability damage that can support on-site fracturing design and reservoir numerical simulation.

[0004] The particle invasion in the artificial fracture belongs to a typical problem of the flow of dispersed solid particles in a porous medium, which widely appears in many engineering fields. For a solid-liquid two-phase flow system with a target suspended particle size greater than 1 μm, deep filtration is the most effective particle separation process. The classical deep filtration model is currently the most commonly used method for simulating the behavior of particles in the filter layer. This model regards the filter layer as a whole, determines the empirical parameters in the model through experimental data, and then predicts the overall filtration efficiency of the filter layer. The particle concentration and pressure distribution in the filter layer are important indicators for evaluating the performance of the filter layer in the model. Due to the advantages of simple form and easy understanding, the classical deep filtration model is widely used in a large number of engineering calculations. However, the classical deep filtration model ignores the particle hydrodynamics behavior in the pore throats, and due to the one-dimensional scale of this model, it cannot characterize the spreading law of the invading particles in the filter layer. Summary of the Invention

[0005] In view of the problem in the prior art that it is difficult to simulate and calculate the dynamic permeability of artificial fractures under long-term mining conditions at the engineering scale, the present invention proposes a method for calculating the long-term permeability damage of artificially filled fractures after hydraulic fracturing. This method can be applied to calculations under engineering scale and long-term mining conditions, meeting the requirements of on-site fracturing design and production capacity prediction.

[0006] The present invention is realized by the following technical solutions: A method for calculating the long-term permeability damage of artificially filled fractures after hydraulic fracturing, comprising the following steps:

[0007] Step A: Calculate the macroscopic fracture flow velocity distribution law to obtain the fluid flow distribution curve in the fracture;

[0008] Step B: Determine the capture coefficient based on the particle trajectory, including:

[0009] Step B1: Construct a deep filtration grid model: Discretize the artificial fracture into a set of interconnected percolation units. Each percolation unit includes an inlet and an outlet in both the horizontal and vertical directions, and determine the geometric structure characteristic parameters of the percolation unit, including the pore radius R p , the throat radius R t , the pore region length L p and the throat region length L t ;

[0010] Step B2: Calculate the particle capture coefficients in the vertical and horizontal directions:

[0011] (B21) Calculate the flow velocity distribution characteristics in each percolation unit, analyze the forces on the particles invading the percolation unit, and determine the velocity sum of the particles in the horizontal and vertical directions and their trajectories in the percolation unit;

[0012] (B22) According to the calculation results of (B21), determine the particle capture coefficients in the horizontal and vertical directions:

[0013]

[0014]

[0015] Among them, λ h is the particle capture coefficient in the horizontal direction, λ v is the particle capture coefficient in the vertical direction, Q in is the flow rate value injected into the percolation unit, φ is the porosity in the unit; V e is the volume of the percolation unit, l vd represents the displacement of the horizontally injected particles in the vertical direction, and l hd represents the displacement of the vertically injected particles in the horizontal direction;

[0016] Step C: Calculate the fracture permeability damage at a single time step;

[0017] Step C1: Determine the amount of particles invading linearly perpendicular to the fracture in each percolation unit, and calculate the total particle concentration entering the percolation unit;

[0018] Step C2: Solve the discrete deep filtration grid model, inject reservoir fluid into the fracture at each time step, and calculate the fracture permeability damage distribution caused by particle plugging;

[0019] Step D: After the end of each time step, process the plugged units, characterize the plugged percolation units as particle source terms; increment the time step by 1, and repeat Step C and Step D until the maximum calculation time is reached to achieve the calculation of permeability damage.

[0020] Further, in Step B2, the flow velocity distribution in the percolation unit is calculated as follows:

[0021]

[0022] In the formula: R0 is the inlet radius, for the pore control region R0 = R p , for the throat control region R0 = R t ; R p is the pore radius, R t is the throat radius, u0 is the fluid flow velocity at the inlet, R(x) is the flow channel radius corresponding to the position x, x is the horizontal position in the percolation unit, and u(r, x) is the flow velocity magnitude at the horizontal position x and flow channel radius r.

[0023] Further, Step C1 is specifically implemented as follows:

[0024] (1) Based on the fluid flow distribution curve in the fracture obtained in Step A, calculate the linear flow rate distribution perpendicular to the fracture under the condition of bi-linear flow;

[0025] (2) Determine the proportion of solid-phase particles in the target reservoir fluid to determine the amount of particle invasion in each percolation unit;

[0026] (3) Further determine the total amount of particles injected into the percolation unit:

[0027] R total = R pi + C h (x - 1, y, t)+ C v (x, y + 1, t)

[0028] In the formula, R piThe amount of particle intrusion entering the unit perpendicular to the crack direction, R total The total particle concentration entering the unit, C h The particle concentration at the horizontal inlet, C v The particle concentration at the vertical inlet.

[0029] Furthermore, the specific implementation of step C2 is as follows:

[0030] (1) For the percolation unit with coordinates (x, y), the horizontal and vertical outlets of the target percolation unit are respectively the inlets of the horizontally adjacent unit and the vertically adjacent unit. The particles in the unit conform to the following mass conservation:

[0031]

[0032] Among them, C v (x,y,t) is the particle concentration at the vertical inlet, C h (x,y,t) is the particle concentration at the horizontal inlet, σ e (x,y,t) is the variation law of the particle concentration with time;

[0033] Apply the finite difference method to discretize and solve each term in the above formula. The discretization method is shown in the following formula:

[0034]

[0035]

[0036]

[0037] (2) Set the initial and boundary conditions of the model, solve the mass conservation equation, and obtain the variation law σ e (x,y,t) of the particle concentration trapped in the percolation unit with time.

[0038] Furthermore, in step D, the blocked unit will be characterized as a particle source term, and its calculation formula is as follows

[0039] C h (x + 1,y,t) = λ rm D c u(x,y,t)

[0040] In the formula, C h (x + 1,y,t) is the particle concentration at the horizontal outlet, λ rm is the re-migration coefficient of the blocked unit, D c is the viscosity coefficient.

[0041] Compared with the prior art, the advantages and positive effects of the present invention are:

[0042] The dynamic permeability damage calculation method for artificial fractures proposed in this scheme can be solved based on the finite difference method, efficiently simulating the crack blockage behavior at an engineering scale. Compared with other technical solutions based on discrete elements, this scheme determines the particle capture efficiency of the unit based on the particle trajectory of the pore space within the fracture, avoiding the generation and trajectory simulation of physical particles. Therefore, fewer computing resources can be used to simulate the fracture damage behavior under long-term mining conditions, making it more suitable for field applications. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Figure 1 A grid decomposition diagram of a reservoir containing artificial fractures according to an embodiment of the present invention;

[0044] Figure 2 is a flow velocity distribution curve in an artificial fracture according to an embodiment of the present invention;

[0045] Figure 3 A schematic diagram of the artificial fracture physical model and its discrete mode according to an embodiment of the present invention;

[0046] Figure 4 This is a flow chart of model calculation according to an embodiment of the present invention;

[0047] Figure 5 This is a curve showing the change of dynamic permeability of artificial fractures with mining time in an embodiment of the present invention. DETAILED DESCRIPTION

[0048] In order to more clearly understand the above-mentioned objects, features and advantages of the present invention, the present invention is further described below with reference to the accompanying drawings and embodiments. In the following description, many specific details are set forth to facilitate a full understanding of the present invention. However, the present invention can also be implemented in other ways than those described herein. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0049] This embodiment discloses a method for predicting dynamic permeability damage of artificial fractures after hydraulic fracturing based on a deep filtration grid model. The method is applicable to artificial fractures formed after various reservoir transformations and can be used in hydraulic fracturing scheme design, repeated fracturing well and layer selection, and post-fracturing production capacity prediction. Figure 4 As shown, the following steps are included:

[0050] Step A: Calculate the macroscopic fracture velocity distribution law to obtain the fluid flow distribution curve in the fracture;

[0051] Step B, determining the capture coefficients in the vertical and horizontal directions based on the particle trajectory;

[0052] Step C: Calculate the total particle concentration entering the filtration unit and calculate the fracture permeability damage for a single time step.

[0053] Step D: Characterize the blocked filtration unit as a particle source term, and repeat Steps C and D until the maximum calculation time is reached to achieve the calculation of permeability damage.

[0054] For a clearer understanding of the solution of the present invention, the present invention will be described in detail below:

[0055] In Step A, when determining the flow velocity distribution law of the macroscopic artificial fracture, the flow velocity distribution in the artificial fracture is obtained through numerical simulation calculation, which can be achieved by using conventional and mature technical means. The specific steps are as follows:

[0056] (1) According to the on-site fracturing data, establish a post-fracture reservoir model with artificial fractures; the modeling parameters include reservoir porosity, permeability, fluid compressibility, fluid viscosity, fluid density, oil well radius, and the fracture geometric model includes fracture length, width, and height. In this embodiment, the fracture width is set to 1 cm, the fracture height is 20 m, and the fracture length is 60 m.

[0057] (2) Mesh the established artificial fracture reservoir model, and densify the mesh around the fracture to obtain the mesh system required for calculation. The numerical simulation mesh is as Figure 1 shown.

[0058] (3) Set the pressure and flow boundary conditions in the model, and use the finite element method to solve the flow velocity distribution in the fracture to obtain the fluid flow distribution curve in the fracture, as Figure 2 shown; in this embodiment, the oil well production is set to 2.5, 5, 7.5, and 10 tons per day respectively.

[0059] In Step B, the key to solving the deep filtration grid model lies in determining the particle capture coefficient in the filtration unit. In this step, the method for determining the capture coefficient based on the particle trajectory model includes the following steps:

[0060] Step B1: Discretize the artificial fracture into a set of interconnected filtration units. Each filtration unit has a horizontal inlet, a horizontal outlet, a vertical inlet, and a vertical outlet, as Figure 3 shown; the geometric structure of the flow-through area of the filtration unit conforms to the pore-throat interlaced characteristics; based on the proppant particle size parameter, determine the geometric characteristic parameters in the artificial fracture filtration unit, including: pore radius R p , throat radius R t , pore area length L p , throat area length L t ;

[0061] Step B2. Calculate the particle capture coefficients in the vertical and horizontal directions:

[0062] (1) Calculate the flow velocity distribution within the infiltration unit:

[0063]

[0064] In the formula: R0 is the inlet radius. For the pore control region, R0 = R p , and for the throat control region, R0 = R t ; u0 is the fluid flow velocity at the inlet, m / s; r is the distance from the target point to the center of the flow channel, m.

[0065] (2) Conduct a force analysis on the particles invading the infiltration unit to determine their velocity and trajectory. The particle capture coefficients in the horizontal and vertical directions are determined by the trajectory of the particles in the infiltration unit. For particles with a particle size range of 1 μm - 1 mm, their motion is dominated by the Stokes force and the gravity-buoyancy resultant force. The Stokes force F d is shown in the following formula:

[0066] F d = D c (v - u)(2)

[0067] In the formula, D c is the viscosity coefficient, D c = 3πμd p ; μ is the fluid viscosity, Pa·s; d p is the particle diameter, m; u is the fluid flow velocity, m / s; v is the particle velocity, m / s.

[0068] The gravity-buoyancy resultant force acting on the particles is shown in the following formula:

[0069]

[0070] In the formula, ρ p is the particle density, kg / m 3 ; ρ f is the fluid density, kg / m 3 ; F v is the gravity-buoyancy resultant force, N; F g is the gravity, N; F b is the buoyancy, N.

[0071] The direction of the gravity-buoyancy resultant force acting on the particles is the vertical direction; the direction of the Stokes force is opposite to the particle motion direction; based on the forces acting on the particles in formulas (2) and (3), calculate the velocities of the particles in the x and y directions according to Newton's law:

[0072]

[0073]

[0074] where m p is the particle mass, in g.

[0075] (3) Based on the particle velocity formulas in Eqs. (4) and (5), combined with the geometric model of the percolation unit, calculate the particle capture coefficients λ h and λ v in the horizontal and vertical directions. The specific process is as follows:

[0076] The particle migration in the horizontal direction conforms to the following equation:

[0077]

[0078] where t1 is the time required for the particle to pass through the entire percolation unit in the horizontal direction, in s; assuming other parameters in Eq. (6) are known, t1 can be calculated. During the time t1, the displacement l vd of the horizontally injected particle in the vertical direction is:

[0079]

[0080] Similarly, the movement of the particle injected from the vertical inlet into the unit conforms to the following law:

[0081]

[0082] where t2 is the time required for the particle to pass through the entire percolation unit in the vertical direction, in s; assuming other parameters in Eq. (8) are known, t2 can be calculated. During the time t2, the displacement l hd [[ID=3�]]of the vertically injected particle in the horizontal direction is:

[0083]

[0084] After determining the above parameters, according to the definition of the capture coefficient, the particle capture efficiency calculation formula can be expressed as follows:

[0085]

[0086] where: σ e is the concentration of the particles captured in the percolation unit; λ h is the particle capture coefficient in the horizontal direction; λ v is the particle capture coefficient in the vertical direction; Q in is the flow rate value injected into the percolation unit, in m 3 / s; φ is the porosity in the unit; V e is the volume of the percolation unit, in m 3 ;

[0087] Based on Equation (10), the expression for the particle capture coefficient can be obtained:

[0088]

[0089]

[0090] In step C for calculating the fracture permeability damage, assume that 1 PV of reservoir fluid is injected into the fracture at each time step. The specific steps are as follows:

[0091] Step C1: Calculate the amount of particles invading linearly perpendicular to the fracture in each percolation unit, and then determine the total particle concentration entering the percolation unit. The flow around the fracture includes linear flow in the fracture and linear flow perpendicular to the fracture direction:

[0092] (1) Based on the fluid distribution curve in the fracture obtained in step A, calculate the linear flow rate distribution perpendicular to the fracture under the condition of double linear flow:

[0093] q in =A(u i+1 -u i ) (13)

[0094] In the formula, u i+1 is the fluid velocity in the fracture at position i + 1, m / s; u i is the fluid velocity in the fracture at position i, m / s; A is the cross-sectional area of the fracture unit, m 3 .

[0095] (2) By means of laboratory experiments and other methods, determine the proportion of solid particles in the target reservoir fluid. Since the amount of particles invading with the double linear flow in the fracture in each percolation unit is proportional to the fluid flow rate entering the unit perpendicular to the fracture direction, the amount of particle invasion in each percolation unit can be determined, that is:

[0096] R pi ∝q in (14)

[0097] In the formula, q in is the linear flow rate perpendicular to the fracture direction, m 3 / s; R pi is the particle concentration entering the unit perpendicular to the fracture direction.

[0098] (3) The calculation formula for the total amount of particles injected into the percolation unit is as follows:

[0099] R total =R pi +C h (x - 1, y, t)+C v (x, y + 1, t) (15)

[0100] wherein, R total is the total particle concentration entering the unit.

[0101] Step C2: Solve the discrete deep filtration grid model, inject reservoir fluid into the fracture at each time step, and calculate the fracture permeability damage distribution caused by particle plugging;

[0102] Discretize the deep filtration grid model in the artificial fracture, and obtain the variation law of the particle concentration trapped in the percolation unit with time σ e (x, y, t), including the following steps:

[0103] (1) The physical model of the deep filtration grid model in the artificial fracture and the particle transfer relationship between units are as Figure 3 shown. For the percolation unit with coordinates (x, y), at time t, the particle concentration at the horizontal inlet is C h (x, y, t), and the particle concentration at the horizontal outlet is C h (x + 1, y, t); the particle concentration at the vertical inlet is C v (x, y, t), and the particle concentration at the vertical outlet is C v (x, y - 1, t); wherein, the horizontal and vertical outlets of the target percolation unit are respectively the inlets of the horizontally adjacent unit and the vertically adjacent unit; the particles in the unit conform to the law of conservation of mass:

[0104]

[0105] According to Figure 3 the particle transfer law between units in the geometric model in, C h (x, y, t), C v (x, y, t), and σ e (x, y, t) can be written in the following forms respectively:

[0106]

[0107]

[0108]

[0109] wherein: λ h is the particle capture coefficient in the horizontal direction; λ v is the particle capture coefficient in the horizontal direction, and C0 is the particle concentration at the inlet of the percolation unit at the fracture end.

[0110] Apply the finite difference method to discretize and solve each term in the particle mass conservation equation (Equation 16) in the unit, and the discretization method is as shown in the following formula:

[0111]

[0112]

[0113]

[0114] (2) The initial and boundary conditions of the mass conservation equation (Equation 16) are as shown in Equations 23 - 25:

[0115] σ e (x, y, 0) = 0 (23)

[0116] C h (0, y, t) = C0 (24)

[0117] C h (x, 0, t) = 0 (25)

[0118] Applying the above initial and boundary conditions, the mass conservation equation (Equation 16) can be solved to obtain the variation law of the particle concentration trapped in the percolation unit with time σ e (x, y, t).

[0119] Step D: After each time step ends, process the blocked units and characterize the blocked percolation units as particle source terms;

[0120] (1) Calculate the porosity and permeability within the unit based on the particle concentration trapped in the unit:

[0121]

[0122]

[0123] In the formula, φ0 is the initial porosity of the percolation unit; φ d is the porosity of the blocked particle accumulation; k is the permeability of the percolation unit; k0 is the initial permeability.

[0124] (2) During the reservoir exploitation process, particle blockage will occur in the fracture units, and the blocked units will be characterized as particle source terms, and its calculation formula is as shown in the following formula:

[0125] C h (x + 1, y, t) = λ rm D c u(x, y, t) (28)

[0126] In the formula, λ rm is the re - migration coefficient of the blocked unit, which is determined by experiments.

[0127] After the above steps are completed, increment the time step by 1 and repeat steps C and D until the maximum calculation time is reached.

[0128] The method of this solution can obtain the dynamic damage characteristics of the fracture permeability during the oil well production process after fracturing. Figure 5 Shown is the curve of the dynamic permeability of the artificial fracture obtained by applying this method with respect to the production time.

[0129] The above are only the preferred embodiments of the present invention and do not limit the present invention in other forms. Any person skilled in the art may use the disclosed technical content to make changes or modifications into equivalent embodiments with equivalent changes and apply them to other fields. However, any simple modification, equivalent change, and modification made to the above embodiments based on the technical essence of the present invention without departing from the technical solution content of the present invention still fall within the protection scope of the technical solution of the present invention.

Claims

1. A calculation method for long-term permeability damage of artificially filled fractures after hydraulic fracturing, characterized in that, It includes the following steps: Step A: Calculate the distribution law of the macroscopic fracture flow velocity to obtain the fluid flow distribution curve in the fracture; Step B: Determine the capture coefficient based on the particle trajectory, including: Step B1. Construct a deep filtration grid model: Discretize the artificial fracture into a set of interconnected percolation units. Each percolation unit includes an inlet and an outlet in both the horizontal and vertical directions, and determine the geometric structure characteristic parameters of the percolation unit, including the pore radius R p , the throat radius R t , the pore region length L p and the throat region length L t ; Step B2: Calculate the particle capture coefficients in the vertical and horizontal directions: (B21) Calculate the flow velocity distribution characteristics in each infiltration unit, analyze the forces acting on the particles invading the infiltration unit, and determine the sum of the velocities of the particles in the horizontal and vertical directions and their trajectories in the infiltration unit; (B22) Determine the particle capture coefficients in the horizontal and vertical directions according to the calculation results of (B21); Among them, λ h is the particle capture coefficient in the horizontal direction, and λ v is the particle capture coefficient in the vertical direction. Q in is the flow rate value injected into the infiltration unit, and φ is the porosity within the unit; V e is the volume of the infiltration unit, and l vd represents the displacement in the vertical direction of the particles injected horizontally, and l hd represents the displacement in the horizontal direction of the particles injected vertically; Step C: Calculate the fracture permeability damage in a single time step; Step C1: Determine the amount of particles invading along the linear flow perpendicular to the fracture in each infiltration unit, and calculate the total particle concentration entering the infiltration unit; Step C2: Solve the discrete deep filtration grid model, inject reservoir fluid into the fracture at each time step, and calculate the fracture permeability damage distribution caused by particle plugging; Step D: After the end of each time step, process the blocked units, characterize the blocked infiltration units as particle source terms; increment the time step by 1, and repeat Steps C and D until the maximum calculation time is reached to achieve the calculation of the permeability damage.

2. The calculation method for long-term permeability damage of artificial filled fractures after hydraulic fracturing according to claim 1, characterized in that: In Step B2, the flow velocity distribution in the infiltration unit is calculated by the following method: Where: R0 is the inlet radius, for the pore control region R0 = R p , for the throat control region R0 = R t ; R p is the pore radius, R t is the throat radius, u0 is the fluid velocity at the inlet, R(x) is the flow channel radius corresponding to the position x, x is the horizontal position within the filtration unit, and u(r,x) is the velocity magnitude at the horizontal position x and the flow channel radius r.

3. The calculation method for long-term permeability damage of artificial filled fractures after hydraulic fracturing according to claim 1, wherein: Step C1 is specifically implemented by the following method: (1) Based on the fluid flow distribution curve in the fracture obtained in Step A, calculate the linear flow rate distribution perpendicular to the fracture direction under the condition of bilinear flow; (2) Determine the proportion of solid-phase particles in the target reservoir fluid to determine the amount of particle invasion in each infiltration unit; (3) Further determine the total amount of particles injected into the infiltration unit: R total = R pi + C h (x - 1, y, t)+ C v (x, y + 1, t) where, R pi is the amount of particle intrusion into the unit in the direction perpendicular to the crack, and R total is the total particle concentration entering the unit, C h is the particle concentration at the horizontal inlet, and C v is the particle concentration at the vertical inlet.

4. The method for calculating the long-term permeability damage of the artificial filled fracture after hydraulic fracturing according to claim 3, characterized in that: Step C2 is specifically implemented by the following method: (1) For the infiltration unit with coordinates (x, y), the horizontal and vertical outlets of the target infiltration unit are the inlets of the horizontally adjacent unit and the vertically adjacent unit respectively, and the particles in the unit satisfy the following mass conservation: Among them, C v (x, y, t) is the particle concentration at the vertical inlet, C h (x, y, t) is the particle concentration at the horizontal inlet, σ e (x, y, t) is the variation law of the particle concentration with time; Apply the finite difference method to discretize and solve each term in the above formula, and the discretization method is shown in the following formula: (2) Set the initial and boundary conditions of the model, solve the mass conservation equation, and obtain the variation law of the particle concentration trapped in the percolation unit with time σ e (x, y, t).

5. The calculation method for long-term permeability damage of artificial filled fractures after hydraulic fracturing according to claim 1, characterized in that: In Step D, the blocked unit is characterized as a particle source term, and its calculation formula is as follows: C h (x + 1, y, t) = λ rm D c u(x, y, t) where C h (x + 1, y, t) is the particle concentration at the horizontal outlet, and λ rm is the re-migration coefficient of the blocked unit, and D c is the viscosity coefficient.