Third-Order Effective Properties for Random-Packing Systems Using Statistical Micromechanics Based on a GPU Parallel Algorithm in Fast Computing n-Point Correlation Functions

Estimating the effective properties of a particulate system is the most direct way to understand its macroscopic performance. In this work, we accurately evaluate the third-order approximations involving the three-point microstructural parameter ζ, which can be calculated from a triple integral involving 1-, 2-, and 3-point correlation functions. A GPU-based parallel algorithm was developed for quickly computing the n-point correlation functions, and the results agree well with analytical solutions. The effective thermal conductivity and diffusion coefficient are calculated by the third-order approximates for the random-packing systems of a super-ellipsoid. By changing the parameters of the super-ellipsoid, the particle-shape effect can be predicted for both the thermal conductivity and diffusion coefficient.


Introduction
The quantitative characterization of heterogeneous materials from their microstructure is a long-studied issue with applications in many scientific fields [1][2][3][4][5][6]. For a particulate system, the most direct way to understand its macroscopic performance is by estimating its effective properties [7]. According to the composite material theory, however, an exact analytical prediction model never exists even for the simplest class of transport processes, such as thermal conduction and diffusion, because an infinite amount of microstructural information is required in the form of n-point correlation functions. Bounds or prediction models involving sole 1-point correlation function, also known as the packing fraction or volume fraction, can only provide a relatively coarse estimation, i.e., the well-known Hashin-Shtrikman (H-S) bounds [8,9]. To improve these bounds, additional information on the statistics of the spatial distribution of the phases is needed.
Brown [10] first brought higher-order parameters to effective properties. Based on H-S bounds, Beran [11] derived third-order bounds for effective thermal conductivity using variational principles in a homogeneous and isotropic system, but the results are complicated due to their involvement of sixfold integrals. Torquato [12] and Milton [13] independently simplified these complex bounds and their formulas are expressed by a volume fraction, φ, and a third-order microstructural parameter, ζ, which is a triple integral involving 1-, 2-, and 3-point correlation functions. In particular, for an isotropic two-phase system, Torquato [14] derived a three-point approximation formula that has shown a good agreement with simulation data [15]. Therefore, with this parameter, ζ, higher-order bounds or approximations can be obtained directly.
Despite decades of progress in computer technology, a direct computation of the microstructural parameter ζ remains a challenge. Most works involving ζ have somehow bypassed the direct computation of this higher-order parameter [16][17][18][19], until Gillman and Matouš [20] proposed adaptive interpolation methods to reduce the complexity of the three-point correlation function, and the results were verified by a penetrable sphere model [21]; later, Gillman et al. [22] extended their work and another three-point parameter η was computed with high accuracy. In their work, 500-7200 CPU computing cores are employed. Hlushkou et al. [23] employed a non-uniform spatial grid for the construction of a sampling template to reduce computation cost because the region of integration near zero makes a substantial contribution to the integral. The results were verified with facecentered cubic (FCC) and body-centered cubic (BCC) sphere arrays by analytical results [24]. Li [25] computed the third-point parameter with different volume fractions of glass spheres, but this yielded the results of simple 2D slices. Therefore, the computation cost of the parameter remains enormous; an ordinary algorithm based on a CPU is not enough to study the effective properties of more complex particulate systems.
In this paper, by employing the non-uniform spatial grid-sampling method [23], a GPU-based parallel algorithm is developed to compute the microstructural parameter, ζ, and an effective thermal conductivity and diffusion coefficient for packing systems of super-ellipsoids is are calculated according to Torquato's third-order approximation [14]. The rest of this paper is organized as follows: Section 2 describes the fundamental theories; Section 2.1 describes the n-point correlation functions, including the standard form of 1-, 2-, and 3-point correlation functions; Section 2.2 details the effective material behavior, including the thermal conductivity and diffusion coefficient; in Section 3, we provide the details of a GPU-based parallel algorithm in computing the microstructural parameter ζ and its verifications; and the results and conclusions are presented in Sections 4 and 5.

Fundamental Theory
This section provides a basic description of n-point correlation functions and an overview of the higher-order statistical micromechanics for calculating the effective thermal conductivity and diffusion coefficient of statistically isotropic microstructures.

