A fluid-solid coupling analysis method based on pore network and CFD-DEM model

Abstract

The present invention provides a fluid-solid coupling analysis method based on a pore network and a CFD-DEM model, including: establishing a porous medium mesoscopic pore network model; dividing the pore channels between the pore balls in the porous medium mesoscopic pore network model; respectively calculating the spatial coordinates of the center points in the divided pore channels; determining the fluid grid area where the pore channels are located according to the position of the center points; calculating the macroscopic permeability tensor of the porous medium based on the permeability of the pore channels in the fluid grid area; establishing a fluid control equation based on the change of the macroscopic permeability tensor; using the fluid control equation and the solid particle equation to establish the fluid-particle interaction control equation, that is, a fluid-solid coupling analysis method based on a pore network and a CFD-DEM model. According to the present invention, by replacing the permeability scalar with the permeability tensor, the flow and transport phenomena in the porous medium can be more accurately described.

Description

Technical Field

[0001] The invention relates to the technical field of numerical simulation, and in particular to a fluid-solid coupling analysis method based on pore network and CFD-DEM model. Background Art

[0002] The discrete element method (DEM) regards the soil as a pile of discrete particles. Given the contact model between particles, the displacement and velocity of the particles are calculated according to Newton's second law, so as to simulate the microscopic response of granular materials under different loading environments. By establishing a fluid-solid coupling method with DEM and computational fluid dynamics (CFD), the migration movement of soil particles at the microscopic scale under the action of liquid is simulated. This simulation method is called the CFD-DEM coupling analysis method. The CFD-DEM coupling analysis method is the basic method of most particle-fluid coupling models. At the same time, because the CFD-DEM analysis method has obvious advantages in numerical calculation consumption and accuracy, it is most widely used in the numerical simulation research of fluid-solid coupling in geotechnical engineering.

[0003] In the CFD-DEM analysis method, the fluid is regarded as a continuum, and the movement of the fluid is calculated by dividing the units. In this method, solid particles are regarded as single particles. Compared with the macro-scale model, the constructed meso-scale model improves the analysis accuracy of solid particle motion in fluid-solid coupling problems, and has obvious advantages in simulating the discontinuous deformation and large deformation flow of rock and soil. More importantly, the simulation of solid meso-scale can understand the meso-mechanism of discontinuous deformation and migration of rock and soil from the particle scale.

[0004] In the CFD-DEM analysis method, the fluid flow behavior around each DEM particle is not solved precisely, but is averaged by volume in the local domain; the CFD-DEM analysis method follows the local homogenization interaction method, ignoring the influence of particle distribution on fluid movement, that is, ignoring the anisotropy of the permeability characteristics of the broken rock mass in the fluid grid, and assuming that the permeability of each fluid grid is a scalar, and the movement direction of the fluid (seepage path) is only related to the flow field boundary and has nothing to do with the distribution of particles. At the same time, it is assumed that the spatial distribution of particles only affects the movement size of the fluid, not the movement direction of the fluid, which is obviously inconsistent with the seepage movement law of groundwater flow in geotechnical engineering. Therefore, the CFD-DEM analysis method cannot accurately describe the flow and transmission phenomena of fluids in porous media. Summary of the invention

[0005] The present invention provides a fluid-solid coupling analysis method based on a pore network and a CFD-DEM model. Based on the pore network model, the permeability tensor is introduced into the fluid-solid coupling analysis, which optimizes the deficiency of the traditional CFD-DEM analysis method that regards the fluid-particle interaction inside the fluid unit as locally averaged and isotropic, and can more accurately describe the flow and transport phenomena in porous media.

[0006] The present invention provides a fluid-solid coupling analysis method based on pore network and CFD-DEM model, which comprises:

[0007] Establish a mesoscopic pore network model for porous media;

[0008] Dividing the pore channels between the pore balls in the porous medium mesoscopic pore network model;

[0009] Calculating the spatial coordinates of the center points in the divided pore channels respectively;

[0010] Determining the fluid grid region where the pore channel is located according to the position of the center point;

[0011] Calculating a macroscopic permeability tensor of the porous medium based on the permeability of the pore channels within the fluid grid region;

[0012] The fluid control equation is established based on the change of the macroscopic permeability tensor;

[0013] The fluid control equation and solid particle equation are used to establish the fluid-particle interaction control equation, that is, the fluid-solid coupling analysis method based on the pore network and CFD-DEM model.

[0014] Preferably, the establishment of a porous medium mesoscopic pore network model comprises:

[0015] Data collection and image processing, and segmenting the processed image to distinguish between porous media and pore phases; extracting pore features from the segmented image to establish a porous media microscopic pore network model.

[0016] Preferably, dividing the pore channels between the pore balls in the porous medium mesoscopic pore network model comprises:

[0017] The Delaunay tetrahedron partitioning method is used to regard all solid particles in the pore network model as representative points of the sphere center or particle mass center to generate tetrahedrons; all representative points in the porous medium microscopic pore network model are partitioned, and the sphericity coefficient is introduced to determine and select the generated tetrahedrons;

[0018] Obtaining a pore sphere within the tetrahedron;

[0019] The distance between adjacent pore balls is used to classify the pore channels.

[0020] Preferably, dividing the pore channels between the pore balls in the porous medium mesoscopic pore network model further includes:

[0021] All representative points in the porous media microscopic pore network model are divided, and the sphericity coefficient is introduced to determine and select the generated tetrahedrons.

[0022] Preferably, the types of the pore channels are:

[0023] Obtain the ID of each tetrahedron and the radius of the four spheres located at the vertices of each tetrahedron and the coordinates of the center of the common tangent sphere in the coordinate tetrahedron, establish a nonlinear equation group of the pore sphere inside the tetrahedron, introduce an iterative algorithm for the related nonlinear equation group, and obtain the coordinates and radius of the center of the largest common tangent sphere inside each tetrahedron;

