Method, device and equipment for simulating migration of proppant and medium

Abstract

The invention discloses a method, a device, equipment and a medium for simulating migration of a proppant, and relates to the technical field of analogue simulation. Simulation is carried out through fluid dynamics and a discrete element analysis model, specifically, fluid parameters are obtained, a calculation model is constructed, and constraints are set; parameters are input for fluid and particle stress analysis, and particle displacement and grid updating are calculated; contact is judged through the particle collision sequence, boundary force is calculated, and the particle motion state is updated; and cyclically feeding back the flow field parameters to particle stress analysis until the number of iterations is met so as to obtain the void ratio, the flow velocity, the position and the momentum of the particles. According to the invention, the motion among particles can be captured and accurately calculated in a high-concentration particle flow or particle dominated collision system.

Description

Technical Field

[0001] The present invention relates to the field of simulation technology, and in particular to a method, device, equipment and medium for simulating proppant migration. Background Art

[0002] Fracturing is a key technology for the effective utilization of conventional and unconventional oil and gas resources. By pumping fracturing fluid, the reservoir rock is fractured to form hydraulic fractures and proppants are injected to prevent the fractures from closing. The sand banks formed by the accumulation of proppants in the fractures form diversion channels for the "outflow" of reservoir resources. Therefore, the accumulation morphology of the sand banks determines the aperture, shape, volumetric conductivity and "sustainability" of the fractures. Therefore, clarifying the migration mechanism of proppants in fractures and optimizing fracturing construction parameters are the keys to improving the fracturing production increase effect.

[0003] Currently, researchers at home and abroad are primarily studying the proppant migration and placement mechanisms through physical experimental simulation, numerical simulation, and theoretical derivation. Numerical simulation technology has attracted considerable attention due to its low cost, high efficiency, and speed. Among numerical simulation technologies, the main models used for solid-liquid two-phase flow are the Euler-Euler model and the Euler-Lagrange model. The Euler-Euler model treats the particle phase as a pseudo-fluid, ignoring the interaction forces between particles and fluid, particles and particles, and particles and walls. The Euler-Lagrange model, on the other hand, treats the particle phase as a number of discrete units and solves the equations of motion in the Lagrangian coordinate system to obtain data such as the position and velocity of the particles at different times, which can more realistically reflect the movement of particles under the action of the fluid.

[0004] Discrete Phase Model (DPM), Lattice Boltzmann Method (LBM), Multiphase Particle-in-Cell Method (MP-PIC), Direct Numerical Simulation-Discrete Element Method (DNS-DEM), Computational Fluid Dynamics-Discrete Element Method (CFD-DEM), etc. Method. Among them, DPM is a sub-function embedded in Fluent, which is suitable for the simulation of sparse particle flow. It ignores the collision between particles and cannot observe the accumulation morphology of particles; LBM is based on the mesoscopic scale and can capture discontinuous effects in the microscopic scale, such as slip flow and rarefied gas flow; the DNS-DEM method is suitable for systems where the particle scale is close to the fluid characteristic scale (microscopic mechanism), with a huge amount of calculation and high difficulty; the MP-PIC method is suitable for large-scale sparse particle flow systems, solid-liquid two-phase flows where fluid dynamics dominates and particle collisions are not intense. This method characterizes the behavior of a large number of particles through particle groups, ignoring most of the interactions between particles. The calculation accuracy is low and it cannot describe the motion trajectory of a single particle and the particle migration mechanism. In other words, the existing technology cannot capture and accurately calculate the motion between particles in systems with high-concentration particle flows or particle-dominated collisions. Summary of the Invention

[0005] Embodiments of the present invention provide a proppant migration simulation method, apparatus, device, and medium, which can address the shortcomings of the prior art in being unable to capture and accurately calculate the motion between particles in systems with high-concentration particle flows or particle-dominated collisions.

[0006] An embodiment of the present invention provides a method for simulating proppant migration, comprising the following steps: Obtain the fluid dynamic parameters of the proppant, including: fluid domain, grid parameters and fluid physical properties; Construct a fluid dynamics model to describe the proppant migration and flow characteristics, and set constraints including boundary parameters, initial conditions, time step, and number of iterations. Input the proppant's fluid dynamics parameters into the model, and use a fluid dynamics CFD solver to obtain proppant fluid information, including flow parameters and mechanical parameter vectors within the flow field. Based on the proppant's fluid dynamics parameters and the number of proppant particles, a discrete element analysis model is used to perform particle force analysis. The discrete element method (DEM) solver is used to solve the particle force analysis results and proppant fluid information to obtain proppant particle information including particle position, particle velocity, and particle momentum. The fluid-solid coupling model in the CFD solver is used to solve the proppant particle information and the proppant fluid information, and the flow field distribution and force distribution during the simulated proppant migration process are updated. The flow field distribution and force distribution are fed back to the DEM solver to determine the error until the preset simulation time is met, completing the simulation of proppant migration.

[0007] Furthermore, the proppant fluid information specifically includes: flow velocity, flow rate and momentum as flow parameters, and pressure, gravity and drag as mechanical parameter vectors in the flow field, which together constitute the proppant fluid information.

