High-volume-fraction particle random distribution generation method based on grid points / grids

Abstract

The invention discloses a grid point / grid-based high-volume-fraction particle random distribution generation method, which comprises the following steps of: firstly, setting input parameters, namely a spherical particle radius R, a ratio delta of the side length of a putting area to the spherical particle radius, and an upper limit value nitr of the total iteration step number, so that the length, the width and the height of the putting area are width = R * delta; assuming that the x-axis of the space rectangular coordinate is horizontally rightward, the y-axis is horizontally backward and the z-axis is vertically upward, each side of the throwing area is respectively parallel to the coordinate axis, and the left lower front angular point is located at the origin of the coordinate; the generation method comprises the following steps: dispersing a delivery area, establishing a fully large three-dimensional background grid to cover the delivery area, uniformly dividing the delivery area into m * m * m grids, enabling the grid width w0 to be equal to width / m, and establishing a deliverable state for all grid points / grids in the delivery area; the method has the beneficial effects that as many balls are thrown into the throwing area as possible through a series of calculation so as to obtain a relatively high volume fraction as much as possible, and the performance of the particle reinforced composite material is researched through the geometric model.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation technology for particle-reinforced composite materials, and in particular to a method for generating high volume fraction particles based on the random distribution of grid points / grids. Background Technology

[0002] Particle-reinforced composite materials have broad application prospects in aerospace, construction and other fields. Mechanical problems exist in the design, processing and application of particle-reinforced composite materials. Researchers often use the finite element method to study their material properties, and the representative volume element (RVE) model is usually used to calculate and analyze the equivalent mechanical properties of the material.

[0003] To establish a suitable representative volume element, in addition to considering the volume fraction of particles, particle size, and the physical properties of each component, it is also necessary to consider the distribution of particles in the matrix. Most existing methods assume that the particle distribution is regular to simplify the microstructural analysis of the representative volume element; however, the actual distribution of particles in the matrix is ​​random. Therefore, the established representative volume element should represent this randomness and inhomogeneity as accurately as possible.

[0004] For particle-reinforced composites, spheres are often used to represent particles. To establish a three-dimensional representative volume element representing randomly distributed particles, spheres need to be randomly placed within the cubic region representing the representative volume element, while satisfying geometric periodicity. That is, when a placed sphere intersects the boundary of the Representation Volume Element (RVE), the spheres located outside the boundary are translated to the other side of the RVE. The volume fraction of the three-dimensional representative volume element is equal to the sum of the volumes of all spheres whose centers are located within the cubic placement region, divided by the volume of the cube.

[0005] The ratio of the width of a representative volumetric unit to the radius of the particle sphere is usually denoted as δ. The larger δ is, the closer the representative volumetric unit is to the macroscopic properties of the material, but the computational cost is also higher. Researchers generally use values ​​between 15 and 50 for δ.

[0006] Random Sequential Adsorption (RSA) is commonly used to establish representative volumetric units with randomly distributed particles. However, the highest volume fraction it can generate is approximately 0.385, making it difficult to generate representative volumetric units with higher volume fractions. To generate representative volumetric units with high volume fractions of randomly distributed particles, researchers have developed several methods. However, each method has limitations; for example, some generate low volume fractions, some generate very few particles, and some produce particles with poor randomness in distribution. Our proposed method achieves a maximum volume fraction exceeding 0.45 when δ is between 8 and 50. At δ of 50, it generates 13429 particles with good randomness in distribution.

[0007] The applicant and related parties filed a patent application with application number "2024104562567" on April 16, 2024, entitled "Method for Generating High Volume Fraction Fiber Random Distribution Based on Grid Points / Grids." However, this method can only generate the fiber distribution on the cross-section of the fiber composite material, which is equivalent to a two-dimensional cross-section. It cannot solve the aforementioned problem: generating randomly distributed reinforcing particles in particle-reinforced composite materials. Particle-reinforced composite materials refer to composite materials composed of a continuous matrix material and dispersed granular phases that provide reinforcement. The matrix material is the continuous phase, whose main functions are to bind the reinforcing particles, transfer and disperse loads, and protect the particles from environmental damage. The matrix can be metal, ceramic, or polymer. The reinforcing particles are the dispersed phase, whose main functions are to bear the main loads and significantly improve the stiffness, strength, hardness, wear resistance, and thermal stability of the matrix material. The particles are usually ceramic materials, such as silicon carbide, alumina, and silicon nitride. Summary of the Invention

[0008] The purpose of this invention is to propose a method for generating high volume fraction particles by random distribution based on grid points / grids, which can generate geometric models with high volume fraction, and then study the performance of particle-reinforced composite materials.

[0009] To achieve the above objectives, the present invention adopts the following technical solution:

[0010] The high volume fraction particle random distribution generation method based on grid points / grids first sets the input parameters: the radius of the sphere particle R, the ratio of the side length of the placement area to the radius of the sphere particle δ, the target volume fraction, and the upper limit of the total number of iterations n_itr. Then the length, width and height of the placement area are all width=R*δ.

[0011] Assume that the x-axis of the spatial rectangular coordinate system is horizontal to the right, the y-axis is horizontal to the back, and the z-axis is vertical to the top. The sides of the projection area are parallel to the coordinate axes, and its lower left front corner is located at the origin.

[0012] The generation method includes the following steps:

[0013] S1. The delivery area is discrete. A sufficiently large three-dimensional background grid is established to cover the delivery area. The delivery area is evenly divided into m×m×m grids. Then the grid width w0 = width / m. The delivery state is established for all grid points / grids in the delivery area and all are set to delivery. The upper right corner of each grid is taken as the grid point. The grid points corresponding to all grid points / grids in the delivery area are called grid points / grid points in the delivery area.

[0014] S2. Let itr = 1, and let i_ball = 0;

[0015] S3. When dropping based on grid points, randomly select a grid point from all the grid points in the dropping area that are dropable as the dropping point; when dropping based on grid, randomly select a grid from all the grid points in the dropping area that are dropable, and randomly select a point in this grid as the dropping point. Let i_ball = i_ball + 1, drop the i_ball ball so that the center of this dropped ball is located at this dropping point, record the coordinates of the center of this dropped ball and its radius, and add it to the dropping ball set;

[0016] S4. Prepare to release new balls, and calculate the repulsion sphere formed by the i_ball-th released ball and the newly released ball.

[0017] S5. When deploying based on grid points, traverse each grid point located inside the repulsion sphere and update the deployable status of the corresponding grid point in the deployment area to undeployable; when deploying based on grids, traverse each grid that intersects with the repulsion sphere and update the deployable status of the corresponding grid point in the deployment area to undeployable.

[0018] S6. Calculate the number of grid points / grids in the delivery area that are in a deliverable state, denoted as n_grid_node;

[0019] S7. Determine if n_grid_node is equal to zero. If it is equal to zero, proceed to the next step; otherwise, jump to step S3.

[0020] S8. Let k_ball be the number of balls that have been placed in the current iteration step, let k_ball = i_ball, and calculate the current volume fraction;

[0021] S9. Determine whether the current volume fraction is greater than or equal to the target volume fraction. If it is, proceed to step S32; if it is less than, proceed to the next step.

[0022] S10. Determine if itr is greater than or equal to n_itr. If it is, jump to step S36; if it is less than, proceed to the next step.

[0023] S11. Let itr = itr + 1;

[0024] S12. Set the occupation status of all grid points / grids within the deployment area to unoccupied;

[0025] S13. Let i_ball = 1;

[0026] S14. When dropping based on grid points, traverse each grid point inside the i_ball dropping ball and update the occupation status of the corresponding grid point in the dropping area to occupied; when dropping based on grid, traverse each grid that intersects with the i_ball dropping ball and update the occupation status of the corresponding grid point in the dropping area to occupied.

[0027] S15. Determine if i_ball is greater than or equal to k_ball. If it is, proceed to the next step; if it is less than, let i_ball = i_ball + 1 and jump to step S14.

[0028] S16. Let i_ball = 1;

[0029] S17. Prepare a reference ball and calculate the repulsion sphere formed by the i-th ball and the reference ball.

[0030] S18. When deploying based on grid points, traverse every grid point inside the repulsion sphere, count the number of grid points inside the repulsion sphere (num2), count the number of grid points inside the repulsion sphere whose corresponding grid point occupies the deployment area (num), calculate the local volume fraction vfl = num / num2 for the i-th deployed ball, and save it to the local volume fraction array; when deploying based on grid, traverse every grid that intersects with the repulsion sphere, count the number of grid points that intersect with the repulsion sphere (num2), count the number of grid points that intersect with the repulsion sphere whose corresponding grid point occupies the deployment area (num), calculate the local volume fraction vfl = num / num2 for the i-th deployed ball, and save it to the local volume fraction array.

[0031] S19. Determine if i_ball is greater than or equal to k_ball. If it is, proceed to the next step; if it is less than, let i_ball = i_ball + 1 and jump to step S17.

[0032] S20. Based on the local volume fractions of all k_ball balls, calculate the maximum value vfl_max and the minimum value vfl_min, and calculate the baseline local volume fraction vfl_base = (1-alpha)*vfl_min + alpha*vfl_max;

[0033] S21. Set the retention status of all thrown balls to non-retention.

[0034] S22, Let i_ball = 1;

[0035] S23. Determine whether the local volume fraction of the i-th ball is greater than or equal to vfl_base. If it is, proceed to the next step; if it is less than, jump to step S25.

[0036] S24. Mark the retention status of the i_ball-th ball as retained;

[0037] S25. Determine if i_ball is greater than or equal to k_ball. If it is, proceed to the next step; if it is less than, let i_ball = i_ball + 1 and jump to step S23.

[0038] S26. Delete the balls whose retention status is not retained from the ball set, update the ball index in the new set, and update the value of k_ball to the number of retained balls;

[0039] S27. Set the drop status of all grid points / grids in the drop area to dropable;

[0040] S28. Let i_ball = 1, and prepare to throw a new ball;

[0041] S29. Calculate the repulsion ball formed by the i-th ball and the newly placed ball;

[0042] S30. When deploying based on grid points, traverse each grid point located inside the repulsion sphere and update the deployable state of the corresponding grid point in the deployment area to be undeployable; when deploying based on grids, traverse each grid that intersects with the repulsion sphere and update the deployable state of the corresponding grid point in the deployment area to be undeployable.

[0043] S31. Determine if i_ball is greater than or equal to k_ball. If it is, jump to step S6; if it is less than, let i_ball = i_ball + 1, and jump to step S29.

[0044] S32. Let i_ball = 1;

[0045] S33. Determine whether the i_th ball intersects with the boundary of the throwing area. If they intersect, generate a supplementary ball for this ball and add it to the supplementary ball set; otherwise, jump to step S34.

[0046] S34. Determine if i_ball is greater than or equal to k_ball. If it is, proceed to the next step; if it is less than, let i_ball = i_ball + 1 and jump to step S33.

[0047] S35. Save and output the set of thrown balls and the set of replenished balls, output a successful throwing flag, and jump to step S37.

[0048] S36, Output delivery failure flag;

[0049] S37, End the deployment.

[0050] Preferably, in steps S4 and S29, the center coordinates of the repulsion ball are equal to the center coordinates of the i-th placed ball, and the radius of the repulsion ball is equal to the radius of the i-th placed ball plus the radius of the newly placed ball.

[0051] Preferably, the traversal and update processes in steps S5 and S30 are the same, both including the following steps:

[0052] S51. Calculate the lower and upper bound coordinates of the cube circumscribed by the repulsion sphere in the z-direction, denoted as z0_out and z1_out respectively. Then z0_out = z0 - R_ex and z1_out = z0 + R_ex. When deploying based on grid points, calculate the lower and upper bound grid numbers in the z-direction, denoted as i0_z and i1_z respectively. Then i0_z = floor(z0_out / w0) + 1 and i1_z = ceil(z1_out / w0) - 1. When deploying based on grid points, calculate the lower and upper bound grid numbers in the z-direction, denoted as i0_z and i1_z respectively. Then i0_z = floor(z0_out / w0) + 1, i1_z = ceil(z1_out / w0), and calculate the layer number and row number of the grid where the center of the repulsion sphere is located, denoted as center_z and center_y respectively, then center_z = floor(z0 / w0) + 1, center_y = floor(y0 / w0) + 1; floor() means round down, ceil() means round up, y0 and z0 are the y and z coordinates of the center of the repulsion sphere, and R_ex is the radius of the repulsion sphere;

[0053] S52. Let i_z = i0_z;

[0054] S53. Determine if i_z is less than or equal to zero. If it is, let i_z_regin = i_z + m, and jump to step S55. If it is not less than or equal to zero, proceed to the next step.

[0055] S54. Determine whether \(i_z\) is greater than \(m\). If it is greater, let \(i_z\_regin = i_z - m\); if it is not greater, let \(i_z\_regin = i_z\).

[0056] S55. When placing based on grid points, calculate the square of the vertical distance from the center of the sphere to the current grid point layer, denoted as \(dz2\). Then \(dz2=(z0 - i_z*w0)\) 2 , calculate the radius of the reference circle formed by the intersection of the current grid point layer and the exclusion sphere, denoted as \(r_i_z\). Then Calculate the lower and upper grid point numbers in the \(y\) direction of the grid points located inside the exclusion sphere in the current grid point layer, denoted as \(i0_y\) and \(i1_y\) respectively. Then \(i0_y = floor((y0 - r_i_z) / w0)+1\), \(i1_y = ceil((y0 + r_(i_z)) / w0)-1\); When placing based on the grid, determine the reference plane of the current grid layer and calculate the \(z\) coordinate \(z_ref\) of the reference plane. If \(i_z < center_z\), the upper surface of the current grid layer is the reference plane of the current grid layer, \(z_ref = i_z*w0\). If \(i_z = center_z\), the horizontal plane where the center of the sphere is located is the reference plane of the current grid layer, \(z_ref = 0\). If \(i_z > center_z\), the lower surface of the current grid layer is the reference plane of the current grid layer, \(z_ref=(i_z - 1)*w0\). Calculate the square of the minimum vertical distance from the center of the sphere to the reference plane of the current grid layer, denoted as \(dz2\), \(dz2=(z0 - z_ref)\) 2 , calculate the radius of the reference circle formed by the intersection of the reference plane of the current grid layer and the exclusion sphere, denoted as \(r_i_z\). Then Calculate the lower and upper grid numbers in the \(y\) direction of the grids intersecting with the exclusion sphere in the current grid layer, denoted as \(i0_y\) and \(i1_y\) respectively. Then \(i0_y = floor((y0 - r_i_z) / w0)+1\), \(i1_y = ceil((y0 + r_i_z) / w0)\); \(floor()\) represents rounding down, \(ceil()\) represents rounding up, \(y0\) and \(z0\) are the \(y\) and \(z\) coordinates of the center of the exclusion sphere, and \(R_ex\) is the radius of the exclusion sphere.

[0057] S56. Let \(i_y = i0_y\).

[0058] S57. Determine whether \(i_y\) is less than or equal to zero. If it is less than or equal to, let \(i_y\_regin = i_y + m\), and jump to step S59; if it is not less than or equal to, proceed to the next step.

[0059] S58. Determine whether \(i_y\) is greater than \(m\). If it is greater, let \(i_y\_regin = i_y - m\); if it is not greater, let \(i_y\_regin = i_y\).

[0060] S59. When placing based on grid points, calculate the square of the vertical distance from the center of the reference circle to the current grid point row, denoted as dy2, then dy2 = (y0 - i_y * w0). 2 , calculate half of the length of the chord formed by the intersection of the line where the current grid point row is located and the reference circle, denoted as dx, then Calculate the lower and upper grid point numbers in the x direction of the grid points located inside the exclusion sphere in the current grid point row of the current grid point layer, denoted as i0_x and i1_x respectively, then i0_x = floor((x0 - dx) / w0)+1, i1_x = ceil((x0 + dx) / w0)-1; When placing based on the grid, determine the reference plane of the current grid row and calculate the y coordinate y_ref of the reference plane. If i_y < center_y, the back surface of the current grid row is the reference plane of the current grid row, y_ref = i_y * w0. If i_y = center_y, the vertical plane where the center of the reference circle is located and parallel to the front and back surfaces of the grid is the reference plane of the current grid row, y_ref = 0. If i_y > center_y, the front surface of the current grid row is the reference plane of the current grid row, y_ref = (i_y - 1) * w0. Calculate the square of the vertical distance from the center of the reference circle to the reference plane of the current grid row, denoted as dy2, dy2 = (y0 - y_ref). 2 , calculate half of the length of the chord formed by the intersection of the intersection line of the reference plane of the current grid layer and the reference plane of the current grid row and the reference circle, denoted as dx, then Calculate the lower and upper grid numbers in the x direction of the grids that intersect with the exclusion sphere in the current grid row of the current grid layer, denoted as i0_y and i1_y respectively, then i0_x = floor((x0 - dx) / w0)+1, i1_x = ceil((x0 + dx) / w0); floor() represents rounding down, ceil() represents rounding up, y0 and x0 are the y and x coordinates of the center of the exclusion sphere;

[0061] S510. Let i_x = i0_x;

[0062] S511. Judge whether i_x is less than or equal to zero. If it is less than or equal to, let i_x_regin = i_x + m, and jump to step S513; If it is not less than or equal to, go to the next step;

[0063] S512. Judge whether i_x is greater than m. If it is greater than, let i_x_regin = i_x - m; If it is not greater than, let i_x_regin = i_x;

[0064] S513. Set the placeable state of the grid point / grid at the i_x_regin-th column, i_y_regin-th row, and i_z_regin-th layer to non-placeable;

[0065] S514. Determine whether i_x is greater than or equal to i1_x. If it is, proceed to the next step; if it is less than, let i_x = i_x + 1 and jump to step S511.

[0066] S515. Determine whether i_y is greater than or equal to i1_y. If it is, proceed to the next step; if it is less than, let i_y = i_y + 1 and jump to step S57.

[0067] S516. Determine whether i_z is greater than or equal to i1_z. If it is, proceed to the next step; if it is less than, let i_z = i_z + 1 and jump to step S53.

[0068] S517, End traversal.

[0069] Preferably, in step S8, the current volume fraction is equal to the sum of the volumes of all the balls released divided by the volume of the release area.

[0070] Preferably, in step S14, the traversal update process includes the following steps:

[0071] S141. Calculate the lower and upper bound coordinates of the circumscribed cube of the deployed sphere in the z-direction, denoted as z0_out and z1_out respectively. Then z0_out = z0 - R and z1_out = z0 + R. When deploying based on grid points, calculate the lower and upper bound grid numbers in the z-direction, denoted as i0_z and i1_z respectively. Then i0_z = floor(z0_out / w0) + 1 and i1_z = ceil(z1_out / w0) - 1. When deploying based on grid points, calculate the lower and upper bound grid numbers in the z-direction, denoted as i0_z and i1_z respectively. Then i0_z = floor(z0_out / w0) + 1, i1_z = ceil(z1_out / w0), and calculate the layer number and row number of the grid where the center of the ball is located, denoted as center_z and center_y respectively, then center_z = floor(z0 / w0) + 1, center_y = floor(y0 / w0) + 1; floor() means round down, ceil() means round up, y0 and z0 are the y and z coordinates of the center of the ball, and R is the radius of the ball.

[0072] S142. Let i_z = i0_z;

[0073] S143. Determine if i_z is less than or equal to zero. If it is, let i_z_regin = i_z + m, and jump to step S145. If it is not less than or equal to zero, proceed to the next step.

[0074] S144. Determine whether \(i_z\) is greater than \(m\). If it is greater, let \(i_z\_regin = i_z - m\); if it is not greater, let \(i_z\_regin = i_z\).

[0075] S145. When placing based on grid points, calculate the square of the vertical distance from the center of the ball to the current grid point layer, denoted as \(dz2\), then \(dz2=(z0 - i_z*w0)\) 2 , calculate the radius of the reference circle formed by the intersection of the current grid point layer and the placed ball, denoted as \(r_i_z\), then Calculate the lower and upper grid point numbers in the \(y\) direction of the grid points located inside the placed ball in the current grid point layer, denoted as \(i0_y\) and \(i1_y\) respectively. Then \(i0_y = floor((y0 - r_i_z) / w0)+1\), \(i1_y = ceil((y0 + r_i_z) / w0)-1\); When placing based on the grid, determine the reference plane of the current grid layer and calculate the \(z\) coordinate \(z\_ref\) of the reference plane. If \(i_z < center\_z\), the upper surface of the current grid layer is the reference plane of the current grid layer, \(z\_ref = i_z*w0\). If \(i_z = center\_z\), the horizontal plane where the center of the ball is located is the reference plane of the current grid layer, \(z\_ref = 0\). If \(i_z > center\_z\), the lower surface of the current grid layer is the reference plane of the current grid layer, \(z\_ref=(i_z - 1)*w0\). Calculate the square of the minimum vertical distance from the center of the ball to the reference plane of the current grid layer, denoted as \(dz2\), \(dz2=(z0 - z\_ref)\) 2 , calculate the radius of the reference circle formed by the intersection of the reference plane of the current grid layer and the placed ball, denoted as \(r_i_z\), then Calculate the lower and upper grid numbers in the \(y\) direction of the grids that intersect with the placed ball in the current grid layer, denoted as \(i0_y\) and \(i1_y\) respectively. Then \(i0_y = floor((y0 - r_i_z) / w0)+1\), \(i1_y = ceil((y0 + r_i_z) / w0)\); \(floor()\) represents rounding down, \(ceil()\) represents rounding up, \(y0\) and \(z0\) are the \(y\) and \(z\) coordinates of the center of the placed ball, and \(R\) is the radius of the placed ball.