[0024] Given two adjacent pore spheres with IDs i and j, the corresponding pore sphere radii are r i and r j , d ij is the distance between the centers of two pore balls. According to the distance between two adjacent pore balls, pore channels are divided into three types:

[0025] Type A pore channel: separated by two pore balls, i.e., d ij ≥r i +r j ;

[0026] Type B pore channel: two pore spheres intersect, i.e. |r i -r j |≤d ij <r i +r j ;

[0027] Type C pore channel: The sphere with smaller radius is contained by the sphere with larger radius, i.e. d ij <|r i -r j |.

[0028] Preferably, the spatial coordinates of the center point in the pore channel after the calculation and division are:

[0029] The unit direction vector of the pore channel and the spatial coordinates of the pore sphere are calculated, and the spatial coordinates of the center point are obtained respectively according to the difference of the pore channels.

[0030] Preferably, the determination of the pore channel and the fluid grid area according to the position of the center point includes:

[0031] When the center point of the pore channel is contained in a certain fluid grid, the corresponding pore channel is completely contained in the fluid grid;

[0032] In the process of determining the pore channels in the fluid grid, the following conditions should be met:

[0033]

[0034]

[0035] Among them, d c is the minimum calculated width of the fluid flow domain, Δx cfd is the calculation length of the fluid grid, and r is the maximum particle size in the fluid grid.

[0036] Preferably, the calculation of the macroscopic permeability tensor of the porous medium based on the permeability of the pore channels in the fluid grid region is calculated as:

[0037] By comparing the macroscopic flow velocity of the fluid grid with the macroscopic Darcy's law, the permeability tensor of the fluid grid model of the macroscopic structure is obtained. Based on the assumption of macroscopic isotropic medium in the local averaging theory, the connection between the permeability in the pore channel and the macroscopic average permeability k is established, and the relationship between the macroscopic permeability tensor and the macroscopic average permeability is further established.

[0038] Preferably, the solid particle control equation is established as follows:

[0039] The movement of particles is governed by Newton's laws, as follows:

[0040]

[0041]

[0042] in, represents the particle movement speed, t represents the numerical calculation time, Represents the particle contact force between particles, including particle-particle contact force and particle-wall contact force, and reflects the interaction between particles through the contact model. Represents the non-contact force between particles (mainly including van der Waals force, electromagnetic force and capillary force), represents the external force, represents the total fluid-particle interaction force acting on the particle, m p represents the particle mass, represents the acceleration due to gravity, represents the particle angular velocity, represents the contact torque on the particle, I p represents the particle moment of inertia.

[0043] Preferably, the fluid-particle interaction governing equation is:

[0044] Interaction force between fluid and particles It consists of three parts: the drag force of the fluid, the force caused by the pressure gradient of the fluid, and the buoyancy of the fluid on the particle, as follows:

[0045]

[0046] According to the coarse-grid fluid-particle coupling numerical analysis system, the fluid drag force is defined as the drag force exerted by the fluid on the particle, always acting about the particle's center of mass, with no rotational torque exerted by the particle, and is calculated as:

[0047]

[0048] in, is the drag force on a single particle, is the average porosity of the fluid grid where the particle is located, χ is an empirical correction factor for the local average porosity, this correction term applies the force to high-porosity and low-porosity systems and a wide range of Reynolds numbers;

[0049]

[0050] Among them, C d is the resistance coefficient, ρ f is the fluid density, r is the particle radius, is the fluid velocity, is the particle velocity;

[0051] Among them, the resistance coefficient C d can be calculated as:

[0052]

[0053] Among them, Re p is the particle Reynolds number;

[0054] In addition, the modified empirical factor χ is defined as:

[0055]

[0056] The Reynolds number Rep is calculated as:

[0057]

[0058] Where: μ f is the dynamic viscosity coefficient of the fluid.

[0059] The fluid-solid coupling analysis method based on pore network and CFD-DEM model provided by the present invention has the following beneficial effects:

[0060] 1. Compared with the existing CFD-DEM fluid-solid coupling analysis method, the PNM-CFD-DEM coupling analysis method of the present invention not only updates the average porosity of each fluid grid in each bidirectional coupling, but also updates the pore network generated by the particle movement, and the anisotropic permeability tensor characterizing the spatial distribution of the pore channels in each fluid grid is also updated accordingly. The CFD seepage field is updated under the fluid grid water pressure / flow velocity boundary and the anisotropic permeability tensor distribution, and then the updated seepage field information is transmitted to the DEM particle movement until the set number of cycles is reached; the permeability tensor established based on the pore network is introduced into the fluid-solid coupling cycle to weaken the shortcomings of the coarse grid fluid-particle coupling numerical analysis system that regards the fluid-particle interaction inside each fluid grid as isotropic and uniform. The direction of the flow velocity is optimized from being controlled by the pressure gradient to being controlled by both the pressure gradient and the pore network of the particles. Therefore, by replacing the permeability scalar with the permeability tensor, the flow and transport phenomena in the porous medium can be more accurately described;

[0061] 2. The present invention also analyzes the sensitivity of fluid grid and particle size in the fluid-particle coupling model of porous media under irregular arrangement, and finds that the PNM-CFD-DEM fluid-solid coupling method reduces the sensitivity of fluid grid and particle size; at the same time, the flow velocity direction in the fluid grid is further corrected by the permeability tensor established based on the pore grid, so that the flow velocity direction of each fluid grid can reflect the microscopic seepage channel of the particles in each fluid grid in real time, making the direction of the interaction drag force between the fluid and the particles more reasonable. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 It is a schematic diagram of the pore channel division of the present invention;

[0063] Figure 2 A schematic diagram of the type A pore channel and pore ball model of the present invention is established;

[0064] Figure 3 This is a schematic diagram of the BOND bond connecting two pore balls of the present invention;

[0065] Figure 4 A schematic diagram of the type B pore channel and pore sphere model of the present invention is established;