[0008] Furthermore, the construction of a fluid dynamics calculation model for describing proppant migration and flow characteristics specifically comprises the following steps: Based on the flow characteristics of proppant migration, the Standard k-ε model is constructed; The Standard k-ε model is constructed based on the eddy viscosity assumption by simultaneously solving the continuity equation, momentum equation, kinetic energy k equation, turbulent dissipation rate ε equation, and turbulent fluid viscosity equation. The velocity field, pressure field, kinetic energy k of the fluid phase and turbulent dissipation rate ε of the Standard k-ε model are solved iteratively through discretization equations, and the Reynolds stress is updated according to the turbulent viscosity to obtain the fluid dynamics calculation model.

[0009] Furthermore, the particle force analysis using the discrete element analysis model specifically includes the following steps: The particles are considered to have two forms of motion: translation and rotation. In the discrete element simulation, Newton's second law is used to obtain the acceleration and velocity of the particles. Use Newton's second law to obtain the equation of motion of the particle at a certain moment: ; ; in, v is the translational velocity of the particle; m is the mass of the particle; t It’s time; F is the net force acting on the particle; I is the moment of inertia; is the angular velocity; M is the resultant moment acting on the particle; Assuming that the force on the particle remains constant within a time step, the velocity and displacement of the particle can be obtained by numerically integrating the acceleration within a time step: ; ; in, Δt is the time step, u is the displacement caused by translation, is the displacement caused by rotation; The obtained particle acceleration, particle velocity and displacement are the displacement increment of each particle in the proppant and the grid parameters of the particle unit are updated according to the increment conditions to complete the force analysis of the particle.

[0010] Furthermore, the discrete element method (DEM) solver is used to solve the particle force analysis results and the proppant fluid information to obtain the proppant particle information including the particle position, particle velocity and particle momentum. The specific steps include: Based on the particle force analysis results and proppant fluid information, the particle collision sequence is obtained using a DEM solver; Conduct contact judgment between particles and obtain boundary contact forces between particles, including: The normal overlap between particles d n The positive or negative value is used to determine whether there is contact collision between particles and walls, or between particles. Obtain the particle contact point, the particle center coordinates and radius, and then obtain the three-dimensional information of the particle. According to the three-dimensional information of the particle, obtain the normal overlap of the particle d n , the formula is: ; Where, R i 、R j For particles i and particles j radius; r i 、r j Particles i and particles j The sphere center position vector; when When , it indicates that the two particles are in contact with each other; when When , the two particles have no contact; When determining whether a particle contacts a wall, the wall is considered to be a stationary spherical particle with infinite radius, and the normal overlap is used. d nThe calculation formula determines whether contact occurs; The boundary contact force between particles is obtained when two particles are in contact, and the normal contact force and tangential contact force are obtained: Obtain the normal contact force through the Hertz contact force model, the normal contact force of the proppant particles F cn,ij is the normal overlap d n The function is expressed as follows: ; in, k n,ij is the normal stiffness coefficient; d n,ij is the normal displacement between the colliding particles; n is the normal unit vector between contact particles; is the normal dissipation coefficient; is the normal component of the relative velocity between colliding particles; is the particle Young's modulus, Pa; is the equivalent particle radius; According to the physical parameters of the particles, the normal stiffness coefficient and normal dissipation coefficient are obtained as follows: ; ; ; ; Where, Y i 、 v i 、 R i and Y j 、 v j 、 R j Particles i Young's modulus, Poisson's ratio and radius of the particles, as well as j Young's modulus, Poisson's ratio and radius; The normal elastic damping of the particle is obtained, and its expression is: ; ; Where, m i 、 m j Particles i , particles j Mass, kg;m* is the equivalent mass, kg; is the normal component of the particle relative velocity; e is the collision restitution coefficient; s n 、 β By the restitution coefficient e Decide; Tangential contact force Obtained through Mindlin-Deresiewicz contact theory, using the tangential overlap d t and tangential stiffness S t , the formula is as follows: ; ; ; Where G is the equivalent shear modulus; The expression of particle tangential elastic damping is: ; in, v t is the tangential component of the relative velocity of the particles; thus, the proppant particles j Acting on proppant particles i The contact force on is expressed as follows: ; in, F t,ij 、 F d t,ij 、 F n,ij and F d n,ij They are tangential elastic force, tangential damping, normal elastic force and normal damping respectively; The setting of the particle motion time step needs to be based on the Rayleigh time step, which refers to the time required for a shear wave to pass through a solid particle, and the maximum time step of the particle motion cannot exceed its Rayleigh time step. , the calculation formula is as follows: ; in, R s is the particle diameter; is the particle density; G s is the particle shear modulus, v s It is dimensionless; The particle velocity and position are updated according to the boundary contact force, and proppant particle information including particle position, particle velocity and particle momentum is obtained.

