Fracture-wellbore system proppant flowback numerical simulation method and apparatus

Abstract

The application discloses a kind of fracture-wellbore system proppant backflow numerical simulation method and device, it is related to petroleum engineering field.The specific implementation scheme includes: determining fracture-wellbore simulation model, the attribute of liquid phase, the attribute of solid phase;Fracture-wellbore simulation model is set simulation condition;According to the attribute of solid phase, determine the time step of discrete element method and the time step of computational fluid dynamics method;Based on liquid phase control equation, solid phase control equation, liquid phase on solid phase force equation, solid phase on liquid phase force equation, the time step of discrete element method and the time step of computational fluid dynamics method, in fracture-wellbore simulation model, according to the attribute of liquid phase and the attribute of solid phase, the coupling calculation of computational fluid dynamics method and discrete element method is carried out, solid phase, liquid phase is simulated and simulated, to realize the proppant backflow numerical simulation of fracture-wellbore system.The present application can improve the accuracy of fracture-wellbore system proppant backflow numerical simulation.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of petroleum engineering, and in particular to a fracture-wellbore system proppant flowback numerical simulation method and device. BACKGROUND

[0002] In the fracture-wellbore proppant system, proppant flowback is a key problem. Proppant flowback refers to the phenomenon that proppants return from the fracture to the wellbore during oil and gas production. This phenomenon can cause the flow conductivity of the fracture to decrease, thereby affecting the oil and gas recovery. By performing proppant flowback numerical simulation, the mechanism of proppant flowback can be deeply understood, and the possibility and degree of flowback can be predicted in advance, so that effective preventive measures can be taken.

[0003] At present, the CFD-DEM coupling algorithm (i.e., the coupling algorithm of Computational Fluid Dynamics (CFD) and Discrete Element Method (DEM)) is combined with simulation software to simulate the proppant flowback value of the fracture-wellbore system.

[0004] However, the accuracy of the simulation of the proppant flowback value of the fracture-wellbore system is poor at present. SUMMARY

[0005] The embodiments of the present application provide a fracture-wellbore system proppant flowback numerical simulation method and device, solve the problem of poor accuracy of numerical simulation in the prior art, and improve the accuracy of fracture-wellbore system proppant flowback numerical simulation.

[0006] In a first aspect, the embodiments of the present application provide a fracture-wellbore system proppant flowback numerical simulation method, including:

[0007] According to the actual fracture-wellbore system, the physical properties of the actual fracturing fluid, and the physical properties of the actual proppant, the fracture-wellbore simulation model, the properties of the liquid phase, and the properties of the solid phase are determined, and the solid phase and the liquid phase are used for simulation in the fracture-wellbore simulation model; the simulation conditions of the fracture-wellbore simulation model are set; according to the physical properties of the solid phase, the time step of the discrete element method and the time step of the computational fluid dynamics method are determined; based on the liquid phase control equation, the solid phase control equation, the liquid phase to solid phase force equation, the solid phase to liquid phase force equation, the time step of the discrete element method, and the time step of the computational fluid dynamics method, in the fracture-wellbore simulation model, according to the properties of the liquid phase and the properties of the solid phase, the coupling calculation of the computational fluid dynamics method and the discrete element method is performed, and the solid phase and the liquid phase are simulated to realize the proppant flowback numerical simulation of the fracture-wellbore system.

[0008] The liquid phase control equations include a continuity equation, a momentum equation and a turbulence model, and the turbulence model includes a turbulent kinetic energy equation and a dissipation rate equation.

[0009] The continuity equation is shown in the following equation (7):

[0010] (7)

[0011] In equation (7), is the volume fraction of the calculation grid occupied by the liquid phase particles; is the density of the liquid phase; is the motion time; is the flow velocity of the liquid phase;

[0012] The momentum equation is shown in the following equation (8):

[0013] (8)

[0014] In equation (8), P is the fluid pressure; is the fluid viscosity; is the gravitational acceleration; is the volume force of the interaction between the solid phase and the liquid phase;

[0015] The turbulent kinetic energy equation is shown in the following equation (11):

[0016] (11)

[0017] In equation (11), is the volume fraction of the liquid phase; is the density of the liquid phase; is the turbulent kinetic energy of the liquid phase; is the dissipation rate of the turbulent kinetic energy; represents the liquid phase velocity vector; is the liquid phase viscosity coefficient; is the dimensionless Prandtl number corresponding to the turbulent kinetic energy; is the source term of the turbulent kinetic energy; is the solid phase exchange coefficient;

[0018] The dissipation rate equation is shown in the following equation (12):

[0019] (12)

[0020] In equation (12), is the continuous phase viscosity coefficient; is the dimensionless Prandtl number corresponding to the turbulent kinetic energy; C 1ε , C 2ε are empirical constants; is the liquid phase exchange coefficient.

[0021] Further, according to the attribute of the solid phase, the time step of the discrete element method and the time step of the computational fluid dynamics method are determined, comprising:

[0022] According to the solid phase particle diameter, the particle Poisson's ratio and the particle shear modulus of the solid phase, the time step of the discrete element method is determined; according to the time step of the discrete element method, the time step of the computational fluid dynamics method is determined.

[0023] Further, the solid phase control equation comprises a solid phase particle motion equation and a solid phase particle contact equation.

[0024] The solid phase particle contact equation is shown in the following formula (13):

[0025] (13)

[0026] In formula (13), is the contact force acting on the solid phase particle j ; i , , , respectively represent the tangential elastic force, the tangential damping, the normal elastic force and the normal damping of the solid phase particle j to the solid phase particle i .

[0027] The normal elastic force of the solid phase particle j to the solid phase particle i is obtained by the following formula (14):

[0028] (14)

[0029] Wherein, , ;

[0030] In formula (14), Y is the Young's modulus of the solid phase particle; R is the equivalent radius of the solid phase particle; delta n is the normal distance between the solid phase particle j and the solid phase particle i ; , respectively are the Young's modulus of the solid phase particle i and the solid phase particle j ; , respectively are the Poisson's ratio of the solid phase particle i and the solid phase particle j ; , ​respectively, are the radii of the solid phase particles i respectively, are the radii of the solid phase particles j .

[0031] respectively, are the radii of the solid phase particles j respectively, are the radii of the solid phase particles i The normal damping of the solid phase particle

[0032] (15)

[0033] where, , , ;

[0034] In equation (15), is the equivalent mass; is the normal component of the relative velocity of the solid phase particle i and the solid phase particle j ; , respectively, are the masses of the solid phase particle i and the solid phase particle j ; is the coefficient of restitution of collision; S n is the normal stiffness.

[0035] The tangential elastic force of the solid phase particle j on the solid phase particle i is obtained by equation (16) as follows:

[0036] (16)

[0037] where, , ;

[0038] In equation (16), denotes the tangential stiffness, G is the equivalent shear modulus, denotes the tangential overlap of the solid phase particle j and the solid phase particle i .

[0039] The tangential damping of the solid phase particle j on the solid phase particle i is obtained by equation (17) as follows:

[0040] (17)

[0041] In equation (17), v t is the tangential component of the relative velocity of the solid phase particle i and the solid phase particle jThe tangential component of the relative velocity.

