A Universal Method for Modeling and Characterizing Non-Circular Packing Systems Based on n-Point Correlation Functions

A universal method for modeling and characterizing non-circular particles is developed. The n-point correlation functions (n = 1, 2 and 3) are efficiently computed with a GPU parallel computing procedure. An algorithm for dynamic packing of impenetrable non-circular particles is developed based on the fast estimation of overlap information using a one-point correlation function. The packing algorithm is independent of particle shape and proved to be reliable by examples of polygons and super-ellipses. In addition, penetrable packings are generated in an efficient and precise way. Using a two-point correlation function, these non-circular packs are accurately characterized and compared in terms of features such as penetrable and impenetrable, packing fraction and particle shape. In addition, three-point correlation functions are also illustrated and discussed.


Introduction
Particulate materials are the most widely used heterogeneous materials in the engineering practice [1]. Typical two-phase composites can be seen as particulate or skeletal material [2], which can be geometrically considered as impenetrable particle model and penetrable particle model. Take cementitious materials, for example. Cement paste, mortar or aggregate can be simulated as two-phase particulate material at different scales [3][4][5][6]. The particulate structure can be either a solid phase [7] (impenetrable model) or a porous phase [8] (penetrable model) accordingly.

Packing of Particulate Model
Impenetrable particle packing is generally called "packing problem", the final particulate structure of which contains no overlap between particles. The target particulate structure in this paper is that the particles are randomly and discretely distributed, surrounded by the other phase. Therefore, the particulate phase itself is not mechanically stable. For this kind of packing problem, the algorithm of random sequential addition [9] (RSA) is the most commonly used method of random static packing based on Monte Carlo sampling mechanism. RSA refers to a process where particles are randomly and sequentially introduced into a system without overlap to the previous ones. It is a fundamental method for particle packing. However, in the final stage, the elapsed time is unpredictable, and it also has a relatively low "saturated" packing fraction, especially for mono-size particles. i.e., the packing fraction of a saturated random packing for spheres was found to be φ = 0.382 [10] and for disks φ = 0.547. To obtain a higher packing fraction, an optimized packing approach is often used in modeling the particulate system [11][12][13][14][15][16][17][18]; most of these works are based on mathematical modeling of the relations between geometric objects and thus reduce the optimal packing problem to a nonlinear programming problem. However,

Overlap Determination for Non-Circular Particles
Particle shape is one of the most important factors that determine the macroscopic properties. In recent years, researchers have extended their attention to non-spherical particles, including 2D and 3D, i.e., ellipses/ellipsoids, super-ellipses/super-quadrics, polygons/polyhedrons, composite particles, discretized particles, etc. The DEM originally developed by Cundall and Strack [19] has been used worldwide to study many different shaped particles. In this method, the overlap between particles is considered to represent the particle deformations, which is used to estimate the elastic, plastic and frictional contact forces between particles. Contact detection and overlap determination for disks/spheres are simple and efficient. However, for non-spherical particles, both the contact detection and the overlap determination are complicated, highly dependent on the specific shape characteristic.
As reviewed in [25,26], composite particles and discretized particles are all able to simulate any non-spherical particles, or any shape theoretically. However, the resolution of these two methods is dependent on the number of units. In practice, it is computationally expensive if we use these two methods to build a well approximated shape. A combination of super-ellipses/super-quadrics and polygons/polyhedrons represents most of the common non-circular/non-spherical particles. Contact detection between particles for these models was well solved by the geometric potential method [27] and the separation axis theorem [28], respectively. Once two particles are confirmed to contact, the following step is to determine the overlap information, including the overlap magnitude, the action point and the normal direction. However, none of these features can be readily defined and obtained by a unified approach due to lack of sound contact theories. Many existing contact models specify the features independently, as we summarized in Table 1. Gradient vector of an inner potential particle Orientation discretization database solution [42] Polygon, Ellipse, Ellipsoid and others.
Average of centers of overlap cells Averaged normal vector of the cell at the surface of the particle