[0076] S146. Let \(i_y = i0_y\).

[0077] S147. Determine whether \(i_y\) is less than or equal to zero. If it is less than or equal to, let \(i_y\_regin = i_y + m\), and jump to step S149; if it is not less than or equal to, proceed to the next step.

[0078] S148. Determine whether \(i_y\) is greater than \(m\). If it is greater, let \(i_y\_regin = i_y - m\); if it is not greater, let \(i_y\_regin = i_y\).

[0079] S149. When placing based on grid points, calculate the square of the vertical distance from the center of the reference circle to the current grid point row, denoted as dy2, then dy2 = (y0 - i_y * w0). 2 , calculate half of the length of the chord formed by the intersection of the straight line where the current grid point row is located and the reference circle, denoted as dx, then Calculate the lower and upper grid point numbers in the x direction of the grid points located inside the placement ball in the current grid point row of the current grid point layer, denoted as i0_x and i1_x respectively, then i0_x = floor((x0 - dx) / w0) + 1, i1_x = ceil((x0 + dx) / w0) - 1; When placing based on the grid, determine the reference plane of the current grid row and calculate the y coordinate y_ref of the reference plane. If i_y < center_y, the back surface of the current grid row is the reference plane of the current grid row, y_ref = i_y * w0. If i_y = center_y, the vertical plane where the center of the reference circle is located and parallel to the front and back surfaces of the grid is the reference plane of the current grid row, y_ref = 0. If i_y > center_y, the front surface of the current grid row is the reference plane of the current grid row, y_ref = (i_y - 1) * w0. Calculate the square of the vertical distance from the center of the reference circle to the reference plane of the current grid row, denoted as dy2, dy2 = (y0 - y_ref). 2 , calculate half of the length of the chord formed by the intersection of the intersection line of the reference plane of the current grid layer and the reference plane of the current grid row and the reference circle, denoted as dx, then Calculate the lower and upper grid numbers in the x direction of the grids intersecting with the placement ball in the current grid row of the current grid layer, denoted as i0_y and i1_y respectively, then i0_x = floor((x0 - dx) / w0) + 1, i1_x = ceil((x0 + dx) / w0); floor() represents rounding down, ceil() represents rounding up, y0 and x0 are the y and x coordinates of the center of the placement ball.

[0080] S1410. Let i_x = i0_x;

[0081] S1411. Judge whether i_x is less than or equal to zero. If it is less than or equal to, let i_x_regin = i_x + m, and jump to step S1413; If it is not less than or equal to, proceed to the next step;

[0082] S1412. Judge whether i_x is greater than m. If it is greater than, let i_x_regin = i_x - m; If it is not greater than, let i_x_regin = i_x;

[0083] S1413. Set the occupancy status of the grid point / grid at the i_x_regin-th column, i_y_regin-th row, and i_z_regin-th layer to occupied;

[0084] S1414. Determine if i_x is greater than or equal to i1_x. If it is, proceed to the next step; if it is less than, let i_x = i_x + 1 and jump to step S1411.

[0085] S1415. Determine whether i_y is greater than or equal to i1_y. If it is, proceed to the next step; if it is less than, let i_y = i_y + 1 and jump to step S147.

[0086] S1416. Determine if i_z is greater than or equal to i1_z. If it is, proceed to the next step; if it is less than, let i_z = i_z + 1 and jump to step S143.

[0087] S1417, End traversal.

[0088] Preferably, in step S17, the center coordinates of the repulsion ball are equal to the center coordinates of the i-th ball, and the radius of the repulsion ball is equal to the radius of the i-th ball plus the radius of the reference ball.

[0089] Preferably, in step S18, the traversal process includes the following steps:

[0090] S181. Calculate the lower and upper bound coordinates of the cube circumscribed by the repulsion sphere in the z-direction, denoted as z0_out and z1_out respectively. Then z0_out = z0 - R_ex and z1_out = z0 + R_ex. When placing the cube based on the grid points, calculate the lower and upper bound grid numbers in the z-direction, denoted as i0_z and i1_z respectively. Then i0_z = floor(z0_out / w0) + 1 and i1_z = ceil(z1_out / w0) - 1. When placing the cube based on the grid, calculate the lower and upper bound grid numbers in the z-direction, denoted as i0_z and i1_ If z is the center of the repulsion sphere, then i0_z = floor(z0_out / w0) + 1, i1_z = ceil(z1_out / w0), and calculate the layer number and row number of the grid where the center of the repulsion sphere is located, denoted as center_z and center_y respectively, then center_z = floor(z0 / w0) + 1, center_y = floor(y0 / w0) + 1; floor() means round down, ceil() means round up, y0 and z0 are the y and z coordinates of the center of the repulsion sphere, and R_ex is the radius of the repulsion sphere;

[0091] S182. Let num = 0, let num2 = 0, let i_z = i0_z;

[0092] S183. Determine whether \(i_z\) is less than or equal to zero. If it is, let \(i_z\_regin = i_z + m\) and jump to step S185; if not, proceed to the next step.

[0093] S184. Determine whether \(i_z\) is greater than \(m\). If it is, let \(i_z\_regin = i_z - m\); if not, let \(i_z\_regin = i_z\).

[0094] S185. When placing based on grid points, calculate the square of the vertical distance from the center of the ball to the current grid point layer, denoted as \(dz2\). Then \(dz2=(z0 - i_z*w0)\) 2 , calculate the radius of the reference circle formed by the intersection of the current grid point layer and the repulsion ball, denoted as \(r_i_z\). Then Calculate the lower and upper grid point numbers in the \(y\) direction of the grid points located inside the repulsion ball in the current grid point layer, denoted as \(i0_y\) and \(i1_y\) respectively. Then \(i0_y = floor((y0 - r_i_z) / w0)+1\), \(i1_y = ceil((y0 + r_i_z) / w0)-1\); When placing based on the grid, determine the reference plane of the current grid layer and calculate the \(z\) coordinate \(z_ref\) of the reference plane. If \(i_z < center_z\), the upper surface of the current grid layer is the reference plane of the current grid layer, \(z_ref = i_z*w0\); if \(i_z = center_z\), the horizontal plane where the center of the ball is located is the reference plane of the current grid layer, \(z_ref = 0\); if \(i_z > center_z\), the lower surface of the current grid layer is the reference plane of the current grid layer, \(z_ref=(i_z - 1)*w0\). Calculate the square of the minimum vertical distance from the center of the ball to the reference plane of the current grid layer, denoted as \(dz2\), \(dz2=(z0 - z_ref)\) 2 , calculate the radius of the reference circle formed by the intersection of the reference plane of the current grid layer and the repulsion ball, denoted as \(r_i_z\). Then Calculate the lower and upper grid numbers in the \(y\) direction of the grids intersecting with the repulsion ball in the current grid layer, denoted as \(i0_y\) and \(i1_y\) respectively. Then \(i0_y = floor((y0 - r_i_z) / w0)+1\), \(i1_y = ceil((y0 + r_i_z) / w0)\); \(floor()\) represents rounding down, \(ceil()\) represents rounding up, \(y0\) and \(z0\) are the \(y\) and \(z\) coordinates of the center of the repulsion ball, and \(R_ex\) is the radius of the repulsion ball.

[0095] S186. Let \(i_y = i0_y\).

[0096] S187. Determine whether \(i_y\) is less than or equal to zero. If it is, let \(i_y\_regin = i_y + m\) and jump to step S189; if not, proceed to the next step.

[0097] S188. Determine whether \(i_y\) is greater than \(m\). If it is greater, let \(i_y\_regin = i_y - m\); if it is not greater, let \(i_y\_regin = i_y\).

[0098] S189. When placing based on grid points, calculate the square of the vertical distance from the center of the reference circle to the current grid point row, denoted as \(dy2\), then \(dy2=(y0 - i_y*w0)\). 2 , calculate half of the length of the chord formed by the intersection of the straight line where the current grid point row is located and the reference circle, denoted as \(dx\), then Calculate the lower and upper grid point numbers in the \(x\) direction of the grid points located inside the exclusion sphere in the current grid point row of the current grid point layer, denoted as \(i0_x\) and \(i1_x\) respectively. Then \(i0_x = floor((x0 - dx) / w0)+1\), \(i1_x = ceil((x0 + dx) / w0)-1\); when placing based on the grid, determine the reference plane of the current grid row and calculate the \(y\) coordinate \(y\_ref\) of the reference plane. If \(i_y < center_y\), the back surface of the current grid row is the reference plane of the current grid row, \(y\_ref = i_y*w0\). If \(i_y = center_y\), the vertical plane where the center of the reference circle is located and parallel to the front and back surfaces of the grid is the reference plane of the current grid row, \(y\_ref = 0\). If \(i_y > center_y\), the front surface of the current grid row is the reference plane of the current grid row, \(y\_ref=(i_y - 1)*w0\). Calculate the square of the vertical distance from the center of the reference circle to the reference plane of the current grid row, denoted as \(dy2\), \(dy2=(y0 - y\_ref)\). 2 , calculate half of the length of the chord formed by the intersection of the intersection line of the reference plane of the current grid layer and the reference plane of the current grid row and the reference circle, denoted as \(dx\), then Calculate the lower and upper grid point numbers in the \(x\) direction of the grids that intersect with the exclusion sphere in the current grid row of the current grid layer, denoted as \(i0_y\) and \(i1_y\) respectively. Then \(i0_x = floor((x0 - dx) / w0)+1\), \(i1_x = ceil((x0 + dx) / w0)\); \(floor()\) represents rounding down, \(ceil()\) represents rounding up, and \(y0\), \(x0\) are the \(y\) and \(x\) coordinates of the center of the exclusion sphere.

[0099] S1810. Let \(i_x = i0_x\).

[0100] S1811. Determine whether \(i_x\) is less than or equal to zero. If it is less than or equal to zero, let \(i_x\_regin = i_x + m\), and jump to step S1813; if it is not less than or equal to zero, proceed to the next step.

[0101] S1812. Determine if i_x is greater than m. If it is, let i_x_regin = i_x - m; if it is not, let i_x_regin = i_x.

[0102] S1813. Let num2 = num2 + 1. If the occupation status of the cell point / cell in the i_z_regin row and i_x_regin column of the i_y_regin layer is occupied, then let num = num + 1; otherwise, num remains unchanged.

[0103] S1814. Determine if i_x is greater than or equal to i1_x. If it is, proceed to the next step; if it is less than, let i_x = i_x + 1 and jump to step S1811.

[0104] S1815. Determine if i_y is greater than or equal to i1_y. If it is, proceed to the next step; if it is less than, let i_y = i_y + 1 and jump to step S187.

[0105] S1816. Determine if i_z is greater than or equal to i1_z. If it is, proceed to the next step; if it is less than, let i_z = i_z + 1 and jump to step S183.

[0106] S1817. Calculate the local volume fraction vfl = num / num2 and save it to the local volume fraction array, then end the traversal.

[0107] Preferably, in step S33, the generation of the supplementary ball includes the following steps:

[0108] S331. Let the x, y, and z coordinates of the center of the i-th ball be x0, y0, and z0, respectively. Let xyz_min be equal to the ball radius and xyz_max be equal to the region side length minus the ball radius. Execute S332, S334, and S336 in parallel.

[0109] S332. If z0 is less than xyz_min, let z1 equal z0 plus the region side length, let z3 equal zero, and proceed to the next step; if z0 is greater than xyz_max, let z1 equal z0 minus the region side length, let z3 equal the region side length, and proceed to the next step; if z0 is greater than or equal to xyz_min and less than or equal to xyz_max, let z3 equal -1, and jump to step S338.

[0110] S333. A new supplementary ball is generated with center coordinates (x0, y0, z1) and radius equal to the radius of the thrown ball. The center coordinates and radius of this supplementary ball are recorded and added to the supplementary ball set. Then, proceed to step S338.

[0111] S334. If y0 is less than xyz_min, let y1 equal y0 plus the region side length, let y3 equal zero, and proceed to the next step; if y0 is greater than xyz_max, let y1 equal y0 minus the region side length, let y3 equal the region side length, and proceed to the next step; if y0 is greater than or equal to xyz_min and less than or equal to xyz_max, let y3 equal -1, and jump to step S338.

[0112] S335. A new supplementary ball is generated with center coordinates (x0, y1, z0) and radius equal to the radius of the thrown ball. The center coordinates and radius of this supplementary ball are recorded and added to the supplementary ball set. Then, proceed to step S338.

[0113] S336. If x0 is less than xyz_min, set x1 to x0 plus the region side length, set x3 to zero, and proceed to the next step; if x0 is greater than xyz_max, set x1 to x0 minus the region side length, set x3 to the region side length, and proceed to the next step; if x0 is greater than or equal to xyz_min and less than or equal to xyz_max, set x3 to -1, and jump to step S338.

[0114] S337. A new supplementary ball is generated with center coordinates (x1, y0, z0) and radius equal to the radius of the thrown ball. The center coordinates and radius of this supplementary ball are recorded and added to the supplementary ball set. Then, proceed to step S338.

[0115] S338, Wait for S332-S337 to complete execution;

[0116] S339, execute S3310, S3314, S3318 and S3322 in parallel;

[0117] S3310. If x3 is greater than or equal to zero and y3 is greater than or equal to zero, proceed to the next step; otherwise, jump to step S3326.

[0118] S3311. Using x3 and y3 as reference points, calculate the distance d30 from the reference point to (x0, y0).

[0119] S3312. If d30 is less than the radius of the ball being thrown, proceed to the next step; if d30 is greater than or equal to the radius of the ball being thrown, jump to step S3326.

[0120] S3313. A new supplementary ball is generated with center coordinates (x1, y1, z0) and radius equal to the radius of the thrown ball. Record the center coordinates and radius of this supplementary ball and add it to the supplementary ball set. Jump to step S3326.

[0121] S3314. If y3 is greater than or equal to zero and z3 is greater than or equal to zero, proceed to the next step; otherwise, jump to step S3326.

[0122] S3315. Using the x and y coordinates of y3 and z3 as reference points, calculate the distance d30 from the reference point to (y0, z0).

[0123] S3316. If d30 is less than the radius of the ball being thrown, proceed to the next step; if d30 is greater than or equal to the radius of the ball being thrown, jump to step S3326.

[0124] S3317. A new supplementary ball is generated with center coordinates (x0, y1, z1) and radius equal to the radius of the thrown ball. Record the center coordinates and radius of this supplementary ball and add it to the supplementary ball set. Jump to step S3326.

[0125] S3318. If z3 is greater than or equal to zero and x3 is greater than or equal to zero, proceed to the next step; otherwise, jump to step S3326.

[0126] S3319. Using x3 and z3 as reference points, calculate the distance d30 from the reference point to (x0, z0).

[0127] S3320. If d30 is less than the radius of the ball being thrown, proceed to the next step; if d30 is greater than or equal to the radius of the ball being thrown, jump to step S3326.

[0128] S3321. A new supplementary ball is generated with center coordinates (x1, y0, z1) and radius equal to the radius of the thrown ball. Record the center coordinates and radius of this supplementary ball and add it to the supplementary ball set. Jump to step S3326.

[0129] S3322. If x3 is greater than or equal to zero, y3 is greater than or equal to zero, and z3 is greater than or equal to zero, proceed to the next step; otherwise, jump to step S3326.

[0130] S3323. Using x3, y3, z3 as reference points, calculate the distance d30 from the reference points to (x0, y0, z0).

[0131] S3324. If d30 is less than the radius of the ball being thrown, proceed to the next step; if d30 is greater than or equal to the radius of the ball being thrown, jump to step S3326.

[0132] S3325. A new supplementary ball is generated with center coordinates (x1, y1, z1) and radius equal to the radius of the thrown ball. Record the center coordinates and radius of this supplementary ball and add it to the supplementary ball set. Jump to step S3326.

[0133] S3326, Wait for S3310-S3325 to complete.

[0134] The beneficial effects of this invention are as follows:

[0135] 1. Through a series of calculations, as many balls as possible are placed in the placement area to obtain a higher volume fraction, and the properties of particle-reinforced composite materials are studied using this geometric model. Attached Figure Description

[0136] Figure 1 This is an overall flowchart of Embodiment 1 of the present invention;

[0137] Figure 2 This is an overall flowchart of Embodiment 2 of the present invention;

[0138] Figure 3 This is a flowchart of step S33 in Embodiments 1 and 2 of the present invention;

[0139] Figure 4 This is a schematic diagram of the placement area in the experimental example;

[0140] Figure 5 This is a schematic diagram of step S1 in the experimental example;

[0141] Figure 6 This is a schematic diagram of the first run step S3 in the experimental example;

[0142] Figure 7 This is a schematic diagram of the first run of step S4 in the experimental example;

[0143] Figure 8 This is a schematic diagram of step S513 during the first run of S5 in the experimental example;

[0144] Figure 9 This is a diagram showing the positional relationship between the first layer of the grid and the repulsion spheres during the first run of S5 in the experimental example.

[0145] Figure 10 This is a diagram showing the positional relationship between the grid cells in the second layer and the repulsion spheres during the first run of S5 in the experimental example.

[0146] Figure 11 This is a diagram showing the positional relationship between the grid cells and the repulsion spheres in the 6th layer during the first run of S5 in the experimental example.

[0147] Figure 12 This is a schematic diagram of all the grids that intersect with the repulsive ball that was first launched at the end of the first run S5 in the experimental example;

[0148] Figure 13 This is a schematic diagram of the second run step S3 in the experimental example;

[0149] Figure 14 This is a schematic diagram of the second run step S4 in the experimental example;

[0150] Figure 15This is a schematic diagram of step S513 during the second run of S5 in the experimental example;

[0151] Figure 16 This is a diagram showing the positional relationship between the grid in the i_y row and 3rd column of the second layer and the repulsion ball during the second run of S5 in the experimental example.

[0152] Figure 17 This is a schematic diagram showing the cells in the second layer of the distribution area that are not dropable during the second run of S5 in the experimental example.

[0153] Figure 18 This is a schematic diagram of the grids in the delivery area that were previously in a non-deliverable state when the second run S5 in the test case ended.

[0154] Figure 19 This is a schematic diagram of the third running step S3 in the experimental example;

[0155] Figure 20 This is a schematic diagram of the fourth run step S3 in the experimental example;

[0156] Figure 21 This is a schematic diagram of step S3, the fifth run in the experimental example;

[0157] Figure 22 This is a schematic diagram of step S1413 during the first run of step S14 in the experimental example.

[0158] Figure 23 This is a diagram showing the relationship between the position of the ball and the grid in the 4th row of the 3rd layer during the first run of step S14 in the experimental example.

[0159] Figure 24 This is a diagram showing the relationship between the positions of the third layer of cells and the ball during the first run of step S14 in the experimental example.

[0160] Figure 25 This is a diagram showing the relationship between the positions of the traversed cells in the 4th layer and the ball during the first run of step S14 in the experimental example.

[0161] Figure 26 This is a schematic diagram of all the grids that intersect with the first ball during the first run of step S14 in the experimental example;

[0162] Figure 27 This is a schematic diagram of the second running step S14 in the experimental example;

[0163] Figure 28 This is a schematic diagram of the fifth run step S14 in the experimental example;

[0164] Figure 29 This is a schematic diagram of the first run step S1813 in the experimental example;

[0165] Figure 30 This is a schematic diagram of the sixth run step S1813 in the experimental example;

[0166] Figure 31 This is a schematic diagram of step S1813, the 253rd run in the test example;

[0167] Figure 32 This is a schematic diagram of step S26 in the experimental example;

[0168] Figure 33 This is a schematic diagram of the first run step S30 in the experimental example;

[0169] Figure 34 This is a schematic diagram of the fourth run step S30 in the experimental example;

[0170] Figure 35 This is a schematic diagram of the sixth run step S3 in the experimental example;

[0171] Figure 36 This is a schematic diagram of the seventh run step S3 in the experimental example;

[0172] Figure 37 This is a schematic diagram of step S335 during the second execution of step S33 in the experimental example;

[0173] Figure 38 This is a schematic diagram of step S333 during the third run of the test example;

