Identify the Micro-Parameters for Optimized Discrete Element Models of Granular Materials in Two Dimensions Using Hexagonal Close-Packed Structures

The widely used simple cubic-centered (SCC) model structure has limitations in handling diagonal loading and accurately representing Poisson’s ratio. Therefore, the objective of this study is to develop a set of modeling procedures for granular material discrete element models (DEM) with high efficiency, low cost, reliable accuracy, and wide application. The new modeling procedures use coarse aggregate templates from an aggregate database to improve simulation accuracy and use geometry information from the random generation method to create virtual specimens. The hexagonal close-packed (HCP) structure, which has advantages in simulating shear failure and Poisson’s ratio, was employed instead of the SCC structure. The corresponding mechanical calculation for contact micro-parameters was then derived and verified through simple stiffness/bond tests and complete indirect tensile (IDT) tests of a set of asphalt mixture specimens. The results showed that (1) a new set of modeling procedures using the hexagonal close-packed (HCP) structure was proposed and was proved to be effective, (2) micro-parameters of the DEM models were transit form material macro-parameters based on a set of equations that were derived based on basic configuration and mechanism of discrete element theories, and (3) that the results from IDT tests prove that the new approach to determining model micro-parameters based on mechanical calculation is reliable. This new approach may enable a wider and deeper application of the HCP structure DEM models in the research of granular material.