Characterization of Packs
The particulate system is typically described by various statistical descriptors. Over the past several decades, various descriptors have been proposed to characterize the structure. Among them, the n-point correlation functions are able to encompass all the details of the particulate structure [43]. The most common correlation function to characterize particulate structure is the one-point correlation function, also known as the packing fraction in a physical perspective, or the area fraction in 2D and the volume fraction in 3D. Although it is the simplest correlation function, many well-known bounds and prediction models contain only this sole information: for example, the Hashin-Shtirkman (H-S) bounds [44,45], the Mori-Tanaka model [46], the Maxwell-Wagner model [47] and so on. Higher-order bounds for effective permittivity with higher-order correlation functions were derived by Beran [48] and simplified by Torquato [49] and Miltion [50], respectively. Compared with the Hashin-Shtirkman bounds, these higher order bounds added a microstructural parameter ζ. Beran and Molneux [51], McCoy [52] and Milton and Phan-Thien [53] derived higher order bounds for effective shear and bulk moduli with an additional microstructural parameter η. Both ζ and η are triple integrals involving one-, two-and three-point correlation functions, which means that these higher-order bounds contain far more microstructural information than lower-order bounds such as H-S bounds, and enable one to characterize the microstructure more accurately [54].

A Universal Method for Modeling and Characterizing
In this paper, one-, two-and three-point correlation functions are computed accurately with the GPU-based parallel method. The one-point correlation function is the probability of finding a given phase at any location. For a particulate system, it can be calculated by judging all points' location. This process is highly parallelizable, so the GPU-based parallel method will increase the computing speed dramatically. With the advantage of a highly efficient calculation of the one-point correlation function, a novel method is introduced to estimate the comprehensive overlap information between two particles in the dynamic packing of impenetrable particles. This method is described in detail in non-circular models, polygons and super-ellipses. The polygons are composed of vertices and have a unique characteristic of shape corners, while the super-ellipse model is described by a quadric curve, |x/a| 2m + |y/b| 2m = 1, where m is the shape parameter characterizing the geometries, a and b are semiaxes. In addition, a penetrable particle model with a precise packing fraction can be obtained with high efficiency. The two-point correlation functions are mainly used to characterize the generated packs because the function curves can be directly compared. The three-point correlation functions are also illustrated in characterization.

Governing Equations
As a well-established dynamic packing method, the DEM has been successfully applied in modeling non-circular/non-spherical particle [25,26,55]. The DEM provides us with an effective way to simulate and analyze the movement of a particulate system; thus, the basic principle of the DEM is utilized in this paper. In the DEM, any particle can have two types of motion: translational and rotational, which are determined by Newton's second law of motion as given below: where m i , I i , v i and ω i are the mass, inertia, translational and angular velocities of particle i, respectively, and k c is the number of particles in collision with the particle. As our targeted particulate structure is that the particles are all randomly and discretely distributed, thus in the dynamic packing process, gravity is not considered. The forces considered in this paper are normal force F n,ij and tangential force F t,ij . The torques acting on particle i by particle j are consisted of two parts: M n,ij , generated by the normal force, and M t,ij , generated by the tangential force that caused particle i to rotate. In addition, the normal force F n,ij are also consisted of two parts: the normal elastic force F cn,ij and the viscous damping force F dn,ij , and the F dn,ij is used to dissipates the energy of the system.
Similarly, the tangential force F t,ij consists of the tangential elastic force F ct,ij and the tangential damping force F dt,ij , According to Coulomb's law, F t,ij = min F ct,ij + F dt,ij , µ F n,ij , where µ is the sliding friction coefficient. The equations used to calculate the particle-particle interaction forces and torques are listed in Table 2. Table 2. Components of forces and torques acting on particle.