[0066] Figure 5 Schematic diagram of two porous balls of the present invention;

[0067] Figure 6 It is a schematic diagram for distinguishing the spatial position of the pore channel of the present invention;

[0068] Figure 7 It is a schematic diagram of the calculation process of the pore network model of the present invention;

[0069] Figure 8 It is a schematic diagram of the PorePy-PFC3D fluid-solid coupling analysis of the present invention;

[0070] Fig. 9 It is a schematic diagram of the numerical model of porous medium loss under uneven arrangement of the present invention;

[0071] Fig.10 Schematic diagram of fluid grids of different sizes according to the present invention;

[0072] Fig.11 It is a schematic diagram of flow field distribution and particle displacement distribution of the CFD-DEM model of the present invention;

[0073] Fig.12 Schematic diagram of flow field distribution and particle displacement distribution of the PNM-CFD-DEM model of the present invention;

[0074] Fig.13 It is a schematic diagram of the maximum velocity change of the fluid grid under different grid sizes of the present invention;

[0075] Fig.14 It is a schematic diagram of the change of the maximum displacement of particles under different grid sizes of the present invention;

[0076] Fig.15 A schematic diagram of the drag force distribution near the inlet boundary of the present invention;

[0077] Fig.16 This is a schematic diagram of the drag force distribution of the CFD-DEM model of the present invention;

[0078] Fig.17 It is a schematic diagram of the drag force distribution of the PNM-CFD-DEM model of the present invention. DETAILED DESCRIPTION

[0079] like Figure 1-17 As shown, the present invention provides a fluid-solid coupling analysis method based on pore network and CFD-DEM model, which includes the following steps:

[0080] Establish a mesoscopic pore network model for porous media;

[0081] Dividing the pore channels between the pore balls in the porous medium mesoscopic pore network model;

[0082] Calculating the spatial coordinates of the center points in the divided pore channels respectively;

[0083] Determining the fluid grid region where the pore channel is located according to the position of the center point;

[0084] Calculating a macroscopic permeability tensor of the porous medium based on the permeability of the pore channels within the fluid grid region;

[0085] The fluid control equation is established based on the change of the macroscopic permeability tensor;

[0086] The fluid control equation and solid particle equation are used to establish the fluid-particle interaction control equation, that is, the fluid-solid coupling analysis method based on the pore network and CFD-DEM model.

[0087] Specifically, the porous media mesoscopic pore network model helps to understand the physical properties of porous media, such as permeability, strength, seepage, etc., as well as their performance in engineering applications, including:

[0088] Data acquisition and image processing, collect sample images or experimental data, such as through CT scans, MRI, microscope images, digital photos, etc.; preprocess the original image by denoising, enhancing contrast, and segmenting multiple pore phases in the image; image segmentation, use image processing algorithms to separate the porous medium and pore phases in the image, which usually involves thresholding, region growing, edge detection and other methods; pore feature extraction, extract pore shape, size, distribution and connectivity, etc., which can be achieved by calculating geometric properties, topological properties or other statistical methods; establish a pore network, based on the extracted pore features, a pore network model can be established, in which pores are regarded as nodes and the relationships between pores are modeled as edges, which can be completed using graph theory methods.

[0089] Specifically, the tetrahedron is generated:

[0090] The Delaunay tetrahedron partitioning method is a commonly used method in the field of computer graphics and computational geometry. It is used to partition a set of discrete points into non-overlapping tetrahedrons for data analysis and processing in three-dimensional space. The Delaunay tetrahedron partitioning method is used to treat all solid particles in the pore network model as representative points of their sphere centers or particle mass centers. All representative points in the entire model are tetrahedron partitioned. In order to avoid the generated tetrahedrons being too deformed, the concept of sphericity coefficient is introduced to screen the tetrahedrons generated on the boundary. The sphericity coefficient of the corresponding tetrahedron is defined as the ratio of the surface area of ​​the tetrahedron to its volume. The sphericity coefficient of the tetrahedron is calculated as:

[0091]

[0092] Where: θ is the tetrahedron sphericity coefficient, V p is the volume of the tetrahedron, S p is the surface area of ​​the tetrahedron;

[0093] Among them, the sphericity coefficient of the sphere is 1, and the sphericity coefficient of the regular tetrahedron is 0.671. According to historical data and empirical values, the threshold of the tetrahedron sphericity coefficient in the pore network model is set to 0.1. It is judged whether the sphericity coefficient of the generated boundary tetrahedron is less than the set threshold. When it is less than the set threshold, the tetrahedron generated at the boundary is judged to be too flat, and it is deleted or re-divided.

[0094] Specifically, the internal pore sphere of the tetrahedron is determined by:

[0095] After the Delaunay tetrahedron is generated, the ID of each tetrahedron and the radius and coordinates of the four spheres located at the vertices of each tetrahedron are recorded, and the index of each tetrahedron and the four adjacent tetrahedrons is established; let the vertices of the tetrahedron with ID i be numbered P i1 , P i2 , P i3 , P i4 The radius of each vertex is r i1 、r i2 、r i3 、r i4 , the coordinates correspond to (P ilx P ily P ilz )、(P i2x P i2y P i2z )、(P i3x P i3y P i3z )、(P i4x P i4y P i4z ), assuming that the coordinates of the center of the tangent sphere in the tetrahedron are X i , Y i , Z i , the radius of the common tangent sphere of the tetrahedron is r i , then the nonlinear equations of the pore sphere inside the tetrahedron are established, and the iterative algorithm of the related nonlinear equations can be introduced to obtain the coordinates and radius of the center of the largest common tangent sphere inside each tetrahedron to calculate the pore part in the pore network structure;

[0096] The nonlinear equations of the tetrahedron internal pore sphere are established as follows:

[0097] (P i1x -X i ) 2 +(P i1y -Y i ) 2 +(P i1z -Z i ) 2 =(r i1 +r i )2