Introduction
The origin of the discrete element method (DEM) can be traced back to the late 1970s, when it was developed by Cundall and Strack to address the complexities associated with granular materials, owing to their inherently discrete nature [1]. Then, the DEM was introduced into the modeling of the asphalt mixture. Buttlar and You utilized it to model and examine the workings and efficiency of asphalt materials [2,3]. Another study simulated the viscoelastic behavior of asphalt mixture using micromechanical parameters obtained from a dynamic shear rheometer in the Simple Performance Test (SPT) [4]. A 3D microstructure-based Discrete Element Method (DEM) model was created by combining multiple 2D models and then was used to calculate the stress-strain behavior during repeated loading conditions [5]. The results of a laboratory test indicated that the 3D model had better agreement with the test results than the associated 2D models [6]. The contact models form the foundational mechanism in the DEM, and its micro-parameter determination is a critical part. Researchers developed several approaches to relate the micro-parameters with the macro-parameters of the asphalt mixture. For the compacted asphalt mixture at room temperature, the dynamic modulus test was used to determine the viscoelastic parameters of the Burgers model [7]. The dynamic modulus could be used to reflect the stress and strain response by specific load directly. The dynamic modulus test is conducted at temperatures of −10 • C, 10 • C, 21 • C, 37 • C, and 54 • C at loading frequencies of 0.1 Hz, 0.5 Hz, 1 Hz, 5 Hz, 10 Hz, and 25 Hz at each temperature and is specified in AASHTO T342. The simulation results were in agreement with the results obtained from laboratory tests. The creep test was used to determine the viscoelastic parameters in models that were based on microstructure [8]. The dynamic shear rheometer test [9] and the constant strain rate uniaxial compression test were conducted to calculate the time-dependent contact stiffness of the Burgers model [10,11]. Similar contact models and parameter calculations were used in the prediction of the mechanical properties of asphalt [12]. The internal forces configuration of the asphalt mixture was evaluated through the use of established DEM models. The parameter determination for the samples in the compaction process is more difficult than the compacted samples due to the high flowability of asphalt at high temperatures. Chen, Huang et al. proposed an indirect approach to predict the viscoelastic parameters at high temperatures [13]. In their study, serval dynamic modulus tests were performed at low temperatures, and the nonlinear regression analysis was used to obtain the mater curve of the asphalt mixture. Then, the viscoelastic parameters at high temperatures were predicted through the asphalt mixture master curve [14,15].
Another crucial aspect of the DEM model is the use of modeling techniques. You and Buttlar were pioneers in utilizing image processing to construct 2D model structures. They utilized grayscale images obtained through optical scanning to establish microfabric DEM models [2,3]. The simple performance test employed similar models to forecast the dynamic modulus and phase angles of the asphalt mixture [4]. The concept of constructing a 3D model by stacking 2D models was developed based on the 2D DEM models [5]. Then, the development of 3D DEM models utilizing ball/clump elements as basic building blocks was carried out. Ball elements are widely adopted in 3D DEM models for simulating the performance of asphalt mixture due to its simplicity and clear visual aid. For example, the uniaxial compression test was simulated by ball-based DEM models [10]. The modified model was demonstrated as having the capability of simulating creep tests. For the purpose of simplifying the DEM models for asphalt mixtures, researchers typically treat the mixture as a two-phase material composed of coarse aggregates and asphalt mastic [16]. The ball elements were also utilized to represent the asphalt mastic, which is composed of fine aggregate, fines, and asphalt. The ball-based models have an obvious disadvantage due to their inability to represent the irregular shapes of aggregates. New modeling approaches were developed to model the aggregates with more realistic shapes. A proposal was made to utilize randomly generated irregular particles to visualize and simulate the micro-scale properties of the asphalt mixture under mechanical loading [17]. The investigation into the effect of aggregate shape on the diffusivity of asphalt mastic utilized random packing models of ellipsoidal and convex polyhedral particles [18]. Additionally, researchers have utilized realistic aggregate shapes to improve the accuracy of asphalt mixture simulations. Techniques such as X-ray CT and image processing were employed to generate DEM models featuring realistic aggregate shapes [19,20]. A more precise method, the individual aggregate reconstruction technology, was proposed to establish DEM models for asphalt mixture [21,22]. Fracture behavior in asphalt concrete laboratory specimens is able to bridge a vital link in the design of asphalt concrete paving mixtures and pavement structures. A two-dimensional particle flow software package (PFC-2D) was used to study the complex crack behavior observed in asphalt concrete fracture tests [23]. A computer simulation using the discrete element method (DEM) is presented in order to understand and visualize how crushing initiates and develops inside a simulated pavement structure [24,25]. Yu et al. [26] studied the effect of aggregate size distribution and angularity distribution on dynamic modulus using a 3D discrete element method (DEM).
The above research primarily focused on the properties of the compacted asphalt mixture. The compaction process of asphalt mixture is characterized by frequent and intense material movement and changes in contact force. As a result, the contact models and modeling techniques for this process are distinct from those used for compacted asphalt mixtures. There are limited studies that focus on the compaction of asphalt mixture. Wang et al. compared the fundamental mechanics of asphalt compaction using both FEM and DEM and emphasized that DEM can simulate aggregate translation and rotation [27]. The DEM models have also been shown to provide valuable theoretical support for intelligent compaction. Chen and Huang et al. utilized the Burgers model to simulate the compaction of asphalt mixture using DEM [28]. In a subsequent study, they simulated gyratory compaction, vibration compaction, and kneading compaction using an open-source code [13]. Gong et al. established shape-based DEM models to simulate Superpave gyratory compactor (SGC) tests, which simulate the field compaction process of asphalt mixture, and reported agreement between the results of laboratory compaction tests and simulation results [29,30]. This research introduced realistic aggregate shapes into the DEM simulation. However, the established models had limitations on the total number of elements, compaction dynamics, and parameters determination.
Yu Liu et al. [7] introduced a set of theoretical calculations for DEM models using the Burgers contact model and cubic-centered cubic (SCC) ball array structures. This theoretical calculation approach has been proved reliable in predicting the dynamic modulus of asphalt mixture. In order to extend the use of DEM to a variety of performance tests, the hexagonal close-packed (HCP) ball array was used as the basic model structure. It is clear that the HCP structure is more complicated than the SCC structure and has a different mechanical structure. Therefore, a new set of theoretical calculations is needed.
This study aims to establish a hexagonal close-packed discrete element model for granular material with the ability to transit diagonal loading and performance based on Poisson's ratio. To achieve this objective, first, new generation procedures of the HCP ball array for the 2D DEM models were proposed; second, a theatrical approach that was used to transition from material macro-properties to contact micro-parameters values in the 2D DEM models were derived based on the basic configuration and mechanism of discrete element theories; and third, the contact stiffness, bond strength, and an-isotropic properties were discussed and verified by comparing IDT results between designed 2D DEM models and laboratory tests.