Forces and Torques Symbols Equations
Normal elastic force F cn,ij −k n,ij δ n,ij n ij Normal damping force Coulomb friction force Torque by normal force Torque by tangential force For symbols in this table, k and γ are the spring coefficient and viscous damping coefficient for the contact model, respectively, and the subscripts of n or t indicate the variables are in the normal or tangential directions. They are assumed to be constants during the collisions of particles.
As summarized above in Table 1, for different models, methods of determining overlap information varied considerably. The orientation discretization database solution [42] almost provides a universal method for the different shaped particles, but the databases take time to build and not suitable for many different shapes. Another approach to describe a particle of arbitrary shape is digitization [56], which will take a huge computational effort if accurate results are ensured. We propose a new way of determining overlap information based on the review above. The key features of overlap information (i.e., magnitude, action point and normal direction) are listed in Table 3, and the corresponding contact models are shown in Figure 1. almost provides a universal method for the different shaped particles, but the databases take time to build and not suitable for many different shapes. Another approach to describe a particle of arbitrary shape is digitization [56], which will take a huge computational effort if accurate results are ensured. We propose a new way of determining overlap information based on the review above.
The key features of overlap information (i.e., magnitude, action point and normal direction) are listed in Table 3, and the corresponding contact models are shown in Figure 1. Table 3. Contact overlap features used in this paper.

Magnitude Action Point
Normal Direction Area (Volume in 3D): Ratio of sampling points Based on the features, we determine the normal elastic force. For disks or spheres, the overlap magnitude is defined as a distance , = + − . But for non-circular/non-spherical particles, it is far more effective when using area or volume (in 3D) [55,57]. Hence, we define the overlap magnitude as a function of the overlap area , which is , = √2 ⁄ . Thus, the normal elastic force , is defined as: According to the real motion of the particle, it is more reasonable to define the normal damping force Based on the features, we determine the normal elastic force. For disks or spheres, the overlap magnitude is defined as a distance δ n,ij = R i + R j − R ij . But for non-circular/nonspherical particles, it is far more effective when using area or volume (in 3D) [55,57]. Hence, we define the overlap magnitude as a function of the overlap area S o , which is δ n,ij = √ 2S o /π. Thus, the normal elastic force F cn,ij is defined as: According to the real motion of the particle, it is more reasonable to define the normal damping force F dn,ij on the basis of the relative velocity v ij in the action point p c , which is calculated according to Figure 2a.
therefore, we obtain the normal damping force, Materials 2022, 15, 5991 6 of 17 which is , = √2 ⁄ . Thus, the normal elastic force , is defined as: According to the real motion of the particle, it is more reasonable to define the normal damping force , on the basis of the relative velocity in the action point , which is calculated according to Figure 2a. In addition, the tangential elastic force F ct,ij is defined on the basis of the relative velocity v ij and the time step ∆t, and the displacement in the tangential direction is v t,ij ·∆t in a single time step. According to the resolution of the relative velocity in Figure 2b, Similarly, the tangential damping force can be defined in the opposite direction,

Determination of Overlap Area
As the one-point correlation function is computed with the GPU parallel method, the overlap area can be estimated with high efficiency. Once two particles are detected to overlap, we set a rectangle box to cover them completely, and uniformly generate sampling points, as can be seen from Figure 3. Then all the sampling points are classified by judging their locations. The portion of the points that land in both particles, as shown in Figure 3b, are employed to estimate the overlap information. In this work, the number of the sampling points is not less than 10 6 to guarantee the precision. Then we execute the kernels on the GPU, when a kernel is launched, a grid of threads that are organized in a 3D hierarchy is generated, with each gird being organized into array of thread blocks, and each block containing up to 1024 threads. Therefore, all the sampling points are allocated and judged in millions of threads (i.e., the number of computing cores N thds = 1024 × 1024 in the case of 2D block and 2D grid). By contrast, each computer has a limited number of CPU cores. To illustrate the parallel speedup ratio of the GPU to the CPU, we tested the speed of the one-point correlation function for the packing of a regular pentagon. An ordinary computer was used, i.e., the CPU is AMD RYZEN 7 3800X, and "GPU" is NVIDIA GeForce RTX 2070. The results are listed in Table 4.
What is worth mentioning is that all points on the "layer" are also classified into the overlap region, to guarantee that, as long as two particles collide theoretically, the overlap will be detected. According to the proportion of points in the overlap region, the overlap area S o is obtained.
For the non-circular/non-spherical model, a reasonable choice of action point should be the centroid. Under normal circumstances, working out the centroid of the overlapping area/volume is computationally expensive. However, in our method, the overlap region consisted of a set of points, the average location of which is exactly the centroid. Therefore, it can be immediately determined once the overlap points are identified.