[0042] Furthermore, the equations of motion for solid particles include translational and rotational equations.

[0043] The translation equation is shown in equation (18) below:

[0044] (18)

[0045] In equation (18), m i solid particles i The quality; solid particles linear velocity; Liquid phase particles i The force; solid particles j Or the wall surface of solid particles i The applied normal elastic force; solid particles j Or the wall surface of solid particles i The applied tangential elastic force; It is the acceleration due to gravity; solid particles i Total number of contacts with other particles and the wall surface.

[0046] The equation of rotation is shown in equation (19) below:

[0047] (19)

[0048] In equation (19), solid particles i Moment of inertia; omega i solid particles i angular velocity; M t,ij Indicates solid particles j For particles i Torque caused by tangential force; M r,ij Indicates solid particles j For particles i Torque caused by rolling friction; z represents solid particles. i Total number of contacts with other solid particles and solid wall surfaces.

[0049] Furthermore, the equations of the interaction force between the liquid phase and the solid phase include the buoyancy model of solid particles under the influence of the liquid phase and the drag model of solid particles under the influence of the liquid phase.

[0050] The solid particles subjected to buoyancy by the liquid phase are shown in the following equation (20):

[0051] (20)

[0052] In formula (20), is the buoyancy of the solid phase particle from the liquid phase, is the particle size of the solid phase particle; is the density of the liquid phase; is the acceleration of gravity.

[0053] Further, the solid phase particle is subjected to the liquid phase drag force model, as shown in the following formula (21):

[0054] (21)

[0055] wherein, , ;

[0056] In formula (21), is the liquid phase drag force on the solid phase particle; is the drag coefficient of the liquid phase to the solid phase; is the volume fraction of the solid phase particle in the calculation grid; is the volume fraction of the liquid phase in the calculation grid; is the velocity of the liquid phase at the position of the solid phase particle; is the velocity of the solid phase particle; is the particle size of the solid phase particle; is the Reynolds number of the solid phase particle.

[0057] Further, the solid phase to the liquid phase force equation includes the liquid phase subjected to the solid phase particle drag force model.

[0058] The liquid phase subjected to the solid phase particle drag force model is shown in the following formula (22):

[0059] (22)

[0060] In formula (22), is the drag force of the liquid phase from the solid phase particle; is the force applied by the liquid phase to the solid phase particle; Δ V is the volume of the calculation grid; n is the number of solid phase particles contained in the calculation grid.

[0061] In a second aspect, the embodiments of the present application provide a fracture-wellbore system proppant flowback numerical simulation device, comprising: a setting module and a calculation simulation module.

[0062] The configuration module is used to determine the properties of the fracture-wellbore simulation model, the liquid phase, and the solid phase based on the actual fracture-wellbore system, the physical properties of the actual fracturing fluid, and the physical properties of the actual proppant. Both the solid and liquid phases are used for simulation in the fracture-wellbore simulation model. The module also sets the simulation conditions for the fracture-wellbore simulation model and determines the time step for the discrete element method and the computational fluid dynamics method based on the properties of the solid phase.

[0063] The computational simulation module is used to perform coupled calculations of computational fluid dynamics and discrete element method in the fracture-wellbore simulation model based on the liquid phase control equation, solid phase control equation, liquid phase interaction force equation, solid phase interaction force equation, time step of discrete element method and time step of computational fluid dynamics method, according to the properties of liquid phase and solid phase. It simulates the solid and liquid phases to realize the numerical simulation of proppant backflow in fracture-wellbore system.

[0064] The liquid phase governing equations include the continuity equation, the momentum equation, and the turbulence model. The turbulence model includes the turbulent kinetic energy equation and the dissipation rate equation.

[0065] The continuity equation is shown below:

[0066]

[0067] In the formula, This represents the volume fraction of the computational grid occupied by liquid phase particles; The density of the liquid phase; For exercise time; The velocity is the liquid phase flow rate.

[0068] The momentum equation is as follows:

[0069]

[0070] In the formula, P For fluid pressure; It is a fluid viscosity; It is the acceleration due to gravity; It is the volume force of the interaction between the solid and liquid phases.

[0071] The turbulent kinetic energy equation is shown below:

[0072]

[0073] In the formula, It represents the liquid volume fraction; The density of the liquid phase; The turbulent kinetic energy of the liquid phase; The dissipation rate of turbulent kinetic energy; Represents the liquid phase velocity vector; is the liquid phase viscosity coefficient; is the non-dimensional Prandtl number corresponding to turbulent kinetic energy; is the source term of turbulent kinetic energy; is the solid phase exchange coefficient.

[0074] The dissipation rate equation is as follows:

[0075]

[0076] In the formula, is the continuous phase viscosity coefficient; is the non-dimensional Prandtl number corresponding to turbulent kinetic energy; C 1ε , C 2ε are all empirical constants; is the liquid phase exchange coefficient.

[0077] In a third aspect, an embodiment of the present application provides a device, the device comprising: a processor; a memory for storing processor-executable instructions; and the processor implements the method of the first aspect or any possible implementation manner of the first aspect when executing the executable instructions.

[0078] In a fourth aspect, an embodiment of the present application provides a nonvolatile computer readable storage medium, the nonvolatile computer readable storage medium comprising a computer program or instructions for storing, when the computer program or instructions are executed, causing the method of the first aspect or any possible implementation manner of the first aspect to be implemented.

[0079] The one or more technical solutions provided in the embodiments of the present application have at least the following technical effects or advantages:

[0080] The embodiments of the present application can accurately simulate the solid phase and the liquid phase by performing coupling calculation of the computational fluid dynamics method and the discrete element method in the fracture-wellbore simulation model obtained according to the actual fracture-wellbore system based on the liquid phase control equation, the solid phase control equation, the liquid phase to solid phase force equation, the solid phase to liquid phase force equation, the time step of the discrete element method and the time step of the computational fluid dynamics method, thereby accurately performing the proppant flowback numerical simulation of the fracture-wellbore system and improving the accuracy of the proppant flowback numerical simulation of the fracture-wellbore system. BRIEF DESCRIPTION OF DRAWINGS

[0081] In order to more clearly illustrate the technical solutions of the embodiments of the present application, the following will briefly introduce the drawings needed to be used in the embodiments of the present application or the prior art description. Obviously, the drawings in the following description are some embodiments of the present application, and those skilled in the art can also obtain other drawings according to these drawings without creative labor.

[0082] Figure 1 A flowchart of a fracture-wellbore system proppant flowback numerical simulation method provided by the embodiments of the present application is shown in FIG. 1.

[0083] Figure 2 A CFD-DEM coupling calculation flowchart is shown in FIG. 2.

[0084] Figure 3 A schematic diagram of real experimental results of the verification scheme is shown in FIG. 3.

[0085] Figure 4 A schematic diagram of simulation results of the verification scheme is shown in FIG. 4.

[0086] Figure 5 A comparison chart of a real flowback experiment and simulation is shown in FIG. 5.

[0087] Figure 6 A comparison chart of another real flowback experiment and simulation is shown in FIG. 6.

