Calculation Method for Grading Size and Grading Bounding Box of Virtual Aggregate Based on DEM

+e grading size, which can be defined as the side length of the smallest square hole through which an aggregate with irregular shape can directly pass, is an important morphology parameter and can be used to calculate the gradation of mixture material.+e grading bounding box, which can be defined as the circumscribed cuboid with central axis length being the grading size, is an important visual tool for observing the size and direction of the aggregate. Virtual test to calculate the grading size of a virtual aggregate is environmentally friendly and efficient, but the result provided by current research is imprecise and the grading bounding box is also rarely mentioned. In this paper, the multilevel complete projection algorithm is proposed to precisely calculate the grading size of a virtual aggregate.+e whole process of the algorithm can be expressed by formula after the operation of sphere discretization by converting the virtual aggregate shell into the discrete aggregate. +en, the discrete aggregate is projected onto a complete series of the plane to form several 2D figures, and then, each 2D figure is projected onto a complete series of the orthogonal axis to form orthogonal segments. +e grading size can finally be obtained by comparing the length of the above orthogonal segments based on the key central axis length principle. +e influencing factors of computational accuracy and efficiency are considered in the algorithm. Finally, the grading bounding box can be built by using the Rodrigues transformation according to the information obtained from the above algorithm.


Introduction
e parameter of aggregate gradation in mixture material, such as asphalt mixture [1], is an important indicator that affects its whole mechanical properties. e sieving size of a single aggregate is an important factor in calculating the gradation of its mixture, which generally refers to the minimum standard side length of the square-hole sieve that can be passed.
A 3D virtual aggregate shell can be simulated by professional model software, such as 3D modeling engine ACIS and Mimics [2][3][4], by combining a series of equidistant 2D images obtained by photographing or X-ray CT tomography. e grading size of the virtual aggregate shell is difficult to be directly obtained because of its irregular shape. ere are many achievements in finding the grading size of virtual aggregates shell in various papers, which are roughly classified as follows: (1) Virtual screening experiment: based on the screening principle, a large number of 3D virtual aggregates are put into the virtual screen model for the virtual screening experiment. Several sieve models, such as multilayer square-hole sieve [5], banana screen [6], horizontal linear screen [7], and double-amplitude screen with variable aperture [8], have been used. e method needs to simulate the sieve model first, and then, the vibration effect is achieved by applying velocity to the sieve model, which is time-consuming (2) Simple measuring method: based on the vernier caliper principle, the virtual aggregates are measured one by one, taking the point from the surface and calculating circularly. e skeleton model is a way to find the point [9,10]. e farthest two points are the long axis, perpendicular to which the farthest two points are the central axis, which is the final size of virtual aggregate [11]. But the long axis of a rectangle aggregate in this method is its diagonal, clearly out of a common sense (3) Graphic fitting method: based on the constant eigenvalue principle, a virtual aggregate is converted into a regular geometric figure, whose characteristic size can be taken as the size of aggregate: the central axis of the fitting ellipsoid or the short axis of the fitting ellipse [12][13][14], the diameter of the minimum circumscribed or maximum inscribed sphere (circle) [15], the diameter of the sphere (circle) based on the equal volume (area) conversion, and the least second moment of figure [16] (4) Grading bounding box method: based on the key central axis length principle [17], the central side length of the grading bounding box is chosen as the grading size. e grading bounding box [2,3] is searched by dynamically changing the three control angles of a circumscribed cuboid of the virtual aggregate, usually rotating a single angle in turn, fixing the other two angles, and making the cuboid reach the optimal value (the minimum central axis length or the minimum volume of the real-time bounding box). In addition, the minimum Ferrette diameter method [18] for determining the size of 2D figures is also included in this method e former three methods have the drawbacks of resultimprecise, especially the wrong principle of the simple measuring method. e grading bounding box method is a relatively precise manner with a complicated searching path, which is hard to understand and easy to miss the best value. In fact, the grading bounding box is better concerning the result, not the cause of the grading size calculation.