[0098] (P i2x -X i ) 2 +(P i2y -Y i ) 2 +(P i2z -Z i ) 2 =(r i2 +r i ) 2

[0099] (P i3x -X i ) 2 +(P i3y -Y i ) 2 +(P i3z -Z i ) 2 =(r i3 +r i ) 2

[0100] (P i4x -X i ) 2 +(P i4y -Y i ) 2 +(P i4z -z i ) 2 =(r i4 +r i ) 2

[0101] Specifically, the pore channels between the pore balls are divided into:

[0102] The pore space inside the Delaunay tetrahedron represents a pore sphere. Each pore sphere generated using the Delaunay tetrahedron method has four adjacent pore spheres. The corresponding pore spheres share the same triangular surface. Given two adjacent pore spheres with IDs i and pore j, the corresponding pore sphere radii are r i and r j , d ij is the distance between the centers of two pore balls. According to the distance between two adjacent pore balls, pore channels are divided into three types:

[0103] Type A pore channels: Figure 1 As shown in a, the two pore balls are separated, that is, d ij ≥r i +r j ;

[0104] Type B pore channels: Figure 1 As shown in b, the two pore balls intersect, that is, |r i -r j |≤d ij <r i +r j ;

[0105] Type C pore channels: Figure 1 As shown in c, the ball with a smaller radius is contained by the ball with a larger radius, that is, d ij <|r i -r j |;

[0106] Among them, in the C-type pore channels, the pore spheres with smaller radius are completely contained by the pore spheres with larger radius, so the C-type pore channels are regarded as non-directional (extension direction) pore channels. It is necessary to merge this type of pore spheres, integrate the pore spheres with smaller radius into the pore spheres with larger radius, recalculate the distance between adjacent pore spheres, and divide the pore channels again until the C-type pore channels no longer appear.

[0107] Specifically, the pore channel direction vector is calculated as:

[0108] In the Delaunay tetrahedron decomposition method, the largest inscribed sphere of each Delaunay tetrahedron represents a pore sphere, and four pore channels extend outward through the four surfaces of the Delaunay tetrahedron. The four pore channels are connected to the four pore spheres respectively, and the direction vector connecting the pore spheres is also the extension direction of the pore channel.

[0109] Among them, the directional unit vector in the spatial parameters of the pore channel can be expressed as:

[0110]

[0111] Specifically, the type A pore channel model calculates:

[0112] By setting convergent-divergent functional derivative bonds between the separated pores in the porous medium to connect the adjacent pores, the microstructure network of the porous medium can be better simulated, such as Figure 2 As shown, the convergent-divergent functionally derived bond morphology of the pore channel satisfies the bipyramidal asymmetric concentric bond, i.e., BACON bond (Biconicalabscissa AsymmetricCONcentric; BACON);

[0113] Among them, Figure 3 As shown in the figure, the BACON bond connects the positive contours of two pore balls. When the radius of pore ball i is r i, at the center point A, ξ=0; when the radius of the pore sphere j is r j , at the center point B, ξ=1; when fluid passes through the pore channel, ξ F and G They are the inlet and outlet of the BACON bond, that is, the inlet and outlet for the fluid to pass through; the functions at both ends of the BACON bond are in ξ t (Point E), this position represents the minimum cross-sectional area of ​​the pore channel, that is, the throat;

[0114] Among them, the dimensionless function r(ξ) at both ends of the BACON bond inlet and outlet can be expressed as:

[0115]

[0116] Among them, n R It is the curvature parameter that determines the longitudinal shape and size of the BACON bond. ξ represents any point in the distance between two pore balls, and its value ranges from 0 to 1.

[0117] The expression of ξ is:

[0118]

[0119] Where x is the distance from the center of the pore sphere to any point between the centers of the adjacent pore spheres;

[0120] Its i and T j are the r(ξ) values ​​of the corresponding pore spheres i and j:

[0121]

[0122] Wherein, ζ is the coordination number on the two-dimensional plane. In the present invention, each pore ball is connected to four adjacent pore balls, that is, the coordination number ζ is 4, π / ζ is the BAF angle, and n R It is the curvature parameter that determines the longitudinal shape and size of the BACON bond;

[0123] Determine the minimum curvature parameter n based on the pore sphere radius at both ends of the BACON bond R , when the pore sphere radius at both ends of the BACON bond is equal, n R = 0; when the pore sphere radius at both ends of the BACON bond is not equal, n R >0;

[0124] According to the pore-sphere spacing d ij For the radius of the pore sphere r i and r j Perform dimensionless processing to obtain dimensionless pore size

[0125]

[0126] Among them, r i(j) is the radius of the pore sphere i or j;

[0127] According to the empirical value, the sphere is required to have a pore radius R that cannot exceed d ij / 2; the dimensionless size of the pore balls at both ends exists and Then the minimum curvature parameter n R can be calculated as:

[0128]

[0129] The volume of the A-type pore channel is the BOND bond volume, which is the total volume V of the functions on both ends of the BOND bond. BB :

[0130] V BB =V BL +V BR

[0131] Among them, V BL Represents the volume of the left end of the BOND bond, such as Figure 3 As shown F ≤ξ≤ξ t , V BR Represents the volume of the right side of the BOND bond, such as Figure 3 As shown t ≤ξ≤ξ G , then the volumes of the left and right ends can be calculated as:

[0132]

[0133]

[0134] Among them, a x is the geometric shape factor, the circular cross section a x =π, square cross section a x =4, a circular cross section in the present invention;

[0135] like Figure 2 As shown in a, the entrance and exit ξ at both ends of the BACON key F and G Respectively expressed as:

[0136]

[0137]

[0138] The throat area of ​​the BACON key t (Point E) can be expressed as:

[0139]

[0140] Then V BL and V BR Respectively expressed as:

[0141]

[0142]

[0143] The total volume of type A pore channels V can be obtained BB ;