n-Point Correlation Functions
The n-point correlation function is a generalized definition of correlation functions. In this work, only standard one-, two-, and three-point correlation functions, which are included in the higher-order parameter ζ, are considered. The n-point correlation function is defined as the probability that n points with locations x 1 , x 2 , . . . , x n will be in the same phase i: it can be expressed as expectation (or average) of the multiplication of indicator functions at the n locations where the angular bracket · · · denotes the expectation or ensemble average, and I (i) represent the indicator function for phase i,

of 12
where V i is the region occupied by phase i with packing fraction φ i . Therefore, one can define one-, two-, and three-point correlation functions when n = 1, 2, and 3, Considering a two-phase composite material (i.e., concrete), the gray phase represents "aggregate" in this work, and the white phase is a matrix. A schematic representation of one-, two-, and three-point correlation functions are shown in Figure 1, where S aa is another form of S (i) 2 (r) when i represents an aggregate phase; in addition, S mm means the phase i of interest is a matrix. There is another function, S am , which means one point lands on the aggregate, and the other point lands in the matrix, it is obvious that: Similarly, S aaa means S (i) 3 (r 1 , r 2 , θ) when phase i of interest is aggregate: In addition, limiting values of S (i) 2 (r) can be obtained when r = 0 and r = ∞: lim r→∞ S (i) Similarly, limiting values of S (i) 3 (r 1 , r 2 , θ) can be obtained when r 1 , r 2 → 0 and r 1 , r 2 → ∞, θ = 0: Materials 2022, 15, x FOR PEER REVIEW 3 of 13 , , 1, 1, 1 , , Considering a two-phase composite material (i.e., concrete), the gray phase represents "aggregate" in this work, and the white phase is a matrix. A schematic representation of one-, two-, and three-point correlation functions are shown in Figure 1, where is another form of when represents an aggregate phase; in addition, means the phase of interest is a matrix. There is another function, , which means one point lands on the aggregate, and the other point lands in the matrix, it is obvious that: 1.
Similarly, means , , when phase of interest is aggregate: In addition, limiting values of can be obtained when 0 and ∞:

Third-Order Models of Effective Material Behavior
For isotropic two-phase materials, the well-known H-S bounds are generally regarded as the best second order bounds on effective properties such as magnetic permeability and electrical or heat conductivity and diffusivity, because of their mathematical equivalents. The H-S lower bound for effective thermal conductivity, κ e , is expressed as follows: The effective property κ e is dependent on the conductivity of matrix κ m , the conductivity of the aggregate phase κ a , and the volume fraction of aggregate φ a . Beran [11] improved the H-S bounds and derived the third-order bounds of κ e , which was then simplified by Torquato [12] and Milton [13] independently, the lower bound can be expressed as: For isotropic two-phase systems, Torquato [14] derived a three-point approximation that has shown good agreement with simulations for packing systems without large clustering. This third-order approximation is given as: where In this work, a strong contrast ratio between particle and matrix phases is assumed for effective thermal conductivity, i.e., κ a κ m ≈ ∞ Then, the Equation (15) yields: In these formulas, the third-order microstructural parameter ζ i (i = a, m) is introduced, involving 1-, 2-, and 3-point correlation functions. For 2D media: For 3D media: Based on the third-order approximate formula, the effective diffusion coefficient D e in the packing system can be evaluated also assuming a strong contrast ratio [19,26], i.e., In addition, for a two-phase material such as concrete, the aggregate phase contributes nothing to diffusion, and the matrix phase of the material (i.e., mortar) must be attributed to Materials 2022, 15, 5799 5 of 12 the whole effective property; therefore, the third-order approximation in terms of diffusion coefficient D e can be expressed as:

GPU-Based Parallel Algorithm
For particulate packing systems, a GPU-based parallel packing algorithm is developed and described in this section. The periodic samplings method is proposed to match the periodic packing system. The efficiency of this algorithm is compared with the results of the CPU. In addition, the results are compared with analytical solutions of a penetrable sphere packing system and the existing results of an impenetrable sphere packing system, which verified the algorithm.