Objectives and Scopes
After synthesizing and sublating the above algorithms for measuring or calculating the grading size of virtual aggregate, and based on the discrete element method and the critical central axis size principle in a practical screening experiment, a multilevel complete projection algorithm is proposed.
Firstly, in order to effectively describe the irregular virtual aggregate shell, the method of sphere discretization is used to convert it into a discrete aggregate, using enough reasonable small particles ranked in a certain order to fill in its inner space.
Secondly, all the particles in one discrete aggregate are projected onto a complete series of verification planes to form 2D projection figures, and then, each 2D projection figure is projected on a complete series of orthogonal axes in the verification plane to form a group of orthogonal projection segments. e grading size of discrete aggregate can be calculated accurately by judging the length relations of all segments. To enforce this rule, necessary regulatory means, such as reasonable distribution of verification planes and their inner coordinate system, should be prepared in advance.
irdly, the influencing factors, such as the sphere discretization method, inclusion standard, diameter of particles, the number of particles, and self-defined angle increment, should be considered to balance the precision and efficiency of calculation.
Lastly, the grading bounding box is constructed to judge the direction and shape of the virtual aggregate according to the finalized controlling angles of the verification plane and the obliquity of its inner orthogonal axes, in which the projection got the final grading size.

Sphere Discretization
Any virtual aggregate can be expressed as a closed outer surface and its internal space. e sphere particles used in the discrete element method (DEM) have accurate characteristics of location and size, which is more motorized than the mesh used in the finite element method (FEM) [19]. Based on this, the paper proposes filling the internal space of the virtual aggregate to be measured with small-diameter particles, forming the discrete aggregate model, and then, the grading size can be calculated according to all inner particles that do not make a difference. e calculation process is also formulated and flowered due to the characteristics of particles and can be applied to all kinds of nondiscrete virtual aggregates (geometric shells) mentioned in the introduction. e discrete aggregate model obtained by the discrete element method can be directly used for subsequent calculation. e aggregate shell is introduced into the DEM software to simulate the size and shape of it with a large number of particles arranged according to certain rules, as shown in Figure 1. e bubble pack method [18] uses a geometric shell to determine the area of the Delaunay triangulation. e cubic method and hexagonal method [20,21] use the geometric shell to cover uniform particles which are not overlapped but arranged in cubic or hexagonal order, and the geometric shell does not participate in the operation. e filling method uses a geometric shell as the boundary of the wall, inside which uniform or uneven particles are generated and balanced. e contour-filling method [22] uses a geometric shell to define the center of the uniform particles to be generated.

Multilevel Complete Projection Algorithm
e grading size of a discrete aggregate model, which is composed of several particles, can be calculated by a multilevel complete projection algorithm.
e principles and procedures are as follows.

Selection of Verification Plane.
e projection of a single discrete aggregate in different verification planes is different; all possibilities should be tested, so it is called the completeness of verification plane selection. Besides, it has the same projection of discrete aggregate on parallel planes with the same or opposite unit normal vector, so the direction of the plane is the only factor to determine the shape and size of the projection. Based on spherical coordinates, all verification planes can be described as the planes passing through a fixed point O(x 0 , y 0 , z 0 ) with its unit normal vector n varying with the angle ϕ and θ, as shown in Figure 2: n � x n , y n , z n � (sin ϕ cos θ, sin ϕ sin θ, cos ϕ), ϕ ∈ 0, 90 ∘ , θ ∈ 0, 360 ∘ ). (1) Given the discrete features of the verification plane, in order to avoid tedious calculation and to ensure the valid result, the verification plane is selected equally spaced: the angle ϕ is selected in every Δ 1 (degree), and the angle θ in every Δ 2 (ϕ) (degree) on the latitudinal circumference determined by ϕ, as shown in Figure 3. Self-defined Δ 1 evenly divides the 1/4 meridian circle of the imaginary sphere whose radius is R ′ into several equidistant arc L ⌢ , so on the latitudinal circle whose radius is R ′ (ϕ) determined by every angle ϕ, the verification plane (ϕ, θ) is taken from θ � 0°at intervals of the equidistant arc L ⌢ , which corresponds to Δ 2 (ϕ) as follows: Particularly, when k � 0, the unique verification plane that is parallel to xOy plane is taken.

Selection of 2D Coordinate
System. It is necessary for the subsequent calculation to establish an appropriate 2D coordinate system in a plane, which passes the point O(x 0 , y 0 , z 0 ) and whose unit normal vector is n � (x n , y n , z n ), describing the point on it.
Let the only line in the plane, which is not only passing through the point O but also parallel to the xOy plane, be the plane's x-axis, which can be written as x 2D and whose direction vector is L 1 � (x L 1 , y L 1 , z L 1 ). Let the only line in the plane, which is not only passing through the point O but also perpendicular to x 2D , be the plane's y-axis, which can be written as y 2D and whose direction vector is L 2 � (x L 2 , y L 2 , z L 2 ). Particularly, choose L 1 � (1, 0, 0) and L 2 � (0, 1, 0) for the verification plane parallel to the xOy plane. e following set of equations can be listed based on the geometric properties:

Mathematical Problems in Engineering
(1) x 2D ⊥n: (3) (3) x 2D � � � �xOy plane: (4) y 2D ⊥n: (6) y 2D ⊥x 2D : By simplifying, the direction vectors of the x-axis and yaxis in the plane, whose unit normal vector is n � (x n , y n , z n ), can be described in the 3D coordinate system: For a certain verification plane (ϕ i , θ j ) chosen in the above section, the inner coordinate system x 2D Oy 2D (ϕ i , θ j ) can be described as 4.3. Calculation of Projection Coordinates. Each particle in the discrete aggregate to be measured is presented for the calculation of the projection coordinates, as shown in Figure 4. Let the coordinates of a particle's center in the 3D coordinate system be A(x A , y A , z A ). en, the values of the coordinates (L 1,A , L 2,A ) of point B, which is the projection point of A in the verification plane (ϕ i , θ j ), can be calculated by the geometric vector rule. As shown in Figure 5, in the established coordinate system x 2D Oy 2D (ϕ i , θ j ), the orthogonal axes L 3 and L 4 for verification are built with the obliquity angle ω k � (L 3 , L 1 ) ∈ [0, 90 ∘ ). e linear equations of the two axes are y 2D � x 2D · tan ω k and y 2D � x 2D · (−1/tan ω k ), respectively. e values of the coordinates (L 3,ωk,A , L 4,ωk,A ) of a certain particle's center A in the coordinate system L 3 OL 4 (ω k ) determined by a specific obliquity angle of ω k is transformed from its projection coordinate values (L 1,A , L 2,A ) in the coordinate system x 2D Oy 2D :

Grading Size Calculation of Projection.
After transforming all the coordinate values of particles in the discrete aggregate on current ω k in current (ϕ i , θ j ), the range length L ωk,L3 or 4 in the orthogonal axes L 3 or L 4 can be found by firstly finding the coordinate values reached most upwards and downwards, respectively, in the axes L 3 and L 4 , and then, taking a subtraction with the consideration of particle's radius, as the bold line in axes shown in Figure 5, e side length L ωk(ϕi,θj) of the smallest square, through which the projection of the discrete aggregate on the obliquity angle of ω k in the verification plane of (ϕ i , θ j ) can pass, can be found by choosing the maximum range length: Figure 4: Plane projection principle. Figure 5: Orthogonal axis projection principle. 4 Mathematical Problems in Engineering where Am and An stand for arbitrary qualified particles inside the discrete aggregate, and r Am and r An stand for the radius of the corresponding particle.

Multilevel Complete Projection.
e specific obliquity angle of ω k is generalized to take the value in the range ω ∈ [0, 90 ∘ ) and calculate all L ωk(ϕi,θj) values in every selfdefined Δ 3 in the specific verification plane (ϕ i , θ j ). en, the side length of the smallest square L ϕi,θj , through which the projection of the discrete aggregate in the specific verification plane of (ϕ i , θ j ) can pass, is calculated: e specific verification plane (ϕ i , θ j ) is generalized to take the value in the range ϕ ∈ [0, 90 ∘ ], θ ∈ [0, 360 ∘ ), and calculate all L ϕi,θj values in every self-defined Δ 1 and Δ 2 (ϕ).
en, the side length of the smallest square L, through which the discrete aggregate can pass, is calculated: e calculation flowchart of grading size for a discrete aggregate is shown in Figure 6.