[0174] Figure 39 This is a schematic diagram of step S335 during the third run of step S33 in the experimental example;

[0175] Figure 40 This is a schematic diagram of step S3317 during the third run of step S33 in the experimental example;

[0176] Figure 41 This is a schematic diagram of step S3326 during the third run of step S33 in the experimental example;

[0177] Figure 42 This is a schematic diagram of the fourth run step S33 in the experimental example;

[0178] Figure 43 This is a schematic diagram of the fifth run step S33 in the experimental example;

[0179] Figure 44 This is a schematic diagram of the sixth run step S33 in the experimental example;

[0180] Figure 45 This is a schematic diagram of all the supplementary balls generated at the end of step S34 in the experimental example;

[0181] Figure 46 This is a schematic diagram of step S35 in the experimental example.

[0182] The accompanying drawings are for illustrative purposes only and should not be construed as limiting the scope of this patent. To better illustrate this embodiment, some components in the drawings may be omitted, enlarged, or reduced, and do not represent the actual dimensions of the product. It is understandable to those skilled in the art that some well-known structures and their descriptions may be omitted in the drawings. Detailed Implementation

[0183] The invention will now be further described with reference to the accompanying drawings.

[0184] Example 1

[0185] like Figure 1 and Figure 3 As shown, the high volume fraction particle random distribution generation method based on grid points in this embodiment of the invention first sets the input parameters: the radius of the sphere particle R, the ratio of the side length of the placement area to the radius of the sphere particle δ, the target volume fraction, and the upper limit of the total number of iterations n_itr. Then the length, width and height of the placement area are all width=R*δ.

[0186] Assume that the x-axis of the spatial rectangular coordinate system is horizontal to the right, the y-axis is horizontal to the back, and the z-axis is vertical to the top. The sides of the projection area are parallel to the coordinate axes, and its lower left front corner is located at the origin.

[0187] The generation method includes the following steps:

[0188] S1. The delivery area is discrete. A sufficiently large three-dimensional background grid is established to cover the delivery area. The delivery area is evenly divided into m×m×m grids. The grid width w0 = width / m. The delivery state is established for all grid points in the delivery area and all are set to delivery. The upper right corner of each grid is taken as the grid point. The grid points corresponding to all grid points in the delivery area are called grid points in the delivery area.

[0189] S2. Let itr = 1, and let i_ball = 0.

[0190] S3. Randomly select a grid point from all the grid points in the drop area that are dropable as the drop point; let i_ball = i_ball + 1, drop the i_ball ball so that the center of the dropped ball is located at the drop point, record the coordinates of the center of the dropped ball and its radius, and add it to the drop ball set.

[0191] S4. Prepare to release new balls and calculate the repulsion sphere formed by the i-th released ball and the new ball; where the center coordinates of the repulsion sphere are equal to the center coordinates of the i-th released ball, and the radius of the repulsion sphere is equal to the radius of the i-th released ball plus the radius of the new ball.

[0192] S5. Traverse each grid point located inside the repulsion sphere and update the throwable status of the corresponding grid point within the throwing area to unthrowable; the traversal and update process includes the following steps:

[0193] S51. Calculate the lower and upper bound coordinates of the cube circumscribed by the repulsion sphere in the z-direction, denoted as z0_out and z1_out respectively. Then z0_out = z0 - R_ex and z1_out = z0 + R_ex. Calculate the grid point numbers of the lower and upper bounds in the z-direction, denoted as i0_z and i1_z respectively. Then i0_z = floor(z0_out / w0) + 1 and i1_z = ceil(z1_out / w0) - 1. Calculate the layer number of the grid where the center of the repulsion sphere is located and the row number of the grid, denoted as center_z and center_y respectively. Then center_z = floor(z0 / w0) + 1 and center_y = floor(y0 / w0) + 1. floor() represents rounding down and ceil() represents rounding up. y0 and z0 are the y and z coordinates of the center of the repulsion sphere, and R_ex is the radius of the repulsion sphere.

[0194] S52. Let i_z = i0_z;

[0195] S53. Determine if i_z is less than or equal to zero. If it is, let i_z_regin = i_z + m, and jump to step S55. If it is not less than or equal to zero, proceed to the next step.

[0196] S54. Determine if i_z is greater than m. If it is, let i_z_regin = i_z - m; if it is not, let i_z_regin = i_z.

[0197] S55. Calculate the square of the vertical distance from the center of the sphere to the current grid point layer, denoted as dz2, then dz2=(z0-i_z*w0). 2 Calculate the radius of the reference circle formed by the intersection of the current lattice point layer and the repulsion sphere, denoted as r_i_z. Calculate the lower and upper bounds of the grid points located inside the repulsion sphere in the current grid point layer in the y-direction, denoted as i0_y and i1_y respectively. Then, i0_y = floor((y0-r_i_z) / w0)+1, i1_y = ceil((y0+r_(i_z)) / w0)-1; floor() represents rounding down, ceil() represents rounding up, y0 and z0 are the y and z coordinates of the center of the repulsion sphere, and R_ex is the radius of the repulsion sphere;

[0198] S56. Let i_y = i0_y;

[0199] S57. Determine if i_y is less than or equal to zero. If it is, let i_y_regin = i_y + m, and jump to step S59. If it is not less than or equal to zero, proceed to the next step.

[0200] S58. Determine if i_y is greater than m. If it is, let i_y_regin = i_y - m; if it is not, let i_y_regin = i_y.

[0201] S59. Calculate the square of the vertical distance from the center of the reference circle to the current grid point row, denoted as dy2, then dy2=(y0-i_y*w0). 2 Calculate half the length of the chord formed by the intersection of the line containing the current grid point and the reference circle, denoted as dx. Calculate the lower and upper bounds of the grid points located inside the repulsion sphere in the current grid row of the current grid point layer in the x-direction, denoted as i0_x and i1_x respectively. Then i0_x = floor((x0-dx) / w0)+1, i1_x = ceil((x0+dx) / w0)-1; floor() means rounding down, ceil() means rounding up, and y0 and x0 are the y and x coordinates of the center of the repulsion sphere;

[0202] S510. Let i_x = i0_x;

[0203] S511. Determine if i_x is less than or equal to zero. If it is, set i_x_regin = i_x + m and jump to step S513. If it is not less than or equal to zero, proceed to the next step.

[0204] S512. Determine whether i_x is greater than m. If it is, let i_x_regin = i_x - m; if it is not, let i_x_regin = i_x.

[0205] S513. Set the projectable state of the cell point in the i_z_regin row and i_x_regin column of the i_y_regin layer to unprojectable.

[0206] S514. Determine whether i_x is greater than or equal to i1_x. If it is, proceed to the next step; if it is less than, let i_x = i_x + 1 and jump to step S511.

[0207] S515. Determine whether i_y is greater than or equal to i1_y. If it is, proceed to the next step; if it is less than, let i_y = i_y + 1 and jump to step S57.

[0208] S516. Determine whether i_z is greater than or equal to i1_z. If it is, proceed to the next step; if it is less than, let i_z = i_z + 1 and jump to step S53.

[0209] S517, End traversal.

[0210] S6. Calculate the number of grid points in the delivery area that are in a deliverable state, denoted as n_grid_node.

[0211] S7. Determine if n_grid_node is equal to zero. If it is equal to zero, proceed to the next step; otherwise, jump to step S3.

[0212] S8. Let k_ball be the number of balls that have been placed in the current iteration step. Let k_ball = i_ball and calculate the current volume fraction. The current volume fraction is equal to the sum of the volumes of all placed balls divided by the volume of the placement area.

[0213] S9. Determine whether the current volume fraction is greater than or equal to the target volume fraction. If it is greater than or equal to the target volume fraction, proceed to step S32; if it is less than the target volume fraction, proceed to the next step.

[0214] S10. Determine if itr is greater than or equal to n_itr. If it is, jump to step S36; if it is less than, proceed to the next step.

[0215] S11. Let itr = itr + 1.

[0216] S12. Set the occupied status of all grid points within the deployment area to unoccupied.

[0217] S13, Let i_ball = 1.

[0218] S14. Traverse each grid point inside the i_ball-th ball being placed, and update the occupancy status of the corresponding grid point within the placement area to "occupied"; the traversal and update process includes the following steps:

[0219] S141. Calculate the lower and upper bound coordinates of the cube circumscribed by the ball in the z-direction, denoted as z0_out and z1_out respectively. Then z0_out = z0 - R and z1_out = z0 + R. Calculate the grid point numbers of the lower and upper bounds in the z-direction, denoted as i0_z and i1_z respectively. Then i0_z = floor(z0_out / w0) + 1 and i1_z = ceil(z1_out / w0) - 1. Calculate the layer number and row number of the grid where the center of the ball is located, denoted as center_z and center_y respectively. Then center_z = floor(z0 / w0) + 1 and center_y = floor(y0 / w0) + 1. floor() represents rounding down and ceil() represents rounding up. y0 and z0 are the y and z coordinates of the center of the ball, and R is the radius of the ball.

[0220] S142. Let i_z = i0_z;

[0221] S143. Determine if i_z is less than or equal to zero. If it is, let i_z_regin = i_z + m, and jump to step S145. If it is not less than or equal to zero, proceed to the next step.

[0222] S144. Determine if i_z is greater than m. If it is, let i_z_regin = i_z - m; if it is not, let i_z_regin = i_z.

[0223] S145. Calculate the square of the vertical distance from the center of the sphere to the current grid point layer, denoted as dz2, then dz2=(z0-i_z*w0). 2 Calculate the radius of the reference circle formed by the intersection of the current grid point layer and the deployed ball, denoted as r_i_z. Calculate the lower and upper bounds of the grid points located inside the ball in the current grid layer in the y-direction, denoted as i0_y and i1_y respectively. Then, i0_y = floor((y0-r_i_z) / w0)+1, i1_y = ceil((y0+r_i_z) / w0)-1; floor() represents rounding down, ceil() represents rounding up, y0 and z0 are the y and z coordinates of the center of the ball, and R is the radius of the ball.

[0224] S146. Let i_y = i0_y;

[0225] S147. Determine if i_y is less than or equal to zero. If it is less than or equal to zero, set i_y_regin = i_y + m and jump to step S149. If it is not less than or equal to zero, proceed to the next step.

[0226] S148. Determine if i_y is greater than m. If it is, let i_y_regin = i_y - m; if it is not, let i_y_regin = i_y.

[0227] S149. Calculate the square of the vertical distance from the center of the reference circle to the current grid point row, denoted as dy2. Then dy2 = (y0 - i_y * w0). 2 Calculate half the length of the chord formed by the intersection of the line containing the current grid point and the reference circle, denoted as dx. Calculate the lower and upper bounds of the grid points located inside the ball in the current grid row of the current grid layer in the x-direction, denoted as i0_x and i1_x respectively. Then i0_x = floor((x0-dx) / w0)+1, i1_x = ceil((x0+dx) / w0)-1; floor() means rounding down, ceil() means rounding up, and y0 and x0 are the y and x coordinates of the center of the ball.

[0228] S1410. Let i_x = i0_x;

[0229] S1411. Determine if i_x is less than or equal to zero. If it is, set i_x_regin = i_x + m and jump to step S1413. If it is not less than or equal to zero, proceed to the next step.

[0230] S1412. Determine if i_x is greater than m. If it is, let i_x_regin = i_x - m; if it is not, let i_x_regin = i_x.

[0231] S1413. Set the occupied state of the cell point in the i_z_regin row and i_x_regin column of the i_y_regin layer to occupied;

[0232] S1414. Determine if i_x is greater than or equal to i1_x. If it is, proceed to the next step; if it is less than, let i_x = i_x + 1 and jump to step S1411.

[0233] S1415. Determine whether i_y is greater than or equal to i1_y. If it is, proceed to the next step; if it is less than, let i_y = i_y + 1 and jump to step S147.

[0234] S1416. Determine if i_z is greater than or equal to i1_z. If it is, proceed to the next step; if it is less than, let i_z = i_z + 1 and jump to step S143.

[0235] S1417, End traversal.

[0236] S15. Determine if i_ball is greater than or equal to k_ball. If it is, proceed to the next step; if it is less than, set i_ball = i_ball + 1 and jump to step S14.

[0237] S16. Let i_ball = 1.

[0238] S17. Prepare a reference sphere and calculate the repulsion sphere formed by the i-th ball and the reference sphere. In this embodiment, for ease of calculation, the size of the reference sphere can be set to be the same as that of the ball. In other embodiments, its size can be set to be different from that of the ball, both of which can meet the calculation requirements of this application. The center coordinates of the repulsion sphere are equal to the center coordinates of the i-th ball, and the radius of the repulsion sphere is equal to the radius of the i-th ball plus the radius of the reference sphere.

[0239] S18. Traverse each grid point inside the repulsion sphere, count the number of grid points inside the repulsion sphere num2, count the number of grid points inside the repulsion sphere whose corresponding grid point in the delivery area is occupied num, calculate the local volume fraction vfl = num / num2 of the i-th delivery ball, and save it to the local volume fraction array; the traversal and update process includes the following steps:

[0240] S181. Calculate the lower and upper bound coordinates of the cube circumscribed by the repulsion sphere in the z-direction, denoted as z0_out and z1_out respectively. Then z0_out = z0 - R_ex and z1_out = z0 + R_ex. Calculate the grid point numbers of the lower and upper bounds in the z-direction, denoted as i0_z and i1_z respectively. Then i0_z = floor(z0_out / w0) + 1 and i1_z = ceil(z1_out / w0) - 1. Calculate the layer number and row number of the grid where the center of the repulsion sphere is located, denoted as center_z and center_y respectively. Then center_z = floor(z0 / w0) + 1 and center_y = floor(y0 / w0) + 1. floor() represents rounding down and ceil() represents rounding up. y0 and z0 are the y and z coordinates of the center of the repulsion sphere, and R_ex is the radius of the repulsion sphere.

[0241] S182. Let num = 0, let num2 = 0, let i_z = i0_z;

[0242] S183. Determine if i_z is less than or equal to zero. If it is, set i_z_regin = i_z + m and jump to step S185. If it is not less than or equal to zero, proceed to the next step.

[0243] S184. Determine if i_z is greater than m. If it is, let i_z_regin = i_z - m; if it is not, let i_z_regin = i_z.

[0244] S185. Calculate the square of the vertical distance from the center of the sphere to the current grid point layer, denoted as dz2, then dz2=(z0-i_z*w0). 2 Calculate the radius of the reference circle formed by the intersection of the current lattice point layer and the repulsion sphere, denoted as r_i_z. Calculate the lower and upper bounds of the grid points located inside the repulsion sphere in the current grid point layer in the y-direction, denoted as i0_y and i1_y respectively. Then, i0_y = floor((y0-r_i_z) / w0)+1, i1_y = ceil((y0+r_i_z) / w0)-1; floor() represents rounding down, ceil() represents rounding up, y0 and z0 are the y and z coordinates of the center of the repulsion sphere, and R_ex is the radius of the repulsion sphere;

[0245] S186. Let i_y = i0_y;

[0246] S187. Determine if i_y is less than or equal to zero. If it is less than or equal to zero, set i_y_regin = i_y + m and jump to step S189. If it is not less than or equal to zero, proceed to the next step.

[0247] S188. Determine if i_y is greater than m. If it is, let i_y_regin = i_y - m; if it is not, let i_y_regin = i_y.

[0248] S189. Calculate the square of the vertical distance from the center of the reference circle to the current grid point row, denoted as dy2. Then dy2 = (y0 - i_y * w0). 2 Calculate half the length of the chord formed by the intersection of the line containing the current grid point and the reference circle, denoted as dx. Calculate the lower and upper bounds of the grid points located inside the repulsion sphere in the current grid row of the current grid point layer in the x-direction, denoted as i0_x and i1_x respectively. Then i0_x = floor((x0-dx) / w0)+1, i1_x = ceil((x0+dx) / w0)-1; floor() means rounding down, ceil() means rounding up, and y0 and x0 are the y and x coordinates of the center of the repulsion sphere;

[0249] S1810, Let i_x = i0_x;

[0250] S1811. Determine if i_x is less than or equal to zero. If it is, set i_x_regin = i_x + m and jump to step S1813. If it is not less than or equal to zero, proceed to the next step.

[0251] S1812. Determine if i_x is greater than m. If it is, let i_x_regin = i_x - m; if it is not, let i_x_regin = i_x.

[0252] S1813. Let num2 = num2 + 1. If the cell in the i_z_regin row and i_x_regin column of the i_y_regin layer is occupied, then let num = num + 1; otherwise, num remains unchanged.

[0253] S1814. Determine if i_x is greater than or equal to i1_x. If it is, proceed to the next step; if it is less than, let i_x = i_x + 1 and jump to step S1811.

[0254] S1815. Determine if i_y is greater than or equal to i1_y. If it is, proceed to the next step; if it is less than, let i_y = i_y + 1 and jump to step S187.

[0255] S1816. Determine if i_z is greater than or equal to i1_z. If it is, proceed to the next step; if it is less than, let i_z = i_z + 1 and jump to step S183.

[0256] S1817. Calculate the local volume fraction vfl = num / num2 and save it to the local volume fraction array, then end the traversal.

[0257] S19. Determine whether i_ball is greater than or equal to k_ball. If it is, proceed to the next step; if it is less than, let i_ball = i_ball + 1 and jump to step S17.

[0258] S20. Based on the local volume fractions of all k_ball balls, calculate the maximum value vfl_max and the minimum value vfl_min, and calculate the baseline local volume fraction vfl_base = (1-alpha)*vfl_min + alpha*vfl_max; where alpha is an adjustment coefficient, and its value ranges from 0.35 to 0.45. In this application, it can be 0.41.

[0259] S21. Set the retention status of all thrown balls to non-retention.

[0260] S22, Let i_ball = 1.

[0261] S23. Determine whether the local volume fraction of the i_ball ball is greater than or equal to vfl_base. If it is greater than or equal to vfl_base, proceed to the next step; if it is less than vfl_base, proceed to step S25.

[0262] S24. Mark the retention status of the i_ball ball as retained.

[0263] S25. Determine if i_ball is greater than or equal to k_ball. If it is, proceed to the next step; if it is less than, let i_ball = i_ball + 1 and jump to step S23.

[0264] S26. Delete the balls whose retention status is not retained from the ball set, update the ball index in the new set, and update the value of k_ball to the number of retained balls.

[0265] S27. Set the drop status of all grid points in the drop area to dropable.

[0266] S28. Let i_ball = 1, and prepare to throw a new ball.

[0267] S29. Calculate the repulsion sphere formed by the i-th ball and the newly placed ball; where the center coordinates of the repulsion sphere are equal to the center coordinates of the i-th ball, and the radius of the repulsion sphere is equal to the radius of the i-th ball plus the radius of the newly placed ball.

[0268] S30. Traverse each grid point located inside the repulsion sphere and update the throwable status of the corresponding grid point in the throwing area to unthrowable; the traversal and update process is the same as step S5, and will not be described here.

[0269] S31. Determine whether i_ball is greater than or equal to k_ball. If it is, jump to step S6; if it is less than, let i_ball = i_ball + 1, and jump to step S29.

[0270] S32, Let i_ball = 1.

[0271] S33. Determine whether the i_th ball intersects with the boundary of the throwing area. If they intersect, generate a replacement ball for this ball and add it to the replacement ball set; if they do not intersect, proceed to step S34. The generation of the replacement ball includes the following steps:

[0272] S331. Let the x, y, and z coordinates of the center of the i-th ball be x0, y0, and z0, respectively. Let xyz_min be equal to the ball radius and xyz_max be equal to the region side length minus the ball radius. Execute S332, S334, and S336 in parallel.

[0273] S332. If z0 is less than xyz_min, let z1 equal z0 plus the region side length, let z3 equal zero, and proceed to the next step; if z0 is greater than xyz_max, let z1 equal z0 minus the region side length, let z3 equal the region side length, and proceed to the next step; if z0 is greater than or equal to xyz_min and less than or equal to xyz_max, let z3 equal -1, and jump to step S338.

[0274] S333. A new supplementary ball is generated with center coordinates (x0, y0, z1) and radius equal to the radius of the thrown ball. The center coordinates and radius of this supplementary ball are recorded and added to the supplementary ball set. Then, proceed to step S338.

[0275] S334. If y0 is less than xyz_min, let y1 equal y0 plus the region side length, let y3 equal zero, and proceed to the next step; if y0 is greater than xyz_max, let y1 equal y0 minus the region side length, let y3 equal the region side length, and proceed to the next step; if y0 is greater than or equal to xyz_min and less than or equal to xyz_max, let y3 equal -1, and jump to step S338.

[0276] S335. A new supplementary ball is generated with center coordinates (x0, y1, z0) and radius equal to the radius of the thrown ball. The center coordinates and radius of this supplementary ball are recorded and added to the supplementary ball set. Then, proceed to step S338.