[0088] Figure 7 A comparison chart of yet another real flowback experiment and simulation is shown in FIG. 7.

[0089] Figure 8 A comparison chart of yet another real flowback experiment and simulation is shown in FIG. 8.

[0090] Figure 9 A composition schematic diagram of a fracture-wellbore system proppant flowback numerical simulation device provided by the embodiments of the present application is shown in FIG. 9. DETAILED DESCRIPTION

[0091] The technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative work fall within the scope of protection of the present application.

[0092] The following descriptions of some technologies related to the embodiments of the present application are provided to help understanding, and should be considered as merely exemplary. Therefore, those of ordinary skill in the art should recognize that various changes and modifications can be made to the embodiments described herein without departing from the scope and spirit of the present application. Also, some descriptions of well-known functions and structures are omitted from the following description for clarity and brevity.

[0093] In the fracture-wellbore proppant system, proppant flowback is a key issue. Proppant flowback refers to the phenomenon that proppants return from the fracture to the wellbore during oil and gas production. This phenomenon can lead to a decrease in the fracture conductivity, thereby affecting the oil and gas recovery. By performing numerical verification of proppant flowback, the mechanism of proppant flowback can be understood in depth, and the possibility and degree of flowback can be predicted in advance, so that effective preventive measures can be taken.

[0094] At present, the numerical simulation of proppant flowback in the fracture-wellbore system is generally performed by combining a CFD-DEM coupling algorithm (i.e., a coupling algorithm of computational fluid dynamics (CFD) and discrete element method (DEM)) with simulation software.

[0095] However, the numerical simulation of proppant flowback in the fracture-wellbore system has poor accuracy at present.

[0096] In this background technology, the present disclosure provides a numerical simulation method of proppant flowback in a fracture-wellbore system, which can improve the accuracy of the numerical simulation of proppant flowback in the fracture-wellbore system.

[0097] The execution subject of the numerical simulation method of proppant flowback in the fracture-wellbore system provided by the embodiments of the present disclosure can be a computer or a server, or can also be other electronic devices with data processing capability; or the execution subject of the method can also be a processor (such as a central processing unit (CPU)) in the above-mentioned electronic devices; or the execution subject of the method can also be an application program (application, APP) installed in the above-mentioned electronic devices and capable of realizing the function of the method; or the execution subject of the method can also be a functional module or unit with the function of the method in the above-mentioned electronic devices, etc. The execution subject of the method is not limited herein.

[0098] In some embodiments, the server can be a single server, or can also be a server cluster composed of multiple servers. In some embodiments, the server cluster can also be a distributed cluster. The specific implementation of the server is not limited by the present disclosure.

[0099] The numerical simulation method of proppant flowback in the fracture-wellbore system will be described below with reference to the accompanying drawings.

[0100] Figure 1 is a flowchart of the numerical simulation method of proppant flowback in the fracture-wellbore system provided by the embodiments of the present disclosure. In the method, Figure 1The execution sequence shown is only one of the embodiments of the present application, and does not represent the only execution sequence of the fracture-wellbore system proppant flowback numerical simulation method. In the case of achieving the final result, Figure 1 The steps shown can be performed in parallel or in reverse. As shown, Figure 1 The method can include:

[0101] S101, according to the actual fracture-wellbore system, the physical properties of the actual fracturing fluid, and the physical properties of the actual proppant, determine the fracture-wellbore simulation model, the properties of the liquid phase, and the properties of the solid phase.

[0102] Exemplarily, the corresponding fracture-wellbore simulation model can be generated according to the actual physical properties (such as wellbore size, position, etc.) of the actual fracture-wellbore system that needs to be studied and verified by using a geometric model processing software.

[0103] Exemplarily, a three-dimensional fracture plate (i.e. a fracture-wellbore simulation model) can be established by using CAD software Solidworks, with a fracture size of 1m x 0.3m x 5mm, 5 inlet ports on the left side simulating perforations, a diameter of 6mm, and a length of 6mm simulating casing wall thickness, a right side surface as a pressure outlet, and other surfaces as walls.

[0104] Exemplarily, the properties of the liquid phase and the solid phase are the simulation properties of the liquid phase and the solid phase during simulation.

[0105] Exemplarily, the physical properties of the actual fracturing fluid can be determined as the properties of the liquid phase, and the physical properties of the actual proppant can be determined as the properties of the solid phase.

[0106] Exemplarily, the physical properties of the actual fracturing fluid can include fracturing fluid density, fracturing fluid viscosity, and fracturing fluid Reynolds number, and the physical properties of the actual proppant can include proppant particle size, proppant Young's modulus, and proppant Poisson's ratio, without limitation.

[0107] Exemplarily, when determining the properties of the liquid phase and the solid phase, appropriate adjustments can also be made according to the physical properties of the actual fracturing fluid and the actual proppant. For example, taking silica sand with a particle size range of 0.25mm to 0.50mm as an example, the solid phase particle size can be approximately determined as the average particle size of the proppant, i.e. 0.375mm.

[0108] S102, set the simulation conditions of the fracture-wellbore simulation model.

[0109] Exemplarily, the actual fracture-wellbore system to be verified can be studied according to needs, and corresponding boundary conditions and other simulation conditions can be set for the fracture-wellbore simulation model. For example, the liquid phase inlet condition can be set as a velocity inlet, the liquid phase enters the fracture at a set injection rate, the fracture wall surface is set as a no-slip boundary condition, and the model outlet is set as a pressure outlet.

[0110] Exemplarily, the turbulent intensity of the liquid phase can be determined by a turbulent intensity determination formula.

[0111] The turbulent intensity determination formula is shown in the following formula (1):

[0112] (1)

[0113] In formula (1), is the turbulent intensity of the liquid phase; is the Reynolds number of the liquid phase.

[0114] The Reynolds number Re can be obtained by first calculating the product of the average flow velocity of the liquid phase and the hydraulic diameter, and then calculating the ratio of the product to the kinematic viscosity of the liquid phase; the hydraulic diameter can be obtained by the following formula (2):

[0115] (2)

[0116] In formula (2), is the hydraulic diameter; is the perimeter of the flow cross section of the circular pipe; is the area of the flow cross section.

[0117] Exemplarily, the tangential velocity and temperature of the solid phase particles at the wall surface can be obtained by the Johnson-Jackson model. The Johnson-Jackson model is shown in the following formula (3) and formula (4):

[0118]

[0119] In formula (3) and formula (4), is the slip velocity of the solid phase particles at the wall surface; is the shear tensor generated by collision between the solid phase particles and the wall surface; is the shear tensor generated by friction between the solid phase particles and the wall surface; is the normal vector pointing to the inside of the wall surface; is the reflection coefficient; is the friction angle of the proppant and the wall surface; is the energy conduction coefficient of the solid phase particles; is the apparent temperature of the solid phase particles; is the packing concentration of the solid phase particles; ​is the maximum packing concentration of the solid phase particles; is the restitution coefficient of the solid phase particles and the wall, and is generally 0.9; is the density of the solid phase particles; is the contact radial distribution function of the solid phase particles.