Model and Methods
The research methodology of this study is shown in Figure 1, where the generation method of the new HCP model is introduced in the modeling procedures of the hexagonal close-packed generation method section. The process involves several steps, starting with the generation of nonoverlap clumps obtained through scanning aggregates of varying grain sizes with a 3D scanner, followed by grain-size expansion of clumps where the clump sizes are progressively increased to reach their intended dimensions. Next, hexagonal close-packed (HCP) balls are generated, and these balls are then grouped based on clump geometry, where the classification of an HCP ball as either coarse aggregate or asphalt mastic depends on whether its center position falls within a clump. Finally, the installation of contact properties is carried out using the contact-bond model, which has been shown to effectively simulate the fracture behavior of asphalt mixtures in previous studies. The corresponding mechanical calculation for contact micro-parameters was then derived and verified through simple stiffness/bond tests in the 2D DEM model and verified with theoretical values. Finally, an indirect tensile (IDT) test in the 3D DEM modeling generated by the HCP model structure and laboratory test results is compared. In previous studies, the simple cubic-centered (SCC) ball array was widely used as the basic model structure, as shown in Figure 2a. This kind of model is proved to be effective in the simulation and in the prediction of the dynamic modulus of the asphalt mixture. The SCC structure can transit loading in vertical and horizontal directions efficiently. However, it lacks the ability to transit diagonal loading and performance based on Poisson's ratio. To make up for this disadvantage, the hexagonal close-packed (HCP) ball array (see Figure 2b) was used. This study aims to carry out a reliable approach to determine the micro-parameters in the DEM models. Based on the basic configuration and mechanism of discrete element theories, the transition from material macro-properties to contact micro-parameters was derived. The contact stiffness, bond strength, and an-isotropic properties were discussed and verified by designed DEM models.

Case of SCC
The basic mechanical unit of SCC can be described as a single contact with two balls (Figure 3). This unit can be treated as a single-spring system. The contact force ( ) and stress ( ) can be expressed as Equations (1) and (2). In previous studies, the simple cubic-centered (SCC) ball array was widely used as the basic model structure, as shown in Figure 2a. This kind of model is proved to be effective in the simulation and in the prediction of the dynamic modulus of the asphalt mixture. The SCC structure can transit loading in vertical and horizontal directions efficiently. However, it lacks the ability to transit diagonal loading and performance based on Poisson's ratio. To make up for this disadvantage, the hexagonal close-packed (HCP) ball array (see Figure 2b) was used. In previous studies, the simple cubic-centered (SCC) ball array was widely used as the basic model structure, as shown in Figure 2a. This kind of model is proved to be effective in the simulation and in the prediction of the dynamic modulus of the asphalt mixture. The SCC structure can transit loading in vertical and horizontal directions efficiently. However, it lacks the ability to transit diagonal loading and performance based on Poisson's ratio. To make up for this disadvantage, the hexagonal close-packed (HCP) ball array (see Figure 2b) was used. This study aims to carry out a reliable approach to determine the micro-parameters in the DEM models. Based on the basic configuration and mechanism of discrete element theories, the transition from material macro-properties to contact micro-parameters was derived. The contact stiffness, bond strength, and an-isotropic properties were discussed and verified by designed DEM models.

Case of SCC
The basic mechanical unit of SCC can be described as a single contact with two balls ( Figure 3). This unit can be treated as a single-spring system. The contact force ( ) and stress ( ) can be expressed as Equations (1) and (2). This study aims to carry out a reliable approach to determine the micro-parameters in the DEM models. Based on the basic configuration and mechanism of discrete element theories, the transition from material macro-properties to contact micro-parameters was derived. The contact stiffness, bond strength, and an-isotropic properties were discussed and verified by designed DEM models.