[0277] S336. If x0 is less than xyz_min, set x1 to x0 plus the region side length, set x3 to zero, and proceed to the next step; if x0 is greater than xyz_max, set x1 to x0 minus the region side length, set x3 to the region side length, and proceed to the next step; if x0 is greater than or equal to xyz_min and less than or equal to xyz_max, set x3 to -1, and jump to step S338.

[0278] S337. A new supplementary ball is generated with center coordinates (x1, y0, z0) and radius equal to the radius of the thrown ball. The center coordinates and radius of this supplementary ball are recorded and added to the supplementary ball set. Then, proceed to step S338.

[0279] S338, Wait for S332-S337 to complete execution;

[0280] S339, execute S3310, S3314, S3318 and S3322 in parallel;

[0281] S3310. If x3 is greater than or equal to zero and y3 is greater than or equal to zero, proceed to the next step; otherwise, jump to step S3326.

[0282] S3311. Using x3 and y3 as reference points, calculate the distance d30 from the reference point to (x0, y0).

[0283] S3312. If d30 is less than the radius of the ball being thrown, proceed to the next step; if d30 is greater than or equal to the radius of the ball being thrown, jump to step S3326.

[0284] S3313. A new supplementary ball is generated with center coordinates (x1, y1, z0) and radius equal to the radius of the thrown ball. Record the center coordinates and radius of this supplementary ball and add it to the supplementary ball set. Jump to step S3326.

[0285] S3314. If y3 is greater than or equal to zero and z3 is greater than or equal to zero, proceed to the next step; otherwise, jump to step S3326.

[0286] S3315. Using the x and y coordinates of y3 and z3 as reference points, calculate the distance d30 from the reference point to (y0, z0).

[0287] S3316. If d30 is less than the radius of the ball being thrown, proceed to the next step; if d30 is greater than or equal to the radius of the ball being thrown, jump to step S3326.

[0288] S3317. A new supplementary ball is generated with center coordinates (x0, y1, z1) and radius equal to the radius of the thrown ball. Record the center coordinates and radius of this supplementary ball and add it to the supplementary ball set. Jump to step S3326.

[0289] S3318. If z3 is greater than or equal to zero and x3 is greater than or equal to zero, proceed to the next step; otherwise, jump to step S3326.

[0290] S3319. Using x3 and z3 as reference points, calculate the distance d30 from the reference point to (x0, z0).

[0291] S3320. If d30 is less than the radius of the ball being thrown, proceed to the next step; if d30 is greater than or equal to the radius of the ball being thrown, jump to step S3326.

[0292] S3321. A new supplementary ball is generated with center coordinates (x1, y0, z1) and radius equal to the radius of the thrown ball. Record the center coordinates and radius of this supplementary ball and add it to the supplementary ball set. Jump to step S3326.

[0293] S3322. If x3 is greater than or equal to zero, y3 is greater than or equal to zero, and z3 is greater than or equal to zero, proceed to the next step; otherwise, jump to step S3326.

[0294] S3323. Using x3, y3, z3 as reference points, calculate the distance d30 from the reference points to (x0, y0, z0).

[0295] S3324. If d30 is less than the radius of the ball being thrown, proceed to the next step; if d30 is greater than or equal to the radius of the ball being thrown, jump to step S3326.

[0296] S3325. A new supplementary ball is generated with center coordinates (x1, y1, z1) and radius equal to the radius of the thrown ball. Record the center coordinates and radius of this supplementary ball and add it to the supplementary ball set. Jump to step S3326.

[0297] S3326, Wait for S3310-S3325 to complete.

[0298] S34. Determine whether i_ball is greater than or equal to k_ball. If it is, proceed to the next step; if it is less than, set i_ball = i_ball + 1 and jump to step S33.

[0299] S35. Save and output the set of balls to be thrown and the set of balls to be replenished, output the successful throwing flag, and jump to step S37.

[0300] S36, Output delivery failure flag.

[0301] S37, End the deployment.

[0302] In this embodiment, during the traversal of grid points in S5 / S14 / S18, all grid points inside the sphere are directly calculated, eliminating the need to determine whether each grid point is inside the sphere, thus greatly improving computational efficiency. The use of "local volume fraction" in steps S21-S25 further improves efficiency. Moreover, in this application, after successful placement, each placed ball in the set is checked for intersection with the boundary, rather than checking for intersection with the boundary after each placed ball is placed, which also improves computational efficiency.

[0303] Example 2

[0304] like Figure 2 and Figure 3 As shown in the embodiment of the present invention, the method for generating a high volume fraction particle random distribution based on a grid is basically the same as that in Embodiment 1, and specifically involves:

[0305] First, set the input parameters: the radius of the ball particle R, the ratio of the side length of the placement area to the radius of the ball particle δ, the target volume fraction, and the upper limit of the total number of iterations n_itr. Then, the length, width, and height of the placement area are all width = R * δ. Assume that the x-axis of the spatial rectangular coordinate system is horizontal to the right, the y-axis is horizontal to the back, and the z-axis is vertical to the up. The sides of the placement area are parallel to the coordinate axes, and its lower left front corner is located at the origin.

[0306] The generation method includes the following steps:

[0307] S1. The delivery area is discrete. A sufficiently large three-dimensional background grid is established to cover the delivery area. The delivery area is evenly divided into m×m×m grids. Then the grid width w0 = width / m. All grids in the delivery area are set to be deliverable. The grids corresponding to all grids in the delivery area are called the grids in the delivery area.

[0308] S2. Let itr = 1, and let i_ball = 0.

[0309] S3. Randomly select a cell from all cells in the drop zone that are dropable. Randomly select a point in this cell as the drop point. Let i_ball = i_ball + 1. Drop the i_ball ball so that the center of the dropped ball is located at this drop point. Record the coordinates of the center and radius of the dropped ball and add it to the drop ball set.

[0310] S4. Prepare to release new balls. Calculate the repulsion sphere formed by the i-th released ball and the new ball. The center coordinates of the repulsion sphere are equal to the center coordinates of the i-th released ball, and the radius of the repulsion sphere is equal to the radius of the i-th released ball plus the radius of the new ball.

[0311] S5. Traverse each cell that intersects with the repulsion ball, and update the throwable status of the corresponding cell within the throwing area to unthrowable. The traversal and update process includes the following steps:

[0312] S51. Calculate the lower and upper coordinates of the circumscribed cube of the exclusion sphere in the z direction, denoted as z0_out and z1_out respectively. Then z0_out = z0 - R_ex, z1_out = z0 + R_ex. Calculate the lower and upper grid numbers in the z direction, denoted as i0_z and i1_z respectively. Then i0_z = floor(z0_out / w0) + 1, i1_z = ceil(z1_out / w0). Also, calculate the layer number and row number of the grid where the center of the exclusion sphere is located, denoted as center_z and center_y respectively. Then center_z = floor(z0 / w0) + 1, center_y = floor(y0 / w0) + 1; floor() represents rounding down, ceil() represents rounding up, y0 and z0 are the y and z coordinates of the center of the exclusion sphere, and R_ex is the radius of the exclusion sphere.

[0313] S52. Let i_z = i0_z;

[0314] S53. Determine whether i_z is less than or equal to zero. If it is, let i_z_regin = i_z + m, and jump to step S55; if it is not, proceed to the next step;

[0315] S54. Determine whether i_z is greater than m. If it is, let i_z_regin = i_z - m; if it is not, let i_z_regin = i_z;

[0316] S55. Determine the reference plane of the current grid layer and calculate the z coordinate z_ref of the reference plane. If i_z < center_z, the upper surface of the current grid layer is the reference plane of the current grid layer, z_ref = i_z * w0. If i_z = center_z, the horizontal plane where the center of the sphere is located is the reference plane of the current grid layer, z_ref = 0. If i_z > center_z, the lower surface of the current grid layer is the reference plane of the current grid layer, z_ref = (i_z - 1) * w0. Calculate the square of the minimum vertical distance from the center of the sphere to the reference plane of the current grid layer, denoted as dz2, dz2 = (z0 - z_ref) 2 , calculate the radius of the reference circle formed by the intersection of the reference plane of the current grid layer and the exclusion sphere, denoted as r_i_z. Then Calculate the lower and upper grid numbers in the y direction of the grids that intersect with the exclusion sphere in the current grid layer, denoted as i0_y and i1_y respectively. Then i0_y = floor((y0 - r_i_z) / w0) + 1, i1_y = ceil((y0 + r_i_z) / w0); floor() represents rounding down, ceil() represents rounding up, y0 and z0 are the y and z coordinates of the center of the exclusion sphere, and R_ex is the radius of the exclusion sphere.

[0317] S56. Let \(i_y = i0_y\);

[0318] S57. Determine whether \(i_y\) is less than or equal to zero. If it is, let \(i_y\_regin = i_y + m\) and jump to step S59; if not, proceed to the next step;

[0319] S58. Determine whether \(i_y\) is greater than \(m\). If it is, let \(i_y\_regin = i_y - m\); if not, let \(i_y\_regin = i_y\);

[0320] S59. Determine the reference plane of the current grid row and calculate the y - coordinate \(y\_ref\) of the reference plane. If \(i_y < center_y\), the back surface of the current grid row is the reference plane of the current grid row, and \(y\_ref = i_y * w0\). If \(i_y = center_y\), the vertical plane passing through the center of the reference circle and parallel to the front and back surfaces of the grid is the reference plane of the current grid row, and \(y\_ref = 0\). If \(i_y > center_y\), the front surface of the current grid row is the reference plane of the current grid row, and \(y\_ref=(i_y - 1)*w0\). Calculate the square of the perpendicular distance from the center of the reference circle to the reference plane of the current grid row, denoted as \(dy2\), \(dy2=(y0 - y\_ref)\) 2 , calculate half of the length of the chord formed by the intersection of the intersection line of the reference plane of the current grid layer and the reference plane of the current grid row with the reference circle, denoted as \(dx\), then Calculate the lower and upper grid numbers in the x - direction of the grids intersecting with the exclusion sphere in the current grid row of the current grid layer, denoted as \(i0_y\) and \(i1_y\) respectively. Then \(i0_x = floor((x0 - dx) / w0)+1\), \(i1_x = ceil((x0 + dx) / w0)\); \(floor()\) represents rounding down, \(ceil()\) represents rounding up; \(y0\) and \(x0\) are the y and x coordinates of the center of the exclusion sphere;

[0321] S510. Let \(i_x = i0_x\);

[0322] S511. Determine whether \(i_x\) is less than or equal to zero. If it is, let \(i_x\_regin = i_x + m\) and jump to step S513; if not, proceed to the next step;

[0323] S512. Determine whether \(i_x\) is greater than \(m\). If it is, let \(i_x\_regin = i_x - m\); if not, let \(i_x\_regin = i_x\);

[0324] S513. Set the castable state of the grid at the \(i_z\_regin\) - th layer, \(i_y\_regin\) - th row, and \(i_x\_regin\) - th column to non - castable;

[0325] S514. Determine whether i_x is greater than or equal to i1_x. If it is, proceed to the next step; if it is less than, let i_x = i_x + 1 and jump to step S511.

[0326] S515. Determine whether i_y is greater than or equal to i1_y. If it is, proceed to the next step; if it is less than, let i_y = i_y + 1 and jump to step S57.

[0327] S516. Determine whether i_z is greater than or equal to i1_z. If it is, proceed to the next step; if it is less than, let i_z = i_z + 1 and jump to step S53.

[0328] S517, End traversal.

[0329] S6. Calculate the number of grids in the drop zone that are dropable, denoted as n_grid_node.

[0330] S7. Determine if n_grid_node is equal to zero. If it is equal to zero, proceed to the next step; otherwise, jump to step S3.

[0331] S8. Let k_ball be the number of balls that have been placed in the current iteration step. Let k_ball = i_ball, and calculate the current volume fraction. The current volume fraction is equal to the sum of the volumes of all placed balls divided by the volume of the placement area.

[0332] S9. Determine whether the current volume fraction is greater than or equal to the target volume fraction. If it is greater than or equal to the target volume fraction, proceed to step S32; if it is less than the target volume fraction, proceed to the next step.

[0333] S10. Determine if itr is greater than or equal to n_itr. If it is, jump to step S36; if it is less than, proceed to the next step.

[0334] S11. Let itr = itr + 1.

[0335] S12. Set the occupied status of all cells in the distribution area to unoccupied.

[0336] S13, Let i_ball = 1.

[0337] S14. Traverse each cell that intersects with the i_ball-th ball, and update the occupancy status of the corresponding cell within the throwing area to "occupied". The traversal and update process includes the following steps:

[0338] S141. Calculate the lower and upper coordinates of the circumscribed cube of the placed ball in the z direction, denoted as z0_out and z1_out respectively. Then z0_out = z0 - R, z1_out = z0 + R. Calculate the lower and upper grid numbers in the z direction, denoted as i0_z and i1_z respectively. Then i0_z = floor(z0_out / w0) + 1, i1_z = ceil(z1_out / w0). Also, calculate the layer number of the grid layer where the center of the placed ball is located and the row number of the grid row, denoted as center_z and center_y respectively. Then center_z = floor(z0 / w0) + 1, center_y = floor(y0 / w0) + 1; floor() represents rounding down, ceil() represents rounding up, y0 and z0 are the y and z coordinates of the center of the placed ball, and R is the radius of the placed ball.

[0339] S142. Let i_z = i0_z;

[0340] S143. Determine whether i_z is less than or equal to zero. If it is, let i_z_regin = i_z + m, and jump to step S145; if it is not, proceed to the next step;

[0341] S144. Determine whether i_z is greater than m. If it is, let i_z_regin = i_z - m; if it is not, let i_z_regin = i_z;

[0342] S145. Determine the reference plane of the current grid layer and calculate the z coordinate z_ref of the reference plane. If i_z < center_z, the upper surface of the current grid layer is the reference plane of the current grid layer, and z_ref = i_z * w0. If i_z = center_z, the horizontal plane where the ball center is located is the reference plane of the current grid layer, and z_ref = 0. If i_z > center_z, the lower surface of the current grid layer is the reference plane of the current grid layer, and z_ref = (i_z - 1) * w0. Calculate the square of the minimum vertical distance from the ball center to the reference plane of the current grid layer, denoted as dz2. Then dz2 = (z0 - z_ref) 2 , calculate the radius of the reference circle formed by the intersection of the reference plane of the current grid layer and the placed ball, denoted as r_i_z. Then Calculate the lower and upper grid numbers in the y direction of the grids intersecting with the placed ball in the current grid layer, denoted as i0_y and i1_y respectively. Then i0_y = floor((y0 - r_i_z) / w0) + 1, i1_y = ceil((y0 + r_i_z) / w0); floor() represents rounding down, ceil() represents rounding up, y0 and z0 are the y and z coordinates of the center of the placed ball, and R is the radius of the placed ball.

[0343] S146. Let \(i_y = i0_y\);

[0344] S147. Determine whether \(i_y\) is less than or equal to zero. If it is, let \(i_y\_regin = i_y + m\) and jump to step S149; if not, proceed to the next step;

[0345] S148. Determine whether \(i_y\) is greater than \(m\). If it is, let \(i_y\_regin = i_y - m\); if not, let \(i_y\_regin = i_y\);

[0346] S149. Determine the reference plane of the current grid row and calculate the y - coordinate \(y\_ref\) of the reference plane. If \(i_y < center_y\), the back surface of the current grid row is the reference plane of the current grid row, and \(y\_ref = i_y * w0\); if \(i_y = center_y\), the vertical plane where the center of the reference circle is located and parallel to the front and back surfaces of the grid is the reference plane of the current grid row, and \(y\_ref = 0\); if \(i_y > center_y\), the front surface of the current grid row is the reference plane of the current grid row, and \(y\_ref=(i_y - 1)*w0\). Calculate the square of the perpendicular distance from the center of the reference circle to the reference plane of the current grid row, denoted as \(dy2\), \(dy2=(y0 - y\_ref)\) 2 , calculate half of the length of the chord formed by the intersection of the intersection line of the reference plane of the current grid layer and the reference plane of the current grid row with the reference circle, denoted as \(dx\), then Calculate the lower and upper grid numbers in the x - direction of the grids that intersect with the dropped ball in the current grid row of the current grid layer, denoted as \(i0_y\) and \(i1_y\) respectively. Then \(i0_x = floor((x0 - dx) / w0)+1\), \(i1_x = ceil((x0 + dx) / w0)\); \(floor()\) represents rounding down, \(ceil()\) represents rounding up; \(y0\) and \(x0\) are the y and x coordinates of the center of the dropped ball;

[0347] S1410. Let \(i_x = i0_x\);

[0348] S1411. Determine whether \(i_x\) is less than or equal to zero. If it is, let \(i_x\_regin = i_x + m\) and jump to step S1413; if not, proceed to the next step;

[0349] S1412. Determine whether \(i_x\) is greater than \(m\). If it is, let \(i_x\_regin = i_x - m\); if not, let \(i_x\_regin = i_x\);

[0350] S1413. Set the occupancy status of the grid at the \(i_z\_regin\) - th layer, \(i_y\_regin\) - th row, and \(i_x\_regin\) - th column to occupied;

[0351] S1414. Determine if i_x is greater than or equal to i1_x. If it is, proceed to the next step; if it is less than, let i_x = i_x + 1 and jump to step S1411.

[0352] S1415. Determine whether i_y is greater than or equal to i1_y. If it is, proceed to the next step; if it is less than, let i_y = i_y + 1 and jump to step S147.

[0353] S1416. Determine if i_z is greater than or equal to i1_z. If it is, proceed to the next step; if it is less than, let i_z = i_z + 1 and jump to step S143.

[0354] S1417, End traversal.

[0355] S15. Determine if i_ball is greater than or equal to k_ball. If it is, proceed to the next step; if it is less than, set i_ball = i_ball + 1 and jump to step S14.

[0356] S16. Let i_ball = 1.

[0357] S17. Prepare a reference sphere and calculate the repulsion sphere formed by the i-th thrown ball and the reference sphere. The coordinates of the center of the repulsion sphere are equal to the coordinates of the center of the i-th thrown ball, and the radius of the repulsion sphere is equal to the radius of the i-th thrown ball plus the radius of the reference sphere.

[0358] S18. Traverse each cell intersecting with the repulsion ball, count the number of cells intersecting with the repulsion ball (num2), count the number of cells intersecting with the repulsion ball and whose corresponding cell within the placement area is occupied (num), calculate the local volume fraction vfl = num / num2 for the i-th placement ball, and save it to the local volume fraction array. The traversal and update process is as follows:

[0359] S181. Calculate the lower and upper coordinates of the circumscribed cube of the exclusion sphere in the z direction, denoted as z0_out and z1_out respectively. Then z0_out = z0 - R_ex, z1_out = z0 + R_ex. Calculate the lower and upper grid numbers in the z direction, denoted as i0_z and i1_z respectively. Then i0_z = floor(z0_out / w0) + 1, i1_z = ceil(z1_out / w0). Also, calculate the layer number of the grid layer where the center of the exclusion sphere is located and the row number of the grid row, denoted as center_z and center_y respectively. Then center_z = floor(z0 / w0) + 1, center_y = floor(y0 / w0) + 1. floor() represents rounding down, ceil() represents rounding up. y0 and z0 are the y and z coordinates of the center of the exclusion sphere, and R_ex is the radius of the exclusion sphere.

[0360] S182. Let num = 0, let num2 = 0, and let i_z = i0_z.

[0361] S183. Determine whether i_z is less than or equal to zero. If it is, let i_z_regin = i_z + m, and jump to step S185. If it is not, proceed to the next step.

[0362] S184. Determine whether i_z is greater than m. If it is, let i_z_regin = i_z - m. If it is not, let i_z_regin = i_z.

[0363] S185. Determine the reference plane of the current grid layer and calculate the z coordinate z_ref of the reference plane. If i_z < center_z, the upper surface of the current grid layer is the reference plane of the current grid layer, and z_ref = i_z * w0. If i_z = center_z, the horizontal plane where the center of the sphere is located is the reference plane of the current grid layer, and z_ref = 0. If i_z > center_z, the lower surface of the current grid layer is the reference plane of the current grid layer, and z_ref = (i_z - 1) * w0. Calculate the square of the minimum vertical distance from the center of the sphere to the reference plane of the current grid layer, denoted as dz2. dz2 = (z0 - z_ref) 2 , calculate the radius of the reference circle formed by the intersection of the reference plane of the current grid layer and the exclusion sphere, denoted as r_i_z. Then Calculate the lower and upper grid numbers in the y - direction of the grids that intersect with the exclusion sphere in the current grid layer, denoted as \(i0\_y\) and \(i1\_y\) respectively. Then \(i0\_y=\lfloor(y0 - r\_i\_z) / w0\rfloor + 1\), \(i1\_y=\lceil(y0 + r\_i\_z) / w0\rceil\); \(\lfloor\)\(\rfloor\) represents floor function (rounding down), \(\lceil\)\(\rceil\) represents ceiling function (rounding up), \(y0\) and \(z0\) are the y and z coordinates of the center of the exclusion sphere, and \(R\_ex\) is the radius of the exclusion sphere.