[0144] By T i(j) and t The distance between the pore throat and the pore ball i is:

[0145]

[0146] Then the unit direction vector of the pore channel and the spatial coordinate of the pore ball i (X i , Y i , z i ) can be obtained, the spatial coordinates of the pore throat E t for:

[0147]

[0148] Specifically, the type B pore channel model calculates:

[0149] When two adjacent pore balls intersect, the flow direction between the pore balls is the intersection between the pore balls as the pore channel, such as Figure 4 As shown in the figure, the 3D model of the pore channel is approximately represented when the pore balls intersect; the cylinders between the pore balls serve as pore channels, as shown in the figure. Figure 4 As shown in a, the radius of the cylinder is the radius of the common cross section between the two pore balls, as shown in Figure 4 As shown in b, the length of the cylinder is 1 / 10 of the radius of the smaller of the two adjacent balls; Figure 5 As shown, this is the side view of two adjacent pore balls intersecting;

[0150] According to Heron's volume formula P = (r i +r j +d ij ) / 2, we can get the intersection radius of the pore sphere:

[0151]

[0152] The volume of the cylinder is:

[0153]

[0154] Then the unit direction vector of the pore channel and the spatial coordinate of the pore ball i (X i , Y i , Z i ) we can get the spatial coordinates of the center point of the cylinder: t for:

[0155]

[0156] Specifically, pore channels and fluid grid areas are distinguished:

[0157] Based on the Particle Centroid Method, the pore channels contained in each fluid grid are judged and calculated, such as Figure 6 As shown in Figure 2, when the center (pore throat) of a pore channel is contained in a fluid grid, the corresponding pore channel is completely contained in this fluid grid;

[0158] When calculating, the particle center body method only needs to determine whether the center of the pore channel (pore throat) is within the fluid grid. Therefore, it has a relatively high calculation efficiency. However, it should be noted that the particle center body method does not distinguish between the volume of the pore channel located within the fluid grid and the volume outside the grid, which will cause errors in the pore volume calculation, especially when the particle size is close to the grid size. This error cannot be ignored. In order to achieve the coupling calculation between particles and fluids and ensure the calculation accuracy of this method, the fluid grid should be divided finely enough to satisfy the calculation of fluid motion. At the same time, a fluid grid must contain a certain number of particles, which should satisfy the following inequality:

[0159]

[0160]

[0161] Among them, d c is the minimum calculated width of the fluid flow domain, Δx cfd is the calculation length of the fluid grid, and r is the maximum particle size in the fluid grid.

[0162] Specifically, the macroscopic permeability tensor of porous media is calculated as:

[0163] Based on the established pore network model, the permeability tensor calculation formula that can characterize the spatial distribution of the pore network is derived. Assuming that the water flow in each pore channel conforms to Darcy's theorem, the flow velocity v of the seepage pipe is pipe It can be expressed as:

[0164]

[0165] Where lp is the length of the pore channel, k is the permeability of the pore channel, the permeabilities of each pore channel are considered to be the same, μ is the dynamic viscosity of water, and P j , P i They represent the pore water pressure at the outflow and inlet of the pore channel respectively. For type A pore channel, the pore channel length l p can be calculated as:

[0166] l p =d ij (ξ F -ξ G )

[0167] The macroscopic seepage velocity in each fluid grid can be defined as:

[0168]

[0169] in, is the flow velocity along channel i. By restricting the flow in the channel within the fluid, this integral can be replaced by the sum of all M channels contained in the fluid grid volume V:

[0170]

[0171] in, is the volume of the pore channel, V BB is the total volume of the pore channels in the fluid grid;

[0172] The flow rate in different channels In addition to being related to the water velocity, it is also related to the unit vector of the pipe, which can be expressed as:

[0173]

[0174] in, Representing the unit vector of channel i, the flow velocity in different pore channels can be calculated as:

[0175]

[0176] Where ΔP is the pressure difference between the pipes;

[0177] Assume that the average pressure gradient within the fluid grid is Applicable to each pressure level within the pipe, the pressure gradient ΔP can then be calculated as:

[0178]

[0179] Macroscopic flow velocity of the fluid grid <v i > can be calculated as:

[0180]

[0181] By comparing this expression with the macroscopic Darcy's law, it can be found that the permeability tensor of the fluid grid model of the macrostructure can be expressed as:

[0182]

[0183] Written in 3×3 symmetric matrix form:

[0184]

[0185] For the average penetration for:

[0186]

[0187] Then the pore channel permeability k can be calculated as:

[0188]

[0189] Then the macroscopic permeability tensor of a single fluid grid can be calculated as:

[0190]

[0191] Specifically, the Darcy governing equation for fluid based on the permeability tensor:

[0192] Based on the fluid control equation of the coarse grid CFD-DEM coupled numerical analysis system, the permeability tensor k is introduced based on the pore network. ij , then the fluid control equation based on Darcy's law and the classical Kozeny-Kalman equation is as follows:

[0193]

[0194]

[0195]

[0196]

[0197] in: represents the macroscopic average permeability, which is a function of the porosity of the fluid grid.

[0198] Specifically, the governing equations for the solid are:

[0199] The movement of particles is governed by Newton's laws, as follows:

[0200]

[0201]

[0202] in, represents the particle movement speed, t represents the numerical calculation time, Represents the particle contact force between particles, including particle-particle contact force and particle-wall contact force, and reflects the interaction between particles through the contact model. Represents the non-contact force between particles (mainly including van der Waals force, electromagnetic force and capillary force), represents the external force, represents the total fluid-particle interaction force acting on the particle, m p represents the particle mass, represents the acceleration due to gravity, represents the particle angular velocity, represents the contact torque on the particle, I p represents the particle moment of inertia;

[0203] The mutual collision motion between particles, that is, the force-displacement motion equation, uses linear contact model, parallel bonding model, rolling friction model, etc. according to needs.

[0204] Specifically, the governing equations for fluid-particle interactions are:

