A parallel computing method for gas-solid simulation

Abstract

The application provides a parallel computing method for gas-solid simulation, which comprises the following steps: forming a geometric boundary as a solid by building a plurality of static geometric shape units according to a multiphase complex boundary condition; setting an initial state of a system to be simulated, and adopting a spherical or ellipsoidal particle model for a gas particle; performing spatial decomposition on the system to be simulated to obtain a plurality of sub-regions, each of which comprises an interior and a boundary, the interior is simulated by a hard sphere, and the boundary is simulated by a pseudo-particle, each of the sub-regions is divided into a plurality of grids, a mapping relationship between the particles and the grids is established, and the geometric boundary in each of the grids is calculated; within a set time step, updating particle information based on event driving and in view of particle collision in the interior of each of the sub-regions and particle collision between the particles and the geometric boundary; updating the particle information based on time driving and in view of particle collision on the boundary of each of the sub-regions and particle collision between the particles and the geometric boundary; and outputting a simulation result.

Description

Technical Field

[0001] This invention relates to the field of parallel computing for discrete particle simulation, and particularly to a parallel computing method for gas-solid simulation. Background Technology

[0002] Gas-solid two-phase flow has received widespread attention as a research hotspot in experiments and simulations, such as the flow of aircraft in hypersonic gases and gas-solid diffusion in fluidized beds. Due to the high cost of experiments, simulation methods such as computational fluid dynamics and direct Monte Carlo simulations provide an effective supplementary means to understand the essential laws governing these processes and to design materials with superior performance. For example, in wind tunnel experiments, by simulating the drag coefficients of solids such as spheres, cones, and planes around gases, the essential laws of flow can be understood, guiding the design of aircraft.

[0003] For gas-solid two-phase flows, more accurate molecular dynamics simulations have made remarkable progress in complementing experiments in the study of gaseous material diffusion and transport processes at the nanoscale and microscale. However, traditional molecular dynamics based on soft-sphere potential models wastes a lot of resources on ineffective calculations and updates when simulating low-density substances such as gases due to the limited number of intermolecular collisions in high Knudsen number systems. Hard-sphere molecular dynamics simulations are more efficient in handling such problems, but hard-sphere models are event-driven, making parallelization difficult. Achieving parallelization through complex methods significantly reduces computational efficiency, negating the simulation advantages of hard-sphere models. As an efficient simulation method, the hard-sphere-pseudo-particle coupling model (HS-PPM) can achieve parallel simulation of hard-sphere models. However, in practical applications, it is no longer possible to simulate simple single-phase fluid flows but requires various complex internal solid boundary shapes. Existing HS-PPM achieves different boundary conditions by stacking frozen particles, which is a relatively coarse method. When simulating flow around a sphere, the event boundary of the geometry is used to describe the sphere, but this simulation is relatively simple and lacks a systematic approach to handling various geometries and complex boundaries.

[0004] Meanwhile, real gas molecules are mostly composed of multiple atoms, while existing HS-PPM systems can only handle spherical gas particles. Using spherical particles to simulate the diffusion and flow behavior of multi-atomic gas molecules often introduces errors. To more accurately simulate gas collision behavior, ellipsoids can be used instead of ordinary spherical particles. By adjusting the lengths of the three semi-axis, high-precision simulations that closely approximate the morphology of real gas molecules can be achieved. Event-driven collision simulations between ellipsoids can determine the positional relationship between the two ellipsoids based on the sign change of the zero-point roots of the characteristic equation formed by the characteristic matrices of the two ellipsoids, and the velocity after the collision can be handled through conservation relationships.

[0005] However, the event-driven ellipsoidal algorithm has a large computational load, and the difficulty in parallelizing the hard sphere model also limits the size of the simulation system. Therefore, parallelizing the hard sphere model of the ellipsoid has become an urgent problem for those skilled in the art to solve. Summary of the Invention

[0006] The purpose of this invention is to overcome the difficulty of parallelizing event-driven geometric boundary and ellipsoidal particle collision algorithms, and to provide a parallel computing method for gas-solid simulation of complex flow simulation of multiphase discrete particle systems that is closer to the real system.

[0007] To achieve the above objectives, this invention proposes a parallel computing method for gas-solid simulation, the method comprising:

[0008] Step 1) Based on the multiphase complex boundary conditions, a geometric boundary is constructed from several static geometric units to form a solid.

[0009] Step 2) Set the initial state of the system to be simulated according to the simulation settings information. The gas particles are modeled as spherical or ellipsoidal particles.

[0010] Step 3) Decompose the system to be simulated into multiple sub-regions. Each sub-region includes the interior and the boundary. Hard sphere simulation is used for the interior and quasi-particle simulation is used for the boundary. Each sub-region is divided into several grids. Establish the mapping relationship between particles and grids and calculate the geometric boundary in each grid.

[0011] Step 4) Within the set time step, based on event-driven updates, update particle information for particle collisions within each sub-region and collisions between particles and geometric boundaries; based on time-driven updates, update particle information for particle collisions at the boundaries of each sub-region and collisions between particles and geometric boundaries.

[0012] Step 5) Determine whether the simulation termination condition has been met. If yes, proceed to step 6); otherwise, proceed to step 4.

[0013] Step 6) Output the simulation results.

[0014] As an improvement to the above method, the multiphase complex boundary includes one or more of the following: conical boundary, triangular surface element boundary, spherical boundary and cylindrical boundary; the geometric shape element includes one or more of the following: sphere, cone, tube, simple plane and triangular surface element; the characteristic parameters of the geometric shape element include shape feature description parameters and geometric shape range parameters.