[0120] Exemplarily, the generation rate of the solid phase particles can be determined according to the physical properties of the actual fracturing fluid and the actual proppant, and the properties of the solid phase and the liquid phase, by using a solid phase particle generation rate determination formula.

[0121] The solid phase particle generation rate determination formula can be shown in the following formula (5):

[0122] (5)

[0123] wherein, .

[0124] In the formula (5), is the generation rate of the solid phase particles; is the injection velocity of the solid phase particles; is the fracture width; is the height of the simulated fracture; is the mass concentration of the solid phase; is the density of the sand-carrying fluid (i.e., the mixture of the fracturing fluid and the proppant); is the radius of the solid phase particles; is the density of the solid phase particles; is the density of the fracturing fluid.

[0125] S103. Determine the time step of the discrete element method and the time step of the computational fluid dynamics method according to the properties of the solid phase.

[0126] It should be noted that the coupling of the computational fluid dynamics method (CFD) and the discrete element method (DEM) requires discretization of the calculation domain, and the flow field changes of the solid-liquid two phases are solved according to the time step, and then the changes of the particle position, velocity and other properties are obtained, so the determination of the time step has an important influence on the calculation efficiency and accuracy.

[0127] In the DEM, its time step is related to the particle size, and the specific calculation value is the percentage of Rayleigh Time, that is, the smaller the particle size, the smaller the time step required to meet the calculation requirements. For particles with fixed particle size, when the time step is set too large, the calculation of the motion properties of the particles is not complete, the particles contact each other and the overlap amount is large, which is manifested as irregular bouncing of the particles, penetration of the model wall, etc., thereby leading to non-convergence of the flow field calculation; if the time step is set too small, although the calculation of the motion properties of the particles after collision is more complete and accurate, the calculation time of the dense particle flow is longer, the calculation efficiency is extremely low, and it does not meet the simulation requirements of large-scale proppant particle flow.

[0128] The coupling of CFD-DEM sets the time step of CFD as an integer multiple of the time step of DEM, which can complete the calculation of multiple time steps of DEM particle motion in one complete CFD time step, and then transfer the motion properties back to CFD for coupled calculation, which makes the calculation more continuous and greatly improves the calculation efficiency.

[0129] In some possible embodiments, determining the time step of the discrete element method and the time step of the computational fluid dynamics method according to the properties of the solid phase can include steps 1 and 2.

[0130] Step 1, determining the time step of the discrete element method according to the solid phase particle diameter, the solid phase particle Poisson's ratio and the solid phase particle shear modulus of the solid phase.

[0131] For example, the Rayleigh Time can be determined according to the solid phase particle diameter, the solid phase particle Poisson's ratio and the solid phase particle shear modulus, and the time step of the discrete element method not greater than the Rayleigh Time can be determined according to the Rayleigh Time and the specific situation of the computing resources.

[0132] For example, the Rayleigh Time can be determined by the following formula (6):

[0133] (6)

[0134] In formula (6), is the Rayleigh Time; is the solid phase particle density; is the solid phase particle shear modulus; is the solid phase particle Poisson's ratio.

[0135] For example, the time step of the discrete element method can be 15.57% of the Rayleigh Time.

[0136] Step 2, determining the time step of the computational fluid dynamics method according to the time step of the discrete element method.

[0137] Exemplarily, according to the specific situation of the computing resource, an integer multiple of the time step of the discrete element method can be determined as the time step of the computational fluid dynamics method under the guarantee of computing efficiency. For example, the time step of the computational fluid dynamics method can be 100 times the time step of the discrete element method. Taking the time step of the discrete element method as 5×10 -7 s, the time step of the computational fluid dynamics method is 5×10 -5 s.

[0138] It should be noted that before S104 is executed, the fracture-wellbore simulation model needs to be discretized and meshed to divide the fracture-wellbore simulation model into a number of calculation grids.

[0139] The fracture-wellbore simulation model preprocessing can use ICEM software, which has powerful CAD model repair capabilities, unique meshing and editing techniques, and extensive solver support capabilities. The fracture-wellbore simulation model is cut using this software, and the O-block structure is used to divide the grid. Viscous fluid is used in the simulation process, and the fluid flowing in the narrow fracture has the maximum velocity gradient near the wall surface, so the inlet, outlet and boundary are encrypted and an expansion layer is added to improve the calculation accuracy.

[0140] Generally speaking, the more dense the calculation grid of the fracture-wellbore simulation model is, the more nodes it has, and the more accurate the calculation result is, so the fracture-wellbore simulation model is usually encrypted, especially at irregular model structures. However, the denser the calculation grid is, the larger the amount of calculation is, and the longer the calculation period is, and the higher the requirements for computing resources such as CPU and memory size of the computer are, otherwise, the calculation result will deviate and the accuracy will decrease. Therefore, the calculation grid of the fracture-wellbore simulation model is verified for independence, the density of the calculation grid is changed, and 420,000, 600,000 and 810,000 calculation grid numbers are divided to test the calculation result and the calculation speed. The preliminary simulation calculation result shows that the calculation results of 600,000 and 810,000 calculation grid numbers are basically the same, and both are more accurate than the calculation result of 420,000 calculation grid numbers. Considering the computer performance and the calculation speed, 620,000 calculation grid numbers are preferred.

[0141] After the calculation grid of the model is divided, the calculation grid quality can be checked. When the calculation grid quality is greater than a preset threshold (such as 0.7), it can be determined that the calculation grid quality check is passed, otherwise, the calculation grid quality check is not passed. After the calculation grid quality check is passed, S104 can be continued.

[0142] S104, based on the liquid phase control equation, the solid phase control equation, the liquid phase to solid phase force equation, the solid phase to liquid phase force equation, the time step of the discrete element method and the time step of the computational fluid dynamics method, in the fracture-wellbore simulation model, according to the properties of the liquid phase and the properties of the solid phase, the coupling calculation of the computational fluid dynamics method and the discrete element method is carried out, and the solid phase and the liquid phase are simulated to realize the proppant backflow numerical simulation of the fracture-wellbore system.

[0143] Specifically, the liquid phase control equation includes a continuity equation, a momentum equation and a turbulence model, and the turbulence model includes a turbulent kinetic energy equation and a dissipation rate equation.

[0144] The continuity equation is shown in the following formula (7):

[0145] (7)

[0146] In formula (7), is the volume fraction of the liquid phase particles in the calculation grid; is the liquid phase density; is the motion time; is the liquid phase flow velocity.

[0147] The momentum equation is shown in the following formula (8):

[0148] (8)

[0149] In formula (8), P is the fluid pressure; is the fluid viscosity; is the gravitational acceleration; is the volume force of the interaction between the solid phase and the liquid phase.

[0150] The volume fraction of the liquid phase particles in the calculation grid is obtained by the following formula (9):

[0151] (9)

[0152] In formula (9), Δ V is the volume of the calculation grid; is the volume of the solid phase particles i in the calculation grid; is the number of solid phase particles contained in the calculation grid.

[0153] The volume force of the interaction between the solid phase and the liquid phase is obtained by the following formula (10):

[0154] (10)

[0155] In formula (10), is the volume of the solid phase particlesi viscous resistance of the liquid phase in the flow process.

