A cross-scale particle flow calculation method based on feature point arrangement

Abstract

The application relates to the field of multiphase numerical simulation, and discloses a cross-scale particle flow calculation method based on feature point arrangement, which comprises the following steps: acquiring the position coordinates, velocity, radius and augmented coefficient of a particle at a current time step, and the position coordinates, volume and velocity of all fluid domain grids; determining the position coordinates of feature points according to the grid volume, particle radius and augmented coefficient; dividing the grid where the feature points are located into an internal grid, a boundary grid and an external grid based on the position relationship between the feature points and the particles, calculating the void fraction of each grid in different ways, and obtaining the background void fraction of the particles; calculating the background flow field velocity of the particles according to the grid void fraction, volume and velocity; and calculating the drag force suffered by the particles according to the background void fraction, velocity and background flow field velocity of the particles, so as to realize numerical calculation of cross-scale particles.

Description

Technical Field

[0001] This application relates to the field of multiphase numerical simulation, and in particular to a method for calculating cross-scale particle flow based on the arrangement of feature points. Background Technology

[0002] Particle-fluid two-phase flow is widely present in natural and engineering fields, and is a complex multi-scale system composed of the interaction between a discrete particle phase and a continuous fluid phase. In deep-sea mining systems, hoisting pipes play a crucial role in transporting ore particles from the seabed to surface vessels. The particulate materials involved in this process exhibit significant multi-scale characteristics, with a wide particle size distribution ranging from a few millimeters to tens of millimeters. Therefore, accurately characterizing the movement mechanism of multi-scale particles within hoisting pipes is of significant scientific and engineering importance for ensuring the efficient and stable operation of the system.

[0003] Currently, numerical simulation methods for this type of problem are mainly divided into two categories: fully analytical methods applicable to large particles and non-analytical methods for fine particles. Fully analytical methods offer high accuracy in handling the interaction between large particles and fluids, but are extremely sensitive to mesh resolution and incur huge computational costs, making them difficult to apply to large-scale multi-particle systems. Non-analytical methods, while efficient in handling small particles, suffer from insufficient flow field analysis and increased deviations in particle force prediction when particle size is comparable to the mesh size. Clearly, existing methods cannot simultaneously address the simulation needs of both large and fine particles within the same framework, limiting their engineering applicability in scenarios with wide particle size distributions. Therefore, a cross-scale computational method applicable to both large and small particles is needed, capable of capturing the perturbation effect of large particles on the flow field while accurately characterizing the in-flow transport characteristics of small particles, thus providing a more reliable numerical simulation tool for the design and operation of deep-sea mining hoisting systems. Summary of the Invention

[0004] To overcome the above problems, this application provides a cross-scale particle flow calculation method based on feature point arrangement. The present invention captures the grid in the augmented region based on the positioning of feature points, and corrects the background porosity and background flow field velocity in the particle drag calculation, thereby realizing the numerical simulation calculation of cross-scale particles.

[0005] According to a first aspect of the embodiments of this application, a method for calculating cross-scale particle flow based on feature point arrangement is provided, comprising:

[0006] S1: Obtain the position coordinates, velocity, radius, and augmentation factor of the particle at the current time step;

[0007] S2: Obtain the position coordinates, volume, and flow velocity of all fluid domain meshes;

[0008] S3: Calculate the coordinates of feature points based on the mesh volume, particle radius and augmentation factor, and then filter out the meshes located in the augmented domain and extract their coordinates, volume and flow velocity.

[0009] S4: Based on the position coordinates of the grid, the volume of the grid, and the position coordinates of the particles within the augmented domain, determine the grid porosity of the grid within the augmented domain in different ways, and calculate the background porosity of the particles;

[0010] S5: Calculate the background flow field velocity of the particles based on the mesh porosity, mesh volume, mesh position coordinates, flow field velocity, particle velocity, and particle position coordinates;

[0011] S6: Calculate the drag force of the particles based on the background porosity of the particles, the velocity of the particles, and the background flow field velocity;

[0012] S7: The drag force of the particle is calculated iteratively to obtain the position coordinates and velocity of the particle in the next time step;