Case of SCC
The basic mechanical unit of SCC can be described as a single contact with two balls ( Figure 3). This unit can be treated as a single-spring system. The contact force (F) and stress (σ) can be expressed as Equations (1) and (2).
where, k n is stiffness of the contact, δ is the displacement at the contact, E is the material modulus, and ε is strain.
where, is stiffness of the contact, is the displacement at the contact, is the material modulus, and is strain.
In regard to unit dimensions, the stress and strain at contact can also be expressed as: where is the area of the contact plane and is the length of contact. Submit Equations (2)-(4) into Equation (1) where is the radius of SCC balls. One important point to mention is that the third dimension L is hidden in the calculation of S.
Eventually, the material modulus (E) of an SCC array can be expressed by contact stiffness (kn):

Case of HCP
The basic unit of HCP is the combination of three closed contact balls, as seen in Figure 4. This unit can be treated as a simple truss system. The contact force and displacement are the combinations of the vertical portion of ′ and ′: In regard to unit dimensions, the stress and strain at contact can also be expressed as: where S is the area of the contact plane and L is the length of contact.
where R is the radius of SCC balls. One important point to mention is that the third dimension L is hidden in the calculation of S. Eventually, the material modulus (E) of an SCC array can be expressed by contact stiffness (k n ):

Case of HCP
The basic unit of HCP is the combination of three closed contact balls, as seen in Figure 4. This unit can be treated as a simple truss system. The contact force F and displacement δ are the combinations of the vertical portion of F and δ : As the angle θ of the truss equals 30 degrees, the contact force F can be expressed as: The length of the truss system is calculated as the vertical portion of the connection between the two balls: As the angle of the truss equals 30 degrees, the contact force can be expressed as: The length of the truss system is calculated as the vertical portion of the connection between the two balls: Eventually, the material modulus (E) of a unidimensional hex array can be related to contact stiffness (kn) as: This section may be divided by subheadings. It should provide a concise and precise description of the experimental results and their interpretation as well as the experimental conclusions that can be drawn.

Validation Example
A set of DEM models was used to verify the reliability of Equations (6) and (12). The ratio of model height versus model width was set as 2.0. Due to hardware and computational power limitations, the validation model dimension is constrained to a limited size, which is worth noting. Consider the impact of model sizes, as shown in Figure 5, where four groups of models were tested with scales ranging from 5 × 10, 10 × 20, and 20 × 40 to 40 × 80. Eventually, the material modulus (E) of a unidimensional hex array can be related to contact stiffness (k n ) as: This section may be divided by subheadings. It should provide a concise and precise description of the experimental results and their interpretation as well as the experimental conclusions that can be drawn.

Validation Example
A set of DEM models was used to verify the reliability of Equations (6) and (12). The ratio of model height versus model width was set as 2.0. Due to hardware and computational power limitations, the validation model dimension is constrained to a limited size, which is worth noting. Consider the impact of model sizes, as shown in Figure 5, where four groups of models were tested with scales ranging from 5 × 10, 10 × 20, and 20 × 40 to 40 × 80.  The contact stiffness was set as 1 × 10 5 for all the tested models. All boundaries were rigid and confined. Vertical displacements of 1% model height per second were applied on the top plane, according to Equations (6) and (22). The theoretical material moduli should be 100 kPa and 129.8 kPa. The obtained material moduli from DEM model are shown in Figure 6. The contact stiffness was set as 1 × 10 5 for all the tested models. All boundaries were rigid and confined. Vertical displacements of 1% model height per second were applied on the top plane, according to Equations (6) and (22). The theoretical material moduli should be 100 kPa and 129.8 kPa. The obtained material moduli from DEM model are shown in Figure 6.
The model scale has a significant impact on the obtained material moduli from the DEM model. HCP was more sensitive to the model scale than was SCC. The obtained material moduli of the SCC group from the DEM model were close to the 100 kPa theoretical value. The HCP group reached 96.22% of the theoretical value (124.9/129.8 kPa) when using the 60 × 120 configuration.
The contact stiffness was set as 1 × 10 5 for all the tested models. All boundaries were rigid and confined. Vertical displacements of 1% model height per second were applied on the top plane, according to Equations (6) and (22). The theoretical material moduli should be 100 kPa and 129.8 kPa. The obtained material moduli from DEM model are shown in Figure 6. The model scale has a significant impact on the obtained material moduli from the DEM model. HCP was more sensitive to the model scale than was SCC. The obtained material moduli of the SCC group from the DEM model were close to the 100 kPa theoretical value. The HCP group reached 96.22% of the theoretical value (124.9/129.8 kPa) when using the 60 × 120 configuration.