[0156] The turbulent kinetic energy equation is shown in the following equation (11):

[0157] (11)

[0158] In equation (11), is the liquid volume fraction; is the liquid density; is the turbulent kinetic energy of the liquid phase; is the dissipation rate of the turbulent kinetic energy; represents the liquid phase velocity vector; is the liquid viscosity coefficient; is the dimensionless Prandtl number corresponding to the turbulent kinetic energy; is the source term of the turbulent kinetic energy; is the liquid exchange coefficient.

[0159] The dissipation rate equation is shown in the following equation (12):

[0160] (12)

[0161] In equation (12), is the continuous phase viscosity coefficient; is the dimensionless Prandtl number corresponding to the turbulent kinetic energy; C 1ε , C 2ε are empirical constants; is the liquid exchange coefficient.

[0162] Specifically, the solid phase control equation includes a solid phase particle motion equation and a solid phase particle contact equation.

[0163] The solid phase particle contact equation is shown in the following equation (13):

[0164] (13)

[0165] In equation (13), is the solid phase particle j contact force acting on the solid phase particle i ; , , , respectively represent the tangential elastic force, the tangential damping, the normal elastic force, and the normal damping of the solid phase particle j on the solid phase particle i .

[0166] The normal elastic force of the solid phase particle j on the solid phase particle i is obtained by the following equation (14):

[0167] (14)

[0168] wherein, , .

[0169] In formula (14), Y is the Young's modulus of the solid-phase particle; R is the equivalent radius of the solid-phase particle; delta n is the solid-phase particle j and the solid-phase particle i ; , are the Young's modulus of the solid-phase particle i and the solid-phase particle j , respectively; , are the Poisson's ratio of the solid-phase particle i and the solid-phase particle j , respectively; , are the particle radius of the solid-phase particle i and the solid-phase particle j , respectively.

[0170] The normal damping of the solid-phase particle j to the solid-phase particle i is obtained by the following formula (15):

[0171] (15)

[0172] wherein, , , .

[0173] In formula (15), is the equivalent mass; is the normal component of the relative velocity of the solid-phase particle i to the solid-phase particle j ; , are the mass of the solid-phase particle i and the solid-phase particle j , respectively; is the collision restitution coefficient; S n is the normal stiffness.

[0174] The tangential elastic force of the solid-phase particle j to the solid-phase particle i is obtained by the following formula (16):

[0175] (16)

[0176] wherein, , .

[0177] In formula (16), represents tangential stiffness, G is equivalent shear modulus, represents the tangential overlap amount of the solid phase particle j and the solid phase particle i .

[0178] The tangential damping of the solid phase particle j to the solid phase particle i is obtained by the following formula (17):

[0179] (17)

[0180] In formula (17), is the relative velocity tangential component of the solid phase particle i and the solid phase particle j .

[0181] In some possible embodiments, the solid phase particle motion equation includes a translation equation and a rotation equation.

[0182] The translation equation is shown in the following formula (18):

[0183] (18)

[0184] In formula (18), m i is the mass of the solid phase particle i ; is the linear velocity of the solid phase particle ; is the force of the liquid phase to the solid phase particle i ; is the normal elastic force of the solid phase particle j or the wall to the solid phase particle i ; is the tangential elastic force of the solid phase particle j or the wall to the solid phase particle i ; is the gravitational acceleration; is the total contact number of the solid phase particle i and other particles and the wall.

[0185] The rotation equation is shown in the following formula (19):

[0186] (19)

[0187] In formula (19), is the moment of inertia of the solid phase particle i . omega i is the angular velocity of the solid phase particle i . M t,ij represents the moment of force of the solid phase particle j on the particle i caused by the tangential force; M r,ij represents the moment of force of the solid phase particle j on the particle i caused by the rolling friction; z is the total contact number of the solid phase particle i with other solid phase particles and solid wall surfaces.

[0188] Specifically, the liquid phase force equation on the solid phase includes a solid phase particle buoyancy model and a solid phase particle drag force model.

[0189] The solid phase particle buoyancy model is shown in the following formula (20):

[0190] (20)

[0191] In formula (20), is the buoyancy of the solid phase particle from the liquid phase, is the particle size of the solid phase particle; is the density of the liquid phase; is the acceleration of gravity.

[0192] Specifically, the solid phase particle drag force model is shown in the following formula (21):

[0193] (21)

[0194] wherein, , .

[0195] In formula (21), is the liquid phase drag force on the solid phase particle; is the liquid phase resistance coefficient on the solid phase; is the volume fraction of the solid phase particle in the calculation grid; is the volume fraction of the liquid phase particle in the calculation grid; is the velocity of the liquid phase at the position of the solid phase particle; is the velocity of the solid phase particle; is the particle size of the solid phase particle; is the Reynolds number of the solid phase particle.

[0196] Specifically, the solid phase force equation on the liquid phase includes a liquid phase resistance model of the solid phase particle.

[0197] The liquid phase resistance to the solid phase particle is shown in the following equation (22):

[0198] (22)

[0199] In equation (22), is the liquid phase resistance to the solid phase particle; is the force applied by the liquid phase to the solid phase particle; Δ V is the volume of the calculation grid; n is the number of solid phase particles contained in the calculation grid.

[0200] Exemplarily, based on the solid phase control equation, the liquid phase control equation, the liquid phase to solid phase force equation, the solid phase to liquid phase force equation, the discrete element method time step and the computational fluid dynamics method time step, the fracture solid-liquid two-phase bidirectional coupling flow simulation is carried out by using the coupling custom model interface in the simulation software Fluent combined with the discrete element simulation software. The basic principle of the solid-liquid two-phase bidirectional coupling is that the interaction between the solid phase particles and the fracture wall is calculated by using the discrete element simulation software, which is combined with CFD. By bidirectional coupling of CFD for calculating continuous fluid and DEM for calculating discrete particles, the proppant transport process is studied by using the CFD particle tracking model, the interaction between the proppants and between the proppants and the wall is studied by using DEM, and the momentum conversion caused by the interaction between the particles and the fluid is used for bidirectional calculation. First, the drag force of the fluid on the particles is calculated in each calculation cell, and then the interaction force between the particles and between the particles and the wall is calculated according to the force-displacement law. Finally, the particle motion is calculated according to Newton's second law. The particle migration causes the change of the porosity in the calculation cell, which further affects the flow field. Between the force-displacement calculation and the motion calculation in DEM, the fluid pressure and velocity are calculated by using CFD, so as to calculate and update the flow field, and the numerical simulation of the proppant flowback in the fracture-wellbore system is completed. The CFD-DEM coupling calculation process is shown in Figure 2 According to the simulation results and the experimental results, the numerical simulation of the proppant flowback in the fracture-wellbore system is realized.

[0201] In this embodiment, based on the liquid phase control equation, the solid phase control equation, the liquid phase to solid phase force equation, the solid phase to liquid phase force equation, the discrete element method time step and the computational fluid dynamics method time step, the coupling calculation of the computational fluid dynamics method and the discrete element method is carried out in the fracture-wellbore simulation model obtained according to the actual fracture-wellbore system, which can accurately simulate the solid phase and the liquid phase, and thus accurately carries out the numerical simulation of the proppant flowback in the fracture-wellbore system, and improves the accuracy of the numerical simulation of the proppant flowback in the fracture-wellbore system.