[0013] S8: Add the drag force as a source term to the momentum equation to iteratively calculate the motion of the fluid and obtain the velocity of the fluid domain grid at the next time step.

[0014] Furthermore, it also includes:

[0015] S9: Determine if the current time step is the last time step. If not, return to S1 to continue execution; if yes, terminate the operation process and obtain the final particle position coordinates, velocity, and flow field velocity information of the fluid computational domain mesh.

[0016] Furthermore, the coordinates of feature points are calculated based on the mesh volume, particle radius, and augmentation factor. Mesh elements located within the augmented domain are then selected, and their coordinates, volume, and flow velocity are extracted simultaneously, including:

[0017] S31: Calculate the minimum grid spacing based on the volume of the fluid domain grid;

[0018] S32: Calculate the total number of layers for feature point arrangement based on the minimum grid spacing, particle radius, and augmentation factor. s ;

[0019] S33: Based on the minimum grid spacing, particle radius, and augmentation coefficient, the feature point is calculated. x The radius of the layer;

[0020] S34: Based on the minimum grid spacing and the... x The radius of the layer is calculated layer by layer from the outside in to determine the number of feature points in each layer;

[0021] S35: According to the x The radius of the layer and the total number of feature points are used to calculate the first layer. x The position of each feature point in the layer;

[0022] S36: If x < s ,but x Increment by one and return S33; otherwise, end the loop to obtain the position coordinates of all feature points.

[0023] S37: Filter out the grid where all feature points are located based on their position coordinates, and further obtain the position coordinates, volume and flow velocity of the grid.

[0024] Furthermore, based on the position coordinates of the mesh within the augmented domain, the mesh volume, and the particle position coordinates, the mesh porosity within the augmented domain is determined using different methods, and the background porosity of the particles is calculated, including:

[0025] S41: Based on the position coordinates of the particle and the position coordinates of the grid, obtain the distance between them;

[0026] S42: Obtain the extension distance based on the distance and the volume of the grid;

[0027] S43: Classify the meshes within the augmented domain using the distance, extension distance, and particle radius, and then calculate the porosity of each type of mesh.

[0028] S44: Calculate the background porosity of the particles based on the porosity and volume of the various types of meshes.

[0029] Further, based on the mesh porosity, mesh volume, mesh position coordinates, flow field velocity, and particle velocity and position coordinates, the background flow field velocity of the particles is calculated, including:

[0030] S51: Determine the weight value based on the distance between the position coordinates of the particle and the position coordinates of each grid in the augmented domain, and obtain the final weight coefficient corresponding to each grid through normalization.

[0031] S52: The background flow velocity of the particles is calculated by weighted averaging based on the mesh porosity, mesh volume, and flow field velocity, combined with the mesh weighting coefficient.

[0032] The technical solutions provided by the embodiments of this application may include the following beneficial effects:

[0033] As can be seen from the above embodiments, this application achieves accurate and efficient positioning of the fluid mesh within the augmented domain by radially constructing multiple layers of feature points. For different types of meshes, corresponding porosity calculation methods are adopted, and the particle background porosity and background flow field velocity are corrected, thereby significantly improving the accuracy of drag force calculation. This method solves the problem that existing numerical calculation methods cannot simultaneously handle the calculation of large and small particles, providing an effective solution for the numerical calculation of cross-scale particle two-phase flow.

[0034] It should be understood that the above general description and the following detailed description are exemplary and explanatory only, and do not limit this application. Attached Figure Description

[0035] Figure 1 This is a flowchart illustrating a cross-scale particle flow calculation method based on feature point arrangement according to an exemplary embodiment.

[0036] Figure 2 This is a schematic diagram (two-dimensional) of the arrangement of multi-layer feature points according to an exemplary embodiment.

[0037] Figure 3 This is a schematic diagram (two-dimensional) illustrating grid determination and classification according to an exemplary embodiment.

[0038] Figure 4 This is a schematic diagram (three-dimensional) illustrating the extension distance principle according to an exemplary embodiment.

[0039] Figure 5 This is a schematic diagram (two-dimensional) of a grid corner point according to an exemplary embodiment. Detailed Implementation