Case of SCC
The bonding condition makes no difference to the SCC ball arrays, since the bonding plane is perpendicular to the vertical direction.

Case of SCC
The bonding condition makes no difference to the SCC ball arrays, since the bonding plane is perpendicular to the vertical direction.

Case of HCP
The bonding condition makes no difference to the SCC ball arrays, since the bonding plane is perpendicular to the vertical direction. The bonding plane in the HCP ball arrays has an angle of θ degrees in the horizontal direction ( Figure 7). Thus, the shear force at the bonding plane contributes a vertical component when applied to vertical loading. The bonding condition makes no difference to the SCC ball arrays, since the bonding plane is perpendicular to the vertical direction. The bonding plane in the HCP ball arrays has an angle of θ degrees in the horizontal direction ( Figure 7). Thus, the shear force at the bonding plane contributes a vertical component when applied to vertical loading. The shear force at contact plane can be expressed as: where is the shear stiffness at bonding plane, is the displacement in shear direction, is the friction coefficient, and is the normal contact force at bonding plane. Then, the total force in vertical direction equals Submit Equations (13) and (14) into (15): The material modulus before slip then equals The shear force at contact plane can be expressed as: where k s is the shear stiffness at bonding plane, δ s is the displacement in shear direction, µ is the friction coefficient, and F n is the normal contact force at bonding plane. Then, the total force in vertical direction equals Submit Equations (13) and (14) into (15): The material modulus before slip then equals

Case of HCP
The configuration of the validation example is the same as the models used in the previous section. Contact stiffness (k n ) without bonding k n and k s were set as 1 × 10 5 . The theoretical value according to Equation (17) was 173.2 kPa. The simulation results are shown in Figure 8. Model scale also has significant impacts on the material moduli of DEM models. As the model scale increased, the obtained material moduli from the DEM model were closer to the theoretical value. When using the 60 × 120 configuration, the obtained material moduli from the DEM model reached 97.29% (168.5/173.2) of the theoretical value. Considering the hardware calculation efficiency and the scale effects influence, a 40 × 80 configuration was used in this study.

Case of SCC
The tensile bond strength can be expressed as: In the case of the unidimensional ball array, refer to Equation (5):

Case of HCP
The tensile bond strength can be expressed as: Model scale also has significant impacts on the material moduli of DEM models. As the model scale increased, the obtained material moduli from the DEM model were closer to the theoretical value. When using the 60 × 120 configuration, the obtained material moduli from the DEM model reached 97.29% (168.5/173.2) of the theoretical value. Considering the hardware calculation efficiency and the scale effects influence, a 40 × 80 configuration was used in this study.

Case of SCC
The tensile bond strength T F can be expressed as: In the case of the unidimensional ball array, refer to Equation (5):

Case of HCP
The tensile bond strength T F can be expressed as: The relationship of contact stiffness in the contact interface is: where k n is the aggregate-mastic interface normal contact stiffness, k n1 is the normal contact stiffness of aggregate, and k n2 is the normal contact stiffness of mastic; k s used the same method.

Validation Example
Setting T F = 0.5 N and R =0.01 m, according to Equations (19) and (21), the theoretical values of SCC and HCP models are 25 Pa and 43.3 Pa, respectively. The obtained value from DEM model are 24.7 Pa (98.8% of theoretical value) and 38.37 Pa (88.6% of theoretical value), respectively.

An-Isotropic of Hexagonal Close-Packed (HCP) Structures
The vertical direction and horizontal direction of the HCP ball array are different. The an-isotropic properties can be written as:  (19) and (21)