[0015] As an improvement to the above method, the gas particles are outside the boundary condition topology and do not overlap with the boundary.

[0016] As an improvement to the above method, step 4) further includes the following:

[0017] Based on the position, velocity, and geometric boundary of the particles within the sub-region, a hard sphere region event list is established for that sub-region. The hard sphere region event list includes events that occur in chronological order and their predicted occurrence times. The events include collision events between particles within the sub-region and the geometric boundary, collision events between particles within the sub-region, and events where particles within the sub-region cross the mesh.

[0018] As an improvement to the above method, when the event is a collision event between a particle and a geometric boundary within a sub-region, the corresponding predicted collision time is calculated by sequentially combining the particle centroid motion equation with the analytical structure equation of each geometric unit, using analytical or numerical iteration methods.

[0019] As an improvement to the above method, when the event is a collision event between particles within a sub-region, the gas particles adopt an ellipsoidal particle model, and an external tangent sphere of the ellipsoidal particle is introduced. If there is no collision between the external tangent spheres, there is also no collision between the ellipsoidal particles. If there is a collision between the external tangent spheres, the possible collision time interval [t1, t2] of the ellipsoidal particles is determined by the time interval between the two tangent spheres. The minimum time t∈[t1, t2] that satisfies the tangency of the two ellipsoidal particles is obtained by the Brent numerical iteration method. This time t is the collision time, and then the velocity and position of the ellipsoidal particles are updated.

[0020] As an improvement to the above method, step 4) updates particle information based on event-driven mechanisms within a set time step, specifically addressing particle collisions within each sub-region and collisions between particles and geometric boundaries; this includes:

[0021] Within a set time step, read the event list of the hard sphere region, obtain the earliest occurring event, and determine whether it is a real event. If not, obtain the next event record. If it is, then based on the type of the event record: if the event is a collision between a particle and the geometric boundary within the sub-region, determine the particle's velocity after the collision based on the analytical geometric expression of the particle's position coordinates and the geometric boundary, and then update the particle's position and velocity; if the event is a collision between particles within the sub-region, calculate the translational and rotational velocities based on the conservation of momentum, angular momentum, and total energy, and then update the particle's position and velocity; if the event is a particle crossing the mesh within the sub-region, update the mapping relationship between the particle and the mesh, and then update the particle's position and velocity.

[0022] Delete the events that have already occurred from the hard ball region event list, and regenerate the hard ball region event list corresponding to the sub-region based on the updated particle position and velocity, combined with the geometric boundary.

[0023] As an improvement to the above method, step 4) is time-driven and employs a pseudo-particle model to update particle information for particle collisions at the boundaries of each sub-region and for collisions between particles and geometric boundaries; specifically, it includes:

[0024] For collisions between particles at the boundaries of sub-regions, a pseudo-particle model is used to update the particle velocity and position. When gas particles are modeled as ellipsoidal particles, the characteristic matrices of the two ellipsoids are obtained and then substituted into the characteristic sub-structure formula to determine whether the two ellipsoids overlap. If the two ellipsoids overlap, the tangent point is determined by the scaling method, and then the particle velocity and position after the collision rebound are obtained.

[0025] For particles and geometric boundaries at the sub-region boundary, the relative positional relationship between them is determined by analytical or numerical iteration methods based on the particle coordinates and the analytical structure equations of each geometric shape element of the corresponding boundary. If they are disjoint, no collision will occur and no action will be taken. If they are tangent or intersecting and tend to approach the interior of the geometric boundary, a collision will occur. The velocity after the collision is calculated, and the position of the particle is updated.

[0026] As an improvement to the above method, if the two ellipsoids overlap, the tangent point is determined by the scaling method, and then the velocity and position of the particles after the collision and rebound are obtained. Specifically, this includes:

[0027] The two ellipsoidal particles are reduced in size using the scaling method. The scaling factor that makes the two ellipsoids tangent is obtained by numerical solution using Newton's iteration method. Then, the velocity and position after the collision are calculated using the reduced tangent ellipsoid.

[0028] As an improvement to the above method, the geometric unit also includes stacked frozen particles.

[0029] Compared with the prior art, the advantages of the present invention are:

[0030] 1. Compared with traditional event-driven geometric boundary conditions, this algorithm provides a simple, easy-to-implement, and highly efficient parallelization method, realizing the parallelization of event-driven geometric boundary conditions that are difficult to parallelize. It can handle boundaries of arbitrary shapes and the boundary conditions can span multiple processes, improving the applicability of event-driven hard sphere simulation and expanding its application scenarios from single-phase simple flow to multi-phase complex boundary flow.

[0031] 2. Compared with the traditional stacked frozen particle boundary conditions, this algorithm has more accurate boundary representation due to the introduction of analytical description of geometric boundaries, and is more efficient for simple boundaries. Furthermore, by combining geometric boundaries with traditional stacked frozen particle boundaries, a more general combined boundary condition is achieved.

[0032] 3. Compared with traditional time-driven geometric boundary conditions, this algorithm only uses time-driven conditions at the boundaries of sub-regions, while using event-driven conditions for most of the internal regions. Since hard sphere models are often used for rarefied gas flows, event-driven conditions are much more efficient than time-driven conditions in such low-density simulations. Therefore, the simulation efficiency of this algorithm is also higher than that of time-driven geometric boundary conditions. Attached Figure Description