Influencing Factors of Calculation Results and Efficiency
An example is shown in Figure 7; a geometric shell of the octahedron is surrounded by eight random planes. e unit normal vector and the distance from the specified point inside the polyhedron to each plane direction are shown in Table 1. e sphere discretization of the aggregate shell is represented by different conditions: (1) Sphere discretization method: cubic method (using ball generate cubic . . . command) or filling method (using ball distribute porosity 0.4 . . . command and cycle 800 calm 100 command), as shown in Figure 8 (2) Inclusion standard (only for cubic method): center inclusion (as long as the sphere center of the particle is inside the polyhedral shell, it can be included) or whole inclusion (as long as the whole particle is located in the polyhedral shell, it can be included) (3) e diameter of particles: 0.6 mm, 1.0 mm, or 1.5 mm By combining, the conditions of discrete aggregate after deleting the wall are formed, as shown in Figure 9. e grading sizes of the above discrete aggregates are calculated by taking Δ 1 � 5°and Δ 3 � 5°, as shown in Table 2.
e smaller the diameter of the particle is, the more accurate the final calculation results are. Condition (a) and condition (e) can be used as more accurate models for the corresponding discretization method. e difference of grading size between center inclusion (condition (b)) and sphere inclusion (condition (c)) is about the diameter of the particle, the former is closer to the value of the condition (a), and all the angle variables are consistent, so sphere center inclusion is a relatively reasonable manner.
When the particle diameter is fixed, the measurement value of the filling method is slightly less than that of the cubic method, which is also because the center inclusion is chosen so that part of the particle is out of the plane, while the wall of the filling method blocks the particles reaching out to some extent, which is shown in Figure 8. Taking into account the fact that the particle number of the cubic method is slightly less and does not need to run a balance, it can be used in the rough calculation.
When the discretization method is fixed, the measurement precision decreases with the increase of diameter. When the diameter reaches 1.5 mm, the measurement value of the filling method (condition (g)) decreases, but that of the cubic method (condition (d)) obviously increases, which is also due to the blocking effect of the wall in the filling method, while the particle diameter has a great influence on the cubic method. Besides, there is a deviation from the angle variable in both methods.
In the rough calculation of grading size, the combination of the cubic method and 1 mm particle diameter can be used with accepted error. In fact, selecting 10% of the estimated size value as the diameter can obtain good results, especially when there are subsequent mechanical calculations, while the shape of the original polyhedron shell is required.
In the fine calculation of grading size, the filling method should be adopted and the smaller diameter should be selected, but this will prolong the calculation time, which is proportional to the following conditions: (1) e number of particles (num p ): the bigger the diameter of particles, the less the total number of particles, but the rougher the shape it will eventually take, so the diameter of particles cannot be increased casually. Given that the inner particles contribute almost nothing to the final result, all the particles whose distance from the geometric shell is less than the diameter of particles can be grouped separately and the others can be deleted. e number of particles whether the entire or grouped is shown in Table 3. (2) e number of verification planes (num v ): num v is inversely proportional to Δ 2 1 . When Δ 1 is taken as 1°, 3°, and 5°, num v is calculated as 20852, 2367, and 869, respectively.
(3) e number of orthogonal axes on each verification plane (num a ): num a is inversely proportional to the increment Δ 3 . When Δ 3 is taken as 1°, 3°, and 5°, num a is calculated as 90, 30, and 18, respectively.
Taking condition (a) and condition (e) as an example, the discrete aggregate is run by taking different angle increments Δ 1 and Δ 3 , as shown in Tables 4 and 5.
Mathematical Problems in Engineering e results still confirm that the filling method is more accurate than the cubic method. Besides, the smaller the angle increment is, the more accurate the measured value is most of the time, but more computation is needed. Due to the discrete nature of particles, occasionally prominent or concave particles at a certain angle make measurement results not smooth, as occurs in Δ 1 � 3°and Δ 3 � 3°in the filling method. e maximum difference of grading size in the two methods is, respectively, 0.031 mm and 0.075 mm, which is much less than the particle diameter, so taking the two angle increment values as 5°can be accepted most of the time for less running time.
For the convenience of observation, the 2D projection of the discrete aggregate in their corresponding verification plane in the increment of Δ 1 � 1°and Δ 3 � 1°is shown in Figure 10.

Start End
Global initialization L = +∞ Grading size L The coordinate system L1 and L2 are built in new (ϕ, θ) The orthogonal axes L3 and L4 are built at new (ω) in current (ϕ, θ)
New particle (radius r) selected to calculate the projection coordinates in the L10L2 coordinate system (ϕ, θ), then converted to coordinates in the L30L4 coordinate system (ϕ, θ, ω)    Table 2: Grading size of discrete aggregate in different sphere discretization conditions.

Grading Bounding Box
e grading bounding box of a discrete aggregate can be defined as the cuboid, whose six surfaces are both circumscribed to the discrete aggregate without any part going beyond, and which has the same grading size as the discrete aggregate. Let the three side lengths of a grading bounding box be H ≥ L ≥ W in descending order; then, its central axis length L is necessarily the well-calculated grading size of its inner discrete aggregate. e grading bounding box can be used to simplify the judgment of the direction and size of the inner discrete aggregate. (1) e normal vector direction of the plane in which the two edges LW are located represents the orientation of the discrete aggregate. (2) When W/H < 0.4, the discrete aggregate can be judged as needle-like aggregate (L << H) or flake-like aggregate (L >> W).
Let a discrete aggregate find its grading size L on the obliquity of ω c in the verification plane (ϕ a , θ b ), the short axis length and the long axis direction of its grading bounding box can be calculated by equations (17) and (18), and the short axis direction, the central axis direction, and the long axis length can be obtained by applying the Rodríguez transform: n LW � ± 1 × sin ϕ a cos θ b , sin ϕ a sin θ b , cos ϕ a . (18)