[0202] Simulation reliability verification

[0203] The simulation of solid phase and liquid phase in the method of the application is verified by numerical simulation of proppant placement and migration in the fracture during the fracturing stage. The verification scheme includes that the viscosity of the fracturing fluid is 2 mPa s, the proppant particle size is 30 / 50 mesh, the proppant density is 2860 kg / m 3 , the proppant volume fraction is 6%, the experimental displacement is 20 L / min, and the simulation speed is 0.2 m / s.

[0204] Figure 3 A schematic diagram of the actual experimental results of the verification scheme is shown in Figure 4 A schematic diagram of the simulation results of the verification scheme is shown in Figure 3 and Figure 4 It can be seen that the error of the sand dam equilibrium height formed by the actual experiment and the simulation is small, and it can be considered that the simulation of the solid phase and the liquid phase is reliable, and can meet the requirements of subsequent numerical verification of proppant flowback in the fracture after fracturing.

[0205] Numerical simulation reliability verification of proppant flowback in the fracture- wellbore system

[0206] The numerical simulation is carried out according to the experimental parameter conversion by establishing a flat fracture model and carrying out grid division.

[0207] Comparison and verification

[0208] Figure 5 A comparison chart of an actual flowback experiment and simulation is shown. Under the conditions of 20 / 40 mesh quartz sand, sand ratio 6%, and flowback fluid viscosity 1 mPa s, as shown in Figure 5 The sand dam shape of the simulation and the actual flowback experiment is basically the same, the actual flowback experiment tests the sand rate, which is the difference between the proportion of the actual initial placement of the sand body cross section to the total cross section and the proportion of the sand body cross section to the total cross section at the actual flowback flow rate of 2.07 cm / s, i.e. 31.6%-27.9%=3.7%, and the simulation result is the difference between the proportion of the simulation initial placement of the sand body cross section to the total cross section and the proportion of the sand body cross section to the total cross section at the simulation flowback flow rate of 2.07 cm / s, i.e. 31.2%-28.5%=3.5%, and the error is about (3.7%-3.5%) / 3.7%=5.4%.

[0209] Figure 6 Another comparison chart of an actual flowback experiment and simulation is shown. Under the conditions of 40 / 70 mesh quartz sand, sand ratio 6%, and flowback fluid viscosity 1 mPa s, as shown in Figure 6As shown in FIG. 6, the sand production rate tested in the actual flowback experiment is the difference between the proportion of the actual initial sand body cross section to the total cross section and the proportion of the sand body cross section at the actual flowback flow rate of 2.07 cm / s to the total cross section, i.e. 28.3%-25.0%=3.3%, and the simulation result is the difference between the proportion of the simulated initial sand body cross section to the total cross section and the proportion of the sand body cross section at the simulated flowback flow rate of 2.07 cm / s to the total cross section, i.e. 28.4%-25.2%=3.2%, and the error is about (3.3%-3.2%) / 3.3%=3%.

[0210] Figure 7 FIG. 7 is another comparison chart of the actual flowback experiment and the simulation. Under the conditions of 0 / 40 mesh quartz sand, sand ratio of 6%, injection flow rate of 10 L / min, and flowback fluid viscosity of 1 mPa·s, the sand dam morphology of the simulation and the actual flowback experiment is basically the same, as shown in FIG. 7. Figure 7 As shown in FIG. 7, the sand production rate tested in the actual flowback experiment is the difference between the proportion of the actual initial sand body cross section to the total cross section and the proportion of the sand body cross section at the actual flowback flow rate of 2.07 cm / s to the total cross section, i.e. 28.6%-25.4%=3.2%, and the simulation result is the difference between the proportion of the simulated initial sand body cross section to the total cross section and the proportion of the sand body cross section at the simulated flowback flow rate of 2.07 cm / s to the total cross section, i.e. 28.1%-25.1%=3.0%, and the error is about (3.2%-3.0%) / 3.2%=6.6%.

[0211] Figure 8 FIG. 8 is another comparison chart of the actual flowback experiment and the simulation. Under the conditions of 40 / 70 mesh quartz sand, sand ratio of 8%, injection flow rate of 15 L / min, and flowback fluid viscosity of 1 mPa·s, the sand dam morphology of the simulation and the actual flowback experiment is basically the same, and the sand production rate tested in the actual flowback experiment is the difference between the proportion of the actual initial sand body cross section to the total cross section and the proportion of the sand body cross section at the actual flowback flow rate of 2.07 cm / s to the total cross section, i.e. 34.3%-29.7%=4.6%, and the simulation result is the difference between the proportion of the simulated initial sand body cross section to the total cross section and the proportion of the sand body cross section at the simulated flowback flow rate of 2.07 cm / s to the total cross section, i.e. 34.1%-29.2%=4.9%, and the error is about (4.9%-4.6%) / 4.6%=6.5%.

[0212] Through the above comparison verification, it can be seen that the error between the simulation result and the actual experiment result is less than 10%, which proves that the method of the present application can accurately simulate the backflow of the proppant in the flowback fracture.

[0213] Although the present application provides method operation steps as embodiments or flowcharts, more or less operation steps can be included based on routine or non-creative labor. The sequence of steps listed in the embodiments is only one of the many ways of executing the steps, and does not represent the only execution sequence. In actual device or client product execution, the method sequence shown in the embodiments or the drawings can be executed in sequence or in parallel (for example, in a parallel processor or multi-thread processing environment).

[0214] As shown in the embodiments of the present application, the present application further provides a proppant flowback numerical simulation device for a fracture-wellbore system, which comprises a setting module 901 and a calculation simulation module 902. Figure 9

[0215] The setting module 901 is configured to determine a fracture-wellbore simulation model, properties of a liquid phase, and properties of a solid phase according to an actual fracture-wellbore system, physical properties of an actual fracturing fluid, and physical properties of an actual proppant, simulate and analyze the solid phase and the liquid phase in the fracture-wellbore simulation model, set simulation conditions of the fracture-wellbore simulation model, and determine a time step of a discrete element method and a time step of computational fluid dynamics according to the properties of the solid phase.

[0216] The calculation simulation module 902 is configured to perform coupled calculation of the computational fluid dynamics and the discrete element method in the fracture-wellbore simulation model according to the properties of the liquid phase and the properties of the solid phase based on a liquid phase control equation, a solid phase control equation, a liquid phase to solid phase force equation, a solid phase to liquid phase force equation, the time step of the discrete element method, and the time step of the computational fluid dynamics, simulate and analyze the solid phase and the liquid phase, and realize proppant flowback numerical simulation of the fracture-wellbore system.

[0217] The liquid phase control equation comprises a continuity equation, a momentum equation, and a turbulent flow model, and the turbulent flow model comprises a turbulent kinetic energy equation and a dissipation rate equation.

[0218] The continuity equation is shown in the following formula:

[0219]

[0220] In the formula, V is a volume fraction of a calculation grid occupied by the liquid phase particle; ρ is a density of the liquid phase; t is a motion time; and u is a flow velocity of the liquid phase.