Dynamic Packing Scheme
The flowchart of the dynamic packing scheme on non-spherical particles in this work is shown in Figure 4a. We set the initial material parameters including: the particle number, the function parameters, the mass, moment of inertia, the particle size and direction angles, etc. In order to eliminate the wall effect, periodic boundary conditions are applied in this paper. First, we generate the initial packings without the overlapping constraints, as shown in Figure 4b. Then direct contact detection between the two particles can be conducted by the geometric potential method and the separation axis theorem for the super-ellipses and  Note: the packing system is closely related to the execution time. In this test, the packing fraction φ a = 0.5, total particle number N p = 254.

Dynamic Packing Scheme
The flowchart of the dynamic packing scheme on non-spherical particles in this work is shown in Figure 4a.  We set the initial material parameters including: the particle number, the function parameters, the mass, moment of inertia, the particle size and direction angles, etc. In order to eliminate the wall effect, periodic boundary conditions are applied in this paper. First, we generate the initial packings without the overlapping constraints, as shown in  We set the initial material parameters including: the particle number, the function parameters, the mass, moment of inertia, the particle size and direction angles, etc. In order to eliminate the wall effect, periodic boundary conditions are applied in this paper. First, we generate the initial packings without the overlapping constraints, as shown in Figure 4b. Then direct contact detection between the two particles can be conducted by the geometric potential method and the separation axis theorem for the super-ellipses and the polygons, respectively. For the determined overlapped particles, we estimate the overlap area (volume in 3D) and centroid with the GPU-based parallel method. Forces and torques can be obtained by the equations list above. The positions are calculated by the equations of motion, which can be integrated by the Verlet scheme [58] as follows: When the particle information is updated, it goes back to the procedure of contact detection, until no overlap is detected. Finally, we store all the positions and direction angles. The results are visualized as Figure 4c.

Penetrable Packing Model
By using the Poisson limit theorem in statistics, the total inclusion rate for a penetrable system has the quantitative relationship between the inclusion rate and the inclusion number [8]. In this paper, these inclusions can be expressed by particles; thus, the packing fraction φ and the particle number N p have the following relationships when N p = ∞.
Therefore, the penetrable models for the various packing fraction can be generated. However, as the particle number is limited, the real packing fraction is only an approximate value, fluctuating around theoretical one. For the continuum percolation models, the pore phase can be simplified to penetrable particles as the porous phase, and by repeatedly generating penetrable models, the percolation threshold of different non-spherical models can be derived. Considering the huge computational cost, the accuracy of the real packing fraction was not considered. In this work, the penetrable packing models are generated efficiently with an accurate packing fraction.
The best approach to acquire the real packing fraction is systematic point sampling, which is equivalent to the one-point correlation function. In our approach, due to the high efficiency of the one-point correlation function, one can repeatedly generate penetrable packings and at the same time calculate the real packing faction, until the one within close tolerance appears. The flowchart of the packing algorithm is shown in Figure 5a.
For example, a penetrable super-ellipse and a regular square with a packing fraction of 0.8 are generated with different precision, the time elapsed in increasing as upgrading precision, as shown in Table 5.

n-Point Correlation Functions
The n-point correlation function is a generalized definition of the correlation functions. The higher-order correlation functions are able to provide more information about the geometric features of a generation, for example, the standard two-point correlation function, the three-point correlation function, and various other two-point correlation functions such as Lineal-Path Function, Chord-Length Density Function, Pore-Size Functions and so on. In this work, only the standard one-, two-and three-point correlation functions, which are included in higher-order parameters and , are concerned.
The n-point correlation function is the probability that the n given points with the locations x 1 , x 2 , … , x will be in the same phase , and it can be expressed as the expectation (or average) of the multiplication of the indicator functions at the n locations, where the angular bracket ⟨⋯ ⟩ denotes the expectation or the ensemble average, ( ) represents indicator function for phase , where is the region occupied by phase with the packing fraction . One can define the one-, two-and three-point correlation functions when n = 1, 2 and 3:  For both generations, the particle number (N p = 819) is the same because the equivalent diameter D eq is used for non-circular particles.