[0033] Figure 1 This is a schematic diagram of domain decomposition and computational grid generation;

[0034] Figure 2 This is a flowchart of the parallel computing method for gas-solid simulation of the present invention;

[0035] Figure 3 This is a schematic diagram of the particles and the cone in the flow around the cone in Example 1;

[0036] Figure 4(a) shows the cone collision partitioning configuration in Example 1;

[0037] Figure 4(b) shows the specific meaning of Figure 4(a);

[0038] Figure 5 This is a schematic diagram showing the positional relationship of the boundary points of each region of the cone in Example 1;

[0039] Figure 6 This is the pseudocode from Example 1 for determining the relationship between the position of particle P and the cone;

[0040] Figure 7 This is a schematic diagram of the rotation algorithm in Example 1;

[0041] Figure 8 This is a graph showing the relationship between the traction coefficient and wall temperature and thermal adaptability coefficient in Example 1;

[0042] Figure 9 This is a schematic diagram of the boundary composition of the triangular facets in Example 2;

[0043] Figure 10(a) is a top plan view of the collision partition of the triangular facet boundary in Example 2;

[0044] Figure 10(b) is a front view of the vertical cross-section of the collision partition of the triangular surface element in Example 2;

[0045] Figure 10(c) shows the specific meaning of Figures 10(a) and 10(b);

[0046] Figure 11 This is a comparison chart of the test and experimental results of the conical boundary and triangular surface element boundary in Example 2;

[0047] Figure 12 This is a schematic diagram of the typical combined geometric boundary of Example 3;

[0048] Figure 13 This is a schematic diagram of the ellipsoidal collision in Example 4;

[0049] Figure 14 This is the pseudocode for the overlapping ellipsoid scaling method in Example 4;

[0050] Figure 15(a) shows the number density distribution of the ellipsoidal particle cone around the flow in Example 4;

[0051] Figure 15(b) shows the dimensionless temperature distribution of the flow number around the ellipsoidal particle cone in Example 4. Detailed Implementation

[0052] The technical solution adopted in this invention includes the following steps:

[0053] Step (1) Based on the multiphase complex boundary conditions to be constructed, select appropriate basic building blocks to build. It is easy to see that even for very complex boundaries, they can be constructed by triangular surface boundary elements or frozen particle combinations.

[0054] Step (2) Parallel task partitioning, generating initial particles, ensuring that the particles are outside the boundary condition topology and do not overlap with the boundary;

[0055] Step (3) Divide the region into a grid, such as Figure 1 As shown, taking two dimensions as an example, the entire simulation region is decomposed into several rectangular sub-regions. Each process processes one sub-region. For each sub-region, a computational mesh is created, and the basic constituent elements of the boundary conditions that may be involved within the mesh are calculated. The boundary mesh of each sub-region is set to a time-driven quasi-particle mode, and the internal mesh is set to an event-driven hard-sphere mode. Particles move and collide according to the mode of their respective mesh. Each process constructs an initial event list for the event-driven hard-sphere mesh within its processed sub-region. In addition to predicting the mesh crossing time and particle collision time, the process also predicts the collision time between the particle and all the basic constituent elements of the geometric boundary within the corresponding mesh. The specific steps are as follows: Figure 2 As shown.

[0056] Step (4) For each pseudo-particle time step, each process updates its own hard sphere region events. In addition to the traditional mesh crossing and particle collision events, it adds the geometric boundary basic constituent unit events recorded in the corresponding mesh for particle collision. Then, it updates all events of the corresponding mesh and related meshes and determines the next earliest occurrence event. This step is repeated until all events in this time step are processed.

[0057] Step (5) updates the particle position status of the pseudo-particle region, handles particle collisions within and between sub-regions, and on this basis, handles the basic element collisions between the particles and the geometric boundaries of the corresponding pseudo-particle region mesh records; then repeats step (4) until the program runs out.

[0058] Preferably, the basic constituent units of the geometric boundary in step (1) include, but are not limited to, conical boundaries, triangular surface boundaries, spherical boundaries, cylindrical boundaries, etc., and complex geometric boundaries are achieved by combining the above basic constituent units to limit the range of action.

[0059] Preferably, the basic constituent unit expressed in step (1) is generally expressed in the form of analytical geometry. Its basic structure and boundary composition can be described in the form of piecewise functions. For example, the boundary of a cone can be described by multiple analytical formulas for the cone surface, the bottom surface and the vertex part respectively. The boundary of a triangular facet can be described by multiple analytical formulas for the plane, the three sides of the plane and the three vertices respectively.

[0060] Preferably, the various basic building blocks described in step (1) contain simple three-dimensional upper and lower bounds to limit the scope of use of each basic structure, thereby achieving the construction of more complex boundaries.

[0061] Preferably, the boundary described in step (1) can be coupled with a stacked frozen particle boundary. The stacked frozen particles are used to process regions that are difficult to construct using basic building blocks that are difficult to express analytically, thereby achieving a more efficient simulation.

[0062] Preferably, the particles in steps (2)-(5) can be spherical particles or ellipsoidal particles.

[0063] Preferably, the mesh division in step (3) can be a cuboid mesh with three unequal directions. Whether a mesh records a constituent unit is determined by the relative relationship between the mesh center and the bounding box of each basic constituent unit on the geometric boundary. The bounding box of each basic constituent unit needs to be extended by at least twice the maximum mesh width based on the upper and lower bounds in the three dimensions.