[0221] The momentum equation is shown in the following formula:

[0222]

[0223] In the formula, p is a fluid pressure. P ​​​​​​ is the fluid viscosity; is the gravitational acceleration; is the volume force of interaction between the solid and liquid phases;

[0224] the turbulent kinetic energy equation is given by:

[0225]

[0226] where, is the liquid volume fraction; is the liquid density; is the turbulent kinetic energy of the liquid phase; is the dissipation rate of turbulent kinetic energy; denotes the liquid phase velocity vector; is the liquid viscosity coefficient; is the dimensionless Prandtl number corresponding to the turbulent kinetic energy; is the source term of the turbulent kinetic energy; is the solid exchange coefficient;

[0227] the dissipation rate equation is given by:

[0228]

[0229] where, is the continuous phase viscosity coefficient; is the dimensionless Prandtl number corresponding to the turbulent kinetic energy; C 1ε , C 2ε are empirical constants; is the liquid exchange coefficient.

[0230] The beneficial effects and specific embodiments of the device can refer to the foregoing method embodiments, and will not be repeated here.

[0231] Some of the modules in the device described in the present application can be described in the general context of computer-executable instructions, such as program modules, which are executed by computers. Generally, program modules include routines, programs, objects, components, data structures, classes, and the like, which perform particular tasks or implement particular abstract data types. The present application can also be practiced in distributed computing environments where tasks are performed by remote processing devices that are connected through a communication network. In a distributed computing environment, program modules can be located in both local and remote computer storage media including storage devices.

[0232] The apparatuses or modules illustrated in the above application examples can be implemented by computer chips or entities, or by products with certain functions. For the convenience of description, the above apparatuses are described as various modules with functions. In the implementation of the application examples, the functions of the modules can be implemented in one or more software and / or hardware. Of course, the modules with certain functions can also be implemented by a combination of multiple sub-modules or sub-units.

[0233] The methods, apparatuses or modules described in the present application can be implemented in a computer readable program code in any appropriate manner, for example, the controller can take the form of, for example, a microprocessor or a processor and a computer readable medium storing computer readable program code (such as software or firmware) executable by the (micro) processor, logic gates, switches, application specific integrated circuits (Application Specific Integrated Circuit; abbreviated as: ASIC), programmable logic controllers and embedded microcontrollers, examples of the controller include but are not limited to the following microcontrollers: ARC 625D, Atmel AT91SAM, Microchip PIC18F26K20 and Silicone Labs C8051F320, the memory controller can also be implemented as part of the control logic of the memory. Those skilled in the art also know that in addition to implementing the controller in a pure computer readable program code, the same function can also be implemented by logically programming the method steps in the form of logic gates, switches, application specific integrated circuits, programmable logic controllers and embedded microcontrollers. Therefore, such a controller can be considered as a hardware component, and the devices included therein for implementing various functions can also be considered as structures within the hardware component. Alternatively, the devices for implementing various functions can also be considered as both software modules implementing the method and structures within the hardware component.

[0234] The embodiments of the present application also provide a device, which comprises: a processor; a memory for storing processor executable instructions; and the processor executes the executable instructions to implement the method as described in the embodiments of the present application.

[0235] The embodiments of the present application also provide a non-volatile computer readable storage medium, which stores a computer program or instructions, and when the computer program or instructions are executed, the method as described in the embodiments of the present application is implemented.

[0236] In addition, the functional modules in each of the embodiments of the present application can be integrated in one processing module, or each module can exist independently, or two or more modules can be integrated in one module.

[0237] The storage medium described above includes, but is not limited to, a random access memory (RAM), a read-only memory (ROM), a cache, a hard disk (HDD), or a memory card. The storage medium can be used to store computer program instructions.

[0238] From the above description of the embodiments, those skilled in the art can clearly understand that the present application can be implemented by means of software and the necessary hardware. Based on such an understanding, the technical solutions of the present application can be embodied in the form of a software product or can be embodied in the form of data migration during implementation. The computer software product can be stored in a storage medium, such as a ROM / RAM, a magnetic disk, an optical disk, etc., and includes a number of instructions for causing a computer device (which can be a personal computer, a mobile terminal, a server, or a network device, etc.) to execute the methods described in the various embodiments or some parts of the embodiments.

[0239] The various embodiments in the specification are described in a progressive manner, and the same or similar parts between the various embodiments can be referred to each other. Each embodiment focuses on the differences from other embodiments. The whole or part of the present application can be used in many general or special computer system environments or configurations. For example: personal computers, server computers, handheld devices or portable devices, tablet devices, mobile communication terminals, multi-processor systems, microprocessor-based systems, programmable electronic devices, network PCs, small computers, large computers, distributed computing environments including any of the above systems or devices, etc.

[0240] The above embodiments are only used to illustrate the technical solutions of the present application, and not to limit the present application; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that they can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacements for some or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the present application.

Claims

1. A method of numerical simulation of proppant flowback in a fracture- wellbore system, characterized by, The method comprises the following steps: determining a fracture-wellbore simulation model, properties of a liquid phase, and properties of a solid phase according to an actual fracture-wellbore system, physical properties of an actual fracturing fluid, and physical properties of an actual proppant, wherein the solid phase and the liquid phase are used for analog simulation in the fracture-wellbore simulation model; setting simulation conditions of the fracture-wellbore simulation model; determining a time step of a discrete element method and a time step of a computational fluid dynamics method according to the properties of the solid phase; performing coupled calculation of the computational fluid dynamics method and the discrete element method in the fracture-wellbore simulation model according to the properties of the liquid phase and the properties of the solid phase based on a liquid phase control equation, a solid phase control equation, a liquid phase to solid phase force equation, a solid phase to liquid phase force equation, the time step of the discrete element method, and the time step of the computational fluid dynamics method, and performing analog simulation on the solid phase and the liquid phase to realize proppant flowback numerical simulation of the fracture-wellbore system; wherein the liquid phase control equation comprises a continuity equation, a momentum equation, and a turbulence model, and the turbulence model comprises a turbulent kinetic energy equation and a dissipation rate equation; the continuity equation is shown in the following formula (7): （7） In formula (7), is the volume fraction of the liquid phase in the computational grid; is the density of the liquid phase; is the time of flight; is the flow velocity of the liquid phase; the momentum equation is shown in the following formula (8): （8） In formula (8), P is the fluid pressure; is the fluid viscosity; is the gravitational acceleration; is the volume force of the interaction of the solid and liquid phases; the turbulent kinetic energy equation is shown in the following formula (11): （11） in formula (11), is the liquid volume fraction; is the liquid density; is the turbulent kinetic energy of the liquid phase; is the dissipation rate of the turbulent kinetic energy; denotes the liquid velocity vector; is the liquid viscosity coefficient; is the dimensionless Prandtl number corresponding to the turbulent kinetic energy; is the source term of the turbulent kinetic energy; is the solid exchange coefficient; the dissipation rate equation is shown in the following formula (12): （12） In formula (12), is the viscosity coefficient of the continuous phase; is the dimensionless Prandtl number of the turbulent energy; C 1ε , C 2ε are empirical constants; is the exchange coefficient of the liquid phase.