Sampling of n-Point (n = 1, 2, and 3) Correlation Functions
The three-point microstructural parameter ζ, for both 2D and 3D systems, can be seen as a triple integral involving 1-, 2-, and 3-point correlation functions. The 1-point correlation function is the volume fraction φ i of the phase i. The 2-point correlation function is computed for each discrete value of r, according to the range of integration; theoretically, it has to be calculated from zero to infinity. The major challenge is the computation of 3-point correlation functions, because not only r 1 and r 2 have the integration range of zero to infinity, but also the angle θ between r 1 and r 2 must be considered θ = [0, π]. It can be observed from the formula that the region of integration near r 1 = 0 and r 2 = 0 makes an enormous contribution to the integral due to the existence of 1/r 1 and 1/r 2 -this will require a fine spatial grid in order to obtain an accurate result. However, a large number of discrete values of 2-and 3-point correlations will result in a huge computation cost.
Following the sampling method by Hlushkou et al. [23], for r, r 1 , and r 2 the nonuniform sampling method is employed to reduce the discrete sampling number: where D eq is the equivalent diameter of a particle; a 0 is a constant near zero, determining the smallest r, r 1 , and r 2 ; and A is a constant slightly larger than one, i.e., in this work, we take A = 1.11. All particle-packing systems considered in this study are employed in periodic boundary conditions. Thus, the periodic sampling method is proposed in this work to match the periodic-packing system. For S (i) 2 (r) and S (i) 3 (r 1 , r 2 , θ), the start points P 1 of a sampling line or sampling triangle are random or uniformly distributed in the system. Therefore, there are chances that the end points of line segment P 2 or other vertices of a triangle P 2 and P 3 will fall outside the box. If the sampling lines or triangles are not allowed to fall outside box, the "boundary effect" will emerge, as shown in Figure 2a. On the contrary, when the sampling line or triangle is much larger than the packing system, they will not characterize the packing system correctly, as shown in Figure 2b. In order to solve these conflicting problems, the end points must be placed in the periodic position, as can be seen in the 2D examples from Figure 3a,b. The combination of the periodic-particle-packing system and the periodic sampling method allows us to characterize approximately infinite systems.
Although the non-uniform sampling methods have been employed, the cost of computation remains enormous, especially for the 3-point correlation function. For example, after a convergence study in this work, the range of discrete values for r 1 or r 2 is chosen from 10 −15 D eq to 15D eq , the constant A = 1.11, θ is divided evenly into 100 parts, resulting in patterns of the sampling triangle N pattern = N r 1 × N r 2 × N θ = 357 × 357 × 100 = 12, 744, 900, and for each pattern, N sample = 2 × 10 6 are needed. It can be observed that the computation progress for each sampling triangle is highly parallelizable; Gillman et al. [22] used 500-7200 CPU cores to perform such challenging work. However, each computer has a limited number of CPU cores, i.e., intel Core i9-12900KS has only 16 CPU Cores; in contrast, Materials 2022, 15, 5799 6 of 12 a GPU is composed of thousands or millions of cores that can manage multiple threads simultaneously. In this study, a GPU-based parallel algorithm is developed to accelerate the computation. Although the non-uniform sampling methods have been employed, the cost of computation remains enormous, especially for the 3-point correlation function. For example, after a convergence study in this work, the range of discrete values for or is chosen from 10 to 15 , the constant 1.11, is divided evenly into 100 parts, resulting in patterns of the sampling triangle 357 357 100 12744900, and for each pattern, 2 10 are needed. It can be observed that the computation progress for each sampling triangle is highly parallelizable; Gillman et al. [22] used 500-7200 CPU cores to perform such challenging work. However, each computer has a limited number of CPU cores, i.e., intel Core i9-12900KS has only 16 CPU Cores; in contrast, a GPU is composed of thousands or millions of cores that can manage multiple threads simultaneously. In this study, a GPU-based parallel algorithm is developed to accelerate the computation. The flowchart of the GPU parallel algorithm for computing the n-point correlation function is shown in Figure 4. Firstly, the particle information in the system is loaded, including the boundary condition (i.e., periodic in this work), the particle coordinate Although the non-uniform sampling methods have been employed, the cost of computation remains enormous, especially for the 3-point correlation function. For example, after a convergence study in this work, the range of discrete values for or is chosen from 10 to 15 , the constant 1.11, is divided evenly into 100 parts, resulting in patterns of the sampling triangle 357 357 100 12744900, and for each pattern, 2 10 are needed. It can be observed that the computation progress for each sampling triangle is highly parallelizable; Gillman et al. [22] used 500-7200 CPU cores to perform such challenging work. However, each computer has a limited number of CPU cores, i.e., intel Core i9-12900KS has only 16 CPU Cores; in contrast, a GPU is composed of thousands or millions of cores that can manage multiple threads simultaneously. In this study, a GPU-based parallel algorithm is developed to accelerate the computation. The flowchart of the GPU parallel algorithm for computing the n-point correlation function is shown in Figure 4. Firstly, the particle information in the system is loaded, including the boundary condition (i.e., periodic in this work), the particle coordinate The flowchart of the GPU parallel algorithm for computing the n-point correlation function is shown in Figure 4. Firstly, the particle information in the system is loaded, including the boundary condition (i.e., periodic in this work), the particle coordinate locations, the particle size, and other parameters (i.e., Euler angles). Then, sampling points/lines/triangles are generated, and the points are randomly or uniformly generated in the system; consequently, the end points of the lines and triangles are generated using the Monte Carlo method following the periodic conditions mentioned above. Then, we execute kernels in parallel with the GPU; when a kernel is launched, a grid of threads that are organized in a 3D hierarchy are generated, each gird is organized into an array of thread blocks, and each block can contain up to 1024 threads. Then, we allocate threads in advance according to the complexity of the computation; for example, for 1-or 2-point correlation functions, a 3D block and 3D grid are more suitable because the computing cores (threads) number in the millions, but for the 3-point correlation function, a 2D block and 2D grid are appropriate. However, in this work we use the 3D block and 3D grid to allocate the threads, because they can be degraded to the 2D block and 2D grid. This information is then input into the GPU kernels and the sampling points/lines/triangles are allocated to threads cores to execute judgement of the location and classification, the 1-point correlation function yields the result directly, and S (i) 2 (r) and S (i) 3 (r 1 , r 2 , θ) will proceed to the next sample line/triangle, until all the patterns of the sample are calculated. Finally, all the discrete data are stored. To show the speedup ratio of the GPU to the CPU, an ordinary computer is used, i.e., the CPU is an AMD RYZEN 7 3800× and the "GPU" is NVIDIA GeForce RTX 2070. The results are listed in Table 1.
locations, the particle size, and other parameters (i.e., Euler angles). Then, sampling points/lines/triangles are generated, and the points are randomly or uniformly generated in the system; consequently, the end points of the lines and triangles are generated using the Monte Carlo method following the periodic conditions mentioned above. Then, we execute kernels in parallel with the GPU; when a kernel is launched, a grid of threads that are organized in a 3D hierarchy are generated, each gird is organized into an array of thread blocks, and each block can contain up to 1024 threads. Then, we allocate threads in advance according to the complexity of the computation; for example, for 1-or 2-point correlation functions, a 3D block and 3D grid are more suitable because the computing cores (threads) number in the millions, but for the 3-point correlation function, a 2D block and 2D grid are appropriate. However, in this work we use the 3D block and 3D grid to allocate the threads, because they can be degraded to the 2D block and 2D grid. This information is then input into the GPU kernels and the sampling points/lines/triangles are allocated to threads cores to execute judgement of the location and classification, the 1point correlation function yields the result directly, and and , , will proceed to the next sample line/triangle, until all the patterns of the sample are calculated. Finally, all the discrete data are stored. To show the speedup ratio of the GPU to the CPU, an ordinary computer is used, i.e., the CPU is an AMD RYZEN 7 3800× and the "GPU" is NVIDIA GeForce RTX 2070. The results are listed in Table 1.   Note: packing system and discrete number of samplings are closely related to the execution time. The parameters in this test are shown as follows: the packing fraction φ a = 0.5 and total particle number N p = 954; the sampling points/lines/triangles of S a /S aa /S aaa : N sampling = 2 × 10 6 ; the discrete values of S aa /S aaa are N r = 100 and N r 1 × N r 2 × N θ = 100 × 100 × 100, respectively.