[0364] S186. Let \(i\_y = i0\_y\);

[0365] S187. Determine whether \(i\_y\) is less than or equal to zero. If it is, let \(i\_y\_regin = i\_y + m\) and jump to step S189; if not, proceed to the next step;

[0366] S188. Determine whether \(i\_y\) is greater than \(m\). If it is, let \(i\_y\_regin = i\_y - m\); if not, let \(i\_y\_regin = i\_y\);

[0367] S189. Determine the reference plane of the current grid row and calculate the y - coordinate \(y\_ref\) of the reference plane. If \(i\_y\lt center\_y\), the back surface of the current grid row is the reference plane of the current grid row, \(y\_ref = i\_y*w0\). If \(i\_y = center\_y\), the vertical plane where the center of the reference circle lies and is parallel to the front and back surfaces of the grid is the reference plane of the current grid row, \(y\_ref = 0\). If \(i\_y\gt center\_y\), the front surface of the current grid row is the reference plane of the current grid row, \(y\_ref=(i\_y - 1)*w0\). Calculate the square of the perpendicular distance from the center of the reference circle to the reference plane of the current grid row, denoted as \(dy2\), \(dy2=(y0 - y\_ref)\) 2 Calculate the half - length of the chord formed by the intersection of the intersection line of the reference plane of the current grid layer and the reference plane of the current grid row with the reference circle, denoted as \(dx\). Then Calculate the lower and upper grid numbers in the x - direction of the grids that intersect with the exclusion sphere in the current grid row of the current grid layer, denoted as \(i0\_x\) and \(i1\_x\) respectively. Then \(i0\_x=\lfloor(x0 - dx) / w0\rfloor+1\), \(i1\_x=\lceil(x0 + dx) / w0\rceil\); \(\lfloor\)\(\rfloor\) represents floor function (rounding down), \(\lceil\)\(\rceil\) represents ceiling function (rounding up), \(y0\) and \(x0\) are the y and x coordinates of the center of the exclusion sphere;

[0368] S1810. Let \(i\_x = i0\_x\);

[0369] S1811. Determine whether \(i\_x\) is less than or equal to zero. If it is, let \(i\_x\_regin = i\_x + m\) and jump to step S1813; if not, proceed to the next step;

[0370] S1812. Determine if i_x is greater than m. If it is, let i_x_regin = i_x - m; if it is not, let i_x_regin = i_x.

[0371] S1813. Let num2 = num2 + 1. If the cell in the i_z_regin row and i_x_regin column of the i_y_regin layer is occupied, then let num = num + 1; otherwise, num remains unchanged.

[0372] S1814. Determine if i_x is greater than or equal to i1_x. If it is, proceed to the next step; if it is less than, let i_x = i_x + 1 and jump to step S1811.

[0373] S1815. Determine if i_y is greater than or equal to i1_y. If it is, proceed to the next step; if it is less than, let i_y = i_y + 1 and jump to step S187.

[0374] S1816. Determine if i_z is greater than or equal to i1_z. If it is, proceed to the next step; if it is less than, let i_z = i_z + 1 and jump to step S183.

[0375] S1817. Calculate the local volume fraction vfl = num / num2 and save it to the local volume fraction array, then end the traversal.

[0376] S19. Determine whether i_ball is greater than or equal to k_ball. If it is, proceed to the next step; if it is less than, let i_ball = i_ball + 1 and jump to step S17.

[0377] S20. Based on the local volume fractions of all k_ball balls, calculate the maximum value vfl_max and the minimum value vfl_min, and calculate the baseline local volume fraction vfl_base = (1-alpha)*vfl_min + alpha*vfl_max.

[0378] S21. Set the retention status of all thrown balls to non-retention.

[0379] S22, Let i_ball = 1.

[0380] S23. Determine whether the local volume fraction of the i_ball ball is greater than or equal to vfl_base. If it is greater than or equal to vfl_base, proceed to the next step; if it is less than vfl_base, proceed to step S25.

[0381] S24. Mark the retention status of the i_ball ball as retained.

[0382] S25. Determine whether i_ball is greater than or equal to k_ball. If it is, proceed to the next step; if it is less than, let i_ball = i_ball + 1 and jump to step S23.

[0383] S26. Delete the balls whose retention status is not retained from the ball set, update the ball index in the new set, and update the value of k_ball to the number of retained balls.

[0384] S27. Set the drop status of all cells in the drop area to dropable.

[0385] S28. Let i_ball = 1, and prepare to throw a new ball.

[0386] S29. Calculate the repulsion sphere formed by the i-th ball and the newly placed ball. The center coordinates of the repulsion sphere are equal to the center coordinates of the i-th ball, and the radius of the repulsion sphere is equal to the radius of the i-th ball plus the radius of the newly placed ball.

[0387] S30. Traverse each cell that intersects with the repulsion ball, and update the throwable state of the corresponding cell in the throwing area to unthrowable. This traversal and update process is the same as S5.

[0388] S31. Determine whether i_ball is greater than or equal to k_ball. If it is, jump to step S6; if it is less than, let i_ball = i_ball + 1, and jump to step S29.

[0389] S32, Let i_ball = 1.

[0390] S33. Determine whether the i_th ball intersects with the boundary of the throwing area. If they intersect, generate a supplementary ball for this ball and add it to the supplementary ball set; otherwise, proceed to step S34. In this embodiment, the steps for generating the supplementary ball are the same as in Embodiment 1, and will not be described in detail here.

[0391] S34. Determine whether i_ball is greater than or equal to k_ball. If it is, proceed to the next step; if it is less than, set i_ball = i_ball + 1 and jump to step S33.

[0392] S35. Save and output the set of balls to be thrown and the set of balls to be replenished, output the successful throwing flag, and jump to step S37.

[0393] S36, Output delivery failure flag.

[0394] S37, End the deployment.

[0395] Test case

[0396] This experimental example is based on the generation results obtained using the method in Example 2:

[0397] This experimental example demonstrates the specific calculation process. The algorithm can be implemented using common software such as C / VB / MATLAB / Python after self-programming. Those skilled in the art should understand that this process is only for illustrating the algorithm and is not intended to limit the scope of this application.

[0398] First, set the input parameters: sphere radius R = 1.6, the ratio of the side length of the placement area to the sphere radius δ = 5, the target volume fraction 0.2, and the upper limit of the total number of iterations n_itr = 5. Then, the length, width, and height of the placement area are all width = R * δ = 8. Assume that the x-axis of the spatial rectangular coordinate system is horizontal to the right, the y-axis is horizontal to the back, and the z-axis is vertical to the top. Each side of the placement area is parallel to the coordinate axes, and its lower left front corner is located at the origin. A schematic diagram of the placement area is shown below. Figure 4 As shown.

[0399] S1. Discretize the delivery area. Establish a sufficiently large 3D background mesh to cover the delivery area, such that m = 8. Divide the delivery area evenly into 8×8×8 meshes, then the mesh width w0 = width / m = 1. Establish a delivery-ready state for all meshes within the delivery area and set them all to delivery-ready. The meshes corresponding to all meshes within the delivery area are called the meshes within the delivery area. The schematic diagram of the result of discretizing the delivery area through the background mesh is shown below. Figure 5 As shown.

[0400] S2. Let itr = 1, and let i_ball = 0.

[0401] In the first run of S3, all cells are in a drop-able state. A cell is randomly selected, for example, the cell in the 5th row, 4th column of the 5th layer. Within this cell, a point is randomly chosen as the drop point, for example, the point with coordinates (3.936, 4.472, 4.096). Let i_ball = i_ball + 1 = 1, and drop the first ball. The result after dropping the ball is illustrated in the diagram below. Figure 6 As shown, the left and right images are 3D views from different perspectives.

[0402] In the first run of S4, preparing to release new balls, calculate the repulsion sphere formed by the i-th released ball and the newly released ball; the repulsion sphere of the first released ball has its center coordinates (3.936, 4.472, 4.096) and radius 1.6 + 1.6 = 3.2; a schematic diagram of this repulsion sphere is shown below. Figure 7 As shown.

[0403] First time running S5.

[0404] S51. Calculate the lower and upper coordinates of the circumscribed cube of the exclusion sphere in the z direction, denoted as z0_out and z1_out respectively. z0_out = z0 - R_ex = 4.096 - 3.2 = 0.896, z1_out = z0 + R_ex = 4.096 + 3.2 = 7.296. Calculate the lower and upper grid numbers in the z direction, denoted as i0_z and i1_z respectively. Then i0_z = floor(z0_out / w0) + 1 = 1, i1_z = ceil(z1_out / w0) = 8. Calculate the layer number of the grid layer where the center of the exclusion sphere is located and the row number of the grid row, denoted as center_z and center_y respectively. center_z = floor(z0 / w0) + 1 = floor(4.096 / 1) + 1 = 5, center_y = floor(y0 / w0) + 1 = floor(4.472 / 1) + 1 = 5.

[0405] S52. Let i_z = i0_z = 1.

[0406] S53. Determine whether i_z is less than or equal to zero. If it is not, proceed to the next step.

[0407] S54. Determine whether i_z is greater than m. If it is not, let i_z_regin = i_z = 1.

[0408] S55. Determine the reference plane of the current grid layer and calculate the z coordinate z_ref of the reference plane. If i_z < center_z, then the upper surface of the current grid layer is the reference plane of the current grid layer, z_ref = i_z * w0 = 1 * 1 = 1. Calculate the square of the minimum vertical distance from the center of the sphere to the reference plane of the current grid layer, denoted as dz2. dz2 = (z0 - z_ref) 2 =(4.096 - 1) 2 =9.585216. Calculate the radius of the circle formed by the intersection of the reference plane of the current grid layer and the exclusion sphere, denoted as r_i_z. Then Calculate the lower and upper grid numbers in the y direction of the grid intersecting with the exclusion sphere in the current grid layer, denoted as i0_y and i1_y respectively. Then i0_y = floor((y0 - r_i_z) / w0) + 1 = floor((4.472 - 0.8092) / 1) + 1 = 4, i1_y = ceil((y0 + r_i_z) / w0) = ceil((4.472 + 0.8092) / 1) = 6.

[0409] S56. Let i_y = i0_y = 4.

[0410] S57. Determine whether i_y is less than or equal to zero. If it is not, proceed to the next step.

[0411] S58. Determine whether \(i_y\) is greater than \(m\). If not, set \(i_y\_regin = i_y = 4\).

[0412] S59. Determine the reference plane of the current grid row and calculate the y - coordinate \(y\_ref\) of the reference plane. If \(i_y < center_y\), then the back surface of the current grid row is the reference plane of the current grid row, and \(y\_ref = i_y * w0 = 4 * 1 = 1\). Calculate the square of the perpendicular distance from the center of the reference circle to the reference plane of the current grid row, denoted as \(dy2\), where \(dy2=(y0 - y\_ref)\) 2 \(=(4.472 - 4)\) 2 \(=0.2227\). Calculate half of the length of the chord formed by the intersection of the intersection line of the reference plane of the current grid layer and the reference plane of the current grid row with the reference circle, denoted as \(dx\). Then Calculate the lower and upper grid numbers in the x - direction of the grids in the current grid row of the current grid layer that intersect with the exclusion sphere, denoted as \(i0_y\) and \(i1_y\) respectively. Then \(i0_x = floor((x0 - dx) / w0)+1 = floor((3.936 - 0.6573) / 1)+1 = 4\), \(i1_x = ceil((x0 + dx) / w0)=ceil((3.936 + 0.6573) / 1)=5\).

[0413] S510. Set \(i_x = i0_x = 4\).

[0414] S511. Determine whether \(i_x\) is less than or equal to zero. If not, proceed to the next step.

[0415] S512. Determine whether \(i_x\) is greater than \(m\). If not, set \(i_x\_regin = i_x = 4\).

[0416] S513. Set the throwable state of the grid at the 4th row, 4th column of the 1st layer to non - throwable. The schematic diagram of the positional relationship between the grid at the 4th row, 4th column of the 1st layer and the exclusion sphere after being trimmed by this layer is as Figure 8 shown.

[0417] S515. Determine whether \(i_y\) is greater than or equal to \(i1_y\). If less, set \(i_y = i_y + 1 = 4+1 = 5\), and jump to step S57.

[0418] ……

[0419] The schematic diagram of the positional relationship between the grids traversed in the 1st layer and the exclusion sphere after being trimmed by this layer is as Figure 9 shown, where the left figure is a three - dimensional view and the right figure is a top - view.

[0420] ……

[0421] The schematic diagram of the positional relationship between the grids traversed in the 2nd layer and the exclusion sphere after being trimmed by this layer is as Figure 10 As shown, the left image is a 3D view, and the right image is a top view.

[0422] ...

[0423] The diagram showing the positional relationship between the grid in the 6th layer and the repulsion sphere after being cut from this layer is as follows: Figure 11 As shown, the left image is a 3D view, and the right image is a top view.

[0424] ...

[0425] S517, End traversal.

[0426] A diagram showing all the squares intersecting with the repulsive ball that was first thrown is shown below. Figure 12 As shown, these cells are all within the drop zone, and their dropable status has changed to non-droppable. Figure 12 This is a diagram showing the grid cells in the delivery area that are currently in a non-deliverable state.

[0427] S6. Calculate the number of grids in the drop zone that are dropable, n_grid_node = 259.

[0428] S7. Determine if n_grid_node is equal to zero. If not, proceed to step S3.

[0429] In the second run of S3, randomly select a cell from all the cells in the drop zone that are drop-able, such as the cell in the 1st row and 5th column of the 5th layer. Then, randomly select a point within this cell as the drop point, such as the point with coordinates (4.032, 0.056, 4.96). Let i_ball = i_ball + 1 = 2, and drop the second ball. The result after the drop is illustrated in the diagram below. Figure 13 As shown, the left and right images are 3D views from different perspectives.

[0430] In the second run S4, new balls are prepared for deployment. The repulsion sphere formed by the i-th deployed ball and the new ball is calculated. The repulsion sphere of the second deployed ball has a center coordinate of (4.032, 0.056, 4.96) and a radius of 1.6 + 1.6 = 3.2. A schematic diagram of this repulsion sphere is shown below. Figure 14 As shown.

[0431] The second run of S5.

[0432] S51. z0_out = z0 - R_ex = 4.96 - 3.2 = 1.76, z1_out = z0 + R_ex = 4.96 + 3.2 = 8.16, i0_z = floor(z0_out / w0) + 1 = 2, i1_z = ceil(z1_out / w0) = 9, center_z = floor(z0 / w0) + 1 = floor(4.96 / 1) + 1 = 5, center_y = floor(y0 / w0) + 1 = floor(0.056 / 1) + 1 = 1.

[0433] S52. Let i_z = i0_z = 2.

[0434] S53. Judge whether i_z is less than or equal to zero. If not, proceed to the next step.

[0435] S54. Judge whether i_z is greater than m. If not, let i_z_regin = i_z = 2.

[0436] S55. If i_z < center_z, the upper surface of the current grid layer is the reference plane of the current grid layer, z_ref = i_z * w0 = 2 * 1 = 2, dz2 = (z0 - z_ref) 2 = (4.96 - 2) 2 = 8.7616, i0_y = floor((y0 - r_i_z) / w0) + 1 = floor((0.056 - 1.2159) / 1) + 1 = -1, i1_y = ceil((y0 + r_i_z) / w0) = ceil((0.056 + 1.2159) / 1) = 2.

[0437] S56. Let i_y = i0_y = -1.

[0438] S57. Judge whether i_y is less than or equal to zero. If so, let i_y_regin = i_y + m = -1 + 8 = 7, and jump to step S59.

[0439] S59. If i_y < center_y, the back surface of the current grid row is the reference plane of the current grid row, y_ref = i_y * w0 = -1 * 1 = -1, dy2 = (0.056 - (-1)) 2 = 1.1151, i0_x = floor((x0 - dx) / w0) + 1 = floor((4.032 - 0.6028) / 1) + 1 = 4, i1_x = ceil((x0 + dx) / w0) = ceil((4.032 + 0.6028) / 1) = 5.

[0440] S510. Let i_x = i0_x = 4.

[0441] S511. Determine if i_x is less than or equal to zero. If it is not less than or equal to zero, proceed to the next step.

[0442] S512. Determine if i_x is greater than m. If not, let i_x_regin = i_x = 4.

[0443] S513. Set the playable state of the cell in the 7th row and 4th column of the 2nd layer to unplayable. This cell had already been set to unplayable before because it intersects with the repulsion ball of the 1st ball. The positional relationship between the cell in the i_yth row and 4th column of the 2nd layer and the repulsion ball after this layer's trimming is shown in the diagram below. Figure 15 As shown, the left and right images are 3D views from different perspectives.

[0444] S514. Determine whether i_x is greater than or equal to i1_x. If it is less than i_x, let i_x = i_x + 1 = 5, and jump to step S511.

[0445] ...

[0446] S514. Determine if i_x is greater than or equal to i1_x. If it is, proceed to the next step.

[0447] S515. Determine whether i_y is greater than or equal to i1_y. If it is less than i_y, set i_y = i_y + 1 = 0 and jump to step S57.

[0448] S57. Determine if i_y is less than or equal to zero. If it is, let i_y_regin = i_y + m = 0 + 8 = 8, and jump to step S59.

[0449] S59, i0_x = 3, i1_x = 6.

[0450] S510. Let i_x = i0_x = 3.

[0451] S511. Determine if i_x is less than or equal to zero. If it is not less than or equal to zero, proceed to the next step.

[0452] S512. Determine if i_x is greater than m. If not, let i_x_regin = i_x = 3.

[0453] S513. Set the dropable state of the cell in the 8th row and 3rd column of the 2nd layer to non-dropable; the positional relationship between the cell in the i_yth row and 3rd column of the 2nd layer and the repulsion ball after this layer is shown in the diagram. Figure 16 As shown, the left and right images are 3D views from different perspectives.

[0454] After setting the dropable state of the cell in the i_y_regin row and 3rd column of the second layer to non-droppable, the diagram showing the cells in the second layer of the drop area that are now non-droppable is as follows. Figure 17 As shown, the left image is a 3D view, and the right image is a top view.

[0455] ...

[0456] After iterating through each cell that intersects with the repulsive ball of the second ball, updating the throwable state of the corresponding cell within the throwing area to unthrowable, the diagram of the cells within the throwing area that are now unthrowable is shown below. Figure 18 As shown.

[0457] S6. Calculate the number of grids in the drop zone that are dropable, n_grid_node = 143.

[0458] S7. Determine if n_grid_node is equal to zero. If not, proceed to step S3.

[0459] In the third run (S3), randomly select a cell from all the cells in the drop zone that are drop-able. For example, select the cell in the 4th row and 8th column of the 8th layer. Then, randomly select a point within this cell as the drop point, such as the point with coordinates (7.872, 3.992, 7.008). Let i_ball = i_ball + 1 = 3, and drop the 3rd ball. The result after the drop is illustrated in the diagram below. Figure 19 As shown, the left and right images are 3D views from different perspectives.

[0460] ...

[0461] S6. Calculate the number of grids in the drop zone that are dropable, n_grid_node = 31.

[0462] S7. Determine if n_grid_node is equal to zero. If not, proceed to step S3.

[0463] In the fourth run (S3), randomly select a cell from all the cells in the drop zone that are drop-able. For example, select the cell in the 8th row and 5th column of the 1st layer. Then, randomly select a point within this cell as the drop point, such as the point with coordinates (4.352, 7.032, 0.896). Let i_ball = i_ball + 1 = 4, and drop the 4th ball. The result after the drop is illustrated in the diagram below. Figure 20 As shown, the left and right images are 3D views from different perspectives.

[0464] ...

[0465] S6. Calculate the number of grids in the drop zone that are dropable, n_grid_node = 11.

[0466] S7. Determine if n_grid_node is equal to zero. If not, proceed to step S3.

[0467] In the fifth run (S3), randomly select a cell from all the cells in the drop zone that are drop-able, such as the cell in the 8th row and 1st column of the 3rd layer. Then, randomly select a point within this cell as the drop point, such as the point with coordinates (0.16, 7.48, 3.008). Let i_ball = i_ball + 1 = 5, and drop the 5th ball. The result after the drop is illustrated in the diagram below. Figure 21 As shown, the left and right images are 3D views from different perspectives.

[0468] ...

[0469] S6. Calculate the number of grids in the drop zone that are dropable, n_grid_node = 0.

[0470] S7. Determine if n_grid_node is equal to zero. If it is equal to zero, proceed to the next step.

[0471] S8. Let k_ball be the number of balls placed in the current iteration step, and let k_ball = i_ball = 5. Calculate the current volume fraction; the current volume fraction is equal to the sum of the volumes of all placed balls divided by the volume of the placement area.