[0064] Preferably, the collision time of the predicted particle with all basic constituent units of the corresponding geometric boundary in step (3) can be obtained by solving the particle's centroid motion equation and the analytical structure equation of each basic constituent unit in sequence. The collision time can be obtained by analytical method or numerical iteration. For particles that overlap with the geometric boundary, if the particle velocity points into the geometric shape, the collision event is executed immediately. If the particle velocity points out of the geometric shape, no processing is performed.

[0065] Preferably, in step (3), the particle collision events with all geometric boundaries and the mesh crossing events within the corresponding mesh are of the same level and jointly participate in the sorting of the earliest occurring events within the corresponding mesh.

[0066] Preferably, in step (4), the collision between the particle and the corresponding geometric boundary must ensure that the particle's velocity after the collision is far away from the boundary condition.

[0067] Preferably, step (4) of determining the next earliest occurrence event can be performed using a complete binary tree to improve computational efficiency.

[0068] Preferably, in step (5), the particles and the corresponding pseudo-particle region grid records of geometric boundaries basically constitute unit collisions. First, the relative positional relationship between the two is determined according to the analytical geometric expression of the particle coordinates and the corresponding boundary. If they are separate, no operation is required. If they are tangent or intersecting and have a tendency to approach the interior of the geometric boundary, then the collision is processed.

[0069] Preferably, the collision mode between the particles and the geometric boundary in step (5) is the same as that in the hard sphere region.

[0070] The technical solution of the present invention will be described in detail below with reference to the accompanying drawings and embodiments.

[0071] Example 1

[0072] This embodiment provides an event-driven parallel method for conical boundary conditions, which can be applied to accurate flow simulation around a cone in wind tunnel experiments. The specific process is as follows:

[0073] A schematic diagram of the conical boundary is shown below. Figure 3 As shown, the equation of the cone is:

[0074] (yb) 2 +(zc) 2 =(xa) 2 tan 2 θ(a<x<a+L)

[0075] In the above formula, L is the cone height, θ is the semi-cone angle, and O(a, b, c) is the cone vertex. The cone height is parallel to the x-axis.

[0076] Because the collision modes of particles differ in different regions of the cone, the area around the cone is divided into five regions, as shown in Figure 4(a). The thick black solid line represents the actual cross-section of the cone's center, the thin black solid line represents the cross-section of the cone extending by one particle radius, the thick gray solid line represents the furthest position that the particle's centroid can reach based on the extended black cone, and the dashed lines are auxiliary lines. The cone is divided into five regions according to the collision modes, and the meaning of each part is shown in Figure 4(b). Specifically, if the particle's centroid is located in region R11, R2, R31, or R4, it will collide with the cone.

[0077] To distinguish between different regions, the coordinates of the boundary points of each region in the x-direction are marked, such as... Figure 5 As shown, the upper half of the cross-sectional view of the cone's center section corresponds to the horizontal coordinates x0 to x10 of each region's boundary. s .

[0078] The pseudo-particle part focuses on determining the relationship between the position of particle P and the cone. If the particle is outside the cone, no processing is required. The pseudocode for determining the positional relationship is as follows: Figure 6 As shown.

[0079] In this algorithm, to determine the positional relationship between particle P and the base of the cone, let point P and the height L of the cone form a plane. The point D closest to point P is the intersection of this plane and the base of the cone. Then, rotate point P along the height of the cone to point RP in the positive y-direction, as follows. Figure 7 As shown, point D is also rotated by the same angle to point RD. The rotation method facilitates solving the relationship between the particle and the cone. The coordinates of point RD can be obtained by rotating in the opposite direction to obtain the coordinates of point D. The expressions for the rotated coordinates and the rotation angle are easily obtained as follows:

[0080]

[0081]

[0082] For the hard ball part, the focus is on predicting the collision time between the particle and the cone in each region. The algorithm sequentially checks whether a particle collides with the cone in regions R5, R10, and R30, and inserts events accordingly. For particles that may traverse multiple regions, the algorithm uses a progressive approach, checking for potential collisions with each region from the outside in.

[0083] It is easy to see that predicting the collision time of particle P with the circular base of the cone in region R4 is relatively simple. The collision with the cone apex in region R1 is also relatively simple, similar to the formula for the collision between two particles. Predicting the collision between the particle and the cone surface in region R2 can be seen as predicting the collision time between the particle's centroid and the outer cone. If a collision occurs, the particle's centroid must pass through the outer cone surface. Therefore, we can obtain a quadratic equation in time t, and take the minimum solution within the range of [x2, x4] for the corresponding collision time. If particle P collides with the bottom ring in region R3, we can construct a quartic equation in time t similar to the collision with the cone surface, solve this equation, and take the minimum solution within the range of [x3, x5] for the corresponding collision time.

[0084] For the coupled algorithm, a grid-marked collision pattern is used to indicate the region where the particle is located. When the particle is located in the inner region grid, the collision between the particle and the cone boundary is automatically updated according to the hard sphere region collision pattern. When the particle is located in the process boundary grid, the collision between the particle and the cone boundary is automatically updated according to the pseudo-particle region collision pattern. It is easy to see that since the cone boundary conditions do not move, the collision between the particle at the process boundary and the cone boundary with the same number in the adjacent process is the same as that in the current process, and no separate processing is required.