[0040] Exemplary embodiments will now be described in detail, examples of which are illustrated in the accompanying drawings. When the following description relates to the drawings, unless otherwise indicated, the same numbers in different drawings denote the same or similar elements. The embodiments described in the following exemplary embodiments do not represent all embodiments consistent with this application. Rather, they are merely examples of apparatuses and methods consistent with some aspects of this application as detailed in the appended claims.

[0041] The terminology used in this application is for the purpose of describing particular embodiments only and is not intended to be limiting of the application. The singular forms “a,” “the,” and “the” used in this application and the appended claims are also intended to include the plural forms unless the context clearly indicates otherwise. It should also be understood that the term “and / or” as used herein refers to and includes any or all possible combinations of one or more of the associated listed items.

[0042] Figure 1This is a flowchart illustrating a cross-scale particle flow calculation method based on feature point arrangement according to an exemplary embodiment, such as... Figure 1 As shown, the method may include the following steps:

[0043] S1: Obtain the position coordinates, velocity, radius, and augmentation factor of the particle at the current time step;

[0044] Specifically, in the process of particle transport in deep-sea mining hoisting pipelines, by calling LIGGGHTS script commands, output files, or secondary development interfaces, parameters such as particle position, velocity, radius, and augmentation coefficient can be obtained in real time. Combining these parameters with corresponding data processing programs (either by calling or writing programs), a large particle population can be filtered and processed to efficiently extract key particle information closely related to the fluid-structure interaction calculations in the hoisting pipeline, ensuring accurate input for subsequent fluid-particle interaction simulations.

[0045] S2: Obtain the position coordinates, volume, and flow velocity of all fluid domain meshes;

[0046] Specifically, using the Open FOAM object library or custom function interface, the computational domain mesh information of the lift pipe fluid is read, including parameters such as mesh location (mesh node coordinates and mesh center coordinates), volume, and flow field velocity. The above data serves as the basic input for particle phase coupling solution and can accurately reflect the characteristics of particle transport and local disturbance within the lift pipe.

[0047] S3: Calculate the coordinates of feature points based on the mesh volume, particle radius, and augmentation factor. Then, filter out meshes located within the augmented domain and extract their coordinates, volume, and flow velocity. (Reference) Figure 2 This step includes the following sub-steps;

[0048] S31: Calculate the minimum grid spacing based on the volume of the fluid domain grid in the lifting pipeline;

[0049] Specifically, based on the discrete grid volume within the flow field computational domain, the volume of each element is calculated, and its cube root is taken as the element's characteristic length. The smallest characteristic length among all elements is selected as the minimum grid spacing, i.e.:

[0050] ;

[0051] In the formula, l Minimum grid spacing m This represents the total number of grid cells. V cell,i For the first iThe volume of each grid cell. Since the grid cell size is not uniform, the feature length of the smallest cell is directly taken as the basis for radial subdivision. This ensures that there will be no missed detections or over-coarsening when distributing feature points, and ensures that the algorithm has universality under both unstructured grids and locally refined grid conditions.

[0052] S32: Calculate the total number of layers for feature point arrangement based on the minimum grid spacing, particle radius, and augmentation coefficient;

[0053] Specifically, based on the particle radius and augmentation factor, combined with the minimum grid spacing, the number of radially augmented layers is calculated, i.e.:

[0054] ;

[0055] In the formula, r Where is the particle radius, k For augmentation coefficients, kr For the radius of the growth domain, s The total number of layers is determined by the ratio of radius to grid scale, ensuring that the number of subdivision layers matches the grid resolution. This avoids both insufficient layers leading to sparse feature points and excessive layers resulting in redundant computation.

[0056] S33: Based on the minimum grid spacing, particle radius, and augmentation coefficient, the feature point is calculated. x The radius of the layer;

[0057] Specifically, starting from the particle surface, the number of layers decreases gradually from the outside in, resulting in the [number of layers]. x The radius of the layer:

[0058] , x= 1, 2, …, s