[0472] S9. Determine if the current volume fraction is greater than or equal to the target volume fraction. If it is less than the target volume fraction, proceed to the next step.

[0473] S10. Determine if itr is greater than or equal to n_itr. If it is less than n_itr, proceed to the next step.

[0474] S11. Let itr = itr + 1 = 1 + 1 = 2.

[0475] S12. Set the occupied status of all cells in the distribution area to unoccupied.

[0476] S13, Let i_ball = 1.

[0477] First run of S14.

[0478] S141. Calculate the lower and upper limit coordinates of the circumscribed cube of the placed ball in the z direction, denoted as z0_out and z1_out respectively. z0_out = z0 - R = 4.096 - 1.6 = 2.496, z1_out = z0 + R = 4.096 + 1.6 = 5.696. Calculate the lower and upper limit grid numbers in the z direction, denoted as i0_z and i1_z respectively. Then i0_z = floor(z0_out / w0) + 1 = 3, i1_z = ceil(z1_out / w0) = 6. And calculate the layer number of the grid layer where the center of the placed ball is located and the row number of the grid row, denoted as center_z and center_y respectively. center_z = floor(z0 / w0) + 1 = floor(4.096 / 1) + 1 = 5, center_y = floor(y0 / w0) + 1 = floor(4.472 / 1) + 1 = 5.

[0479] S142. Let i_z = i0_z = 3.

[0480] S143. Determine whether i_z is less than or equal to zero. If not, proceed to the next step.

[0481] S144. Determine whether i_z is greater than m. If not, let i_z_regin = i_z = 3.

[0482] S145. Determine the reference plane of the current grid layer and calculate the z coordinate z_ref of the reference plane. If i_z < center_z, then the upper surface of the current grid layer is the reference plane of the current grid layer, z_ref = i_z * w0 = 3 * 1 = 3. Calculate the square of the minimum vertical distance from the center of the ball to the reference plane of the current grid layer, denoted as dz2. dz2 = (z0 - z_ref) 2 =(4.096 - 3) 2 =1.201216. Calculate the radius of the circle formed by the intersection of the reference plane of the current grid layer and the placed ball, denoted as r_i_z. Then Calculate the lower and upper limit grid numbers in the y direction of the grids intersecting with the placed ball in the current grid layer, denoted as i0_y and i1_y respectively. Then i0_y = floor((y0 - r_i_z) / w0) + 1 = floor((4.472 - 1.1657) / 1) + 1 = 4, i1_y = ceil((y0 + r_i_z) / w0) = ceil((4.472 + 1.1657) / 1) = 6.

[0483] S146. Let i_y = i0_y = 4.

[0484] S147. Determine whether i_y is less than or equal to zero. If not, proceed to the next step.

[0485] S148. Determine whether \(i_y\) is greater than \(m\). If not, set \(i_y\_regin = i_y = 4\).

[0486] S149. Determine the reference plane of the current grid row and calculate the y - coordinate \(y\_ref\) of the reference plane. If \(i_y < center_y\), then the back surface of the current grid row is the reference plane of the current grid row, and \(y\_ref = i_y * w0 = 4 * 1 = 4\). Calculate the square of the vertical distance from the center of the reference circle to the reference plane of the current grid row, denoted as \(dy2\), \(dy2=(y0 - y\_ref)\) 2 =(4.472 - 4) 2 =0.2227. Calculate half of the length of the chord formed by the intersection of the intersection line of the reference plane of the current grid layer and the reference plane of the current grid row with the reference circle, denoted as \(dx\). Then Calculate the lower and upper grid numbers in the x - direction of the grids intersecting with the dropped ball in the current grid row of the current grid layer, denoted as \(i0_y\) and \(i1_y\) respectively. Then \(i0_x = floor((x0 - dx) / w0)+1 = floor((3.936 - 1.066) / 1)+1 = 3\), \(i1_x = ceil((x0 + dx) / w0)=ceil((3.936 + 1.066) / 1)=6\).

[0487] S1410. Set \(i_x = i0_x = 3\).

[0488] S1411. Determine whether \(i_x\) is less than or equal to zero. If not, proceed to the next step.

[0489] S1412. Determine whether \(i_x\) is greater than \(m\). If not, set \(i_x\_regin = i_x = 3\).

[0490] S1413. Set the occupancy status of the grid at the 3rd layer, 4th row, and 3rd column to occupied; the schematic diagram of the positional relationship between the grid at the 3rd layer, 4th row, and 3rd column and the dropped ball after being cut by this layer is as Figure 22 shown.

[0491] ......

[0492] The schematic diagram of the positional relationship between the grids traversed in the 4th row of the 3rd layer and the dropped ball after being cut by this layer is as Figure 23 shown, where the left figure is a 3 - D view and the right figure is a top view.

[0493] The schematic diagram of the positional relationship between the grids in the 3rd layer traversed and the dropped ball after being cut by this layer is as Figure 24 shown, where the left figure is a 3 - D view and the right figure is a top view.

[0494] ......

[0495] The diagram showing the relationship between the grid of the 4th layer and the position of the ball after being cut from this layer is as follows: Figure 25 As shown, the left image is a 3D view, and the right image is a top view.

[0496] ...

[0497] A diagram showing all the squares intersecting with the first ball is shown below. Figure 26 As shown, the left and right images are 3D views from different perspectives. These grids are all within the projection area, and their occupancy status has changed to occupied. Figure 26 This is a diagram showing the occupied grids within the current deployment area.

[0498] S15. Determine if i_ball is greater than or equal to k_ball. If it is less than k_ball, let i_ball = i_ball + 1 = 2, and jump to step S14.

[0499] In the second run of S14, each cell intersecting with the second ball is traversed, and the occupancy status of the corresponding cell within the throwing area is updated to "occupied". After the update, the diagram showing the occupied cells within the throwing area is as follows. Figure 27 As shown, the left and right images are 3D views from different perspectives.

[0500] S15. Determine if i_ball is greater than or equal to k_ball. If it is less than k_ball, let i_ball = i_ball + 1 = 3, and jump to step S14.

[0501] ...

[0502] S15. Determine if i_ball is greater than or equal to k_ball. If it is less than k_ball, let i_ball = i_ball + 1 = 5, and jump to step S14.

[0503] In the fifth run of S14, each cell intersecting with the fifth ball is traversed, and the occupancy status of the corresponding cell within the throwing area is updated to "occupied". After the update, the diagram showing the occupied cells within the throwing area is as follows. Figure 28 As shown, the left and right images are 3D views from different perspectives.

[0504] S15. Determine if i_ball is greater than or equal to k_ball. If it is, proceed to the next step.

[0505] S16. Let i_ball = 1.

[0506] S17. Prepare a reference ball and calculate the repulsion sphere formed by the i-th ball and the reference ball. The repulsion sphere of the first ball has the center coordinates (3.936, 4.472, 4.096) and the radius 1.6 + 1.6 = 3.2.

[0507] S181. Calculate the lower and upper coordinates of the circumscribed cube of the exclusion sphere in the z direction, denoted as z0_out and z1_out respectively. z0_out = z0 - R_ex = 4.096 - 3.2 = 0.896, z1_out = z0 + R_ex = 4.096 + 3.2 = 7.296. Calculate the lower and upper grid numbers in the z direction, denoted as i0_z and i1_z respectively. Then i0_z = floor(z0_out / w0) + 1 = 1, i1_z = ceil(z1_out / w0) = 8. And calculate the layer number of the grid layer where the center of the exclusion sphere is located and the row number of the grid row, denoted as center_z and center_y respectively. center_z = floor(z0 / w0) + 1 = floor(4.096 / 1) + 1 = 5, center_y = floor(y0 / w0) + 1 = floor(4.472 / 1) + 1 = 5.

[0508] S182. Let num = 0, let num2 = 0, and let i_z = i0_z = 1.

[0509] S183. Judge whether i_z is less than or equal to zero. If not, proceed to the next step.

[0510] S184. Judge whether i_z is greater than m. If not, let i_z_regin = i_z = 1.

[0511] S185. Determine the reference plane of the current grid layer and calculate the z coordinate z_ref of the reference plane. If i_z < center_z, then the upper surface of the current grid layer is the reference plane of the current grid layer, z_ref = i_z * w0 = 1 * 1 = 1. Calculate the square of the minimum vertical distance from the center of the sphere to the reference plane of the current grid layer, denoted as dz2. dz2 = (z0 - z_ref) 2 =(4.096 - 1) 2 =9.585216. Calculate the radius of the circle formed by the intersection of the reference plane of the current grid layer and the exclusion sphere, denoted as r_i_z. Then Calculate the lower and upper grid numbers in the y direction of the grids intersecting with the exclusion sphere in the current grid layer, denoted as i0_y and i1_y respectively. Then i0_y = floor((y0 - r_i_z) / w0) + 1 = floor((4.472 - 0.8092) / 1) + 1 = 4, i1_y = ceil((y0 + r_i_z) / w0) = ceil((4.472 + 0.8092) / 1) = 6.

[0512] S186. Let i_y = i0_y = 4.

[0513] S187. Determine whether \(i_y\) is less than or equal to zero. If not, proceed to the next step.

[0514] S188. Determine whether \(i_y\) is greater than \(m\). If not, set \(i_y\_regin = i_y = 4\).

[0515] S189. Determine the reference plane of the current grid row and calculate the y - coordinate \(y\_ref\) of the reference plane. If \(i_y < center_y\), then the back surface of the current grid row is the reference plane of the current grid row, and \(y\_ref = i_y * w0 = 4 * 1 = 1\). Calculate the square of the perpendicular distance from the center of the reference circle to the reference plane of the current grid row, denoted as \(dy2\), where \(dy2=(y0 - y\_ref)\) 2 =(4.472 - 4) 2 =0.2227. Calculate half of the length of the chord formed by the intersection of the intersection line of the reference plane of the current grid layer and the reference plane of the current grid row with the reference circle, denoted as \(dx\). Then Calculate the lower and upper grid numbers in the x - direction of the grids that intersect the exclusion sphere in the current grid row of the current grid layer, denoted as \(i0_y\) and \(i1_y\) respectively. Then \(i0_x = floor((x0 - dx) / w0)+1 = floor((3.936 - 0.6573) / 1)+1 = 4\), \(i1_x = ceil((x0 + dx) / w0)=ceil((3.936 + 0.6573) / 1)=5\).

[0516] S1810. Set \(i_x = i0_x = 4\).

[0517] S1811. Determine whether \(i_x\) is less than or equal to zero. If not, proceed to the next step.

[0518] S1812. Determine whether \(i_x\) is greater than \(m\). If not, set \(i_x\_regin = i_x = 4\).

[0519] For the first run of S1813: Set \(num2 = num2 + 1 = 0+1 = 1\). The occupancy status of the grid at the 4th row, 4th column, and 1st layer is unoccupied, and \(num\) remains unchanged at 0. The schematic diagram of the positional relationship between the grid at the 4th row, 4th column, and 1st layer and the exclusion sphere after being trimmed by this layer is as shown in the left figure below, and the schematic diagram of the occupancy status of the 1st layer is as shown in the right figure of 29. Figure 29 As shown in the left figure below, and the schematic diagram of the occupancy status of the 1st layer is as shown in the right figure of 29.

[0520] ……

[0521] S1813. Set \(num2 = num2 + 1 = 3+1 = 4\). The occupancy status of the grid at the 5th row, 5th column, and 1st layer is unoccupied, and \(num\) remains unchanged at 0.

[0522] S1814. Determine if i_x is greater than or equal to i1_x. If it is, proceed to the next step.

[0523] S1815. Determine if i_y is greater than or equal to i1_y. If it is less than i1_y, let i_y = i1_y. y =i y +1 = 6, jump to step S187.

[0524] S187. Determine if i_y is less than or equal to zero. If not, proceed to the next step.

[0525] S188. Determine if i_y is greater than m. If not, set i_y_regin = i_y = 6.

[0526] S189, if i_y > center_y, then the front surface of the current cell row is the reference plane of the current cell row, y_ref = (i_y-1)*w0 = (6-1)*1 = 5; dy2 = (y0-y_ref) 2 =(4.472-5) 2 =0.278784, i0_x=floor((x0-dx) / w0)+1=floor((3.936-0.6132) / 1)+1=4, i1_x=ceil((x0+dx) / w0)=ceil((3.936+0.6132) / 1)=5.

[0527] S1810, Let i_x = i0_x = 4.

[0528] S1811. Determine if i_x is less than or equal to zero. If it is not less than or equal to zero, proceed to the next step.

[0529] S1812. Determine if i_x is greater than m. If not, let i_x_regin = i_x = 4.

[0530] S1813. Let num2 = num2 + 1 = 4 + 1 = 5. The cell in the 6th row and 4th column of the 1st layer is occupied, so num = num + 1 = 0 + 1 = 1.

[0531] ...

[0532] In the sixth run of S1813, let num2 = num2 + 1 = 5 + 1 = 6. The cell in the 6th row and 5th column of the first layer is occupied, and num = num + 1 = 1 + 1 = 2. The diagram showing the positional relationship between the cells in the first layer and the repulsion spheres after this layer is shown below. Figure 30 The diagram in the middle left shows the occupancy status of the first layer. Figure 30 As shown in the middle right figure.

[0533] ...

[0534] In the 253rd run S1813, let num2 = num2 + 1 = 252 + 1 = 253. The cell in the 6th row and 6th column of the 8th layer is occupied, num = num + 1 = 132 + 1 = 133. The top view diagram shows the positional relationship between the traversed cells of the 8th layer and the repulsion ball after this layer's trimming. Figure 31 The top view diagram of the occupied state of the 8th layer is shown in the middle left image. Figure 31 As shown in the middle right figure.

[0535] S1814. Determine if i_x is greater than or equal to i1_x. If it is, proceed to the next step.

[0536] S1815. Determine if i_y is greater than or equal to i1_y. If it is, proceed to the next step.

[0537] S1816. Determine if i_z is greater than or equal to i1_z. If it is, proceed to the next step.

[0538] S1817. Calculate the local volume fraction vfl = num / num2 = 133 / 253 = 0.5257, and save it to the local volume fraction array, then end the traversal.

[0539] S19. Determine if i_ball is greater than or equal to k_ball. If it is less than k_ball, let i_ball = i_ball + 1 = 2, and jump to step S17.

[0540] S17. Prepare a reference ball and calculate the repulsion sphere formed by the i-th ball and the reference ball; the repulsion sphere of the second ball has the center coordinates (4.032, 0.056, 4.96) and the radius 1.6 + 1.6 = 3.2.

[0541] S18. Traverse each cell that intersects with the repulsion ball, count the number of cells that intersect with the repulsion ball (num2 = 264), count the number of cells that intersect with the repulsion ball and whose corresponding cells in the placement area are occupied (num = 153), calculate the local volume fraction of the second placement ball (vfl = num / num2 = 0.5795), and save it to the local volume fraction array.

[0542] ...

[0543] S18. Traverse each cell that intersects with the repulsion ball, count the number of cells that intersect with the repulsion ball num2 = 269, count the number of cells that intersect with the repulsion ball and whose corresponding cells in the placement area are occupied num = 91, calculate the local volume fraction vfl = num / num2 = 0.3383 of the third placement ball, and save it to the local volume fraction array.

[0544] ...

[0545] S18. Traverse each cell that intersects with the repulsion ball, count the number of cells that intersect with the repulsion ball num2 = 258, count the number of cells that intersect with the repulsion ball and whose corresponding cells in the placement area are occupied num = 137, calculate the local volume fraction vfl = num / num2 = 0.5310 of the fourth placement ball, and save it to the local volume fraction array.

[0546] ...

[0547] S18. Traverse each cell that intersects with the repulsion ball, count the number of cells that intersect with the repulsion ball (num2 = 250), count the number of cells that intersect with the repulsion ball and whose corresponding cells in the placement area are occupied (num = 116), calculate the local volume fraction of the 5th placement ball (vfl = num / num2 = 0.4640), and save it to the local volume fraction array.

[0548] S19. Determine if i_ball is greater than or equal to k_ball. If it is, proceed to the next step.

[0549] S20. Based on the local volume fractions of all k_ball balls, calculate the maximum value vfl_max = 0.5795 and the minimum value vfl_min = 0.3383, and calculate the baseline local volume fraction vfl_base = (1-alpha)*vfl_min + alpha*vfl_max = (1-0.41)*vfl_min + 0.41*vfl_max = 0.4372.

[0550] S21. Set the retention status of all thrown balls to non-retention.

[0551] S22, Let i_ball = 1.

[0552] S23. Determine whether the local volume fraction of the first ball is greater than or equal to vfl_base. If it is, proceed to the next step.

[0553] S24. Mark the state of the first ball as "reserved".

[0554] S25. Determine if i_ball is greater than or equal to k_ball. If it is less than k_ball, let i_ball = i_ball + 1 = 2, and jump to step S23.

[0555] S23. Determine whether the local volume fraction of the second ball is greater than or equal to vfl_base. If it is, proceed to the next step.

[0556] S24. Mark the status of the second ball as "reserved".

[0557] S25. Determine if i_ball is greater than or equal to k_ball. If it is less than k_ball, let i_ball = i_ball + 1 = 3, and jump to step S23.

[0558] S23. Determine whether the local volume fraction of the third ball is greater than or equal to vfl_base. If it is less than vfl_base, proceed to step S25.

[0559] S25. Determine if i_ball is greater than or equal to k_ball. If it is less than k_ball, let i_ball = i_ball + 1 = 4, and jump to step S23.

[0560] S23. Determine whether the local volume fraction of the fourth ball is greater than or equal to vfl_base. If it is, proceed to the next step.

[0561] S24. Mark the status of the 4th ball as "reserved".

[0562] S25. Determine if i_ball is greater than or equal to k_ball. If it is less than k_ball, let i_ball = i_ball + 1 = 5, and jump to step S23.

[0563] S23. Determine whether the local volume fraction of the 5th ball is greater than or equal to vfl_base. If it is, proceed to the next step.

[0564] S24. Mark the 5th ball's reserved status as reserved.

[0565] S25. Determine if i_ball is greater than or equal to k_ball. If it is, proceed to the next step.

[0566] S26. Delete the balls whose retention status is "not retained" from the ball-dropping set, update the ball indices in the new set, and update the value of k_ball to the number of retained balls. The ball-dropping set now contains 4 balls. After updating, k_ball = 4. After updating the ball indices in the new set, the center coordinates of ball 1 are (3.936, 4.472, 4.096), ball 2 is (4.032, 0.056, 4.96), ball 3 is (4.352, 7.032, 0.896), and ball 4 is (0.16, 7.48, 3.008), with a radius of 1.6. The positions of the 4 balls are shown in the diagram below. Figure 32 As shown, the left and right images are 3D views from different perspectives.

[0567] S27. Set the drop status of all cells in the drop area to dropable.

[0568] S28. Let i_ball = 1, and prepare to throw a new ball.

[0569] S29. Calculate the repulsion sphere formed by the first ball and the newly placed ball; the repulsion sphere of the first ball has the center coordinates (3.936, 4.472, 4.096) and the radius 1.6 + 1.6 = 3.2.

[0570] The first run of S30 iterates through each cell intersecting with the repulsion ball, updating the throwable state of the corresponding cell within the throwing area to unthrowable. After the update, a diagram showing the cells within the throwing area that were previously throwable but are now unthrowable is shown below. Figure 33 As shown.

[0571] ...

[0572] S31. Determine if i_ball is greater than or equal to k_ball. If it is less than k_ball, let i_ball = i_ball + 1 = 4, and jump to step S29.

[0573] S29. Calculate the repulsion sphere formed by the fourth ball and the newly placed ball; the repulsion sphere of the fourth ball has the coordinates of (0.16, 7.48, 3.008) and the radius of 1.6 + 1.6 = 3.2.

[0574] In the fourth run of S30, each cell intersecting with the repulsion ball is traversed, and the throwable state of the corresponding cell within the throwing area is updated to unthrowable. After the update, the diagram showing the cells within the throwing area that were previously throwable but are now unthrowable is as follows: Figure 34 As shown, the left and right images are 3D views from different perspectives.

[0575] S31. Determine if i_ball is greater than or equal to k_ball. If it is, proceed to step S6.

[0576] S6. Calculate the number of grids in the drop zone that are dropable, n_grid_node = 40.

[0577] S7. Determine if n_grid_node is equal to zero. If not, proceed to step S3.

[0578] In the sixth run (S3), randomly select a cell from all the cells in the drop zone that are drop-able, such as the cell in the 5th row and 1st column of the 7th layer. Then, randomly select a point within this cell as the drop point, such as the point with coordinates (0.704, 4.28, 6.88). Let i_ball = i_ball + 1 = 5, and drop the 5th ball. The result after the drop is illustrated in the diagram below. Figure 35 As shown, the left and right images are 3D views from different perspectives.