n-Point Correlation Functions
The n-point correlation function is a generalized definition of the correlation functions. The higher-order correlation functions are able to provide more information about the geometric features of a generation, for example, the standard two-point correlation function, the three-point correlation function, and various other two-point correlation functions such as Lineal-Path Function, Chord-Length Density Function, Pore-Size Functions and so on. In this work, only the standard one-, two-and three-point correlation functions, which are included in higher-order parameters ζ and η, are concerned.
The n-point correlation function is the probability that the n given points with the locations x 1 , x 2 , . . . , x n will be in the same phase i, and it can be expressed as the expectation (or average) of the multiplication of the indicator functions at the n locations, where the angular bracket · · · denotes the expectation or the ensemble average, I (i) represents indicator function for phase i, where V i is the region occupied by phase i with the packing fraction φ i . One can define the one-, two-and three-point correlation functions when n = 1, 2 and 3: For example, for a two-phase composite material, the gray phase represents the "aggregate", while the white phase is the matrix. Schematic representation of the one-, twoand three-point correlation functions are shown in Figure 6a. S aa is another form of S (i) 2 (r) when i represents the aggregate phase. In addition, S mm means that phase i of interest is the matrix, but there exists another function S am , which means that one point lands on the aggregate, and the other point lands in the matrix. Therefore, it is obvious that: Materials 2022, 15, x FOR PEER REVIEW 10 of 17 For example, for a two-phase composite material, the gray phase represents the "aggregate", while the white phase is the matrix. Schematic representation of the one-, twoand three-point correlation functions are shown in Figure 6a.
is another form of 2 ( ) ( ) when represents the aggregate phase. In addition, means that phase of interest is the matrix, but there exists another function , which means that one point lands on the aggregate, and the other point lands in the matrix. Therefore, it is obvious that: Similarly, means 3 ( ) ( 1 , 2 , ) when the phase of interest is the aggregate. In this work, particles in the realizations are randomly distributed and the periodic boundary conditions are considered. It is reasonable to assume all generations are statistically homogeneous and isotropic; therefore, these functions can be obtained by the random sampling technique. These correlation functions can be calculated as follows.
The one-point correlation function can be calculated by systematic point sampling method, and the points can be generated as both randomly distributed or uniformly arrayed. Considering the overlap area calculation that was mentioned above and was shown in Figure 3b, the uniformly arrayed points are used as the sampling points. Moreover, the number of points in each sampling is 10 6 − 10 8 accordingly.
(a) (b) (c) The two-point correlation function ( ) is the probability of the two points at the phase . The result of ( ) is not a single number, because there is a probability value for each length . Theoretically = [0, ∞], for the purpose of comparison, only the selected values of are calculated. For example, if one wants to compute and , the value often started at an extremely small value [54]. In our work, the start point of the segment is uniformly arrayed, and the end point is determined according to and a random angle. The number of sampling lines is chosen 10 5 − 10 7 accordingly. When = 0 and = ∞, the limiting values are obtained: The three-point correlation function ( 1 , 2 , ) is not commonly used because it is computationally expensive. For example, when 1 , 2 and are all divided into 100 values, there will be 100 3 different patterns of the sampling triangle; for each pattern, the Similarly, S aaa means S (i) 3 (r 1 , r 2 , θ) when the phase i of interest is the aggregate. In this work, particles in the realizations are randomly distributed and the periodic boundary conditions are considered. It is reasonable to assume all generations are statistically homogeneous and isotropic; therefore, these functions can be obtained by the random sampling technique. These correlation functions can be calculated as follows.
The one-point correlation function S a can be calculated by systematic point sampling method, and the points can be generated as both randomly distributed or uniformly arrayed. Considering the overlap area calculation that was mentioned above and was shown in Figure 3b, the uniformly arrayed points are used as the sampling points. Moreover, the number of points in each sampling is 10 6 -10 8 accordingly.
The two-point correlation function S aa (r) is the probability of the two points at the phase a. The result of S aa (r) is not a single number, because there is a probability value for each length r. Theoretically r = [0, ∞], for the purpose of comparison, only the selected values of r are calculated. For example, if one wants to compute ζ and η, the r value often started at an extremely small value [54]. In our work, the start point of the segment is uniformly arrayed, and the end point is determined according to r and a random angle. The number of sampling lines N line is chosen 10 5 -10 7 accordingly. When r = 0 and r = ∞, the limiting values are obtained: The three-point correlation function S aaa (r 1 , r 2 , θ) is the probability of the three points at the aggregate phase. The sampling triangle can be determined by r 1 , r 2 and θ. Theoretically r 1 , r 2 = [0, ∞], θ = [0, π]. S aaa (r 1 , r 2 , θ) is not commonly used because it is computationally expensive. For example, when r 1 , r 2 and θ are all divided into 100 values, there will be 100 3 different patterns of the sampling triangle; for each pattern, the number of the sampling number N 3pt in this work is 10 5 -10 7 accordingly. In addition, the limiting values are obtained when r 1 , r 2 → 0 and r 1 , r 2 → ∞, θ = 0: lim