[0059] In the formula, x Assign a layer number to the current feature point, with an initial value of 1. R x For the first x The radius of the feature points in the layer layout ensures that the interval between each layer and the next is equal to the minimum grid step size. The radial layering is uniform, and the coverage is complete.

[0060] S34: Based on the minimum grid spacing and the... x The radius of the layer is calculated layer by layer from the outside in to determine the number of feature points in each layer;

[0061] Specifically, based on the minimum grid spacing, the number of feature points is determined layer by layer from the outside in. x The formula for calculating the number of feature points in a layer is as follows:

[0062] ;

[0063] In the formula, N This represents the total number of feature points in this layer. The ratio of area to spacing ensures that the distribution density of feature points in each layer adaptively adjusts with changes in radius. A higher number of points in the outer layer and a lower number in the inner layer conforms to geometric principles and avoids wasting computational resources.

[0064] S35: According to the x The radius of the layer and the total number of feature points are used to calculate the first layer. x The position of each feature point in the layer;

[0065] Specifically, no. x The feature point distribution of the layer is as follows:

[0066] ;

[0067] In the formula, x n , y n , z n For the first n The components of the position vector of each feature point The golden ratio, Desirable The golden angle distribution ensures uniform sampling on the sphere, avoiding clustering at the poles. It can quickly generate a near-uniform distribution without complex iterative algorithms, offering high efficiency and numerical stability.

[0068] S36: If x < s ,but x Increment by one and return S33; otherwise, end the loop to obtain the position coordinates of all feature points.

[0069] S37: Based on the position coordinates of all the feature points, obtain the position coordinates, volume, and velocity of the grid where each feature point is located;

[0070] Specifically, based on the obtained coordinates of all feature points, the C++ library provided by Open FOAM is invoked to further process these feature points. This library enables the mapping of feature points to computational grids, thereby extracting the position coordinates, volume information, and local velocity field data of the grid cell containing each feature point, providing fundamental data support for subsequent flow characteristic analysis and result verification.

[0071] S4: Based on the grid position coordinates, grid volume, and particle position coordinates within the augmented domain, determine the grid porosity of the grid within the augmented domain using different methods, and calculate the background porosity of the particles; Figure 3This is a two-dimensional diagram illustrating grid determination and classification according to an exemplary embodiment, which will be explained through the following sub-steps:

[0072] S41: Based on the position coordinates of the particle and the position coordinates of the grid, obtain the distance between them;

[0073] Specifically, it calculates the distance between the particle center and the center of the currently traversed mesh. d 1. The calculation formula is as follows:

[0074] ;

[0075] in, d 1 is the first j The particle and the first i The distance between the coordinates of each grid position x m,i For the first i The position vector of each grid, x p,j For the first j The position vector of each particle;

[0076] S42: Obtain the extension distance based on the distance and the volume of the grid;

[0077] Specifically, based on the distance d 1, and further calculations were performed using the mesh volume to obtain... d 1 extension distance d 2. Figure 4 This is a schematic diagram (3D) illustrating the principle of extension distance according to an exemplary embodiment, and the corresponding calculation formula is as follows:

[0078]

[0079] in, V cell,i For the first i The volume of each grid cell;

[0080] S43: Classify the grids within the augmented domain using the distance, extension distance, and particle radius, and then calculate the porosity of each type of grid separately;

[0081] Specifically, through the distance d 1. Extension distance d 2. The mesh within the augmented domain is classified based on the particle radius, as follows:

[0082] like d 1 + d 2 ≤ r If so, then the grid is defined as an internal grid;

[0083] like d 1 - d 2 ≥ r If so, then the grid is defined as an external grid;

[0084] like d 1 + d 2 > r and d 1 - d 2 < r If so, then the grid is defined as a boundary grid;

[0085] The porosity of the three types of meshes was calculated and assigned values ​​in different ways:

[0086] The initial values ​​for grid porosity are all 0;

[0087] A. Determining the internal grid porosity;

[0088] Based on the specific calculation example, the internal mesh porosity is assigned a value according to the given values;

[0089] B. Determining the porosity of the external grid;

[0090] The porosity of the external mesh is not processed;

[0091] C. Determining the porosity of the boundary grid;