Verification
Berryman [27] derived an analytical expression for the n-point correlation functions of phase i using a fully penetrable sphere-packing model, which is one of the few packing systems for which the n-point correlation functions are defined analytically.
For 2-point correlation function of phase i, where ρ is the number density of the spheres, and V n represents the union volume of n spheres, when n = 2, where R is the sphere radius. Figure 5a shows the comparison of the analytical and estimated 2-point correlation function of the matrix phase for the penetrable sphere packing. It is evident that the algorithm is able to calculate the 2-point correlation function with a high accuracy. In addition, the 2-point correlation function holds the relationships S mm (0) = 0.5 = φ m and lim r→∞ S mm (r) = 0.25 = φ 2 m , thereby satisfying the limiting value of S (i) 2 (r). Figure 5b shows the results of the ζ m of the matrix phase for penetrable sphere packing and said results compared with other results in the literature [20,21]. This verified our algorithm in penetrable packing systems. As the packing of an impenetrable sphere model has been studied extensively [16][17][18][19][20], to better verify our algorithm, an impenetrable packing system of a monodisperse sphere was generated using a DEM-based method [28] to ensure the statistical isotropy; after a convergence study, the container side was set to be 20 times the length of the particle diameter (i.e., L/D = 20), and the results of the impenetrable sphere packing system are shown in Figure 6a. Figure 6b is also shown to ensure that there is not any overlap between the particles. Impenetrable packing spheres with a packing fraction of φ a = 0.1 − 0.6 were calculated using our approach, and the results were compared with existing ones, as shown in Figure 6c. An excellent agreement is shown between our results and the latest results by Gillman [20] and the simulation results by Chan and Torquato [15]. Therefore, the results of both penetrable and impenetrable packing systems verified our approach.