[0085] Figure 8One of the simulation results is presented. Wind tunnel experiments of supersonic rarefied flow often use cold walls. During testing, it was found that for cold walls, the drag coefficient does not change monotonically with the thermal adaptability coefficient. For a cone flow test with a Reynolds number of 126 and a Knudsen number of 0.278, the ratio T of the wall temperature to the infinite inflow temperature was varied. w / T ∞ The coefficients were set to 0.01, 0.1, and 1. The thermal adaptability coefficient was changed from 0.7 to 1.0, and the traction coefficient was tested to show changes with wall temperature and thermal adaptability coefficient. It can be seen that T... w / T ∞ When the value is 0.01, the traction coefficient is the largest when the thermal adaptability coefficient is 0.9. As the wall temperature increases, the traction coefficient gradually transitions to a monotonic change with the change of the thermal adaptability coefficient. When the wall temperature is equal to the infinite flow temperature, the maximum traction coefficient is obtained when the thermal adaptability coefficient is 1.0 (complete diffuse reflection).

[0086] Example 2

[0087] This embodiment provides an event-driven parallel method for triangular element boundaries, applicable to the simulation of aircraft flow in hypersonic rarefied gases. The specific process is as follows:

[0088] Triangular facets simplify complex walls into easily describable, regular walls. The basic elements of a triangular facet are described as follows: Figure 9 As shown.

[0089] Each triangular element contains three vertices n. k i The incenter of the triangular element is c. k The normal vector is no k o Then the three edge vectors are:

[0090]

[0091]

[0092]

[0093] The perpendicular vector from the incenter of the triangular element to the three sides is:

[0094]

[0095] The center of mass is located at c s And the radius is r s The relative vector between the sphere and the surface element is:

[0096] v C =c S -c K

[0097] To facilitate differentiation, the triangular facet is divided into 7 parts, as shown in Figure 10. Figure 10(a) is a top view of the plane containing the triangular facet, and Figure 10(b) is a cross-sectional view in the vertical direction corresponding to the thick gray dashed line in the left figure, which is equivalent to the front view of Figure 10(a). In both figures, the thick black solid lines represent the actual triangular facet, and the thin black solid lines represent the area extending outward from the facet with a radius r equal to one particle radius. s The finite large triangular prism is constructed, with the thick gray solid lines representing the positions where the centroid of the particles does not collide with the face elements on the corresponding cross-sectional view, and the dashed lines being auxiliary lines. Some areas marked in the two figures overlap. The meaning of each part is shown in Figure 10(c), where the units digit x takes values ​​of 0, 1, and 2, representing the numbers of the three sides or three vertices of the triangular face element, respectively. When dividing the region, the entire space is sequentially cut and divided from top to bottom in the table, like a sculpture; that is, each row in the table is operated on after taking the complement of the regions of all preceding rows within the space.

[0098] For pseudo-particle regions, it is necessary to determine the region where the particles are located and then calculate the normal vector.

[0099] The distance l between the particle's centroid and the triangular element is shown in the following formula:

[0100] a = v C (no K ) o

[0101] l = |a| - r S

[0102] It is easy to see that if l > 0, then the particle is located in region R7. The particle does not overlap with the triangular element.

[0103] If l < 0, then calculate the projection x of the particle on the plane. c As shown in the formula below:

[0104] x C =c S -a(no K ) o

[0105] Then the relative vector of this vector with respect to the incenter of the triangular element is vx. c for:

[0106] vx C =x C -c K

[0107] Then, it is necessary to determine the particle projection x. c Which edge of the triangular element forms the subregion? This is determined by calculating vx. c With OG k The included angle acos(b) iThe smallest included angle is the region where the particle is located, and b i The calculation is shown in the following expression:

[0108]

[0109] The region where the particle is located is determined by the smallest included angle, iSub. This is achieved through vx. c In this area og k The projection on the surface is used to determine whether it overlaps with the region. If the projection satisfies the following formula, then the particle is located in region R1, and the particle overlaps with the triangular element.

[0110]

[0111] If the projection satisfies the following formula, then the particle is located in the R6x region. The particle does not overlap with the triangular element.

[0112]

[0113] In other cases, the edges and corners need to be further subdivided. First, calculate vx. c The projection vector x on the corresponding edge of the region c(new) for:

[0114]

[0115] relative vector vx c(new) for:

[0116] vx C(new) =x C(new) -c K

[0117] The center of the particle and the projection of the particle onto the edge x c(new) The distance between them is tnl:

[0118] tnl=|x C(new) -c S |

[0119] If tnl>r s If the particle is in the R3x region, then the particle does not overlap with the triangular element.

[0120] definition:

[0121]

[0122]

[0123] like

[0124]

[0125] The particle is located in the R2x region and overlaps with the triangular element.

[0126] Finally, calculate the distance between the particle and its two adjacent vertices. If the distance is less than r... s If it is located in the R4x region, the particle overlaps with the triangular element.

[0127] Otherwise, the particle is located in the R5x region and does not overlap with the triangular element.

[0128] For the hard spherical region, similar to the flow around a cone, a progressive algorithm is still used. However, unlike the flow around a cone, since direct judgment is too complex, if a particle traverses multiple regions, the particle is moved sequentially until it collides with the boundary of the triangular element.

[0129] If the particle is located in region R7 outside the plane, then the collision time between the particle and the plane needs to be determined, and the particle should be moved to a distance r from the plane. s By considering the position of the particle and the collision of the surface element with respect to time t, we can obtain a quadratic equation. Solving this equation, the smallest root greater than 0 is the collision time.

[0130] If the particle is located in the R6x or R3x region, i.e., the sub-region outside the edge, then the particle needs to be moved to a distance r from the edge first. s The position is similar to the collision with the surface, and the particle motion and the collision with the edge are combined into a quadratic equation with respect to time t.