[0092] The porosity of the mesh is calculated based on the geometric relationship between all corner points of the mesh and the particle boundaries;

[0093] For boundary meshes, the solid volume fraction of the mesh is calculated based on the geometric relationship between all corner points of the mesh and the particle boundaries. Figure 5 This is a schematic diagram (two-dimensional) of a grid corner point according to an exemplary embodiment. The grid corner points can be divided into three categories: (I) the corner point and the grid center are both located inside the particle, such as corner points A1 and A2; (II) the corner point is located outside the particle and the grid center is located inside the particle, such as corner points B1 and B2; (III) the corner point is located inside the particle and the grid center is located outside the particle, such as corner point C1.

[0094] The initial value of the solid volume fraction of the boundary mesh is 0, that is... α i =0, if the mesh exists at a Class I corner, then the solid volume fraction is:

[0095] ;

[0096] In the formula, NI The number of Class I corner points N cor This represents the total number of grid corner points. α i For the first i Solid volume fraction of each grid.

[0097] If the grid exists at type II corner points, such as corner points B1 and B2, then the solid volume fraction is:

[0098] ;

[0099] In the formula, x m For point m Coordinates and context x e1 The symbols are similar, representing the coordinates of the corresponding points;

[0100] If the mesh exists at a Class III corner point, such as corner point C1, then the solid volume fraction is:

[0101] ;

[0102] The resulting mesh porosity is:

[0103] ;

[0104] Inside the lifting pipe, particles are often located near the pipe wall or in clustered areas, leading to complex boundary mesh generation. The classification method described above, based on the relationship between corner points and particle boundaries, can more accurately distinguish between internal, external, and boundary meshes, thus obtaining a reasonable porosity distribution. This not only improves the accuracy of fluid-structure interaction calculations in the boundary region but also ensures numerical stability when large and small particles interact within the lifting pipe.

[0105] S44: Calculate the background porosity of the particles based on the porosity and volume of the various types of meshes;

[0106] Specifically, the final determination of the first [number] mesh is based on the mesh porosity and volume using the following formula. j Background porosity of individual particles

[0107] ;

[0108] In the formula, m This represents the total number of grid cells. r Let be the radius of the particle. ε i For the first i The porosity of each grid, ε j For the first jBackground porosity of individual particles V cell,i For the first i The volume of each grid cell is considered; in the lifting pipe, the particle population is unevenly distributed, with significant differences between locally dense and sparse regions. Calculating the background porosity using the above method can more realistically reflect the stress environment of particles in different flow regions, thereby improving the reliability of drag prediction and transport simulation.

[0109] S5: Calculate the background flow field velocity of the particles based on the mesh porosity, mesh volume, mesh position coordinates, flow field velocity, particle velocity, and particle position coordinates;

[0110] S51: Determine the weight value based on the distance between the position coordinates of the particle and the position coordinates of each grid in the augmented domain, and obtain the final weight coefficient corresponding to each grid through normalization processing. w i ;

[0111] ;

[0112] In the formula, x p The coordinates of the particle center are, x i Here, represents the grid position coordinates, and m represents the total number of particles. Within the lift pipe, particles at different locations are significantly affected by the local flow field. Distance-weighted normalization highlights the influence of neighboring grids and avoids interference from distant grids, thus more accurately reflecting the flow field around the particles and improving computational accuracy.

[0113] S52: The background flow velocity of the particles is calculated by weighted averaging based on the mesh porosity, mesh volume, and flow field velocity, combined with the mesh weighting coefficient.

[0114] Specifically, based on the porosity, volume, and flow field velocity of the mesh within the augmented range, the background flow field velocity of the particles is calculated using a weighting function, as shown in the following formula:

[0115] ;

[0116] In the formula, For the first j The background flow velocity of each particle ε i For the first i The porosity of each grid, V cell,i For the first i The volume of each grid, U f,i For the first iThe flow field velocity is calculated for each grid cell. In the lift pipe, local velocity differences are significant. Weighted averaging integrates grid porosity, volume, and velocity, avoiding single-point errors and more accurately reflecting the flow field around the particles. Therefore, this method improves the stability and accuracy of background flow field velocity calculations.