[0011] An embodiment of the present invention provides a proppant migration simulation device, comprising: Parameter setting module, used to obtain the fluid dynamic parameters of the proppant, including: fluid domain, grid parameters and fluid physical properties parameters; The migration simulation module is used to construct a fluid dynamics calculation model to describe the migration and flow characteristics of the proppant, and set constraints including boundary parameters, initial conditions, time step, and number of iterations. The fluid dynamics parameters of the proppant are input into the fluid dynamics calculation model, and a fluid dynamics CFD solver is used to obtain proppant fluid information including flow parameters and mechanical parameter vectors within the flow field. Based on the proppant fluid dynamics parameters and the number of proppant particles, a discrete element analysis model is used to perform particle force analysis. The discrete element method (DEM) solver is used to solve the particle force analysis results and proppant fluid information to obtain proppant particle information including particle position, particle velocity, and particle momentum. The iterative feedback module is used to solve the proppant particle information and proppant fluid information using the fluid-solid coupling model in the CFD solver, and update the flow field distribution and force distribution during the simulated proppant migration process; the flow field distribution and force distribution are fed back to the DEM solver to determine the error until the preset simulation time is met, completing the simulation of proppant migration.

[0012] An embodiment of the present invention provides a computer device, comprising: a memory and a processor; the memory stores a computer program, and the processor implements the above-mentioned proppant migration simulation method when executing the computer program.

[0013] An embodiment of the present invention provides a computer-readable storage medium storing a computer program. When the computer program is executed by a processor, the computer program implements the above-mentioned method for simulating proppant migration.

[0014] The embodiments of the present invention provide a proppant migration simulation method, device, equipment, and medium. Compared with the prior art, the beneficial effects thereof are as follows: The present invention accurately simulates the proppant migration process through the coupling of fluid dynamics CFD and discrete element method DEM. Through the specific analysis of the fluid and particles of the proppant, the flow field distribution and force distribution during the proppant migration process are fully considered, and the information of each particle in the proppant can be accurately captured. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Figure 1 A schematic structural diagram of an embodiment of the present invention. DETAILED DESCRIPTION

[0016] To make the above-mentioned objects, features, and advantages of the present invention more readily apparent, specific embodiments of the present invention are described in detail below with reference to the accompanying drawings. The following description sets forth numerous specific details to facilitate a full understanding of the present invention. However, the present invention can be implemented in many other ways than those described herein, and those skilled in the art may make similar modifications without departing from the scope of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0017] See also Figure 1 , an embodiment of the present invention provides a method for simulating proppant migration, comprising the following steps: Step 1: Obtain the fluid dynamic parameters of the proppant, including: fluid domain, grid parameters and fluid physical properties parameters.

[0018] Step 2: Construct a fluid dynamics model to describe the proppant migration and flow characteristics, and set constraints including boundary parameters, initial conditions, time step, and number of iterations. Input the proppant's fluid dynamics parameters into the model, and use a fluid dynamics CFD solver to obtain proppant fluid information, including flow parameters and flow field mechanical parameter vectors. This proppant fluid information includes flow parameters such as velocity, flow rate, and momentum, and flow field mechanical parameter vectors such as pressure, gravity, and drag.

[0019] Among them, the Standard k-ε model is constructed based on the flow characteristics of proppant migration; the Standard k-ε model is constructed by simultaneously solving the continuity equation, momentum equation, kinetic energy k equation, turbulent dissipation rate ε equation and turbulent fluid viscosity equation, and is based on the eddy viscosity assumption; the velocity field, pressure field, kinetic energy k and turbulent dissipation rate ε of the Standard k-ε model are iteratively solved by discretizing the equations, and the Reynolds stress is updated according to the turbulent viscosity to obtain the fluid dynamics calculation model.

[0020] Step 3: Based on the proppant fluid dynamics parameters and the number of proppant particles, a discrete element analysis model is used to perform particle force analysis. The discrete element method (DEM) solver is used to solve the particle force analysis results and proppant fluid information to obtain proppant particle information including particle position, particle velocity, and particle momentum.

[0021] Among them, the particle stress analysis using discrete element analysis model includes: The particles are considered to have two forms of motion: translation and rotation. In the discrete element simulation, Newton's second law is used to obtain the acceleration and velocity of the particles.

[0022] Use Newton's second law to obtain the equation of motion of the particle at a certain moment: .

[0023] .

[0024] in, v is the translational velocity of the particle; m is the mass of the particle; t It’s time; F is the net force acting on the particle; I is the moment of inertia; is the angular velocity; M is the resultant moment acting on the particle.

[0025] Assuming that the force on the particle remains constant within a time step, the velocity and displacement of the particle can be obtained by numerically integrating the acceleration within a time step: .

[0026] .

[0027] in, Δt is the time step, u is the displacement caused by translation, and is the displacement caused by rotation.

[0028] The obtained particle acceleration, particle velocity and displacement are the displacement increment of each particle in the proppant and the grid parameters of the particle unit are updated according to the increment conditions to complete the force analysis of the particle.

[0029] Step 4: Use the fluid-solid coupling model in the CFD solver to solve the proppant particle information and proppant fluid information, and update the flow field distribution and force distribution during the simulated proppant migration process; feed the flow field distribution and force distribution back to the DEM solver to determine the error until the preset simulation time is met, completing the simulation of proppant migration.

[0030] 1. Liquid phase governing equation: The governing equation for liquid phase flow is the solid-liquid two-phase coupled Navier-Stokes equation, and the volume of the particle phase needs to be considered when calculating the liquid phase volume fraction.

[0031] (1) Mass conservation equation (continuity equation).