[0579] ...

[0580] S6. Calculate the number of grids in the drop zone that are dropable, n_grid_node = 1.

[0581] S7. Determine if n_grid_node is equal to zero. If not, proceed to step S3.

[0582] In the seventh run (S3), randomly select a cell from all the cells in the drop zone that are drop-able. Only the cell in the 4th row and 8th column of the 3rd layer can be selected. Then, randomly select a point within this cell as the drop point. For example, if the point with coordinates (7.264, 3.896, 2.4) is selected, let i_ball = i_ball + 1 = 6, and drop the 6th ball. The result after the drop is illustrated in the diagram below. Figure 36 As shown, the left and right images are 3D views from different perspectives.

[0583] ...

[0584] S6. Calculate the number of grids in the drop zone that are dropable, n_grid_node = 0.

[0585] S7. Determine if n_grid_node is equal to zero. If it is equal to zero, proceed to the next step.

[0586] S8. Let k_ball be the number of balls that have been placed in the current iteration step. Let k_ball = i_ball = 6, and calculate the current volume fraction. The current volume fraction is equal to the sum of the volumes of all placed balls divided by the volume of the placement area, i.e.

[0587] S9. Determine whether the current volume fraction is greater than or equal to the target volume fraction. If it is, proceed to step S32.

[0588] S32, Let i_ball = 1.

[0589] First run of S33.

[0590] S331. Let the x, y, and z coordinates of the center of the i-th ball be x0, y0, and z0, respectively. Let xyz_min equal the ball radius and xyz_max equal the area side length minus the ball radius. Then x0 = 3.936, y0 = 4.472, z0 = 4.096, xyz_min = 1.6, and xyz_max = 8 - 1.6 = 6.4. Execute S332, S334, and S336 in parallel.

[0591] S332, z0 is greater than or equal to xyz_min and less than or equal to xyz_max, set z3 to -1, and jump to step S338.

[0592] S334, y0 is greater than or equal to xyz_min and less than or equal to xyz_max, set y3 to -1, and jump to step S338.

[0593] S336, x0 is greater than or equal to xyz_min and less than or equal to xyz_max, set x3 to -1, and jump to step S338.

[0594] S338, Wait for S332-S337 to complete.

[0595] S339, execute S3310, S3314, S3318 and S3322 in parallel.

[0596] S3310, x3 is less than zero, jump to step S3326.

[0597] If S3314 and y3 are less than zero, proceed to step S3326.

[0598] If S3318 and z3 are less than zero, proceed to step S3326.

[0599] S3322, x3 is less than zero, jump to step S3326.

[0600] S3326, Wait for S3310-S3325 to complete.

[0601] S34. Determine if i_ball is greater than or equal to k_ball. If it is less than k_ball, let i_ball = i_ball + 1 = 2, and jump to step S33.

[0602] The second run of S33.

[0603] S331. Let the x, y, and z coordinates of the center of the i-th ball be x0, y0, and z0, respectively. Let xyz_min equal the ball radius and xyz_max equal the area side length minus the ball radius. Then x0 = 4.032, y0 = 0.056, z0 = 4.96, xyz_min = 1.6, and xyz_max = 8 - 1.6 = 6.4. Execute S332, S334, and S336 in parallel.

[0604] S332, z0 is greater than or equal to xyz_min and less than or equal to xyz_max, set z3 to -1, and jump to step S338.

[0605] S334, y0 is less than xyz_min, let y1 be equal to y0 plus the side length of the region, then y1 = 0.056 + 8 = 8.056, let y3 be equal to zero, and proceed to the next step.

[0606] S335. A new supplementary ball is generated with center coordinates (4.032, 8.056, 4.96) and radius equal to the radius of the thrown ball (1.6). Record the center coordinates and radius of this supplementary ball, add it to the supplementary ball set, and jump to step S338. A schematic diagram of the second thrown ball and its supplementary ball is shown below. Figure 37 As shown, the left and right images are 3D views from different perspectives. The ball with the coarse grid is the thrown ball, and the ball with the fine grid is the supplementary ball.

[0607] S336, x0 is greater than or equal to xyz_min and less than or equal to xyz_max, set x3 to -1, and jump to step S338.

[0608] S338, Wait for S332-S337 to complete.

[0609] S339, execute S3310, S3314, S3318 and S3322 in parallel.

[0610] S3310, x3 is less than zero, jump to step S3326.

[0611] S3314, y3 equals zero and z3 is less than zero, jump to step S3326.

[0612] If S3318 and z3 are less than zero, proceed to step S3326.

[0613] S3322, x3 is less than zero, jump to step S3326.

[0614] S3326, Wait for S3310-S3325 to complete.

[0615] S34. Determine if i_ball is greater than or equal to k_ball. If it is less than k_ball, let i_ball = i_ball + 2 + 1 = 3, and jump to step S33.

[0616] The third run of S33.

[0617] S331. Let the x, y, and z coordinates of the center of the i-th ball be x0, y0, and z0, respectively. Let xyz_min equal the ball radius and xyz_max equal the area side length minus the ball radius. Then x0 = 4.352, y0 = 7.032, z0 = 0.896, xyz_min = 1.6, and xyz_max = 8 - 1.6 = 6.4. Execute S332, S334, and S336 in parallel.

[0618] S332, z0 is less than xyz_min, let z1 be equal to z0 plus the side length of the region, then z1 = 0.896 + 8 = 8.896, let z3 be equal to zero, and proceed to the next step.

[0619] S333. A new supplementary ball is generated with center coordinates (4.352, 7.032, 8.896) and radius equal to the radius of the thrown ball (1.6). Record the center coordinates and radius of this supplementary ball and add it to the supplementary ball set. Jump to step S338. A schematic diagram of the third thrown ball and its supplementary ball is shown below. Figure 38 As shown, the balls with coarse grids are the balls that are thrown, and the balls with fine grids are the balls that are replenished.

[0620] S334, y0 is greater than xyz_max, let y1 be equal to y0 minus the side length of the region, then y1 = 7.032 - 8 = -0.968, let y3 be equal to the side length of the region, that is, y3 = 8, and proceed to the next step.

[0621] S335. A new supplementary ball is generated with center coordinates (4.352, -0.968, 0.896) and radius equal to the radius of the thrown ball (1.6). Record the center coordinates and radius of this supplementary ball and add it to the supplementary ball set. Jump to step S338. A schematic diagram of the third thrown ball and its supplementary ball is shown below. Figure 39 As shown, the balls with coarse grids are the balls that are thrown, and the balls with fine grids are the balls that are replenished.

[0622] S336, x0 is greater than or equal to xyz_min and less than or equal to xyz_max, set x3 to -1, and jump to step S338.

[0623] S338, Wait for S332-S337 to complete.

[0624] S339, execute S3310, S3314, S3318 and S3322 in parallel.

[0625] S3310, x3 is less than zero, jump to step S3326.

[0626] S3314, y3 is greater than zero and z3 is equal to zero, proceed to the next step; otherwise, jump to step S3326.

[0627] S3315. Using y3 and z3 as reference points, calculate the distance d30 from the reference point to (y0, z0); that is, calculate the distance d30 from (8, 0) to (7.032, 0.896), d30 = 1.3190.

[0628] S3316 and d30 are less than the radius of the ball being thrown, proceed to the next step.

[0629] S3317. A new supplementary ball is generated with center coordinates (4.352, -0.968, 8.896) and radius equal to the radius of the thrown ball. Record the center coordinates and radius of this supplementary ball and add it to the supplementary ball set. Jump to step S3326. A schematic diagram of the third thrown ball and its supplementary ball is shown below. Figure 40As shown, the balls with coarse grids are the balls that are thrown, and the balls with fine grids are the balls that are replenished.

[0630] S3318, z3 equals zero and x3 is less than zero, jump to step S3326.

[0631] S3322, x3 is less than zero, jump to step S3326.

[0632] S3326, Waiting for S3310-S3325 to complete, the diagram of the third ball and its three replacement balls is as follows. Figure 41 As shown, the left and right images are 3D views from different perspectives. The ball with the coarse grid is the thrown ball, and the ball with the fine grid is the supplementary ball.

[0633] ...

[0634] In the fourth run of S33, it determines whether the fourth ball intersects with the boundary of the throwing area. If they intersect, a replacement ball is generated and added to the replacement ball set. A schematic diagram of the fourth ball and its three replacement balls is shown below. Figure 42 As shown, the left and right images are 3D views from different perspectives. The ball with the coarse grid is the thrown ball, and the ball with the fine grid is the supplementary ball.

[0635] ...

[0636] In the fifth run of S33, it determines whether the fifth ball intersects with the boundary of the throwing area. If they intersect, a replacement ball is generated and added to the replacement ball set. A schematic diagram of the fifth ball and its three replacement balls is shown below. Figure 43 As shown, the left and right images are 3D views from different perspectives. The ball with the coarse grid is the thrown ball, and the ball with the fine grid is the supplementary ball.

[0637] ...

[0638] In the sixth run of S33, it determines whether the sixth ball intersects with the boundary of the throwing area. If they intersect, a replacement ball is generated and added to the replacement ball set. A schematic diagram of the sixth ball and its one replacement ball is shown below. Figure 44 As shown, the left and right images are 3D views from different perspectives. The ball with the coarse grid is the thrown ball, and the ball with the fine grid is the supplementary ball.

[0639] S34. Determine if i_ball is greater than or equal to k_ball. If it is, proceed to the next step. A diagram showing all 11 replacement balls is shown below. Figure 45 As shown, the left and right images are 3D views from different perspectives.

[0640] S35. Save and output the set of thrown balls and the set of replenished balls, output a successful throwing flag, and jump to step S37. The set of thrown balls contains 6 balls with center coordinates of (3.936, 4.472, 4.096), (4.032, 0.056, 4.96), (4.352, 7.032, 0.896), (0.16, 7.48, 3.008), (0.704, 4.28, 6.88), and (7.264, 3.896, 2.4). The set of replenished balls contains 11 balls with center coordinates of (4.032, 8.056, 4.096). 96), (4.352,7.032,8.896), (4.352,-0.968,0.896), (4.352,-0.968,8.896), (8.16,7.48,3.008), (0.16,-0.52,3.008), (8.16,-0.52,3.008), (8.704,4.28,6.88), (0.704,4.28,-1.12), (8.704,4.28,-1.12), (-0.736,3.896,2.4), a schematic diagram of all balls is shown below. Figure 46 As shown, the left and right images are 3D views from different perspectives. The spheres with coarse grids are the thrown spheres, and the spheres with fine grids are the supplementary spheres. All spheres that extend beyond the throwing area are hidden.

[0641] S37, End the deployment.

[0642] The results obtained in this experimental example are only under the condition of a relatively small target volume fraction. They are intended to demonstrate that our method has the ability to increase the volume fraction with increasing iteration steps. For clarity, the grid division in this experimental example is relatively sparse. Those skilled in the art will understand that in some experimental examples, to obtain a larger volume fraction, the grid division can be denser and the upper limit of the total number of iterations can be set larger; as the number of iterations increases, the maximum volume fraction can ultimately be obtained.

[0643] The above embodiments are not intended to limit the shape, material, structure, etc. of the present invention in any way. Any simple modifications, equivalent changes, and alterations made to the above embodiments based on the technical essence of the present invention shall fall within the protection scope of the present invention.

[0644] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. However, these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for generating high volume fraction particles based on grid points / grids through random distribution, characterized in that, First, set the input parameters: the radius of the ball particle R, the ratio of the side length of the placement area to the radius of the ball particle δ, the target volume fraction, and the upper limit of the total number of iterations n_itr. Then the length, width, and height of the placement area are all width = R * δ. Assume that the x-axis of the spatial rectangular coordinate system is horizontal to the right, the y-axis is horizontal to the back, and the z-axis is vertical to the top. The sides of the projection area are parallel to the coordinate axes, and its lower left front corner is located at the origin. The generation method includes the following steps: S1. The delivery area is discrete. A sufficiently large three-dimensional background grid is established to cover the delivery area. The delivery area is evenly divided into m×m×m grids. Then the grid width w0 = width / m. The delivery state is established for all grid points / grids in the delivery area and all are set to delivery. The upper right corner of each grid is taken as the grid point. The grid points corresponding to all grid points / grids in the delivery area are called grid points / grid points in the delivery area. S2. Let itr = 1, and let i_ball = 0; S3. When dropping based on grid points, randomly select a grid point from all the grid points in the dropping area that are dropable as the dropping point; when dropping based on grid, randomly select a grid from all the grid points in the dropping area that are dropable, and randomly select a point in this grid as the dropping point. Let i_ball = i_ball + 1, drop the i_ball ball so that the center of this dropped ball is located at this dropping point, record the coordinates of the center of this dropped ball and its radius, and add it to the dropping ball set; S4. Prepare to release new balls, and calculate the repulsion sphere formed by the i_ball-th released ball and the newly released ball. S5. When deploying based on grid points, traverse each grid point located inside the repulsion sphere and update the deployable status of the corresponding grid point in the deployment area to undeployable; when deploying based on grids, traverse each grid that intersects with the repulsion sphere and update the deployable status of the corresponding grid point in the deployment area to undeployable. S6. Calculate the number of grid points / grids in the delivery area that are in a deliverable state, denoted as n_grid_node; S7. Determine if n_grid_node is equal to zero. If it is equal to zero, proceed to the next step; otherwise, jump to step S3. S8. Let k_ball be the number of balls that have been placed in the current iteration step, let k_ball = i_ball, and calculate the current volume fraction; S9. Determine whether the current volume fraction is greater than or equal to the target volume fraction. If it is, proceed to step S32; if it is less than, proceed to the next step. S10. Determine if itr is greater than or equal to n_itr. If it is, jump to step S36; if it is less than, proceed to the next step. S11. Let itr = itr + 1; S12. Set the occupation status of all grid points / grids within the deployment area to unoccupied; S13. Let i_ball = 1; S14. When dropping based on grid points, traverse each grid point inside the i_ball dropping ball and update the occupation status of the corresponding grid point in the dropping area to occupied; when dropping based on grid, traverse each grid that intersects with the i_ball dropping ball and update the occupation status of the corresponding grid point in the dropping area to occupied. S15. Determine if i_ball is greater than or equal to k_ball. If it is, proceed to the next step; if it is less than, set i_ball = i_ball + 1 and jump to step S14. S16. Let i_ball = 1; S17. Prepare a reference ball and calculate the repulsion sphere formed by the i-th ball and the reference ball. S18. When deploying based on grid points, traverse every grid point inside the repulsion sphere, count the number of grid points inside the repulsion sphere (num2), count the number of grid points inside the repulsion sphere whose corresponding grid point occupies the deployment area (num), calculate the local volume fraction vfl = num / num2 for the i-th deployed ball, and save it to the local volume fraction array; when deploying based on grid, traverse every grid that intersects with the repulsion sphere, count the number of grid points that intersect with the repulsion sphere (num2), count the number of grid points that intersect with the repulsion sphere whose corresponding grid point occupies the deployment area (num), calculate the local volume fraction vfl = num / num2 for the i-th deployed ball, and save it to the local volume fraction array. S19. Determine if i_ball is greater than or equal to k_ball. If it is, proceed to the next step; if it is less than, let i_ball = i_ball + 1 and jump to step S17. S20. Based on the local volume fractions of all k_ball balls, calculate the maximum value vfl_max and the minimum value vfl_min, and calculate the baseline local volume fraction vfl_base = (1-alpha)*vfl_min + alpha*vfl_max; S21. Set the retention status of all thrown balls to non-retention. S22. Let i_ball = 1; S23. Determine whether the local volume fraction of the i-th ball is greater than or equal to vfl_base. If it is, proceed to the next step; if it is less than, jump to step S25. S24. Mark the retention status of the i_ball-th ball as retained; S25. Determine if i_ball is greater than or equal to k_ball. If it is, proceed to the next step; if it is less than, let i_ball = i_ball + 1 and jump to step S23. S26. Delete the balls whose retention status is not retained from the ball set, update the ball index in the new set, and update the value of k_ball to the number of retained balls; S27. Set the drop status of all grid points / grids in the drop area to dropable; S28. Let i_ball = 1, and prepare to throw a new ball; S29. Calculate the repulsion ball formed by the i-th ball and the newly placed ball; S30. When deploying based on grid points, traverse each grid point located inside the repulsion sphere and update the deployable state of the corresponding grid point in the deployment area to be undeployable; when deploying based on grids, traverse each grid that intersects with the repulsion sphere and update the deployable state of the corresponding grid point in the deployment area to be undeployable. S31. Determine if i_ball is greater than or equal to k_ball. If it is, jump to step S6; if it is less than, let i_ball = i_ball + 1, and jump to step S29. S32. Let i_ball = 1; S33. Determine whether the i_th ball intersects with the boundary of the throwing area. If they intersect, generate a supplementary ball for this ball and add it to the supplementary ball set; otherwise, jump to step S34. S34. Determine if i_ball is greater than or equal to k_ball. If it is, proceed to the next step; if it is less than, let i_ball = i_ball + 1 and jump to step S33. S35. Save and output the set of thrown balls and the set of replenished balls, output a successful throwing flag, and jump to step S37. S36, Output delivery failure flag; S37, End the deployment.

2. The method for generating high volume fraction particles based on grid points / grids according to claim 1, characterized in that, In steps S4 and S29, the center coordinates of the repulsion ball are equal to the center coordinates of the i-th ball, and the radius of the repulsion ball is equal to the radius of the i-th ball plus the radius of the newly placed ball.

3. The method for generating high volume fraction particles based on grid points / grids according to claim 1, characterized in that, The traversal and update processes in steps S5 and S30 are the same, both including the following steps: S51. Calculate the lower and upper bound coordinates of the cube circumscribed by the repulsion sphere in the z-direction, denoted as z0_out and z1_out respectively. Then z0_out = z0 - R_ex and z1_out = z0 + R_ex. When deploying based on grid points, calculate the lower and upper bound grid numbers in the z-direction, denoted as i0_z and i1_z respectively. Then i0_z = floor(z0_out / w0) + 1 and i1_z = ceil(z1_out / w0) - 1. When deploying based on grid points, calculate the lower and upper bound grid numbers in the z-direction, denoted as i0_z and i1_z respectively. Then i0_z = floor(z0_out / w0) + 1, i1_z = ceil(z1_out / w0), and calculate the layer number and row number of the grid where the center of the repulsion sphere is located, denoted as center_z and center_y respectively, then center_z = floor(z0 / w0) + 1, center_y = floor(y0 / w0) + 1; floor() means round down, ceil() means round up, y0 and z0 are the y and z coordinates of the center of the repulsion sphere, and R_ex is the radius of the repulsion sphere; S52. Let i_z = i0_z; S53. Determine if i_z is less than or equal to zero. If it is, let i_z_regin = i_z + m, and jump to step S55. If it is not less than or equal to zero, proceed to the next step. S54. Determine if i_z is greater than m. If it is, let i_z_regin = i_z - m; if it is not, let i_z_regin = i_z. When placing based on grid points, calculate the square of the vertical distance from the center of the sphere to the current grid point layer, denoted as dz2, then dz2 = (z0 - i_z * w0) 2 , calculate the radius of the reference circle formed by the intersection of the current grid point layer and the exclusion sphere, denoted as r_i_z, then Calculate the lower and upper grid point numbers in the y direction of the grid points located inside the exclusion sphere in the current grid point layer, denoted as i0_y and i1_y respectively. Then i0_y = floor((y0 - r_i_z) / w0) + 1, i1_y = ceil((y0 + r_(i_z)) / w0) - 1; When placing based on the grid, determine the reference plane of the current grid layer and calculate the z coordinate z_ref of the reference plane. If i_z < center_z, the upper surface of the current grid layer is the reference plane of the current grid layer, z_ref = i_z * w0. If i_z = center_z, the horizontal plane where the center of the sphere is located is the reference plane of the current grid layer, z_ref = 0. If i_z > center_z, the lower surface of the current grid layer is the reference plane of the current grid layer, z_ref = (i_z - 1) * w0. Calculate the square of the minimum vertical distance from the center of the sphere to the reference plane of the current grid layer, denoted as dz2, dz2 = (z0 - z_ref) 2 , calculate the radius of the reference circle formed by the intersection of the reference plane of the current grid layer and the exclusion sphere, denoted as r_i_z, then Calculate the lower and upper grid numbers in the y direction of the grids that intersect with the exclusion sphere in the current grid layer, denoted as i0_y and i1_y respectively. Then i0_y = floor((y0 - r_i_z) / w0) + 1, i1_y = ceil((y0 + r_i_z) / w0); floor() represents rounding down, ceil() represents rounding up, y0 and z0 are the y and z coordinates of the center of the exclusion sphere, and R_ex is the radius of the exclusion sphere. S56. Let i_y = i0_y; S57. Determine if i_y is less than or equal to zero. If it is, let i_y_regin = i_y + m, and jump to step S59. If it is not less than or equal to zero, proceed to the next step. S58. Determine if i_y is greater than m. If it is, let i_y_regin = i_y - m; if it is not, let i_y_regin = i_y. When placing based on grid points, calculate the square of the vertical distance from the center of the reference circle to the current grid point row, denoted as dy2, then dy2 = (y0 - i_y * w0) 2 , calculate half of the length of the chord formed by the intersection of the straight line where the current grid point row is located and the reference circle, denoted as dx, then Calculate the lower and upper grid point numbers in the x direction of the grid points located inside the exclusion sphere in the current grid point row of the current grid point layer, denoted as i0_x and i1_x respectively, then i0_x = floor((x0 - dx) / w0) + 1, i1_x = ceil((x0 + dx) / w0) - 1; When placing based on the grid, determine the reference plane of the current grid row and calculate the y coordinate y_ref of the reference plane. If i_y < center_y, the back surface of the current grid row is the reference plane of the current grid row, y_ref = i_y * w0. If i_y = center_y, the vertical plane passing through the center of the reference circle and parallel to the front and back surfaces of the grid is the reference plane of the current grid row, y_ref = 0. If i_y > center_y, the front surface of the current grid row is the reference plane of the current grid row, y_ref = (i_y - 1) * w0. Calculate the square of the vertical distance from the center of the reference circle to the reference plane of the current grid row, denoted as dy2, dy2 = (y0 - y_ref) 2 , calculate half of the length of the chord formed by the intersection of the intersection line of the reference plane of the current grid layer and the reference plane of the current grid row and the reference circle, denoted as dx, then Calculate the lower and upper grid numbers in the x direction of the grids intersecting with the exclusion sphere in the current grid row of the current grid layer, denoted as i0_y and i1_y respectively, then i0_x = floor((x0 - dx) / w0) + 1, i1_x = ceil((x0 + dx) / w0); floor() represents rounding down, ceil() represents rounding up, y0 and x0 are the y and x coordinates of the center of the exclusion sphere S510. Let i_x = i0_x; S511. Determine if i_x is less than or equal to zero. If it is, set i_x_regin = i_x + m and jump to step S513. If it is not less than or equal to zero, proceed to the next step. S512. Determine whether i_x is greater than m. If it is, let i_x_regin = i_x - m; if it is not, let i_x_regin = i_x. S513. Set the projectable state of the grid point / grid in the i_y_regin row and i_x_regin column of the i_z_regin layer to unprojectable; S514. Determine whether i_x is greater than or equal to i1_x. If it is, proceed to the next step; if it is less than, let i_x = i_x + 1 and jump to step S511. S515. Determine whether i_y is greater than or equal to i1_y. If it is, proceed to the next step; if it is less than, let i_y = i_y + 1 and jump to step S57. S516. Determine whether i_z is greater than or equal to i1_z. If it is, proceed to the next step; if it is less than, let i_z = i_z + 1 and jump to step S53. S517, End traversal.