[0117] S6: Calculate the drag force of the particles based on the background porosity of the particles, the velocity of the particles, and the background flow field velocity;

[0118] Specifically, the Gidaspow drag model is used. The Gidaspow drag model calculates the drag force on the particles based on the Gidaspow correlation coefficient, taking into account porosity, background flow field velocity, and particle state. The formula is as follows:

[0119] ;

[0120] ;

[0121] ;

[0122] ;

[0123] ;

[0124] In the formula, V p For particle volume, d p Particle size, F d,j For the first j The drag force of each particle β d The momentum exchange coefficient, ε i For the first i The porosity of each grid, ε j For the first j Background porosity of individual particles For the first j The background flow velocity of each particle U p For particle velocity, ρ f For fluid density, C d Re is the drag coefficient. p The particle Reynolds number, μ f This refers to dynamic viscosity.

[0125] In deep-sea mining hoisting pipelines, large particles significantly disturb the surrounding flow field and alter the drag force distribution due to locally high Reynolds numbers. The Gidaspow drag force model can uniformly describe particle forces under different porosities, applicable to high-speed transport in sparse regions and covering group motion in dense regions. This accurately reflects the force state of particles at multiple scales, improves the accuracy of drag force calculations, and provides support for pipeline transport analysis and equipment design.

[0126] S7: The drag force of the particle is calculated iteratively to obtain the position coordinates and velocity of the particle in the next time step;

[0127] Specifically, LIGGGHTS uses Newton's second law to integrally solve for particle motion, considering both inter-particle collisions and impacts with the pipe wall. In deep-sea mining hoisting pipes, this approach accurately reflects the dynamic characteristics of large particles as they ascend vertically, experiencing friction, collisions, and mass accumulation against the pipe wall. During this process, fluid forces and other external forces act on the particles, dynamically updating their motion state. This allows for synchronized updates of particle position and velocity at each time step, ensuring consistency between particle dynamic evolution and changes in the flow field within the hoisting pipe, thus providing reliable support for accurate calculations using the CFD–DEM coupled model.

[0128] S8: Add the drag force as a source term to the momentum equation to iteratively calculate the motion of the fluid and obtain the velocity of the fluid domain grid in the next time step;

[0129] Specifically, in Open FOAM, the force exerted by particles on the fluid is uniformly transformed into a momentum source term and incorporated into the governing equations to maintain overall momentum conservation. Through coupled solution with the fluid momentum equation, a two-way interaction between particles and fluid is achieved. During numerical iteration, the modified momentum equation is used for calculation to obtain an updated fluid velocity field, which is then used to correct the mesh information of the entire computational domain. This step ensures the synchronous evolution of the fluid field and particle trajectories within the lifting pipe, providing flow field data support for subsequent transport prediction and analysis.

[0130] To ensure the accuracy and consistency of the results obtained from iterative iterations across multiple time steps, this method also includes:

[0131] S9: Determine if the current time step is the last time step. If not, return to S1 to continue execution; if yes, terminate the operation process and obtain the final particle position coordinates, velocity, and flow field velocity information of the fluid computational domain mesh.

[0132] Specifically, after each round of calculation, the system increments the time step count. When the accumulated time step count has not reached the preset total number of steps, it returns to step S1 to continue updating and coupling the particle information and fluid domain mesh information; when the accumulated time step count reaches the preset total number of steps, the loop ends and the final state data of the particles and the mesh results of the fluid computational domain are output.

[0133] As demonstrated by the above embodiments, this invention achieves accurate and efficient positioning of the fluid mesh within the augmented domain by constructing feature points through radial layering. For different mesh types, corresponding porosity processing methods are adopted, and the particle background porosity and background flow field velocity are corrected, thereby effectively improving the accuracy of drag force calculation. This method solves the problem that existing numerical calculation methods cannot simultaneously handle the calculation of large and small particles, realizing effective numerical simulation of cross-scale particle two-phase flow.