[0131] If the particle is located in the R5x region, i.e., the sub-region inside the outer corner line, then the particle needs to be moved to a distance of angle r first. s Given the location, it's easy to predict the collision time between a particle and an angle, which is similar to predicting the collision time between two particles, making it relatively simple.

[0132] For the coupled algorithm, a mesh-marked collision pattern is used to indicate the region where the particle is located. When the particle is located in the internal mesh region, the collision pattern between the particle and the triangular element is automatically updated according to the hard sphere region collision pattern. When the particle is located in the mesh at the process boundary, the collision pattern between the particle and the triangular element is automatically updated according to the pseudo-particle region collision pattern. It is easy to see that since the boundary conditions of the triangular element do not move, the collision between the particle at the process boundary and the triangular element with the same number in the adjacent process is the same as that in the current process, and no separate processing is required.

[0133] like Figure 11 As shown, identical cones were constructed using both the conical geometric boundary parallelization method and the triangular element boundary parallelization method. The test system with the minimum computational cost under the same Knudsen number and Mach number as the wind tunnel experiment was tested, and the drag coefficient was compared with the Knudsen number. All particle-cone collisions were performed under completely diffuse reflection conditions, i.e., the thermal adaptation coefficient was set to 1, and the Mach number was 7.6. The simulation results were within a reasonable range, and this simulation verified the accuracy of the conical geometric boundary and the triangular element boundary.

[0134] Example 3

[0135] This embodiment provides a method for designing complex boundaries based on basic geometric boundary units such as spheres, cones, circular pipes, simple planes, and triangular facets, as well as combinations of accumulated frozen particles, and for parallelizing event-driven simulations. This method is suitable for simulating gas-solid two-phase flows with complex boundaries in general. Complex boundaries include, but are not limited to: multi-sphere, multi-cone, funnel-shaped boundaries, spherical-conical boundaries, cross-shaped and T-shaped pipes, complex circular pipe boundaries, etc. Figure 12 As shown, the black areas represent frozen particles, while the gray areas represent other analytical geometric shapes. The specific process is as follows:

[0136] The composite boundary is composed of basic geometric boundary units, which are sequentially numbered. A grid is used to record the basic unit numbers involved in the region. Since multiple basic units may exist in the same grid, when a particle moves to the grid in the hard sphere region, it is necessary to predict whether the particle will collide with the basic geometric boundary unit, and to sort the events according to the time of collision. In the pseudo-particle region, after each time step, after the particle moves, it is necessary to check whether it collides with the basic geometric boundary unit in the grid, and to process the collisions in sequence.

[0137] For funnel-shaped boundaries, spherical-conical boundaries, and circular-tube boundaries, the constituent geometries need to be limited in scope. For each basic geometry, this can be achieved by specifying the height and depth boundaries in the x, y, and z directions, thus confining it within a cuboid-shaped enclosure. For example, in the spherical-conical region, the cuboid enclosure restricts the sphere to use its left side while the cone uses its right side. In the hard sphere region, when a collision event occurs with a basic constituent unit, after moving to the collision point, it first checks whether the collision point coordinates are within the cuboid enclosure. If they are, the collision occurs normally; otherwise, the event is invalidated, and the events within the mesh are updated. In the pseudo-granular region, when a collision with a basic constituent unit is detected, it is also necessary to determine whether the coordinates are within the cuboid enclosure corresponding to the geometry. If not, the collision is not processed.

[0138] In special cases, the method of stacking frozen particles can be used to supplement the lack of geometric boundaries. Since the algorithm of moving particle collision freezing particles is relatively simple and not much different from ordinary particle collision, it has little impact on the algorithm process. Therefore, stacking frozen particles can be used in areas with few particles but complex configurations, such as pipe connections.

[0139] Example 4

[0140] This embodiment provides an event-driven parallel design method for ellipsoidal hard sphere programs, where ellipsoidal particles can provide shapes closer to real gas particles. The parallel program connects process boundaries using pseudo-particles and determines a unique tangent point using a scaling method. The specific process is as follows:

[0141] The motion of an ellipsoid includes translational motion and rotation, such as Figure 13 As shown, let the lengths of the three semi-axes of the ellipsoid be a, b, and c, respectively, and let the position of the centroid at time t be r. t (x0, y0, z0), with initial rotation quaternion q t (q0, q1, q2, q3), translational velocity is v(v x ,v y ,v z ), with a rotational angular velocity of w(w x ,w y ,w z The translational position update of the center of mass is similar to that of a spherical particle, and the rotational phase change is represented by quaternions. After t+Δt, the position of ellipsoid i and the rotational quaternion are respectively:

[0142]

[0143]

[0144]

[0145] The fundamental equation of an ellipsoid centered at the origin is:

[0146]

[0147] In three-dimensional Euclidean space, the matrix form of the ellipsoid equation can be expressed in two ways: a three-dimensional characteristic matrix and a four-dimensional characteristic matrix. This paper adopts the four-dimensional matrix form. For the four-dimensional matrix form, the equation of the ellipsoid centered at the origin can be expressed as:

[0148] XSX T =0

[0149] in,

[0150]

[0151] X = [xyz 1]

[0152] The characteristic matrix of the ellipsoidal translational motion is:

[0153]

[0154] The characteristic matrix of the ellipsoidal rotational motion is:

[0155]

[0156] The general form of the ellipsoid equation is:

[0157] XAX T =0

[0158] in,

[0159] A = TRSR T T T

[0160] Consider two ellipsoids A and B, whose matrix forms are XAX respectively. T =0, XBX T =0. Define the characteristic polynomial as f = det(λA + B) / det(A). The phase relationship between the two ellipsoids depends on the sign of the solution to this characteristic polynomial f = 0. It is easy to see that the highest term of the equation with respect to λ in the expansion of this formula is fourth degree, and its expansion form is:

[0161] f(λ)=λ 4 +aλ 3 +bλ 2 +cλ+d

[0162] Since the positional relationships of the ellipsoids (separated, tangent, coincident) depend only on the sign of the solution to this quartic equation, determining the positional relationships does not require solving the specific equation; it only requires simplification using the characteristic sub-formula to determine the sign change of the equation's roots. The characteristic sub-formula for this equation is:

[0163] sr0 = -192cd 2 a+256d 3 +144c 2 db+b 2 a 2 c 2 -6c 2 da 2 +18c 3 ba+144ba 2 d 2 -4b 3 a 2 d+16b 4 d-4c 3 a 3 -128d 2 b 2 -4b 3 c 2 -27a 4 d 2 -80cb 2 ad+18cba 3 d-27c 4

[0164] sr1(λ)=(-6a 3 c+2a 2 b 2 -12a 2 d+28abc-8b3 -36c 2 +32bd)λ-9a 3 d+a 2 bc+3ac 2 +32abd-4b 2 c-48cd

[0165] sr2(λ)=(-8b+3a 2 )λ 2 +(2ab-12c)λ+ac-16d

[0166] Define Var(1, m, n, p, q) as the number of sign changes of the coefficient sequence.

[0167] If the two ellipsoids are separated, then Var(1, a, b, c, d) = 2 and

[0168] (a)sr 22 >0, sr 11 >0, sr0>0; or

[0169] (b)sr 22 >0, sr 11 >0, sr 10 >0, sr0=0.

[0170] If the two ellipsoids are tangent, then:

[0171] (a)sr 22 >0, sr 11 >0, sr 10 <0, sr0=0; or

[0172] (b)sr 22 >0, sr 20 >0, sr 11 =0, sr0=0.

[0173] In other cases, the two ellipsoids overlap.

[0174] In the pseudo-particle part, at any time t, the characteristic matrices of the two ellipsoids are obtained, and then substituted into the characteristic sub-node relationship to determine whether the two ellipsoids overlap. If they overlap, the collision bounce update speed is performed.

[0175] It is worth noting that for quasi-particles, due to overlap, the two positive roots λ0, which should be identical, will differ slightly, leading to errors in the calculation of the tangent point. Therefore, for overlapping quasi-particles, the three semi-axes of the two ellipsoids are reduced by the same multiple using a scaling method until the two particles are tangent. A constructor is then used to determine the positional relationship of the ellipsoids, such as... Figure 14As shown in the pseudocode, the input scaling factor is used to scale down two ellipsoids. If the scaled ellipsoids overlap, it returns 1; if they are tangent, it returns 0; and if they have no intersection, it returns -1. Although the function is piecewise and singular at zero, it is easy to see that within the scaling range, the function monotonically changes from 1 to -1 and has only one zero. The scaling factor that makes the two ellipsoids tangent can be easily obtained by numerically solving the problem using Newton's iteration method. Then, the velocity after the collision can be calculated using the scaled tangent ellipsoid.

[0176] For the hard sphere part, assume two ellipsoids A and B are initially separated. To save computation, an external tangent sphere is introduced. If the external tangent sphere cannot collide, then the ellipsoids definitely cannot collide. If the external tangent sphere can collide, the possible collision time interval [t1, t2] of the ellipsoids is determined by the time interval between the two tangent sphere collisions. Due to the introduction of the time variable t, the characteristic sub-equations sr are no longer constants, but variables with respect to time t. From the formulas in the previous section, the characteristic sub-equations sr0(t) and sr0(t) corresponding to any given time t can be obtained. 11 (t), from which sr0 and sr can be constructed. 11 For a function of time t, using the Brent numerical iteration method, we only need to find the minimum t∈[t1,t2] that satisfies the tangency of the two ellipsoids, and this moment is the collision moment.

[0177] A conical static wall was constructed using a method of stacked frozen particles. By setting the mass of the wall particles to infinity, the velocity update formula after gas particles collide with the wall particles was calculated, simulating a flow around an ellipsoidal particle cone with a half-axis ratio of 0.5:0.25:0.25. The density and temperature field distributions at the central section of the system after stabilization are shown in Figures 15(a) and (b). Both the density and temperature distributions are reasonable. This test demonstrates the potential of using ellipsoids, which are closer to the shape of real gas molecules, instead of spherical particles for flow simulation.

[0178] In summary, through the examples of event-driven parallelization of conical geometry boundaries, event-driven parallelization of triangular surface boundary boundaries, and event-driven parallelization of ellipsoidal collisions in the embodiments, the feasibility of the parallel computing method for gas-solid simulation of the present invention and its broad application prospects are demonstrated.

[0179] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to the embodiments, those skilled in the art should understand that modifications or equivalent substitutions to the technical solutions of the present invention do not depart from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A parallel computing method for gas-solid simulation, the method comprising: Step 1) Based on the multiphase complex boundary conditions, a geometric boundary is constructed from several static geometric units to form a solid. Step 2) Set the initial state of the system to be simulated according to the simulation settings information. The gas particles are modeled as spherical or ellipsoidal particles. Step 3) Decompose the system to be simulated into multiple sub-regions. Each sub-region includes the interior and the boundary. Hard sphere simulation is used for the interior and quasi-particle simulation is used for the boundary. Each sub-region is divided into several grids. Establish the mapping relationship between particles and grids and calculate the geometric boundary in each grid. Step 4) Within the set time step, based on event-driven updates, update particle information for particle collisions within each sub-region and collisions between particles and geometric boundaries; based on time-driven updates, update particle information for particle collisions at the boundaries of each sub-region and collisions between particles and geometric boundaries. Step 5) Determine whether the simulation termination condition has been met. If yes, proceed to step 6); otherwise, proceed to step 4. Step 6) Output the simulation results.