[0032] The increase in mass of a fluid element per unit time is equal to the net mass flowing into the element during the same time interval.

[0033] .

[0034] Where, is the liquid density, kg / m³; is the liquid volume fraction, %; t is the movement time, s; is the liquid flow velocity, m / s.

[0035] (2) Momentum conservation equation.

[0036] The change in momentum of the fluid in a microelement per unit time is equal to the superposition of the forces exerted on the microelement during the flow process, which is essentially Newton's second law.

[0037] .

[0038] .

[0039] Where, P is the fluid pressure, Pa; is the liquid phase viscous stress tensor; is the fluid viscosity, Pa·s; g is the acceleration due to gravity, m / s 2 ; is the solid-liquid interaction force, N.

[0040] (3) Liquid-solid phase volume fraction: The volume fraction of the particle phase in the computational grid unit is obtained using the falling point method. The volume fraction of the solid phase and the volume fraction of the fluid are and solid-liquid interaction volume forces It can be obtained by the following formula:

[0041] .

[0042] .

[0043] .

[0044] .

[0045] Where N is the total number of points in the grid unit, n p is the total number of sample points of particles in the calculation grid unit; n is the total number of solid phase particles in the calculation unit body; To calculate the volume of the control unit, m³; is the length of the cell in the x, y, and z directions of the three-dimensional space, m; is the volume of solid particles in the calculation unit, m 3 ; is the viscous resistance of the fracturing fluid on solid particle i during the flow process, N.

[0046] 2. Solid phase motion equation: (1) Translational equation: .

[0047] (2) Rotation equation: .

[0048] Where, v i is the linear velocity of particle i, m / s; ω i represents the angular velocity of particle i, rad / s; z is the total number of times particle i contacts other particles and solid walls; F pf is the force exerted by the fluid on particle i, N; F n is the normal elastic force exerted by particle j or the wall on particle i, N; F t is the tangential elastic force exerted by particle j or the wall on particle i, N; m i g is the gravity of particle i, N; M t,ij 、M r,ij are the moments caused by tangential force and rolling friction, respectively, in N·m.