Rodríguez Transform.
e Rodríguez transform [23][24][25] constructs a cross-product matrix K and rotation matrix P, and then, the vector n 1 is rotated around the axis of k to the vector n 2 by the degree of δ according to the right-hand rule, as shown in Figure 11: where k � (k x , k y , k z ) T is a unit vector and I is a third-order unit matrix. As shown in Figure 12, when L 1 and L 2 in equations (10) and (11) are rotated around n LW by degree of ω c according to the Rodríguez transform, L 3 and L 4 can be obtained in the Table 4: Grading size of discrete aggregate using the cubic method in different angle increments of projection. able 5: Grading size of discrete aggregate using the filling method in different angle increments of projection.   Mathematical Problems in Engineering 3D coordinate system by taking k � n LW , δ � ω c , and n 1 � L 1 and L 2 , respectively: (20) After being unitized, L 3 and L 4 can be taken as the unit normal vector of the plane in which the two edges HW or HL are located: if the central axis length L is obtained on the axis L 4 , then the short axis direction n HL � ± 1 × L 3 /|L 3 |, the central axis direction n HW � ± 1 × L 4 /|L 4 |, and vice versa. e long axis length can be obtained by projecting algorithm, choosing the higher value range in the axes, either n HL or n HW , combined with n LW , while the lower one is the known L or W for verification.

Application.
According to the process, the grading bounding box of two discrete aggregate models mentioned above can be constructed: three controlling angles ϕ � 30°, θ � 104°, and ω � 5°in the cubic method and ϕ � 33°, θ � 102.82°, and ω � 7°in the filling method are used to obtain two groups of grading bounding box parameters, as shown in Table 6.
With the help of the view tool, a targeted observation can be made: Figures 13 and 14 show the position relationship between the grading bounding box and the discrete aggregate is observed from the original angle of view and from the long axis direction, respectively. e former objectively reflects the whole orientation of the discrete aggregate in the space, and the latter is corresponding to the 2D projection image obtained in Figure 10, which verifies the rationality of the whole process. Figure 15 shows a three-view diagram of the polyhedral shell enclosed in its grading bounding box obtained by the    filling method. Because the discrete aggregate as the medium inherits the shape, size, and orientation of the polyhedral shell to a large extent, the grading bounding box is basically suitable for the polyhedral shell.

Summary and Conclusions
(1) Sphere discretization is an effective means to tackle a virtual aggregate with irregular shapes. After transferring the whole geometric shell into discrete aggregate, which consists of lots of particles ranging in a certain order, the grading size and the grading bounding box can be calculated by tracing each particle's exact position. e discrete aggregate inherits the shape, size, and direction of the original virtual aggregate, and the whole process can be formulated and flowed (2) Multilevel complete projection algorithm can be used to calculate the grading size of discrete aggregate with a relatively accurate result. After a complete projection of the three dimensions, the grading size can be obtained by easily comparing the lengths of a series of line segments. Meanwhile, the following additional results are achieved: equal arc angle equation and the method of establishing coordinate system on any plane (3) e filling method is more accurate to calculate the grading size than the cubic method but less efficient. A more accurate result can also be obtained by decreasing the particle diameter and the angle increment Δ 1 and Δ 3 , but both of them are time-consuming. When there is a follow-up mechanical calculation or low precision requirement, it is recommended that the particle diameter be 10% of the estimated size with a relatively precise shape and fewer particles (4) Grading bounding box can be built to find the long axis direction and axis length of the virtual aggregate by using the data generated from the above algorithm and by realizing the vector rotation according to the Rodríguez transform

Data Availability
Some or all data, models, or codes that support the findings of this study are available from the corresponding author upon reasonable request.

Disclosure
Any opinions, findings, and conclusions or recommendations expressed in this material are those of the writers and do not necessarily reflect the views of the National Science Foundation of China and the Hunan Transportation Science and Technology Program.

Conflicts of Interest
e authors declare that they have no conflicts of interest.