2. The parallel computing method for gas-solid simulation according to claim 1, characterized in that, The multiphase complex boundary includes one or more of the following: conical boundary, triangular surface boundary, spherical boundary, and cylindrical boundary; the geometric shape unit includes one or more of the following: sphere, cone, tube, simple plane, and triangular surface; the characteristic parameters of the geometric shape unit include shape feature description parameters and geometric shape range parameters.

3. The parallel computing method for gas-solid simulation according to claim 2, characterized in that, The gas particles are outside the boundary condition topology and do not overlap with the boundary.

4. The parallel computing method for gas-solid simulation according to claim 3, characterized in that, Before step 4), the following also applies: Based on the position, velocity, and geometric boundary of the particles within the sub-region, a hard sphere region event list is established for that sub-region. The hard sphere region event list includes events that occur in chronological order and their predicted occurrence times. The events include collision events between particles within the sub-region and the geometric boundary, collision events between particles within the sub-region, and events where particles within the sub-region cross the mesh.

5. The parallel computing method for gas-solid simulation according to claim 4, characterized in that, When the event is a collision between a particle and a geometric boundary within a sub-region, the corresponding predicted collision time is calculated by sequentially combining the particle centroid motion equations with the analytical structure equations of each geometric unit, using analytical or numerical iteration methods.

6. The parallel computing method for gas-solid simulation according to claim 5, characterized in that, When the event is a collision event between particles within a sub-region, the gas particles adopt an ellipsoidal particle model. An external tangent sphere of the ellipsoidal particle is introduced. If there is no collision between the external tangent spheres, there is also no collision between the ellipsoidal particles. If there is a collision between the external tangent spheres, the possible collision time interval [t1, t2] of the ellipsoidal particles is determined by the time interval between the two tangent spheres. The minimum time t∈[t1, t2] that satisfies the tangency of the two ellipsoidal particles is obtained by the Brent numerical iteration method. This time t is the collision time, and then the velocity and position of the ellipsoidal particles are updated.

7. The parallel computing method for gas-solid simulation according to claim 6, characterized in that, Step 4) involves updating particle information based on event-driven events within a set time step, specifically addressing particle collisions within each sub-region and collisions between particles and geometric boundaries; this includes: Within a set time step, read the event list of the hard sphere region, obtain the earliest occurring event, and determine whether it is a real event. If not, obtain the next event record. If it is, then based on the type of the event record: if the event is a collision between a particle and the geometric boundary within the sub-region, determine the particle's velocity after the collision based on the analytical geometric expression of the particle's position coordinates and the geometric boundary, and then update the particle's position and velocity; if the event is a collision between particles within the sub-region, calculate the translational and rotational velocities based on the conservation of momentum, angular momentum, and total energy, and then update the particle's position and velocity; if the event is a particle crossing the mesh within the sub-region, update the mapping relationship between the particle and the mesh, and then update the particle's position and velocity. Delete the events that have already occurred from the hard ball region event list, and regenerate the hard ball region event list corresponding to the sub-region based on the updated particle position and velocity, combined with the geometric boundary.

8. The parallel computing method for gas-solid simulation according to claim 7, characterized in that, Step 4) is time-driven and uses a pseudo-particle model to update particle information for particle collisions at the boundaries of each sub-region and for collisions between particles and geometric boundaries; specifically, it includes: For collisions between particles at the boundaries of sub-regions, a pseudo-particle model is used to update the particle velocity and position. When gas particles are modeled as ellipsoidal particles, the characteristic matrices of the two ellipsoids are obtained and then substituted into the characteristic sub-structure formula to determine whether the two ellipsoids overlap. If the two ellipsoids overlap, the tangent point is determined by the scaling method, and then the particle velocity and position after the collision rebound are obtained. For particles and geometric boundaries at the sub-region boundary, the relative positional relationship between them is determined by analytical or numerical iteration methods based on the particle coordinates and the analytical structure equations of each geometric shape element of the corresponding boundary. If they are disjoint, no collision will occur and no action will be taken. If they are tangent or intersecting and tend to approach the interior of the geometric boundary, a collision will occur. The velocity after the collision is calculated, and the position of the particle is updated.

9. The parallel computing method for gas-solid simulation according to claim 8, characterized in that, If the two ellipsoids overlap, the tangent point is determined by the scaling method, thereby obtaining the particle velocity and position after the collision and rebound. Specifically, this includes: The two ellipsoidal particles are reduced in size using the scaling method. The scaling factor that makes the two ellipsoids tangent is obtained by numerical solution using Newton's iteration method. Then, the velocity and position after the collision are calculated using the reduced tangent ellipsoid.

10. The parallel computing method for gas-solid simulation according to claim 2, characterized in that, The geometric unit also includes stacked frozen particles.

Images