[0205] Interaction force between fluid and particles It consists of three parts: the drag force of the fluid, the force caused by the pressure gradient of the fluid, and the buoyancy of the fluid on the particle, as follows:

[0206]

[0207] According to the coarse-grid fluid-particle coupling numerical analysis system, the fluid drag force is defined as the drag force exerted by the fluid on the particle, always acting about the particle's center of mass, with no rotational torque exerted by the particle, and is calculated as:

[0208]

[0209] in, is the drag force on a single particle, is the average porosity of the fluid grid where the particle is located, χ is an empirical correction factor for the local average porosity, this correction term applies the force to high-porosity and low-porosity systems and a wide range of Reynolds numbers;

[0210]

[0211] Among them, C d is the resistance coefficient, ρ f is the fluid density, r is the particle radius, is the fluid velocity, is the particle velocity;

[0212] Among them, the resistance coefficient C d can be calculated as:

[0213]

[0214] Among them, Re p is the particle Reynolds number;

[0215] In addition, the modified empirical factor χ is defined as:

[0216]

[0217] Reynolds number p Calculated as:

[0218]

[0219] Where: μ f is the dynamic viscosity coefficient of the fluid.

[0220] Specifically, the technical solution of the present invention is further explained in combination with the calculation case provided in the PFC3D software user manual document. A PNM-CFD-DEM fluid-solid coupling analysis method based on a pore network model and a traditional CFD-DEM model is used. The mesoscopic pore network model is introduced to improve the existing coarse grid CFD-DEM coupling numerical analysis method based on the local averaging theory. First, the pore space in each fluid grid is simplified into a combination of pore balls and pore channels (pore throats). The pore balls between particles are calculated by the Delaunay tetrahedron method. Based on the distance between the pore balls and the radius of the pore balls, the pore channels between different pore balls are calculated to obtain the spatial calculation parameters of the pore channels. The spatial calculation parameters include: direction vector, throat position coordinates, and volume. It is assumed that the flow velocity in each pore channel is laminar flow in accordance with Darcy's law. The permeability in the mesoscopic pore channel is k, and the macroscopic permeability tensor k in the fluid grid is established. ii Based on the assumption of macroscopically isotropic medium in the local averaging theory, the relationship between the permeability k in the pore channel and the macroscopic average permeability is established. The macroscopic permeability tensor k is further established. ij Macro-average permeability The relationship expression between the macroscopic permeability tensor k in each fluid grid ij The calculation process is as follows Figure 7 shown.

[0221] For the PNM-CFD-DEM coupling analysis method, the calculation of the DEM particle motion field is completed using the commercial software PFC3D, and the seepage field calculation of the rock and soil body is completed using the open source porous media seepage simulation open source software Porepy. The pore network model is established by the Delaunay database in the Python computer language combined with secondary development, and then the Taichi language that supports GPU parallel computing engine technology that can run in the Python language environment is used to complete the identification of pore channels in each fluid grid and the calculation of the permeability tensor.

[0222] In order to facilitate the information interaction between fluid motion and solid particle motion in the PNM-CFD-DEM coupling analysis method and to promote the application of the traditional CFD-DEM coupling analysis method, the advantages of Porepy and PFC3D software supporting both Python computer language and Window computer operating system are utilized, and the TCPsocket communication transmission control protocol under the itasca module in the Python computer language is used to realize the direct interaction of fluid motion and solid particle motion in Porepy and PFC3D software, and establish the Porepy-PFC3D fluid-solid coupling analysis framework to improve the information interaction efficiency of the fluid-solid bidirectional coupling process, without building an intermediate communication platform and virtual machine, so as to avoid the distortion of program calculation caused by the blockage of communication channels in the information interaction of multiple computer languages ​​and multiple computer operating systems; the PNM-CFD-DEM coupling analysis method is applied to the Porepy-PFC3D fluid-solid coupling analysis framework such as Figure 8 shown.

[0223] Based on the pore network, the permeability scalar of each fluid grid is converted into a permeability tensor that can reflect the anisotropic seepage characteristics of particles in the fluid grid. In order to prove that the PNM-CFD-DEM fluid-solid coupling analysis method can reduce the sensitivity of the ratio of fluid grid to particle size and improve its applicability in simulations with different fluid grid to particle size ratios, two calculation methods are used to demonstrate and analyze the existing coarse grid fluid-particle coupling calculation method and the coarse grid fluid-particle coupling calculation method based on pore network optimization through the PorePy-PFC3D fluid-solid coupling framework; Fig. 9 The numerical model of porous media under uneven arrangement is shown in the figure. The total number of particles is 10587. According to the initial porosity of 0.4, the particles of two sizes (0.0025m and 0.00375m) are randomly distributed and generated. The mass proportion of particles with a size of 0.0025m is 60%, and the mass proportion of particles with a size of 0.00375m is 40%. Other boundary conditions, particles, inter-particle contacts and fluid parameters are as follows: Fig. 9 As shown in the model, the physical calculation time is set to 1s, and every 5×10 -2 s coupling once, a total of 20 couplings; each fluid grid of different sizes is as follows Fig.10 As shown, in order to ensure that the fluid boundaries of each model are the same, the grid size in the Y direction is adjusted, where Figure a represents a 10×20×10 fluid grid, Figure b represents a 10×15×10 fluid grid, Figure c represents a 10×10×10 fluid grid, and Figure d represents a 10×5×10 fluid grid.