Results and Discussion
With the 3-point microstructural parameter accurately estimated, the effective properties of thermal conductivity, , and the diffusion coefficient can be calculated for almost any packing system. In this work, super-ellipsoid was employed due to its morphological diversity, it is a joint name used to define a variety of 3D regular-shaped particles including ellipsoids, cuboids, etc. In mathematics, it can be expressed by adding a deformation coefficient to the standard formula of an ellipsoid, where , , and denote the three semi-axes in coordinates. Table 2 displays a fraction of the super-ellipsoid family with different values of the aspect ratio ⁄ and deformation coefficient , in which and 1. By changing the aspect ratio and deformation coefficient, we can study the effect of particle shape in two directions.

Results and Discussion
With the 3-point microstructural parameter ζ accurately estimated, the effective properties of thermal conductivity, κ e , and the diffusion coefficient D e can be calculated for almost any packing system. In this work, super-ellipsoid was employed due to its morphological diversity, it is a joint name used to define a variety of 3D regular-shaped particles including ellipsoids, cuboids, etc. In mathematics, it can be expressed by adding a deformation coefficient m to the standard formula of an ellipsoid, where a, b, and c denote the three semi-axes in coordinates. Table 2 displays a fraction of the super-ellipsoid family with different values of the aspect ratio a/b and deformation coefficient m, in which b = c and m ≥ 1. By changing the aspect ratio and deformation coefficient, we can study the effect of particle shape in two directions.   Figure 7a presents the third-order approximations of thermal conductivity using Equation (18) for different packing fractions. First, the effect of the aspect ratio is explored. It can be observed that as the packing fracture increases, the effect of the particle shape becomes more pronounced; for example, the normalized effective thermal conductivity ( ⁄ ) for spheres ( ) and ellipsoids ( ⁄ 2.5, 1, ) at packing fraction 0.1 are basically the same, while at the packing fraction 0.5, the effective thermal conductivity for an ellipsoid ( ⁄ 2.5, 1, ) is almost two times as high as the value for spheres. In general, "long" particles (i.e., a high aspect ratio ellipsoid) lead to a greater thermo conductivity. These results are also compared to the H-S lower bound,   Figure 7a presents the third-order approximations of thermal conductivity using Equation (18) for different packing fractions. First, the effect of the aspect ratio is explored. It can be observed that as the packing fracture increases, the effect of the particle shape becomes more pronounced; for example, the normalized effective thermal conductivity ( ⁄ ) for spheres ( ) and ellipsoids ( ⁄ 2.5, 1, ) at packing fraction 0.1 are basically the same, while at the packing fraction 0.5, the effective thermal conductivity for an ellipsoid ( ⁄ 2.5, 1, ) is almost two times as high as the value for spheres. In general, "long" particles (i.e., a high aspect ratio ellipsoid) lead   Figure 7a presents the third-order approximations of thermal conductivity using Equation (18) for different packing fractions. First, the effect of the aspect ratio is explored. It can be observed that as the packing fracture increases, the effect of the particle shape becomes more pronounced; for example, the normalized effective thermal conductivity ( ⁄ ) for spheres ( ) and ellipsoids ( ⁄ 2.5, 1, ) at packing fraction 0.1 are basically the same, while at the packing fraction 0.5, the effective thermal conductivity for an ellipsoid ( ⁄ 2.5, 1, ) is almost two times as high as the value for spheres. In general, "long" particles (i.e., a high aspect ratio ellipsoid) lead   Figure 7a presents the third-order approximations of thermal conductivity using Equation (18) for different packing fractions. First, the effect of the aspect ratio is explored. It can be observed that as the packing fracture increases, the effect of the particle shape becomes more pronounced; for example, the normalized effective thermal conductivity ( ⁄ ) for spheres ( ) and ellipsoids ( ⁄ 2.5, 1, ) at packing fraction 0.1 are basically the same, while at the packing fraction 0.5, the effective thermal conductivity for an ellipsoid ( ⁄ 2.5, 1, ) is almost two times as high as the value for spheres. In general, "long" particles (i.e., a high aspect ratio ellipsoid) lead   Figure 7a presents the third-order approximations of thermal conductivity using Equation (18) for different packing fractions. First, the effect of the aspect ratio is explored. It can be observed that as the packing fracture increases, the effect of the particle shape becomes more pronounced; for example, the normalized effective thermal conductivity ( ⁄ ) for spheres ( ) and ellipsoids ( ⁄ 2.5, 1, ) at packing fraction 0.1 are basically the same, while at the packing fraction 0.5, the effective thermal conductivity for an ellipsoid ( ⁄ 2.5, 1, ) is almost two times as high as the value for spheres. In general, "long" particles (i.e., a high aspect ratio ellipsoid) lead to a greater thermo conductivity. These results are also compared to the H-S lower bound, which is unable to characterize any shape effect. Then, we studied the effect of the defor-   Figure 7a presents the third-order approximations of thermal conductivity using Equation (18) for different packing fractions. First, the effect of the aspect ratio is explored. It can be observed that as the packing fracture increases, the effect of the particle shape becomes more pronounced; for example, the normalized effective thermal conductivity ( ⁄ ) for spheres ( ) and ellipsoids ( ⁄ 2.5, 1, ) at packing fraction 0.1 are basically the same, while at the packing fraction 0.5, the effective thermal conductivity for an ellipsoid ( ⁄ 2.5, 1, ) is almost two times as high as the value for spheres. In general, "long" particles (i.e., a high aspect ratio ellipsoid) lead to a greater thermo conductivity. These results are also compared to the H-S lower bound, which is unable to characterize any shape effect. Then, we studied the effect of the defor-   Figure 7a presents the third-order approximations of thermal conductivity using Equation (18) for different packing fractions. First, the effect of the aspect ratio is explored. It can be observed that as the packing fracture increases, the effect of the particle shape becomes more pronounced; for example, the normalized effective thermal conductivity ( ⁄ ) for spheres ( ) and ellipsoids ( ⁄ 2.5, 1, ) at packing fraction 0.1 are basically the same, while at the packing fraction 0.5, the effective thermal conductivity for an ellipsoid ( ⁄ 2.5, 1, ) is almost two times as high as the value for spheres. In general, "long" particles (i.e., a high aspect ratio ellipsoid) lead to a greater thermo conductivity. These results are also compared to the H-S lower bound, which is unable to characterize any shape effect. Then, we studied the effect of the defor- Figure 7a presents the third-order approximations of thermal conductivity κ e using Equation (18) for different packing fractions. First, the effect of the aspect ratio is explored. It can be observed that as the packing fracture φ a increases, the effect of the particle shape becomes more pronounced; for example, the normalized effective thermal conductivity (κ e /κ m ) for spheres (a = b = c) and ellipsoids (a/b = 2.5, m = 1, b = c) at packing fraction φ a = 0.1 are basically the same, while at the packing fraction φ a = 0.5, the effective thermal conductivity for an ellipsoid (a/b = 2.5, m = 1, b = c) is almost two times as high as the value for spheres. In general, "long" particles (i.e., a high aspect ratio ellipsoid) lead to a greater thermo conductivity. These results are also compared to the H-S lower bound, which is unable to characterize any shape effect. Then, we studied the effect of the deformation coefficient on the thermo conductivity; when m is changes from one to infinity, the super-ellipsoid (i.e., a = b = c) transforms from a sphere into a regular hexahedron, as shown in Figure 7b. With the increase of the packing fraction, the shape effect becomes more and more obvious, but the overall magnitude is less pronounced.