Discussion
Based on the packing algorithm, the various non-circular packing systems are generated. The n-point correlation functions describe the probabilities of the different phase encounters and other geometric features and aim to encompass all details of the packing system, i.e., the two-point correlation function has been extensively used in the characterization of short-range information. In what follows, we characterize the packing systems with the different packing fraction and the particle shape using the two-point correlation functions. In addition, the three-point correlation function is illustrated and discussed.

Packing Fraction
We generated the realizations with the packing fraction φ a = 0.1 to φ a = 0.8 using the packing method described above in this work. Figure 7a illustrates examples of the generations for the super-ellipse (a/b, m = 2) with the packing fraction φ a = 0.1, 0.4 and 0.8. Figure 7b correspondingly presents the penetrable ones.

Discussion
Based on the packing algorithm, the various non-circular packing systems are generated. The n-point correlation functions describe the probabilities of the different phase encounters and other geometric features and aim to encompass all details of the packing system, i.e., the two-point correlation function has been extensively used in the characterization of short-range information. In what follows, we characterize the packing systems with the different packing fraction and the particle shape using the two-point correlation functions. In addition, the three-point correlation function is illustrated and discussed.

Packing Fraction
We generated the realizations with the packing fraction = 0.1 to = 0.8 using the packing method described above in this work. Figure 7a  The two-point correlation function S aa (r) of the particulate phase for the super-ellipse with the packing fraction φ a = 0.1 ∼ 0.8 is shown in Figure 7c,d. When r = 0, only two types of the function exist: S aa and S mm , at this time S aa + S mm = 1 which is actually a one-point correlation function. As r grows, S aa is decreased due to the emergence of S am . By comparing Figure 7c,d, we can clearly see that for impenetrable models, the decay of S aa is faster, and obvious oscillations begin to show up when r ≥ D eq . At that point, the decay corresponds to their long-ranger values, while the penetrable results did not show clear oscillations. It is worth noting that when the packing fraction is very low, i.e., φ a = 0.1, S aa results are basically the same for the impenetrable and the penetrable packing models, because the generations are similar, as shown in Figure 7a. In addition, the speed at which the fluctuations level off is different for S aa in the impenetrable system.