4. The method for generating high volume fraction particles based on grid points / grids according to claim 1, characterized in that, In step S8, the current volume fraction is equal to the sum of the volumes of all the balls placed divided by the volume of the placement area.

5. The method for generating high volume fraction particles based on grid points / grids according to claim 1, characterized in that, In step S14, the traversal update process includes the following steps: S141. Calculate the lower and upper bound coordinates of the circumscribed cube of the deployed sphere in the z-direction, denoted as z0_out and z1_out respectively. Then z0_out = z0 - R and z1_out = z0 + R. When deploying based on grid points, calculate the lower and upper bound grid numbers in the z-direction, denoted as i0_z and i1_z respectively. Then i0_z = floor(z0_out / w0) + 1 and i1_z = ceil(z1_out / w0) - 1. When deploying based on grid points, calculate the lower and upper bound grid numbers in the z-direction, denoted as i0_z and i1_z respectively. Then i0_z = floor(z0_out / w0) + 1, i1_z = ceil(z1_out / w0), and calculate the layer number and row number of the grid where the center of the ball is located, denoted as center_z and center_y respectively, then center_z = floor(z0 / w0) + 1, center_y = floor(y0 / w0) + 1; floor() means round down, ceil() means round up, y0 and z0 are the y and z coordinates of the center of the ball, and R is the radius of the ball. S142. Let i_z = i0_z; S143. Determine if i_z is less than or equal to zero. If it is, let i_z_regin = i_z + m, and jump to step S145. If it is not less than or equal to zero, proceed to the next step. S144. Determine if i_z is greater than m. If it is, let i_z_regin = i_z - m; if it is not, let i_z_regin = i_z. When placing based on grid points, calculate the square of the vertical distance from the center of the sphere to the current grid point layer, denoted as dz2, then dz2 = (z0 - i_z * w0) 2 , calculate the radius of the reference circle formed by the intersection of the current grid point layer and the placed sphere, denoted as r_i_z, then Calculate the lower and upper grid point numbers in the y direction of the grid points located inside the placed sphere in the current grid point layer, denoted as i0_y and i1_y respectively. Then i0_y = floor((y0 - r_i_z) / w0) + 1, i1_y = ceil((y0 + r_i_z) / w0) - 1; When placing based on the grid, determine the reference plane of the current grid layer and calculate the z coordinate z_ref of the reference plane. If i_z < center_z, the upper surface of the current grid layer is the reference plane of the current grid layer, z_ref = i_z * w0. If i_z = center_z, the horizontal plane where the center of the sphere is located is the reference plane of the current grid layer, z_ref = 0. If i_z > center_z, the lower surface of the current grid layer is the reference plane of the current grid layer, z_ref = (i_z - 1) * w0. Calculate the square of the minimum vertical distance from the center of the sphere to the reference plane of the current grid layer, denoted as dz2, dz2 = (z0 - z_ref) 2 , calculate the radius of the reference circle formed by the intersection of the reference plane of the current grid layer and the placed sphere, denoted as r_i_z, then Calculate the lower and upper grid numbers in the y direction of the grids intersecting with the placed sphere in the current grid layer, denoted as i0_y and i1_y respectively. Then i0_y = floor((y0 - r_i_z) / w0) + 1, i1_y = ceil((y0 + r_i_z) / w0); floor() represents rounding down, ceil() represents rounding up. y0 and z0 are the y and z coordinates of the center of the placed sphere, and R is the radius of the placed sphere S146. Let i_y = i0_y; S147. Determine if i_y is less than or equal to zero. If it is less than or equal to zero, set i_y_regin = i_y + m and jump to step S149. If it is not less than or equal to zero, proceed to the next step. S148. Determine if i_y is greater than m. If it is, let i_y_regin = i_y - m; if it is not, let i_y_regin = i_y. When placing based on grid points, calculate the square of the vertical distance from the center of the reference circle to the current grid point row, denoted as dy2, then dy2 = (y0 - i_y * w0) 2 , calculate half of the length of the chord formed by the intersection of the line where the current grid point row is located and the reference circle, denoted as dx, then Calculate the lower and upper grid point numbers in the x direction of the grid points located inside the placement ball in the current grid point row of the current grid point layer, denoted as i0_x and i1_x respectively. Then i0_x = floor((x0 - dx) / w0) + 1, i1_x = ceil((x0 + dx) / w0) - 1; When placing based on the grid, determine the reference plane of the current grid row and calculate the y coordinate y_ref of the reference plane. If i_y < center_y, the back surface of the current grid row is the reference plane of the current grid row, y_ref = i_y * w0. If i_y = center_y, the vertical plane where the center of the reference circle is located and parallel to the front and back surfaces of the grid is the reference plane of the current grid row, y_ref = 0. If i_y > center_y, the front surface of the current grid row is the reference plane of the current grid row, y_ref = (i_y - 1) * w0. Calculate the square of the vertical distance from the center of the reference circle to the reference plane of the current grid row, denoted as dy2, dy2 = (y0 - y_ref) 2 , calculate half of the length of the chord formed by the intersection of the intersection line of the reference plane of the current grid layer and the reference plane of the current grid row and the reference circle, denoted as dx, then Calculate the lower and upper grid numbers in the x direction of the grids intersecting with the placement ball in the current grid row of the current grid layer, denoted as i0_y and i1_y respectively. Then i0_x = floor((x0 - dx) / w0) + 1, i1_x = ceil((x0 + dx) / w0); floor() represents rounding down, ceil() represents rounding up, y0 and x0 are the y and x coordinates of the center of the placement ball S1410. Let i_x = i0_x; S1411. Determine if i_x is less than or equal to zero. If it is, set i_x_regin = i_x + m and jump to step S1413. If it is not less than or equal to zero, proceed to the next step. S1412. Determine if i_x is greater than m. If it is, let i_x_regin = i_x - m; if it is not, let i_x_regin = i_x. S1413. Set the occupied status of the cell point / cell in the i_z_regin row and i_x_regin column of the i_y_regin layer to occupied; S1414. Determine if i_x is greater than or equal to i1_x. If it is, proceed to the next step; if it is less than, let i_x = i_x + 1 and jump to step S1411. S1415. Determine if i_y is greater than or equal to i1_y. If it is, proceed to the next step; if it is less than, let i_y = i_y + 1 and jump to step S147. S1416. Determine if i_z is greater than or equal to i1_z. If it is, proceed to the next step; if it is less than, let i_z = i_z + 1 and jump to step S143. S1417, End traversal.

6. The method for generating high volume fraction particles based on grid points / grids according to claim 1, characterized in that, In step S17, the center coordinates of the repulsion ball are equal to the center coordinates of the i-th ball, and the radius of the repulsion ball is equal to the radius of the i-th ball plus the radius of the reference ball.

7. The method for generating high volume fraction particles based on grid points / grids according to claim 1, characterized in that, In step S18, the traversal process includes the following steps: S181. Calculate the lower and upper bound coordinates of the cube circumscribed by the repulsion sphere in the z-direction, denoted as z0_out and z1_out respectively. Then z0_out = z0 - R_ex and z1_out = z0 + R_ex. When placing the cube based on the grid points, calculate the lower and upper bound grid numbers in the z-direction, denoted as i0_z and i1_z respectively. Then i0_z = floor(z0_out / w0) + 1 and i1_z = ceil(z1_out / w0) - 1. When placing the cube based on the grid, calculate the lower and upper bound grid numbers in the z-direction, denoted as i0_z and i1_ If z is the center of the repulsion sphere, then i0_z = floor(z0_out / w0) + 1, i1_z = ceil(z1_out / w0), and calculate the layer number and row number of the grid where the center of the repulsion sphere is located, denoted as center_z and center_y respectively, then center_z = floor(z0 / w0) + 1, center_y = floor(y0 / w0) + 1; floor() means round down, ceil() means round up, y0 and z0 are the y and z coordinates of the center of the repulsion sphere, and R_ex is the radius of the repulsion sphere; S182. Let num = 0, let num2 = 0, let i_z = i0_z; S183. Determine if i_z is less than or equal to zero. If it is, set i_z_regin = i_z + m and jump to step S185. If it is not less than or equal to zero, proceed to the next step. S184. Determine if i_z is greater than m. If it is, let i_z_regin = i_z - m; if it is not, let i_z_regin = i_z. When placing based on grid points, calculate the square of the vertical distance from the center of the sphere to the current grid point layer, denoted as dz2, then dz2 = (z0 - i_z * w0). 2 , calculate the radius of the reference circle formed by the intersection of the current grid point layer and the repulsion sphere, denoted as r_i_z, then Calculate the lower and upper grid point numbers in the y direction of the grid points located inside the repulsion sphere in the current grid point layer, denoted as i0_y and i1_y respectively. Then i0_y = floor((y0 - r_i_z) / w0) + 1, i1_y = ceil((y0 + r_i_z) / w0) - 1; When placing based on the grid, determine the reference plane of the current grid layer and calculate the z coordinate z_ref of the reference plane. If i_z < center_z, the upper surface of the current grid layer is the reference plane of the current grid layer, z_ref = i_z * w0. If i_z = center_z, the horizontal plane where the center of the sphere is located is the reference plane of the current grid layer, z_ref = 0. If i_z > center_z, the lower surface of the current grid layer is the reference plane of the current grid layer, z_ref = (i_z - 1) * w0. Calculate the square of the minimum vertical distance from the center of the sphere to the reference plane of the current grid layer, denoted as dz2, dz2 = (z0 - z_ref). 2 , calculate the radius of the reference circle formed by the intersection of the reference plane of the current grid layer and the repulsion sphere, denoted as r_i_z, then Calculate the lower and upper grid numbers in the y direction of the grids that intersect with the repulsion sphere in the current grid layer, denoted as i0_y and i1_y respectively. Then i0_y = floor((y0 - r_i_z) / w0) + 1, i1_y = ceil((y0 + r_i_z) / w0); floor() represents rounding down, ceil() represents rounding up. y0 and z0 are the y and z coordinates of the center of the repulsion sphere, and R_ex is the radius of the repulsion sphere. S186. Let i_y = i0_y; S187. Determine if i_y is less than or equal to zero. If it is less than or equal to zero, set i_y_regin = i_y + m and jump to step S189. If it is not less than or equal to zero, proceed to the next step. S188. Determine if i_y is greater than m. If it is, let i_y_regin = i_y - m; if it is not, let i_y_regin = i_y. When placing based on grid points, calculate the square of the vertical distance from the center of the reference circle to the current grid point row, denoted as dy2, then dy2 = (y0 - i_y * w0) 2 , calculate half of the length of the chord formed by the intersection of the line where the current grid point row is located and the reference circle, denoted as dx, then Calculate the lower and upper grid point numbers in the x direction of the grid points within the exclusion sphere in the current grid point row of the current grid point layer, denoted as i0_x and i1_x respectively. Then i0_x = floor((x0 - dx) / w0) + 1, i1_x = ceil((x0 + dx) / w0) - 1; When placing based on the grid, determine the reference plane of the current grid row and calculate the y coordinate y_ref of the reference plane. If i_y < center_y, the back surface of the current grid row is the reference plane of the current grid row, y_ref = i_y * w0. If i_y = center_y, the vertical plane passing through the center of the reference circle and parallel to the front and back surfaces of the grid is the reference plane of the current grid row, y_ref = 0. If i_y > center_y, the front surface of the current grid row is the reference plane of the current grid row, y_ref = (i_y - 1) * w0. Calculate the square of the vertical distance from the center of the reference circle to the reference plane of the current grid row, denoted as dy2, dy2 = (y0 - y_ref) 2 , calculate half of the length of the chord formed by the intersection of the intersection line of the reference plane of the current grid layer and the reference plane of the current grid row and the reference circle, denoted as dx, then Calculate the lower and upper grid point numbers in the x direction of the grids intersecting with the exclusion sphere in the current grid row of the current grid point layer, denoted as i0_y and i1_y respectively. Then i0_x = floor((x0 - dx) / w0) + 1, i1_x = ceil((x0 + dx) / w0); floor() represents rounding down, ceil() represents rounding up, y0 and x0 are the y and x coordinates of the center of the exclusion sphere S1810, Let i_x = i0_x; S1811. Determine if i_x is less than or equal to zero. If it is, set i_x_regin = i_x + m and jump to step S1813. If it is not less than or equal to zero, proceed to the next step. S1812. Determine if i_x is greater than m. If it is, let i_x_regin = i_x - m; if it is not, let i_x_regin = i_x. S1813. Let num2 = num2 + 1. If the occupation status of the cell point / cell in the i_z_regin row and i_x_regin column of the i_y_regin layer is occupied, then let num = num + 1; otherwise, num remains unchanged. S1814. Determine if i_x is greater than or equal to i1_x. If it is, proceed to the next step; if it is less than, let i_x = i_x + 1 and jump to step S1811. S1815. Determine if i_y is greater than or equal to i1_y. If it is, proceed to the next step; if it is less than, let i_y = i_y + 1 and jump to step S187. S1816. Determine if i_z is greater than or equal to i1_z. If it is, proceed to the next step; if it is less than, let i_z = i_z + 1 and jump to step S183. S1817. Calculate the local volume fraction vfl = num / num2 and save it to the local volume fraction array, then end the traversal.

8. The method for generating high volume fraction particles based on grid points / grids according to claim 1, characterized in that, In step S33, the generation of the supplementary ball includes the following steps: S331. Let the x, y, and z coordinates of the center of the i-th ball be x0, y0, and z0, respectively. Let xyz_min be equal to the ball radius and xyz_max be equal to the region side length minus the ball radius. Execute S332, S334, and S336 in parallel. S332. If z0 is less than xyz_min, let z1 equal z0 plus the region side length, let z3 equal zero, and proceed to the next step; if z0 is greater than xyz_max, let z1 equal z0 minus the region side length, let z3 equal the region side length, and proceed to the next step; if z0 is greater than or equal to xyz_min and less than or equal to xyz_max, let z3 equal -1, and jump to step S338. S333. A new supplementary ball is generated with center coordinates (x0, y0, z1) and radius equal to the radius of the thrown ball. The center coordinates and radius of this supplementary ball are recorded and added to the supplementary ball set. Then, proceed to step S338. S334. If y0 is less than xyz_min, let y1 equal to y0 plus the region side length, let y3 equal to zero, and proceed to the next step; if y0 is greater than xyz_max, let y1 equal to y0 minus the region side length, let y3 equal to the region side length, and proceed to the next step; if y0 is greater than or equal to xyz_min and less than or equal to xyz_max, let y3 equal to -1, and jump to step S338. S335. A new supplementary ball is generated with center coordinates (x0, y1, z0) and radius equal to the radius of the thrown ball. The center coordinates and radius of this supplementary ball are recorded and added to the supplementary ball set. Then, proceed to step S338. S336. If x0 is less than xyz_min, set x1 to x0 plus the region side length, set x3 to zero, and proceed to the next step; if x0 is greater than xyz_max, set x1 to x0 minus the region side length, set x3 to the region side length, and proceed to the next step; if x0 is greater than or equal to xyz_min and less than or equal to xyz_max, set x3 to -1, and jump to step S338. S337. A new supplementary ball is generated with center coordinates (x1, y0, z0) and radius equal to the radius of the thrown ball. The center coordinates and radius of this supplementary ball are recorded and added to the supplementary ball set. Then, proceed to step S338. S338, Wait for S332-S337 to complete execution; S339, execute S3310, S3314, S3318 and S3322 in parallel; S3310. If x3 is greater than or equal to zero and y3 is greater than or equal to zero, proceed to the next step; otherwise, jump to step S3326. S3311. Using x3 and y3 as reference points, calculate the distance d30 from the reference point to (x0, y0). S3312. If d30 is less than the radius of the ball being thrown, proceed to the next step; if d30 is greater than or equal to the radius of the ball being thrown, jump to step S3326. S3313. A new supplementary ball is generated with center coordinates (x1, y1, z0) and radius equal to the radius of the thrown ball. Record the center coordinates and radius of this supplementary ball and add it to the supplementary ball set. Jump to step S3326. S3314. If y3 is greater than or equal to zero and z3 is greater than or equal to zero, proceed to the next step; otherwise, jump to step S3326. S3315. Using the x and y coordinates of y3 and z3 as reference points, calculate the distance d30 from the reference point to (y0, z0). S3316. If d30 is less than the radius of the ball being thrown, proceed to the next step; if d30 is greater than or equal to the radius of the ball being thrown, jump to step S3326. S3317. A new supplementary ball is generated with center coordinates (x0, y1, z1) and radius equal to the radius of the thrown ball. Record the center coordinates and radius of this supplementary ball and add it to the supplementary ball set. Jump to step S3326. S3318. If z3 is greater than or equal to zero and x3 is greater than or equal to zero, proceed to the next step; otherwise, jump to step S3326. S3319. Using x3 and z3 as reference points, calculate the distance d30 from the reference point to (x0, z0). S3320. If d30 is less than the radius of the ball being thrown, proceed to the next step; if d30 is greater than or equal to the radius of the ball being thrown, jump to step S3326. S3321. A new supplementary ball is generated with center coordinates (x1, y0, z1) and radius equal to the radius of the thrown ball. Record the center coordinates and radius of this supplementary ball and add it to the supplementary ball set. Jump to step S3326. S3322. If x3 is greater than or equal to zero, y3 is greater than or equal to zero, and z3 is greater than or equal to zero, proceed to the next step; otherwise, jump to step S3326. S3323. Using x3, y3, z3 as reference points, calculate the distance d30 from the reference points to (x0, y0, z0). S3324. If d30 is less than the radius of the ball being thrown, proceed to the next step; if d30 is greater than or equal to the radius of the ball being thrown, jump to step S3326. S3325. A new supplementary ball is generated with center coordinates (x1, y1, z1) and radius equal to the radius of the thrown ball. Record the center coordinates and radius of this supplementary ball and add it to the supplementary ball set. Jump to step S3326. S3326, Wait for S3310-S3325 to complete.

Images