[0224] like Fig.11 The flow field distribution and particle displacement distribution of different fluid grid sizes under the CFD-DEM model are shown, where the particle displacement field is a cross-sectional view of X=0.05, where Figure a represents a 10×20×10 fluid grid, Figure b represents a 10×15×10 fluid grid, Figure c represents a 10×10×10 fluid grid, and Figure d represents a 10×5×10 fluid grid; Fig.12 The flow field distribution and particle displacement distribution of different fluid grid sizes under PNM-CFD-DEM are shown in the figure, where Figure a represents a 10×20×10 fluid grid, Figure b represents a 10×15×10 fluid grid, Figure c represents a 10×10×10 fluid grid, and Figure d represents a 10×5×10 fluid grid; from the grid size division, it can be seen that when the number of grids in the Y direction is 20, the volume of each fluid grid is the smallest, and they are all regular hexahedrons, and the calculation results are more accurate; when the number of grids in the Y direction is 20, the maximum velocity and maximum displacement of the CFD-DEM flow field are taken as the benchmark values, and the relative error changes of the maximum flow field velocity and particle displacement of the two fluid-solid coupling models under different grid sizes are compared, as shown Fig.13 and 14 as shown in .

[0225] like Fig.13 The maximum velocity error change of the fluid grid shown in Figure 1 shows that the correction effect of the permeability based on the pore grid is significant. Although the number of fluid grids is decreasing, the error of the corrected velocity change is relatively insignificant. In particular, when the grid size in the Y direction of the fluid grid is 5, the velocity in the fluid grid decreases sharply, while the corrected velocity increases. The relative change is not very obvious. Fig.14 It can be found from the change of the maximum displacement error of the particles shown that the change of the corrected particle displacement error is also relatively gentle.

[0226] Among them, comparison Fig.11 and Fig.12From the velocity field distribution, it can be found that the velocity direction of each unit of the optimized velocity has changed, especially near the inlet flow boundary, it is no longer uniformly distributed. In the CFD-DEM fluid-solid coupling model, the velocity direction of the fluid is only affected by the overall pressure gradient direction, while in the PNM-CFD-DEM fluid-solid coupling model, the velocity direction in each fluid grid is not only affected by the overall pressure gradient direction, but also by the seepage channel formed by the particle accumulation state inside the fluid grid, that is, the pore network. By comparing the particle displacement field, it can be found that the particle displacement near the inlet flow boundary is affected by Fig.11 The plane in becomes Fig.12 For the arc surface in the figure, when the grid size in the Y direction is 20, the drag force on the particles in the range of Y = 0 to 0.05 near the inlet flow boundary before and after optimization is analyzed, as shown in Figure 2. Fig.15 As shown in the figure, Figure a represents before optimization, and Figure b represents after optimization. The permeability tensor calculated based on the pore network can reflect the stacking state in the fluid grid in real time, so that the flow velocity in each fluid grid seeps along different directions. At the same time, the overall flow velocity at the critical boundary is relatively low, and the flow velocity near the center is generally higher. This law is closer to the actual flow field distribution, so the drag force on the particles is greater at the center than at the periphery, which further makes the displacement center of the particles greater than the periphery, and an arc-shaped displacement field is generated.

[0227] In order to further analyze the influence of the establishment of pore network structure on the drag force in fluid-solid coupling force, the drag force at different times in the range of Y = 0 ~ 0.05, X = 0 ~ 0.05, and Z = 0 ~ 0.05 near the inlet flow boundary is selected for analysis, such as Fig.16 As shown in the figure, when the particles in the fluid grid migrate, the average porosity increases, and the corresponding permeability increases accordingly, which eventually leads to an increase in flow velocity and drag force. However, as the particles migrate, the stacking state of the particles also changes accordingly, and the average porosity cannot reflect the change in the stacking state of the particles, so the drag force only increases in value. Fig.17 It can be seen from the figure that the green model is the pore pipe distribution, and the khaki model is the pore ball distribution. With the migration of particles and the increase of seepage erosion time, the particle pore network (particle stacking state characteristics) is also changing accordingly, especially the pore channel in the middle of the entire flow field becomes larger, and the corresponding changes are also reflected based on the permeability tensor calculated based on the pore channel. The size and direction of the seepage velocity are adjusted, which further leads to changes in the size and direction of the drag force.

[0228] Compared with the existing CFD-DEM fluid-solid coupling analysis method, the PNM-CFD-DEM coupling analysis method not only updates the average porosity of each fluid grid in each bidirectional coupling, but also updates the pore network generated by particle movement. The anisotropic permeability tensor that characterizes the spatial distribution of pore channels in each fluid grid is also updated accordingly. The CFD seepage field is updated under the fluid grid water pressure / velocity boundary and anisotropic permeability tensor distribution, and then the updated seepage field information is transmitted to the DEM particle movement until the set number of cycles is reached; the permeability tensor established based on the pore network is introduced into the fluid-solid coupling cycle to weaken the shortcomings of the coarse grid fluid-particle coupling numerical analysis system that regards the fluid-particle interaction within each fluid grid as isotropic and uniform. The direction of the flow velocity is optimized from being controlled by the pressure gradient to being controlled by both the pressure gradient and the pore network of the particles. Therefore, by replacing the permeability scalar with the permeability tensor, the flow and transport phenomena in porous media can be described more accurately.

[0229] By analyzing the sensitivity of fluid grid and particle size in the fluid-particle coupling model of porous media under irregular arrangement, it is found that the PNM-CFD-DEM fluid-solid coupling method reduces the sensitivity of fluid grid and particle size; at the same time, the flow velocity direction in the fluid grid is further corrected by the permeability tensor established based on the pore grid, so that the flow velocity direction of each fluid grid can reflect the microscopic seepage channel of the particles in each fluid grid in real time, making the direction of the interaction drag force between the fluid and the particles more reasonable.

[0230] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principle of the present invention shall be included in the protection scope of the present invention. The above description is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and variations can be made without departing from the technical principles of the present invention. These improvements and variations should also be regarded as the protection scope of the present invention.

Claims