Particle Shape
To study the effect of the particle shape of S aa (r), we generated a series of packing models with the same packing fraction φ a = 0.5. For simplicity without loss of generality, we used the special cases of the super-ellipse. First, we let m = 1, super-ellipse is equivalent to ellipse, we compared S aa (r) for different aspect ratio κ. Then we let κ = 1, and computed S aa (r) for the different shape parameter m. The results are shown in Figure 8a, a growing κ led to a faster decay, and the fluctuations were smaller. There will be an obvious difference when κ is changing from 1 to 5, while when κ is fixed, as we see in Figure 8b, the overall pattern of S aa (r) is similar and the slope before first trough of wave is almost identical, but a visible difference appears near the first peak of wave.
Then we compared S aa (r) of the different regular polygons: triangle, cube and hexagon. As we can observe in Figure 8c, S aa (r) of the regular cube has the largest trough, because in this paper, the equivalent circular diameter D eq is used for non-circular particles, and the regular square has the longest line segment of these three shapes. Moreover, the regular square has the fastest decay. We also computed S aa (r) in penetrable models, and an example is shown in Figure 8d. A relatively smaller difference can be observed compared with the impenetrable ones.
Roundness, also called circularity, is a 2D measure of how closely the shape of an object approaches that of a mathematically perfect circle, which is defined as the perimeter ratio of the particle and the circle with the same area [59]. The roundness of the super-ellipse could be any value only by changing the aspect ratio κ, as shown in Table 6. Therefore, it is reasonable to compare the two particles with the same roundness, i.e., the regular square and the ellipse with an aspect ratio κ = 2.283670. S aa (r) are calculated for both impenetrable and penetrable models, and the resulting comparisons are shown in Figure 9a,b. We can see a clear difference in the impenetrable model and a visible distinction in the penetrable model, which means the two-point correlation functions are able to describe more morphological features besides roundness.
we used the special cases of the super-ellipse. First, we let = 1, super-ellipse is eq lent to ellipse, we compared ( ) for different aspect ratio . Then we let = 1 computed ( ) for the different shape parameter . The results are shown in F 8a, a growing led to a faster decay, and the fluctuations were smaller. There will b obvious difference when is changing from 1 to 5, while when is fixed, as we s Figure 8b, the overall pattern of ( ) is similar and the slope before first trough of w is almost identical, but a visible difference appears near the first peak of wave.  Then we compared ( ) of the different regular polygons: triangle, cube and hexagon. As we can observe in Figure 8c, ( ) of the regular cube has the largest trough, because in this paper, the equivalent circular diameter is used for non-circular particles, and the regular square has the longest line segment of these three shapes. Moreover, the regular square has the fastest decay. We also computed ( ) in penetrable models, and an example is shown in Figure 8d. A relatively smaller difference can be observed compared with the impenetrable ones.
Roundness, also called circularity, is a 2D measure of how closely the shape of an object approaches that of a mathematically perfect circle, which is defined as the perimeter ratio of the particle and the circle with the same area [59]. The roundness of the superellipse could be any value only by changing the aspect ratio , as shown in Table 6. Therefore, it is reasonable to compare the two particles with the same roundness, i.e., the regular square and the ellipse with an aspect ratio = 2.283670.
( ) are calculated for both impenetrable and penetrable models, and the resulting comparisons are shown in Figure 9a,b. We can see a clear difference in the impenetrable model and a visible distinction in the penetrable model, which means the two-point correlation functions are able to describe more morphological features besides roundness.
(a) (b) Figure 9. ( ) of different particle packing systems with the same circularity and the same packing fraction = 0.5: (a) Impenetrable packing system of regular square and ellipse; (b) Penetrable packing system of regular square and ellipse.

Three-Point Correlation Function
The three-point correlation functions ( 1 , 2 , ) for the generations of regular squares with the packing fraction = 0.7 are presented in Figure 10. We only show some selected values of for illustration. It can be observed that the Figure 9. S aa (r) of different particle packing systems with the same circularity and the same packing fraction φ a = 0.5: (a) Impenetrable packing system of regular square and ellipse; (b) Penetrable packing system of regular square and ellipse.