2. The method of claim 1, wherein, the determination of the time step of the discrete element method and the time step of the computational fluid dynamics method according to the properties of the solid phase comprises: determining the time step of the discrete element method according to a solid phase particle diameter, a particle Poisson ratio, and a particle shear modulus of the solid phase; determining the time step of the computational fluid dynamics method according to the time step of the discrete element method.

3. The method of claim 1, wherein, the solid phase control equation comprises a solid phase particle motion equation and a solid phase particle contact equation; the solid phase particle contact equation is shown in the following formula (13): （13） In formula (13), solid phase particles j acting on the solid phase particles i ; , , , represent the tangential elastic force, the tangential damping, the normal elastic force, and the normal damping of the solid phase particles j to the solid phase particles i , respectively. Solid phase particles j The normal elastic force of the solid phase particles i is obtained by the following equation (14): （14） wherein , ; In equation (14), Y The Young's modulus of the solid particles; R The effective radius of the solid particles is denoted as . the solid phase particle motion equation comprises a translation equation and a rotation equation; n solid particles j With solid particles i The normal spacing; , They are solid particles i solid particles j Young's modulus; , They are solid particles i solid particles j Poisson's ratio; , They are solid particles i solid particles j The particle radius; Solid phase particles j The normal damping of the solid phase particles i is obtained by the following equation (15): （15） wherein , , ; In formula (15), is the equivalent mass; is the solid phase particle i is the normal component of the relative velocity of the solid phase particle j ; , is the mass of the solid phase particle i is the mass of the solid phase particle j ; is the restitution coefficient of the collision; S n is the normal stiffness; Solid phase particles j The tangential elastic force of the solid phase particles i is obtained by the following equation (16): （16） wherein , ; In formula (16), represents the tangential stiffness, G is the equivalent shear modulus, represents the solid phase particles j and the tangential overlap amount of the solid phase particles i ; Solid phase particles j The tangential damping of the solid phase particles i is obtained by the following equation (17): （17） In formula (17), v t are solid phase particles i with the solid phase particles j a tangential component of the relative velocity of the solid phase particles 4. The method of claim 3, wherein, the translation equation is shown in the following formula (18): the rotation equation is shown in the following formula (19): （18） In formula (18), m i is the mass of the solid phase particles i ; is the linear velocity of the solid phase particles ; is the force of the liquid phase on the solid phase particles i ; is the normal elastic force of the solid phase particles j or the wall surface on the solid phase particles i ; is the tangential elastic force of the solid phase particles j or the wall surface on the solid phase particles i ; is the acceleration of gravity; is the total contact times of the solid phase particles i with other particles and the wall surface; the liquid phase to solid phase force equation comprises a solid phase particle liquid phase buoyancy model and a solid phase particle liquid phase drag model; （19） In formula (19), is the moment of inertia of the solid particles i is the moment of inertia of the solid particles the solid phase particle liquid phase buoyancy model is shown in the following formula (20): i angular velocity of the solid phase particles i of the solid phase particles M t,ij torque on the particles due to the tangential force j of the solid phase particles i torque due to the tangential force M r,ij torque on the particles due to the rolling friction j of the solid phase particles i torque due to the rolling friction i total number of contacts of the solid phase particles with other solid phase particles and solid walls 5. The method of claim 1, wherein, the solid phase particle liquid phase drag model is shown in the following formula (21): the solid phase to liquid phase force equation comprises a liquid phase solid phase particle resistance model; （20） in formula (20), is the buoyancy force on the solid phase particles from the liquid phase, is the particle size of the solid phase particles; is the density of the liquid phase; is the acceleration due to gravity.

6. The method of claim 5, wherein, the liquid phase solid phase particle resistance model is shown in the following formula (22): （21） wherein , ; In formula (21), is the drag force of the liquid phase on the solid phase particle; is the drag coefficient of the liquid phase on the solid phase; is the volume fraction of the solid phase particle in the calculation grid; is the volume fraction of the liquid phase particle in the calculation grid; is the velocity of the liquid phase at the position of the solid phase particle; is the velocity of the solid phase particle; is the particle size of the solid phase particle; is the Reynolds number of the solid phase particle.

7. The method of claim 1, wherein, The method comprises the following steps: a setting module is configured to determine a fracture-wellbore simulation model, properties of a liquid phase, and properties of a solid phase according to an actual fracture-wellbore system, physical properties of an actual fracturing fluid, and physical properties of an actual proppant, wherein the solid phase and the liquid phase are used for analog simulation in the fracture-wellbore simulation model; （22） In formula (22), is the resistance of the liquid phase to the solid phase particles; is the force exerted by the liquid phase on the solid phase particles; Δ V to calculate the volume of the grid; n to calculate the number of solid particles contained in the grid.

8. A fracture-wellbore system proppant flowback numerical simulation apparatus, characterized in that, setting simulation conditions of the fracture-wellbore simulation model; and determining a time step of a discrete element method and a time step of a computational fluid dynamics method according to the properties of the solid phase. ​ ​ The computing simulation module is configured to perform coupled computation of the computational fluid dynamics method and the discrete element method in the fracture-wellbore simulation model based on a liquid phase control equation, a solid phase control equation, a liquid phase to solid phase force equation, a solid phase to liquid phase force equation, a time step of the discrete element method, and a time step of the computational fluid dynamics method, simulate the solid phase and the liquid phase according to properties of the liquid phase and properties of the solid phase, and achieve numerical simulation of proppant flowback of a fracture-wellbore system. The liquid phase control equation comprises a continuity equation, a momentum equation, and a turbulence model, and the turbulence model comprises a turbulent kinetic energy equation and a dissipation rate equation. The continuity equation is as follows: wherein is the volume fraction of the computational grid occupied by the liquid phase; is the density of the liquid phase; is the transit time; is the flow velocity of the liquid phase; The momentum equation is as follows: wherein P is the fluid pressure; is the fluid viscosity; is the gravitational acceleration; is the volume force of the interaction of the solid and liquid phases; The turbulent kinetic energy equation is as follows: wherein is the liquid volume fraction; is the liquid density; is the turbulent kinetic energy of the liquid phase; is the dissipation rate of the turbulent kinetic energy; denotes the liquid phase velocity vector; is the liquid viscosity coefficient; is the dimensionless Prandtl number corresponding to the turbulent kinetic energy; is the source term of the turbulent kinetic energy; is the solid phase exchange coefficient; The dissipation rate equation is as follows: wherein is the viscosity coefficient of the continuous phase; is the dimensionless Prandtl number of the turbulent kinetic energy; C 1ε , C 2ε are empirical constants; is the liquid phase exchange coefficient.

9. An apparatus for performing a fracture-wellbore system proppant flowback numerical simulation method, comprising: The method comprises the following steps: a processor; a memory for storing processor-executable instructions; The processor executes the executable instructions to implement the method of any one of claims 1 to 7.

10. A non-transitory computer readable storage medium, comprising: The computer program or instructions are stored in a computer readable storage medium, and when the computer program or instructions are executed, the method of any one of claims 1 to 7 is implemented.

Images