An-Isotropic of Hexagonal Close-Packed (HCP) Structures
The vertical direction and horizontal direction of the HCP ball array are different. The an-isotropic properties can be written as: in which, ϕ is among 0-30 degree, θ equals to 30 degrees. In the case of material modulus, the modulus constant in the horizontal direction is | | = √ , and in the vertical direction it is | | = √ . Then, the an-isotropic properties of the material are plotted in Figure 9.

Modeling Procedures of Modified Random Generation Method with Realistic Coarse Aggregate Shapes
To enhance the efficiency of the model, a modified random generation method with realistic coarse aggregate shapes was introduced in this study. The new method employed

Modeling Procedures of Modified Random Generation Method with Realistic Coarse Aggregate Shapes
To enhance the efficiency of the model, a modified random generation method with realistic coarse aggregate shapes was introduced in this study. The new method employed the cross-section of 3D models as the 2D model geometries rather than the directed generation of 2D models. The new modeling procedures are described in the following steps: •

Generation of Nonoverlap Clumps
The clump geometries were obtained through scanning aggregates of varying grain sizes using a 3D scanner. The methods for generating clumps and determining grain size have been described in prior studies [21,31]. The clump grain sizes were determined based on the mixture design. The clumps were generated within a 100 × 63 mm cylinder container, with each clump being generated at 70% of its intended grain size to ensure successful generation. The coarse aggregates of grain sizes G2, G3, and G4 were generated in succession, as depicted in Figure 10(1). It is worth mentioning that the coarse aggregate can also be directly introduced via a compacted model through the compaction process. calculated by iterating through the entire clump set, this method would add unnecessary computational strain to the computer. As an alternative, the maximum overlap ratio could be estimated by monitoring the leading contact force. To limit clump movement and enhance efficiency, a high damping ratio of 0.7 was assigned to all clumps.

•
Generation of Hexagonal Close-Packed (HCP) Balls The two lattice structures that result in the highest density for equal-diameter ball arrangements are the cubic-centered cubic (SCC) and the hexagonal close-packed (HCP). In this study, the hexagonal arrangement was chosen. For ease of ball labeling, the balls were generated within a cubic space, and then any balls outside the cylindrical container boundary were removed, as illustrated in Figure 10(3).

•
Grouping HCP Balls Based on Clump Geometry (Objective Search Efficiency Improved Algorism) The classification of an HCP ball into either coarse aggregate or asphalt mastic depends on whether its center position falls within a clump. The number of HCP balls representing rubber particles and voids was determined based on the mixture design. These two groups of HCP balls were then randomly selected from within the mastic. The final grouping results are displayed in Figure 10(4). The most time-consuming step in this section is the objective search of overlap detection, which determines the group properties of HCP balls. Thus, here we proposed an improved objective algorism. The original algorism needs to loop the ball list and clump list (pebble list) from beginning to end, as shown in Figure 11. The required steps for a model with 88,489 balls and 142,266 pebbles are 12.6 billion.  •

Grain-Size Expansion of Clumps
The clump grain sizes were increased until they reached their target dimensions, as depicted in Figure 10(2). The expansion procedure involved several iterations to prevent excessive overlapping in a single expansion. In the example shown, the clump grain sizes were expanded 10 times with an expansion factor of approximately 1.03631121 for each step, calculated as (1/0.7) (1/10). To minimize overlap between clumps, the model was run until the maximum overlap ratio dropped below 0.1%. While the overlap ratio could be calculated by iterating through the entire clump set, this method would add unnecessary computational strain to the computer. As an alternative, the maximum overlap ratio could be estimated by monitoring the leading contact force. To limit clump movement and enhance efficiency, a high damping ratio of 0.7 was assigned to all clumps.

•
Generation of Hexagonal Close-Packed (HCP) Balls The two lattice structures that result in the highest density for equal-diameter ball arrangements are the cubic-centered cubic (SCC) and the hexagonal close-packed (HCP). In this study, the hexagonal arrangement was chosen. For ease of ball labeling, the balls were generated within a cubic space, and then any balls outside the cylindrical container boundary were removed, as illustrated in Figure 10(3).