1. A fluid-solid coupling analysis method based on pore network and CFD-DEM model, characterized in that: include: Establish a mesoscopic pore network model for porous media; Dividing the pore channels between the pore balls in the porous medium mesoscopic pore network model; Calculating the spatial coordinates of the center points in the divided pore channels respectively; Determining the fluid grid region where the pore channel is located according to the position of the center point; Calculating a macroscopic permeability tensor of the porous medium based on the permeability of the pore channels within the fluid grid region; The fluid control equation is established based on the change of the macroscopic permeability tensor; The fluid control equation and solid particle equation are used to establish the fluid-particle interaction control equation, that is, the fluid-solid coupling analysis method based on the pore network and CFD-DEM model; The calculation of the macroscopic permeability tensor of porous media based on the permeability of pore channels in the fluid grid region is: By comparing the macroscopic flow velocity of the fluid grid with the macroscopic Darcy's law, the permeability tensor of the fluid grid model of the macroscopic structure is obtained. Based on the assumption of macroscopic isotropic medium in the local averaging theory, the connection between the permeability in the pore channel and the macroscopic average permeability k is established, and the relationship between the macroscopic permeability tensor and the macroscopic average permeability is further established.

2. The fluid-solid coupling analysis method based on pore network and CFD-DEM model according to claim 1 is characterized in that: The method of establishing a porous medium microscopic pore network model comprises: Data collection and image processing, and segmenting the processed image to distinguish between porous media and pore phases; extracting pore features from the segmented image to establish a porous media microscopic pore network model.

3. The fluid-solid coupling analysis method based on pore network and CFD-DEM model according to claim 1, characterized in that: The pore channels between the pore balls in the porous medium mesoscopic pore network model are divided, including: The Delaunay tetrahedron partitioning method is used to regard all solid particles in the pore network model as representative points of the sphere center or particle mass center to generate tetrahedrons; all representative points in the porous medium microscopic pore network model are partitioned, and the sphericity coefficient is introduced to determine and select the generated tetrahedrons; Obtaining a pore sphere within the tetrahedron; The distance between adjacent pore balls is used to classify the pore channels.

4. The fluid-solid coupling analysis method based on pore network and CFD-DEM model according to claim 3 is characterized in that: The pore channels between the pore balls in the porous medium mesoscopic pore network model are divided, and further comprising: All representative points in the porous media microscopic pore network model are divided, and the sphericity coefficient is introduced to determine and select the generated tetrahedrons.

5. The fluid-solid coupling analysis method based on pore network and CFD-DEM model according to claim 3, characterized in that: The categories of the pore channels are: Obtain the ID of each tetrahedron and the radius of the four spheres located at the vertices of each tetrahedron and the coordinates of the center of the common tangent sphere in the coordinate tetrahedron, establish a nonlinear equation group of the pore sphere inside the tetrahedron, introduce an iterative algorithm for the related nonlinear equation group, and obtain the coordinates and radius of the center of the largest common tangent sphere inside each tetrahedron; Given two adjacent pore spheres with IDs i and j, the corresponding pore sphere radii are r i and r j , d ij is the distance between the centers of two pore balls. According to the distance between two adjacent pore balls, pore channels are divided into three types: Type A pore channel: separated by two pore balls, i.e., d ij ≥r i +r j ; Type B pore channel: two pore spheres intersect, i.e. |r i -r j |≤d ij <r i +r j ; Type C pore channel: The sphere with smaller radius is contained by the sphere with larger radius, i.e. d ij <|r i -r j |.

6. The fluid-solid coupling analysis method based on pore network and CFD-DEM model according to claim 3, characterized in that: The spatial coordinates of the center point in the pore channel after the calculation and division are: The unit direction vector of the pore channel and the spatial coordinates of the pore sphere are calculated, and the spatial coordinates of the center point are obtained respectively according to the difference of the pore channels.

7. The fluid-solid coupling analysis method based on pore network and CFD-DEM model according to claim 1, characterized in that: The determination of the pore channel and the fluid grid area according to the position of the center point includes: When the center point of the pore channel is contained in a certain fluid grid, the corresponding pore channel is completely contained in the fluid grid; In the process of determining the pore channels in the fluid grid, the following conditions should be met: Among them, d c is the minimum calculated width of the fluid flow domain, Δx cfd is the calculation length of the fluid grid, and r is the maximum particle size in the fluid grid.

8. The fluid-solid coupling analysis method based on pore network and CFD-DEM model according to claim 1, characterized in that: The solid particle control equation is established as follows: The movement of particles is governed by Newton's laws, as follows: in, represents the particle movement speed, t represents the numerical calculation time, Represents the particle contact force between particles, including particle-particle contact force and particle-wall contact force, and reflects the interaction between particles through the contact model. Represents the non-contact force between particles (mainly including van der Waals force, electromagnetic force and capillary force), represents the external force, represents the total fluid-particle interaction force acting on the particle, m p represents the particle mass, represents the acceleration due to gravity, represents the particle angular velocity, represents the contact torque on the particle, I p represents the particle moment of inertia.

9. The fluid-solid coupling analysis method based on pore network and CFD-DEM model according to claim 1, characterized in that: The governing equation for the fluid-particle interaction is: Interaction force between fluid and particles It consists of three parts: the drag force of the fluid, the force caused by the pressure gradient of the fluid, and the buoyancy of the fluid on the particle, as follows: According to the coarse-grid fluid-particle coupling numerical analysis system, the fluid drag force is defined as the drag force exerted by the fluid on the particle, always acting about the particle's center of mass, with no rotational torque exerted by the particle, and is calculated as: in, is the drag force on a single particle, is the average porosity of the fluid grid where the particle is located, χ is an empirical correction factor for the local average porosity, this correction term applies the force to high-porosity and low-porosity systems and a wide range of Reynolds numbers; Among them, C d is the resistance coefficient, ρ f is the fluid density, r is the particle radius, is the fluid velocity, is the particle velocity; Among them, the resistance coefficient C d can be calculated as: Among them, Re p is the particle Reynolds number; In addition, the modified empirical factor χ is defined as: Reynolds number p Calculated as: Where: μ f is the dynamic viscosity coefficient of the fluid.

Images