Conclusion
In this work, third-order effective properties including thermal conductivity and the diffusion coefficient were estimated with a high accuracy and efficiency. The GPU-based parallel algorithm was developed to calculate the 3-point parameter. Verifications were conducted for both the penetrable and impenetrable sphere packing systems. The results of the 2-point correlation function and the 3-point parameter were verified by analytical expressions and recent simulation results. The effective thermo conductivity of the impenetrable sphere packing system was compared to the existing results and showed a good agreement.
With our GPU-based parallel algorithm, we calculated third-order approximations of thermal conductivity and a diffusion coefficient for a differently shaped super-ellipsoids packing system with a different packing fraction. A significant particle shape effect was predicted for both the thermal conductivity and diffusion coefficient, especially when increasing the aspect ratio of the particle. In addition, a clear shape effect could be observed when the particle transformed from a sphere to a regular hexahedron, but the magnitude was less pronounced. With these techniques, effective properties for a wide range of particle-packing systems can be predicted with high accuracy and speed. The third-order approximations of the diffusion coefficient D e were calculated for different packing systems using equation (22). Figure 7c shows the normalized D e as a function of the packing fraction. Unlike the effective thermal conductivity, the effective diffusion coefficient in all the packing systems decreases with the increase of the packing fraction. This is because the increasing particle phase will obstruct the diffusion. For the aspect ratio's effect, a clear shape effect was also discovered, which becomes most obvious when the packing fraction φ a is around 0.3-0.4. Overall, the packing of the particles with a high aspect ratio will more rapidly decrease the D e . These results are also compared with the H-S lower bounds; it should be noted that some packing systems (i.e., a/b = 2.5, a/b = 2.0, b = c) in a range of packing fractions (i.e., φ a = 0 ∼ 0.3) will exceed the H-S lower bounds. Figure 7d represents the effect of the deformation coefficient; it shows a similar tendency, but with a smaller magnitude.

Conclusions
In this work, third-order effective properties including thermal conductivity and the diffusion coefficient were estimated with a high accuracy and efficiency. The GPUbased parallel algorithm was developed to calculate the 3-point parameter. Verifications were conducted for both the penetrable and impenetrable sphere packing systems. The results of the 2-point correlation function and the 3-point parameter were verified by analytical expressions and recent simulation results. The effective thermo conductivity of the impenetrable sphere packing system was compared to the existing results and showed a good agreement.
With our GPU-based parallel algorithm, we calculated third-order approximations of thermal conductivity and a diffusion coefficient for a differently shaped super-ellipsoids packing system with a different packing fraction. A significant particle shape effect was predicted for both the thermal conductivity and diffusion coefficient, especially when increasing the aspect ratio of the particle. In addition, a clear shape effect could be observed when the particle transformed from a sphere to a regular hexahedron, but the magnitude was less pronounced. With these techniques, effective properties for a wide range of particle-packing systems can be predicted with high accuracy and speed.