[0049] 3. Turbulent flow model Turbulence numerical models include: single equation Spalart-Allmaras model; k- Model (standard k- Model, pressure correction k- model, Reynolds pressure model, large eddy simulation model); k- Model (Standard k- Model, RNGk- Model, Realizable (with swirl correction) k- The most commonly used standard k-ε model is selected in this simulation. It is a typical empirical model. Due to its high convergence and calculation accuracy, its applicability to high Reynolds numbers, its wide range of applicability, and its economic rationality, combined with the high fluid inlet velocity and drastic flow field changes in this paper, the standard k-ε model is adopted as the turbulent flow model.

[0050] (1) Turbulent kinetic energy k equation: .

[0051] (2) Turbulent dissipation rate equation: .

[0052] Where k is the kinetic energy of the fluid phase, m 2 / s 2 ; Turbulence dissipation rate, m 2 / s 3 ; G kis the generation term of turbulent kinetic energy, kg / m·s 3 ; k is the Prandtl number corresponding to the turbulent kinetic energy, dimensionless, and is taken as 1.0; is the Prandtl number corresponding to the turbulent dissipation rate, dimensionless, and is taken as 1.3; S k 、S is the flow exchange term between liquid and solid phases, kg / m·s 3 ; 、 is an empirical constant, dimensionless, and is taken 、 .

[0053] 4. Particle contact model: (1) Particle contact determination.

[0054] The normal overlap between particles δ n The positive or negative value is used to determine whether there is contact collision between particles and walls, or between particles.

[0055] Normal overlap δ n The calculation formula is: .

[0056] Where R i 、R j is the radius of particles i and j, m; r i 、r j are the spherical center position vectors of particles i and j respectively. When , it indicates that the two particles are in contact with each other; when When , the two particles are not in contact. When determining whether a particle is in contact with a wall, the wall is considered to be a stationary spherical particle with an infinite radius, and then the above formula is substituted to determine whether contact occurs.

[0057] (2) Particle contact model.

[0058] Based on the soft-sphere model, which can handle three-body and higher particle collisions, particle contact models include the Hertz-Mindlin no-slip contact model, the Hertz-Mindlin adhesive contact model, and the linear elastic contact model. This paper uses the Hertz-Mindlin no-slip contact model to calculate particle-particle and particle-wall contact forces. It considers the damping coefficients of the normal and tangential forces and simulates particle collisions by introducing springs, dampers, and sliders.

[0059] When particles i and j, with radii Ri and Rj, collide, the contact point is C, and the spring and damper represent the normal contact force. Here, δn,ij represents the normal displacement between the colliding particles, and δt,ij represents the tangential displacement between the colliding particles, represented by the spring, damper, and slider. When particles collide, the normal and tangential contact forces are first calculated and then superimposed to obtain the net force and net moment of the particle contact, as shown in the equation below. When particles i and k are in contact simultaneously, the inter-particle forces and torques acting on them are vectorially superimposed.

[0060] .

[0061] Where, is the contact force between particles i and j, N; is the normal contact force between particles i and j, N; is the tangential contact force between particles i and j, N.

[0062] The normal contact force of the proppant particles is calculated using the Hertz contact force model. The normal contact force is a function of the normal overlap δn: .

[0063] Where, is the normal stiffness coefficient; is the normal displacement between colliding particles, m; is the normal unit vector between contact particles; is the normal dissipation coefficient; is the normal component of the relative velocity between colliding particles; is the particle Young's modulus, Pa; is the equivalent particle radius, m.

[0064] The normal stiffness coefficient and normal dissipation coefficient are calculated based on the physical parameters of the particles, such as Young's modulus and Poisson's ratio, as shown in the following formula: .

[0065] .

[0066] .

[0067] .

[0068] Where, Y i 、 v i 、 R i and Y j 、v j 、 R j Particles i , particles j Young's modulus, Poisson's ratio (dimensionless) and radius. The particle normal elastic damping is calculated as:

[0069] .

[0070] .

[0071] Where, m i 、 m j Particles i , particles j Mass, kg; m* is the equivalent mass, kg; is the normal component of the particle relative velocity; e is the collision restitution coefficient; s n 、 β By the restitution coefficient e (dimensionless) decision.

[0072] Tangential contact force Calculated by Mindlin-Deresiewicz contact theory, depending on the amount of tangential overlap d t and tangential stiffness S t : .

[0073] .

[0074] .

[0075] Where G is the equivalent shear modulus, Pa.

[0076] The expression of particle tangential elastic damping is: .

[0077] Where, v t is the tangential component of the particle relative velocity. j Acting on proppant particles i The contact force on can be expressed as:

[0078] .

[0079] Where,F t,ij 、 F d t,ij 、 F n,ij 、 F d n,ij are the tangential elastic force, tangential damping, normal elastic force and normal damping.

[0080] The setting of the particle motion time step needs to be based on the Rayleigh time step, which refers to the time required for a shear wave to pass through a solid particle, and the maximum time step of the particle motion cannot exceed its Rayleigh time step. , the calculation formula is as follows: .

[0081] in, R s is the particle diameter; is the particle density; G s is the particle shear modulus, v s It is dimensionless.

[0082] 5. Solid-liquid two-phase coupling equation.

[0083] During solid-liquid two-phase coupled flow, particles are subject to various liquid forces such as drag, buoyancy, pressure gradient force, gravity, and virtual mass force.

[0084] (1) The buoyancy of the fluid acting on the particle is.

[0085] .

[0086] Where, d p is the particle size, m; r f is the liquid density, kg / m 3 .

[0087] (2) Drag force.

[0088] The drag force on proppant particles when flowing with fluid is composed of pressure difference resistance and friction resistance. It is the main reason for the momentum exchange between the fluid and particles. Currently, common drag models include: Di-Felice model, Wen & Yu model and Gidsapow model. You can choose according to your needs. The drag force can be uniformly expressed as follows: .

[0089] .

[0090] Where, D p is the particle drag coefficient, dimensionless; S g is the specific gravity of the particle, dimensionless; C d is the drag coefficient, dimensionless; u f is the fluid velocity, m / s; u p is the particle velocity, m / s; r is the particle radius, m.

[0091] (3) The force exerted by particles on fluid: Based on the CFD-DEM coupling and the fluid drag force on the particles, it is necessary to calculate the momentum exchange source term between the particles and the fluid F p In order to achieve the coupling between the two, the momentum exchange source term represents the resistance of the particle to the fluid (equivalent to the reaction force of the fluid on the particle drag), and its specific expression is: .

[0092] Where, F d,i is the force exerted by the fluid on the particle, N; To calculate the unit cell volume, m3; Δx, Δy, Δz are the lengths of the calculation grid in the x, y, and z directions, m; n is the number of particles contained in the calculation grid, which is dimensionless.

[0093] (4) Pressure gradient force: If the pressure gradient along the flow direction is P / l, the pressure gradient force acting on a single particle is: .

[0094] 6.CFD-DEM coupling process.

[0095] The CFD solver calculates the velocity, flow rate, momentum and force information of the liquid phase fluid at the initial inlet according to the initial boundary conditions, and transmits the fluid information to the DEM solver in real time. The DEM solver solves the particle velocity, force and position change information according to the fluid force and flow rate and particle characteristics, and transmits the particle information to the CFD solver in real time. The CFD solver uses the physical properties, velocity, position and other parameters of the particles to solve the porosity, velocity, position and momentum based on the fluid-solid coupling model and returns it to the DEM solver to judge the error, thus forming an iterative closed loop until the calculation is completed.

[0096] The present invention provides a proppant migration simulation method based on CFD-DEM coupling, which is mainly divided into the following steps.

[0097] 1. Determine the crack model, construct the calculation domain, and establish the "entrance, exit, and wall".

[0098] 2. Define the computational domain boundary conditions (inlet, outlet, wall) in the CFD solver.

[0099] 3. Construct particles in the DEM solver and set the physical and mechanical properties of the particles according to actual conditions.

[0100] 4. Couple the CFD solver with the DEM solver and start the proppant transport simulation based on CFD-DEM coupling until the preset simulation time is reached.

[0101] This invention, based on bidirectional CFD-DEM coupling, fully considers the interactions between particles, between particles and walls, and between particles and fluids, enabling realistic simulation of proppant particle migration within fractures under the influence of fluids. Furthermore, the invention discloses methods for constructing different fracture models and applying them to different working conditions. Compared to DPM, LBM, MP-PIC, and DNS-DEM, CFD-DEM coupling offers the following advantages:

[0102] 1. The fluid-particle bidirectional coupling is fully considered.

[0103] DPM: It is generally assumed that the feedback force of particles on the fluid is small (one-way coupling), which cannot accurately describe the strong feedback effect of particles on the fluid. LBM: Although LBM has high accuracy in fluid simulation, its coupling with particles is complex. Especially in areas with high particle concentration, the computational efficiency of fluid-particle coupling is low. MP-PIC: It uses a statistical model to handle the feedback force of particles on the fluid, which is less accurate and cannot capture the microscopic behavior of particles. DNS-DEM: It can accurately describe fluid-particle coupling, but the computational cost is extremely high, making it difficult to apply to large-scale proppant migration simulations. CFD-DEM coupling can accurately describe the forces acting on particles by the fluid (drag, buoyancy, pressure gradient, etc.) as well as the feedback forces of particles on the fluid (such as the momentum source term).

[0104] 2. The interaction between proppant particles and the particle deposition behavior are described more accurately.

[0105] DPM: Direct interactions between particles are usually ignored, or processed through simple statistical models, which cannot accurately describe the collision and friction between particles. LBM: LBM itself does not deal with the interaction between particles, and requires the additional introduction of DEM or other methods, which increases the computational complexity. MP-PIC: The use of statistical models to deal with the interaction between particles has low accuracy and cannot capture the microscopic behavior of particles. DNS-DEM: It can accurately describe the interaction between particles, but the computational cost is extremely high. The CFD-DEM coupling can calculate the contact forces (such as collision force, friction, and rolling resistance) between particles pair by pair, accurately describe the microscopic dynamic behavior of particles, and is suitable for simulating particle accumulation, deposition, and blockage behavior in high-concentration proppant flows.

[0106] An embodiment of the present invention provides a proppant migration simulation device, comprising: The parameter setting module is used to obtain the fluid dynamic parameters of the proppant, including: fluid domain, grid parameters and fluid physical properties parameters.

[0107] The migration simulation module is used to construct a fluid dynamics calculation model to describe the migration and flow characteristics of the proppant, and set constraints including boundary parameters, initial conditions, time step and number of iterations; the fluid dynamics parameters of the proppant are input into the fluid dynamics calculation model, and the fluid dynamics CFD solver is used to obtain the proppant fluid information including flow parameters and mechanical parameter vectors in the flow field; based on the fluid dynamics parameters of the proppant and the number of proppant particles, the discrete element analysis model is used to perform particle force analysis; the discrete element method DEM solver is used to solve the particle force analysis results and the proppant fluid information to obtain the proppant particle information including particle position, particle velocity and particle momentum.

[0108] The iterative feedback module is used to solve the proppant particle information and proppant fluid information using the fluid-solid coupling model in the CFD solver, and update the flow field distribution and force distribution during the simulated proppant migration process; the flow field distribution and force distribution are fed back to the DEM solver to determine the error until the preset simulation time is met, completing the simulation of proppant migration.

[0109] An embodiment of the present invention provides a computer device, comprising: a memory and a processor; the memory stores a computer program, and the processor implements the steps of a proppant migration simulation method when executing the computer program.

[0110] An embodiment of the present invention provides a computer-readable storage medium storing a computer program. When the computer program is executed by a processor, the computer program implements the steps of a proppant migration simulation method.

[0111] A specific embodiment is as follows: This embodiment discloses a method for simulating proppant migration, and the specific steps are as follows: In this example, Space Claim is used as the modeling software, Meshing is used as the meshing software, Fluent is used as the CFD solver, and Rocky DEM is used as the DEM solver.

[0112] S1. Crack model construction: Open ANSYS Workbench, select the fluid flow module (Fluent), right-click the geometry structure and select Space Claim to model and design different crack structures (first draw a two-dimensional model of the crack with a straight line, and then stretch it to draw the desired crack structure).

[0113] S2. Meshing: Close SpaceClaim, click Meshing software to import the established model, perform meshing and check the mesh quality.

[0114] S3. Export the divided mesh model and import it into Fluent (CFD solver). Set the corresponding gravity direction, turbulence model, fluid properties, crack inlet and outlet, crack wall conditions, time step and other parameters in Fluent. After setting, export the .cas file.

[0115] S4. Open the discrete element software Rocky DEM, also known as the DEM solver. First, compare the proppant stacking state obtained from the cylinder simulation experiment in the "virtual test" method of stacking angle with the proppant stacking state obtained from the indoor cylinder experiment, mainly comparing the repose angle of the two. By continuously adjusting the three physical parameters of the collision recovery coefficient, static friction coefficient, and kinetic friction coefficient within the proppant particles and between the proppant particles and the cylinder wall, the experimental results in the cylinder laboratory and the simulation results are more closely matched, and the optimal parameters are calibrated.

[0116] S5. Calibrate the collision recovery coefficient, static friction coefficient, and kinetic friction coefficient of the proppant particles based on the indoor cylinder test and the virtual cylinder test.

[0117] S6. After parameter calibration, open RockyDEM again and set the gravity direction, normal force model, and tangential force model in Physics. This case is based on the Hertzian Spring Dashpot model and the Mindlin Deresiewicz model. After setting these, select "CFD Coupled Particle Statistics" in Modules to detect coupled data. After the calculation is complete, perform post-processing and import the data in "Geometries."cas file, select and delete inlet and outlet, then right-click "Geometries" and select "Createinlet" to create a particle inlet (multiple inlets can be set at the same time to simulate multiple injection hole conditions) click "Materials" to set the density, elastic modulus, and Poisson's ratio of the fracture and proppant; then click "MaterialInteractions" to input the calibrated collision recovery coefficient, static friction coefficient, and dynamic friction coefficient; right-click "Particle" and select "CreateParticle" to set the required proppant particle shape in "shape and size" (you can also scan the proppant particles by laser and then import the model), particle size, and particle ratio (multiple particles can be mixed); right-click "Inputs" and select "CreateContinuousInjection" (if multiple perforations are simulated, the number of ContinuousInjection creations must be the same as the number of perforations (inlets created in Geometries)), select the corresponding inlet in "EntryPoint" in ContinuousInjection, select the corresponding particles, set the appropriate "MassFlowRate" according to the proppant mass fraction (it can be said to be sand ratio, proppant concentration, etc.), and in "Ti Set the particle injection duration in "Me" and the proppant particle injection rate in "Entry". Click "CFDCoupling", select 2-way-Fluent in "CouplingMode", click "2-way-Fluent" below "CFDCoupling", click "particle" in the "Interactions" sub-option, set "DragLaw" (Drag force; the Gidaspow model is recommended for proppant transport simulation) and "VirtualMassLaw". Then click the "Fluent" sub-option and change "Serial" in "ExecutionMode" to "LocalParallel". Then, according to your computer configuration, set the "SolverProcess" used by Fluent during coupling and select whether to keep all data (keepallfiles). Click "Solver", set "SimulationDuration" and "OutputFrequencies", and check "ReleaseParticleswithoutOverlapCheck" to disable particle overlap when releasing particles to prevent particle overlap from causing calculation problems.

[0118] S7. After the above conditions are set, click "START" to start the CFD-DEM bidirectional coupling calculation.

[0119] The above-described embodiments merely illustrate several implementations of the present invention, and while their descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the patent. It should be noted that a person skilled in the art would be able to make numerous variations and improvements without departing from the spirit of the present invention, all of which fall within the scope of protection of the present invention. Therefore, the scope of protection of the patent for this invention shall be determined by the appended claims.

Claims

1. A method for simulating proppant migration, characterized in that: The following steps are involved: Obtain the fluid dynamic parameters of the proppant, including: fluid domain, grid parameters and fluid physical properties; Construct a fluid dynamics model to describe the proppant migration and flow characteristics, and set constraints including boundary parameters, initial conditions, time step, and number of iterations. Input the proppant's fluid dynamics parameters into the model, and use a fluid dynamics CFD solver to obtain proppant fluid information, including flow parameters and mechanical parameter vectors within the flow field. Based on the proppant's fluid dynamics parameters and the number of proppant particles, a discrete element analysis model is used to perform particle force analysis. The discrete element method (DEM) solver is used to solve the particle force analysis results and proppant fluid information to obtain proppant particle information including particle position, particle velocity, and particle momentum. The fluid-solid coupling model in the CFD solver is used to solve the proppant particle information and the proppant fluid information, and the flow field distribution and force distribution during the simulated proppant migration process are updated. The flow field distribution and force distribution are fed back to the DEM solver to determine the error until the preset simulation time is met, completing the simulation of proppant migration.

2. A method for simulating proppant migration according to claim 1, characterized in that: The proppant fluid information specifically includes: Velocity, flow rate and momentum as flow parameters, and pressure, gravity and drag as mechanical parameter vectors in the flow field, together constitute the proppant fluid information.

3. A method for simulating proppant migration according to claim 1, characterized in that: The construction of a fluid dynamics calculation model for describing proppant migration and flow characteristics specifically comprises the following steps: Based on the flow characteristics of proppant migration, the Standard k-ε model is constructed; The Standard k-ε model is constructed based on the eddy viscosity assumption by simultaneously solving the continuity equation, momentum equation, kinetic energy k equation, turbulent dissipation rate ε equation, and turbulent fluid viscosity equation. The velocity field, pressure field, kinetic energy k of the fluid phase and turbulent dissipation rate ε of the Standard k-ε model are solved iteratively through discretization equations, and the Reynolds stress is updated according to the turbulent viscosity to obtain the fluid dynamics calculation model.

4. A method for simulating proppant migration according to claim 1, characterized in that: The particle force analysis using the discrete element analysis model specifically includes the following steps: The particles are considered to have two forms of motion: translation and rotation. In the discrete element simulation, Newton's second law is used to obtain the acceleration and velocity of the particles. Use Newton's second law to obtain the equation of motion of the particle at a certain moment: ； ； in, v is the translational velocity of the particle; m is the mass of the particle; t It’s time; F is the net force acting on the particle; I is the moment of inertia; is the angular velocity; M is the resultant moment acting on the particle; Assuming that the force on the particle remains constant within a time step, the velocity and displacement of the particle can be obtained by numerically integrating the acceleration within a time step: ； ； in, Δt is the time step, u is the displacement caused by translation, is the displacement caused by rotation; The obtained particle acceleration, particle velocity and displacement are the displacement increment of each particle in the proppant and the grid parameters of the particle unit are updated according to the increment conditions to complete the force analysis of the particle.

5. A method for simulating proppant migration according to claim 1, characterized in that: The discrete element method (DEM) solver is used to solve the particle force analysis results and proppant fluid information to obtain proppant particle information including particle position, particle velocity, and particle momentum. The specific steps include: Based on the particle force analysis results and proppant fluid information, the particle collision sequence is obtained using a DEM solver; Conduct contact judgment between particles and obtain boundary contact forces between particles, including: The normal overlap between particles δ n The positive or negative value is used to determine whether there is contact collision between particles and walls, or between particles. Obtain the particle contact point, the particle center coordinates and radius, and then obtain the three-dimensional information of the particle. According to the three-dimensional information of the particle, obtain the normal overlap of the particle δ n , the formula is: ； Where, R i 、R j For particles i and particles j radius; r i 、r j Particles i and particles j The sphere center position vector; when When , it indicates that the two particles are in contact with each other; when When , the two particles have no contact; When determining whether a particle contacts a wall, the wall is considered to be a stationary spherical particle with infinite radius, and the normal overlap is used. δ n The calculation formula determines whether contact occurs; The boundary contact force between particles is obtained when two particles are in contact, and the normal contact force and tangential contact force are obtained: Obtain the normal contact force through the Hertz contact force model, the normal contact force of the proppant particles F cn,ij is the normal overlap δ n The function is expressed as follows: ； in, k n,ij is the normal stiffness coefficient; δ n,ij is the normal displacement between the colliding particles; n is the normal unit vector between contact particles; is the normal dissipation coefficient; is the normal component of the relative velocity between colliding particles; is the particle Young's modulus, Pa; is the equivalent particle radius; According to the physical parameters of the particles, the normal stiffness coefficient and normal dissipation coefficient are obtained as follows: ； ； ； ； Where, Y i 、 v i 、 R i and Y j 、 v j 、 R j Particles i Young's modulus, Poisson's ratio and radius of the particles, as well as j Young's modulus, Poisson's ratio and radius; The normal elastic damping of the particle is obtained, and its expression is: ； ； Where, m i 、 m j Particles i , particles j Mass, kg; m* is the equivalent mass, kg; is the normal component of the particle relative velocity; e is the collision restitution coefficient; s n 、 β By the restitution coefficient e Decide; Tangential contact force Obtained through Mindlin-Deresiewicz contact theory, using the tangential overlap δ t and tangential stiffness S t , the formula is as follows: ； ； ； Where G is the equivalent shear modulus; The expression of particle tangential elastic damping is: ； in, v t is the tangential component of the relative velocity of the particles; thus, the proppant particles j Acting on proppant particles i The contact force on is expressed as follows: ； in, F t,ij 、 F d t,ij 、 F n,ij and F d n,ij They are tangential elastic force, tangential damping, normal elastic force and normal damping respectively; The setting of the particle motion time step needs to be based on the Rayleigh time step, which refers to the time required for a shear wave to pass through a solid particle, and the maximum time step of the particle motion cannot exceed its Rayleigh time step. , the calculation formula is as follows: ； in, R s is the particle diameter; is the particle density; G s is the particle shear modulus, v s It is dimensionless; The particle velocity and position are updated according to the boundary contact force, and proppant particle information including particle position, particle velocity and particle momentum is obtained.

6. A proppant migration simulation device, characterized in that: include: Parameter setting module, used to obtain the fluid dynamic parameters of the proppant, including: fluid domain, grid parameters and fluid physical properties parameters; The migration simulation module is used to construct a fluid dynamics calculation model to describe the migration and flow characteristics of the proppant, and set constraints including boundary parameters, initial conditions, time step, and number of iterations. The fluid dynamics parameters of the proppant are input into the fluid dynamics calculation model, and a fluid dynamics CFD solver is used to obtain proppant fluid information including flow parameters and mechanical parameter vectors within the flow field. Based on the proppant fluid dynamics parameters and the number of proppant particles, a discrete element analysis model is used to perform particle force analysis. The discrete element method (DEM) solver is used to solve the particle force analysis results and proppant fluid information to obtain proppant particle information including particle position, particle velocity, and particle momentum. The iterative feedback module is used to solve the proppant particle information and proppant fluid information using the fluid-solid coupling model in the CFD solver, and update the flow field distribution and force distribution during the simulated proppant migration process; the flow field distribution and force distribution are fed back to the DEM solver to determine the error until the preset simulation time is met, completing the simulation of proppant migration.

7. A computer device comprising: memory and processor; The memory stores a computer program, wherein the processor implements a proppant migration simulation method according to any one of claims 1 to 5 when executing the computer program.

8. A computer-readable storage medium storing a computer program, characterized in that: When the computer program is executed by a processor, the method for simulating proppant migration according to any one of claims 1 to 5 is implemented.

Images