[0134] Other embodiments of this application will readily occur to those skilled in the art upon consideration of the specification and practice of the disclosure herein. This application is intended to cover any variations, uses, or adaptations of this application that follow the general principles of this application and include common knowledge or customary techniques in the art not disclosed herein. The specification and embodiments are to be considered exemplary only, and the true scope and spirit of this application are indicated by the claims.

[0135] It should be understood that this application is not limited to the precise structure described above and shown in the accompanying drawings, and various modifications and changes can be made without departing from its scope. The scope of this application is limited only by the appended claims.

Claims

1. A method for calculating cross-scale particle flow based on feature point arrangement, characterized in that, include: S1: Obtain the position coordinates, velocity, radius, and augmentation factor of the particle at the current time step; S2: Obtain the position coordinates, volume, and flow velocity of all fluid domain meshes; S3: Calculate the coordinates of feature points based on the mesh volume, particle radius and augmentation factor, and then filter out the meshes located in the augmented domain and extract their coordinates, volume and flow velocity. S4: Based on the position coordinates of the grid, the volume of the grid, and the position coordinates of the particles within the augmented domain, determine the grid porosity of the grid within the augmented domain using different methods, and calculate the background porosity of the particles; specifically including: S41: Obtain the distance between the position coordinates of the particles and the position coordinates of the grid; S42: Obtain the extension distance based on the distance and the volume of the grid; S43: Classify the grids within the augmented domain using the distance, the extension distance, and the radius of the particles, and then calculate the porosity of each type of grid; S44: Calculate the background porosity of the particles based on the porosity and volume of each type of grid; S5: Calculate the background flow field velocity of the particles based on the mesh porosity, mesh volume, mesh position coordinates, flow field velocity, particle velocity, and particle position coordinates; S6: Calculate the drag force of the particles based on the background porosity of the particles, the velocity of the particles, and the background flow field velocity; S7: The drag force of the particle is calculated iteratively to obtain the position coordinates and velocity of the particle in the next time step; S8: Add the drag force as a source term to the momentum equation to iteratively calculate the motion of the fluid and obtain the velocity of the fluid domain grid at the next time step.

2. The method for calculating cross-scale particle flow based on feature point arrangement according to claim 1, characterized in that, Also includes: S9: Determine if the current time step is the last time step. If not, return to S1 and continue execution. If so, the operation process is terminated, and the final particle position coordinates, velocity, and flow field velocity information of the fluid computational domain mesh are obtained.

3. The method for calculating cross-scale particle flow based on feature point arrangement according to claim 1, characterized in that, The coordinates of feature points are calculated based on the mesh volume, particle radius, and augmentation factor. Mesh elements located within the augmented domain are then selected, and their coordinates, volume, and flow velocity are extracted simultaneously, including: S31: Calculate the minimum grid spacing based on the volume of the fluid domain grid; S32: Calculate the total number of layers for feature point arrangement based on the minimum grid spacing, particle radius, and augmentation factor. s ; S33: Based on the minimum grid spacing, particle radius, and augmentation coefficient, the feature point is calculated. x The radius of the layer; S34: Based on the minimum grid spacing and the... x The radius of the layer is calculated layer by layer from the outside in to determine the number of feature points in each layer; S35: According to the... x The radius of the layer and the total number of feature points are used to calculate the first layer. x The position of each feature point in the layer; S36: If x < s ,but x Increment by one and return S33; otherwise, end the loop to obtain the position coordinates of all feature points. S37: Filter out the grid where all feature points are located based on their position coordinates, and further obtain the position coordinates, volume and flow velocity of the grid.

4. The method for calculating cross-scale particle flow based on feature point arrangement according to claim 1, characterized in that, Based on the mesh porosity, mesh volume, mesh position coordinates, flow field velocity, and particle velocity and position coordinates, the background flow field velocity of the particles is calculated, including: S51: Determine the weight value based on the distance between the position coordinates of the particle and the position coordinates of each grid in the augmented domain, and obtain the final weight coefficient corresponding to each grid through normalization. S52: The background flow velocity of the particles is calculated by weighted averaging based on the mesh porosity, mesh volume, and flow field velocity, combined with the mesh weighting coefficient.

Images