•
Grouping HCP Balls Based on Clump Geometry (Objective Search Efficiency Improved Algorism) The classification of an HCP ball into either coarse aggregate or asphalt mastic depends on whether its center position falls within a clump. The number of HCP balls representing rubber particles and voids was determined based on the mixture design. These two groups of HCP balls were then randomly selected from within the mastic. The final grouping results are displayed in Figure 10(4). The most time-consuming step in this section is the objective search of overlap detection, which determines the group properties of HCP balls. Thus, here we proposed an improved objective algorism. The original algorism needs to loop the ball list and clump list (pebble list) from beginning to end, as shown in Figure 11. The required steps for a model with 88,489 balls and 142,266 pebbles are 12.6 billion. The improved algorism decreases the calculation steps by narrowing down the search area, as shown in Figure 12.
• First, find the location and diameter of the current pebble. • Second, calculate the extended coverage area. • Third, determine if the ball is within the pebble area.
By estimation, the improved objective algorism requires 21.4 million steps to finish the calculation, which saves about 98.9% of calculation time. •

Installation of Contact Properties
The contact-bond model was selected because it has been demonstrated to effectively The improved algorism decreases the calculation steps by narrowing down the search area, as shown in Figure 12. The improved algorism decreases the calculation steps by narrowing down the search area, as shown in Figure 12.

•
First, find the location and diameter of the current pebble. • Second, calculate the extended coverage area. • Third, determine if the ball is within the pebble area.
By estimation, the improved objective algorism requires 21.4 million steps to finish the calculation, which saves about 98.9% of calculation time. •

Installation of Contact Properties
The contact-bond model was selected because it has been demonstrated to effectively simulate the fracture behavior of asphalt mixtures. Although nearly all aggregates and rubber particles are covered by asphalt, there is bond strength between directly connected aggregates and rubber particles. Furthermore, the linear contact model was designated as the default model for all subsequent contacts (following fracture), and the contact properties would be derived from the parent particles. By estimation, the improved objective algorism requires 21.4 million steps to finish the calculation, which saves about 98.9% of calculation time.

Installation of Contact Properties
The contact-bond model was selected because it has been demonstrated to effectively simulate the fracture behavior of asphalt mixtures. Although nearly all aggregates and rubber particles are covered by asphalt, there is bond strength between directly connected aggregates and rubber particles. Furthermore, the linear contact model was designated as the default model for all subsequent contacts (following fracture), and the contact properties would be derived from the parent particles.

Validation Example with Indirect Tension (IDT) Tests
The IDT test is an effective method for evaluating the low-temperature cracking performance of asphalt mixture [32]. This section designed a group of indirect tension (IDT) tests in laboratory to verify the reliability of the proposed mechanical parameters transition.

Mixture Design and DEM models
Three mixture designs were selected; see Table 1. The IDT test setup is shown in Figure 13. The test speed is 50 mm/s, and the load and displacements during the test are recorded and compared with the DEM model.

Mixture Design and DEM models
Three mixture designs were selected; see Table 1. The IDT test setup is shown in Figure 13. The test speed is 50 mm/s, and the load and displacements during the test are recorded and compared with the DEM model.  The models with three mixture designs are shown in Figure 14.

Calculation of Model Micro-Parameters
The micro-parameters were calculated based on the equations derived from the previous section, and the contact model parameters were calculated based on materials' macro-properties [7,[33][34][35], as shown in Table 2.

IDT Results and Discussion
The IDT results from laboratory tests and DEM simulation are shown in Figure 15. As shown in Figure 15a, a total of 9 laboratory specimens belonging to 3 groups were tested. For Mix1, Mix2, and Mix3, the average peak forces were 6.25 kN, 10.12 kN, and

Calculation of Model Micro-Parameters
The micro-parameters were calculated based on the equations derived from the previous section, and the contact model parameters were calculated based on materials' macro-properties [7,[33][34][35], as shown in Table 2.