Three-Point Correlation Function
The three-point correlation functions S aaa (r 1 , r 2 , θ) for the generations of regular squares with the packing fraction φ a = 0.7 are presented in Figure 10.
the is doubled, which explains the fast decay of and the zig zag line afterward. When 1 = 0 or 2 = 0, it degenerates to the standard two-point correlation function ( ). It can be observed that the two-point correlation function is only a small section of the three-point correlation function ( 1 , 2 , ), even for the one selected . When calculating the function values for 3 ( ) ( 1 , 2 , ), 3 random sampling triangles are used. Each sampling triangle has 1 × 2 × shapes. The accuracies are often verified by these their degeneracies because the fundamentals of computing the one-, twoand three-point correlation functions are the same, that is to determine a single point's location. By comparing ( ) to ( ), as shown in Figure 8c, the three-point correlation characterizes the structure more comprehensively, while the two-point correlation can describe the obvious features, i.e., particle shape.

Conclusions
The n-point correlation functions, including the one-, two-and three-point correlations in this paper, are utilized for modeling and characterizing a non-spherical packing system.
A dynamic packing algorithm for impenetrable non-circular particles was developed based on the DEM. In this algorithm, a novel method of determining overlap information (i.e., overlap area and centroid) for non-circular particles was developed by means of efficiently calculating the one-point correlation function. In addition, a penetrable non- We only show some selected values of θ for illustration. It can be observed that the one-and the two-point correlation function are contained within the three-point correlation function, and the function can be degenerated. When r 1 = r 2 = 0, S aaa (r 1 , r 2 , θ) degenerates to the one-point correlation function S a . We let θ = 0, r 1 = r 2 = 0, as shown in Figure 10. The zig zag diagonal line is a special case of two-point correlation function, and it has the steepest slop. In this function, the two end points of r 1 and r 2 coincide and have to be judged twice. According to Equation (21), it is obvious that for the different phase of interest, S aaa (r 1 , r 2 , θ) in two-phase materials, S aaa + S aam + S amm + S mmm = 1, (26) and in this situation, when it degenerates to S aa (r), the S am is doubled, which explains the fast decay of S aa and the zig zag line afterward. When r 1 = 0 or r 2 = 0, it degenerates to the standard two-point correlation function S aa (r). It can be observed that the two-point correlation function S aa is only a small section of the three-point correlation function S aaa (r 1 , r 2 , θ), even for the one selected θ. When calculating the function values for S (i) 3 (r 1 , r 2 , θ), N 3pt random sampling triangles are used. Each sampling triangle has N r 1 × N r 2 × N θ shapes. The accuracies are often verified by their degeneracies because the fundamentals of computing the one-, twoand three-point correlation functions are the same, that is to determine a single point's location. By comparing S aaa (r) to S aa (r), as shown in Figure 8c, the three-point correlation characterizes the structure more comprehensively, while the two-point correlation can describe the obvious features, i.e., particle shape.

Conclusions
The n-point correlation functions, including the one-, two-and three-point correlations in this paper, are utilized for modeling and characterizing a non-spherical packing system.
A dynamic packing algorithm for impenetrable non-circular particles was developed based on the DEM. In this algorithm, a novel method of determining overlap information (i.e., overlap area and centroid) for non-circular particles was developed by means of efficiently calculating the one-point correlation function. In addition, a penetrable nonspherical packing algorithm with high precision and efficiency was proposed via the one-point correlation function, and it is applicable to arbitrary shape in principle.
With the packing algorithms, packs of non-circular particles, both impenetrable and penetrable, with a different packing fraction (φ = 0.1 − 0.8) and various geometry shapes (regular polygons and super-ellipses), have been generated. The two-point correlation function are chosen as a statistical descriptor in this work for characterizing the packs. For the impenetrable models, clear differences can be observed with different packing fractions or various geometry shapes, even if the roundness is the same. The differences in corresponding penetrable models are less pronounced. The three-point correlation function is illustrated in three selected θ. It characterizes far more details than two-point correlation function. Moreover, with the efficient computation of the one-, two-and three-point correlations, an effective material behavior of the particulate systems can be predicted by third-order bounds, which will provide us a straightforward route to quantitative characterization.