IDT Results and Discussion
The IDT results from laboratory tests and DEM simulation are shown in Figure 15. As shown in Figure 15a, a total of 9 laboratory specimens belonging to 3 groups were tested. For Mix1, Mix2, and Mix3, the average peak forces were 6.25 kN, 10.12 kN, and 10.81 kN, respectively. In general, the mixture type has the largest coarse aggregate grain size (Mxi3) and presented the highest peak force (tensile strength). Accordingly, Mix3 showed the steepest increasing rate (material moduli) and minimum displacement at the force peak (ultimate strain). After specimen failure, the decrease curves were relatively gentle compared to the increase curve. In this test, the standard deviations caused by the differences in specimens and test errors were relatively large but were still in a reasonable range.   Figure 15b shows the results of the DEM simulation for Mix1, Mix2, and Mix3. The average peak forces were 6.21 kN, 9.96 kN, and 10.56 kN, respectively. The results were close to that of laboratory tests, with relative errors ranging from 0.64% to 2.31%. The standard deviations of IDT results were at the same level as the comparison group. With zigzag data point curves, the results were much "rougher" than that of the laboratory control group. This is caused by the limited model scale. Actually, in the authors' other studies there are models with millions of elements (less than 10,000 elements were used in this study) that could present more smooth curves, especially in 2D. The authors' intention is to showcase the ability and reliability of their models by utilizing limited scales. After specimen failure, sharp drops were observed. There were two major reasons. First, the limited model scale led to large jumps at each of the failures, and second, the 2D models had less freedom than reality in which cracks could develop in lateral directions.
To compare the results, one force/displacement curve (whose test value is in the middle) for each mixture type was chosen, as shown in Figure 15c. The peak value of DEM simulations causes more displacement than the in lab because the initial loading stage of the DEM needs a process to "compact" the model into a denser status to achieve better loading transfer efficiency. However, prior to specimen failure, the peaking value and other parts of loading curves have good consistency with lab results. Considering that all the parameters are based on theoretical calculation without adjustment and with limitations on minimum element size, the DEM simulation results are reasonable, and the parameter calculation is reliable.
Compared to other models that utilize the discrete element method, which often exhibit a relative error exceeding 10%, the error in the results of this study is relatively small. It is noteworthy that the parameters in this study are derived using formulas rather than iteratively fitting them based on simulation results, as is commonly practiced in general studies. Additionally, this study employs a minimum-cost two-dimensional model with limited dimensions and scale, and the calculation of contact parameters is based on laboratory experiments that may have some degree of error fluctuation. Given these constraints, simulating loading curves of three different graded mixtures with consistent trends is still a challenging task, despite slight differences in the curves.
The crux of this paper lies in utilizing a theoretical calculation method to derive contact parameters for the discrete element method when simulating particulate matter instead of relying on iterative back-calculation fitting, as in typical research. Starting from a theoretical level, this method is more logical, resource-efficient, and reproducible.

Conclusions
In this study, a new approach for modeling procedures and determining parameters was proposed to enhance the integration of discrete element models (DEMs) in asphalt simulation. The following conclusions were drawn: (1) A new approach for modeling procedures using the hexagonal close-packed (HCP) structure was proposed. This method, which employs realistic coarse aggregate morphology from 3D scanning, was demonstrated to be effective and can help save time and resources by reducing the need for laboratory samples. An objective searchefficiency improvement algorism is developed in this process. (2) Micro-parameters of the DEM models were transformed from material macro-parameters using a set of equations that were derived based on the basic configuration and mechanism of discrete element theories. The effectiveness of the DEM models in simulating the indirect tensile strength for asphalt mixtures was demonstrated. The results were close to that of laboratory tests, with relative errors ranging from 0.64% to 2.31%.
The key contribution of this research is the use of a reliable approach for determining model micro-parameters through mechanical calculation instead of a radical and inefficient iteration method for model parameter fitting. This new approach has the potential to expand and deepen the application of HCP structure DEM models in granular material research.