Experimental and numerical analyses of crushing resistance of unbound road materials

ABSTRACT

Aggregate breakage in unbound pavement layers can lead to pavement distresses that affect their functionality and service life. Thus understanding the mechanics and clarifying the factors affecting materials breakage resistance are important for ensuring adequate performance of these layers. In this study, aggregate breakage in unbound granular materials (UGM) is investigated experimentally and numerically. Experimentally, aggregate breakage under uniaxial compression is examined for two UGMs prepared with the same aggregate type but different gradations. To capture the experimentally observed influence of gradation and load magnitude on aggregate breakage, a Discrete Element Method (DEM) model was developed, based on granular mechanics particle contact and failure laws. A simple procedure to identify the contact and failure law parameters from experiments is proposed. With those parameters, the model’s capability of capturing the effect of gradation and loading on the aggregate breakage in UGM is evaluated. Based on comparison with experimental findings, it is shown that the model can capture macro-scale properties of UGM, such as its deformation response under uniaxial compression, as well as the amount of aggregate breakage in the material.

KEYWORDS:

• Unbound granular materials (UGM)
• aggregate breakage
• discrete element method (DEM)
• gradation
• contact mechanics

1 Introduction

Aggregate degradation through crushing and abrasion may profoundly affect unbound layers performance and, accordingly, pavements service life and its functional properties. As reported by Saeed et al. (2001) aggregate degradation in unbound layers may contribute to rutting, fatigue and surface unevenness of the asphalt pavement layers on top. Thus understanding the mechanics of aggregate damage in unbound road materials and layers is crucial for optimal material selection, pavements structural design and performance prediction. The present study aims to investigate this important topic by examining experimentally and numerically aggregate crushing in unbound road materials.

To characterise crushing resistance of aggregates, numerous empirical methods have been proposed, such as Los Angeles (LA) abrasion test, aggregate impact test and aggregate value tests, as described in Bjarnason et al. (2000) and Saeed et al. (2001). These test methods allow to measure empirical indices capturing the strength of the aggregates and provide thus a practical way for ranking the materials with respect to their crushing resistance and establishing minimum acceptability requirements. As being pointed out in several studies, however, the relationship between the empirical indices and the aggregates performance in the field is not necessarily strong, cf. e.g. Wu et al. (1998). Furthermore, the empirical test methods may not capture the effects of loading conditions and material parameters, such as gradation and aggregate toughness, on crushing resistance of unbound road materials.

To establish a quantitative link between crushing performance of unbound road layers and the parameters and loading conditions of unbound granular material (UGM) an analysis framework based on granular mechanics is needed (Cole and Peters 2008, Celma Cervera et al. 2017). Discrete element method (DEM) is a particularly promising tool in this context, as it allows linking the local deformations and damage arising at stone-to-stone contact points with the materials performance on the macro scale, cf. Cole and Peters (2008) and Celma Cervera et al. (2017). In a number of recent studies, DEM has been applied to model macro-mechanical behaviour of road materials, e.g. Chen et al. (2020), Wang et al. (2021) and railway ballast materials, e.g. McDowell and Lu (2010), Huang and Tutumluer (2011). Triaxial response of UGMs used in railway ballast was studied numerically by McDowell and Lu (2010), using 3D DEM, with particular emphasise on the influence of particle interlock on the macro-scale response of the material. McDowel and Lu (2010) showed that, provided that adequate contact and particle breakage law parameters are chosen, the model can capture the materials response under static and cyclic triaxial tests reasonably well and may be used to gain micromechanical insight into the material behaviour. Furthermore, Huang and Tutumluer (2011) showed that adjusting contact law parameters in DEM modelling allows to capture the effect of fine content on ballast performance in shear box tests. Chen et al. (2020) employed DEM to study particle segregation in UGMs and uncompacted asphalt mixtures at vibratory loading. They showed that DEM can accurately capture the effects of vibratory loading magnitude, bitumen content and aggregate gradation on the amount of segregation induced in the material. DEM models developed by Olsson et al. (2019b) and Wang et al. (2021) were shown to capture successfully the asphalt mixture densification process during compaction in the laboratory (Olsson et al. 2019b) and in the field (Wang et al. 2021).

DEM modelling of aggregate crushing in UGMs has also received some attention in the literature. Bono and McDowell (2016) evaluated several DEM particle breakage criteria with respect to their applicability for evaluating silica sand particle breakage in uniaxial compression tests on single size UGM. They showed that breakage criteria based on maximum contact force and maximum octahedral shear stress resulted in simulated macroscopic behaviour in agreement with the experiments. Bisht and Das (2021) used 2D DEM to evaluate the effect of assumed particle shape evolution after breakage on the resulting macroscopic response of UGM. They showed that characteristics of particle fragmentation may have a profound effect on the stress–strain response of the material at macro scale.

The existing DEM studies of aggregate crushing rely on contact and particle failure laws, calibrated with the experimental observations of UGM behaviour on the macro level. The empirical nature of those particle interaction laws precludes quantitative analysis of the influence of UGM’s parameters, such as gradation, aggregate stiffness and toughness, on the materials performance, as the empirical particle interaction laws may not capture accurately the aggregate crushing under conditions, which are significantly different from those used for calibration of models. This issue is of particular importance with respect to the implementation of novel materials for unbound granular layers, such as, for example, using marginal types for UGMs. The intention of this study is to address this gap, by establishing a new DEM-based model for aggregate crushing in UGMs, relying on granular mechanics-based contact and damage laws.

In this study, a DEM model is developed and implemented in commercial DEM software PFC3D^TM (Itasca Consulting Group Inc. 2019) to investigate aggregate crushing in UGMs under uniaxial monotonic compression. The new model builds up on a modelling approach developed by Olsson et al. (2019a, 2019b) which incorporates particle contact and damage laws based on granular mechanics. Accordingly, aggregate contact and failure laws proposed by Olsson and Jelagin (2019a, 2019b) based on elastic contact mechanics and the experimental results by Celma Cervera et al. (2017) are implemented in PFC3D and incorporated into the model. In this study, the existing model approach is further extended with a simplified aggregate breakage simulation allowing to account for the effect of aggregate breakage on contact force redistribution in the material. The developed model is used to investigate the effects of UGM gradations and load levels on aggregate crushing. The aggregate crushing under uniaxial monotonic compression is also examined experimentally for two gradations and three maximum applied loads. The capability of the model to quantitatively capture the effects of gradation and load magnitude on aggregate crushing is validated through comparison of experimental and numerical findings. Feasibility of using the proposed DEM model for UGMs material design optimisation and for gaining further insight into mechanical aggregate degradation in unbound road layers is discussed.

2. Problem formulation and methodology

In UGMs, the external forces are transferred through the network of aggregate contacts. In addition to aggregate properties, aggregate propensity to breakage in unbound layers will thus be controlled by the magnitude and distribution of contact forces at stone-to-stone contact points. Aggregates in unbound layers of asphalt pavements have to withstand a large number of load cycles during service life as well as repeated high loading during construction. The magnitude of pressure applied to unbound layers may be quite high, reaching approximately 1.6 MPa in case of heavy construction traffic travelling directly on unbound layers, cf. Bjarnason et al. (2000). Furthermore, high-pressure pulses are applied to UGM during compaction (Holtz and Kovacs 1981).

To obtain the best computational insight into the aggregate crushing process, models should explicitly account for stone-to-stone interactions, particle rearrangements and address stone breakage as a function of fundamental materials properties. The DEM developed by Cundall and Strack (1979) is a promising tool for investigating this problem, as it allows to capture the relationship between the macroscopic load experienced by UGM and the local forces arising at aggregate contact points. Accordingly, the main goal of the present paper is to develop a modelling framework for evaluating aggregate crushing in UGMs based on DEM. This study aims to establish a general analysis framework for evaluating aggregate breakage at arbitrary loading conditions. As a first step towards this aim, the aggregate crushing is evaluated experimentally and numerically in uniaxial compression test. The DEM model for uniaxial compression test is implemented in a commercial software, PFC3D^TM. The user subroutines for the contact law and the damage law, detailed in following subsections, are implemented into DEM software through C++ and fish codes respectively.

The main analysis steps of this study are summarised in Figure 1. As may be seen there, an UGM, characterised by its gradation, porosity and aggregate properties is idealised in the DEM model as collection of spherical particles, with their size distribution following the UGM gradation. The external load is applied incrementally in DEM, and at each increment the normal and tangential contact forces are calculated for each particle, following the contact laws detailed in Section 2.2 and used as an input to particle failure laws for evaluating aggregate breakage by employing the particle failure law as detailed in Section 2.3. Once, a given particle breaks, the effect of its fragmentation on the macroscopic response of UGM is accounted for through a fragmentation model, cf. Section 2.3. In the following, the main elements of the developed modelling procedure are summarised in detail. The materials and methods used in the experiments are presented in Section 3, while the DEM model of the uniaxial compression experiment is presented in Section 4.

Figure 1. Flowchart of the DEM modelling framework (pink particles represents fractured particles).

2.1. DEM modelling

In DEM, the individual aggregates are represented as spherical particles for computational efficiency. The rotational degrees of freedom of these particles are constrained to simulate aggregate interlock. The model geometry and coordinate system used in this study are shown in Figure 2. The UGM is placed into a stationary rigid cylinder with one open-end and a rigid horizontal moving loading head is applied on top of the specimen. The models in Figure 2(a, b) are illustrated for the two gradations investigated: a model consisting of single size particles as shown in Figure 2(a) and a model consisting of densely graded particles as shown in Figure 2(b).

Figure 2. (a) Single size gradation model and (b) denser gradation model.

The movement of the particles is governed by explicit numerical integration of Newton’s laws of motion. The forces acting on the particles are gravity and contact forces arising from particles interaction with other particles or the rigid wall elements. To capture the physics of the system, an accurate model for the contact behaviour is of utmost importance and will be discussed in the following section. In addition, local damping, $\zeta$, as implemented by PFC3D^TM is used, so that the model can reach steady state quickly during the initial UGM packing due to gravity. It should be noted that it affects the acceleration of the particles and does not affect the particles when they are in steady state. After the particles fall to the bottom of the confining cylinder and compressed by the loading head, they are constrained and will not be able to accelerate. Hence, the contact forces between the particles will not be affected by the damping. The effect of the local damping is given as (1) $F-F_{\zeta }=\mathrm{ma}$(1) where $F$ is the force applied on the particle, $m$ is the mass of the particle and $a$ is the acceleration of the particle. $F_{\zeta }$ is the damping force and shown as (2) $F_{\zeta }=\mathrm{\zeta F}$(2)

2.2. Contact law

The local contact behaviour between stone particles was investigated by Olsson et al. (2019a) based on spherical indentation experimental results of Cervera et al. (2017). To transfer the results from the indentation experiments where a rigid spherical indenter is pressed into a flat piece of rock material to the problem with two spherical particles locally in contact, the effective contact radius: (3) $R_0={\left({\displaystyle \frac{1}{R_1}+{\displaystyle \frac{1}{R_2}}}\right)}^{-1}$(3) is introduced with quantities defined in Figure 3.

Figure 3. Illustration of two spheres in contact and variables to define contact –overlap relation.

The normal contact force model is divided into three different parts to account for different physical processes observed experimentally during loading. Initially up to an indentation depth, $h_s$, a linear response is seen experimentally, interpreted as crushing of surface asperities, cf. Cervera et al. (2017). This relation can be stated as (4) $F_N={\displaystyle \frac{F_s}{h_s}h}$(4) where $h$ is the indentation depth and $F_s$ is the normal force at the indentation depth $h_s$. In the indentation experiments by Celma Cervera et al (2017), $F_s=100\mathrm{N}$ and $h_s=0.02\mathrm{mm}$ were found using an indenter with a radius of 6.25 mm. This value is scaled presently using the effective contact radius $R_0$ by (5) $F_s=100\mathrm{N}\sqrt{{\displaystyle \frac{R_0}{6.25\mathrm{mm}}}}$(5) To make continuous transition, the force is scaled with $\sqrt{R_0}$. When the indentation depth $h$ has exceeded $h_s$, i.e. when the surface asperities are crushed, the contact force follows Hertz elastic contact model: (6) $F_N=k_p(h-h_1{)}^{1.5}$(6) where $h_1$ is introduced to ensure a continuous contact force as function of $h$ and is given by (7) $h_1=h_s-{\left({\displaystyle \frac{F_s}{k_p}}\right)}^{2/3}$(7) and the contact stiffness $k_p$ is given by contact theory from Hertz (1881): (8) $k_p={\displaystyle \frac{4E_{\mathrm{eff}}\sqrt{R_0}}{3}}$(8) where the effective Young’s modulus, $E_{\mathrm{eff}}$, of the two bodies in contact is defined as (9) $E_{\mathrm{eff}}={\left({\displaystyle \frac{1-v_1^2}{E_1}+{\displaystyle \frac{1-v_2^2}{E_2}}}\right)}^{-1}$(9) $E_1$, $E_2$, $v_1$, $v_2$ represent Young’s moduli and Poisson’s ratios of the interacting bodies. At large contact forces, a deviation from the Hertzian relationship occurs experimentally due to shear driven damage that can be modelled as plasticity (Olsson et al 2019a). In the present work, this effect is modelled by a linear force–displacement relationship when the indentation depth $h$ exceeds a value $h_L$: initiation of plastic deformation. To ensure a continuous relationship and to have a continuous derivative, as seen experimentally, the contact force is given by (10) $F_N=k_p(h_L-h_1{)}^{1.5}+1.5k_p\sqrt{h_L-h_1}(h-h_L)$(10) where $h_L$ is determined to be $0.08\mathrm{mm}$ from experimental observations (Celma Cervera et al. 2017). In the simulations, $h_L$ is scaled with particles’ radius (Olsson et al. 2019a): (11) $h_L=0.08\mathrm{mm}{\displaystyle \frac{R_0}{6.25\mathrm{mm}}+h_s}$(11) The complete $F_N\left(h\right)$ relation is illustrated in Figure 4 with initial linear asperity crushing regime (dashed line), the elastic Hertzian regime (dotted line) and the linear regime with plastic deformations (straight line).

Figure 4. The force–overlap relationship.

In the DEM model, the tangential contact law determines frictional sliding between the particles. This slipping condition is determined by the equilibrium between the tangential forces ($F_t$) between particles and the shear strength ($F_t^{\mu }$) as shown in Equation (14). The calculation of $F_t$ is based on Agnolin and Roux (2007) and it can be expressed as (12) $F_t=\{\begin{array}{c}{k}_s\delta ,\left|\left|F_t\right|\right|\le F_t^{\mu }\mathrm{no}\mathrm{slip}\\ F_t^{\mu },\left|\left|F_t\right|\right|>F_t^{\mu }\mathrm{slip}\mathrm{occurs}\end{array}$(12) where $\delta$ is the tangential displacement in the contact point between particles, $k_s$ is the tangent shear stiffness and represented as (13) $k_s={\displaystyle \frac{2\left(1-v\right)}{2-v}1.5\left(k_p{)}^{2/3}\right(F_N{)}^{1/3}}$(13) Models with a non-constant tangential stiffness are cumbersome as this leads to spurious increases of elastic energy at some contact displacement histories. This is solved by scaling the tangential force as suggested by Agnolin and Roux (2007). The relation between shear strength and normal force is given according to Coulomb friction law, as shown as (14) $F_t^{\mu }=\mu F_N$(14) where $\mu$ is the friction coefficient. When the tangential force reaches the contact shear strength, slipping occurs.

2.3. Failure law

To simulate fracture of the aggregates in DEM, an accurate failure law is essential for a computational framework that can provide more accurate predictions at different loading conditions and configurations than the framework was calibrated for. Contact induced damage in rock materials was studied in detail experimentally in indentation experiments of Celma Cervera et al. (2017) and numerically in Olsson et al. (2019a). These results will be used to formulate the aggregate failure law.

Fracture of the aggregate is assumed to occur if the normal contact force $F_N$ exceeds a critical fracture force $F_C$ calculated from the contact failure stress ${\sigma }_F$ by (15) $F_C={\sigma }_F{R}_0^2$(15) The experiments by Celma Cervera et al. (2017) show a significant scatter in critical contact forces and, hence, a probabilistic description of the failure stress is relied upon. As the probability of finding a critical defect causing fracture is larger in a large loaded volume than in a small loaded volume, the Weibull weakest-link model is used for describing the probability distribution $S$ of the failure stress according to (16) $S=1-\exp \left[-{\left({\displaystyle \frac{{\sigma }_F}{{\sigma }_w}}\right)}^m{\displaystyle \frac{R_0^3}{V_{\mathrm{ref}}}}\right]$(16) where ${\sigma }_w$ is the Weibull stress, controlling the median strength of the particles, and $m$ is controlling the scatter of the distribution. The reference volume $V_{\mathrm{ref}}=(6.25\mathrm{mm}{)}^3$ is based on the radius of the indenter in Celma Cervera et al. (2017) and is only introduced in the description for dimensional consistency.

Once a contact force between particles reaches to the critical force value of one of these particles, that particle is designated as fractured. To obtain a realistic response of the sample, the post-fracture behaviour of the aggregate has to be modelled appropriately. One popular route in the literature is to model each aggregate as an agglomerate of many smaller particles cf . Cheng et al. (2003). The benefit is that the fracture process and the post-fracture behaviour can be captured in detail but comes with large computational costs as the number of particles in the simulation is by magnitudes larger than the number of aggregates. Another route is to replace the fractured particle with smaller bonded or unbounded particles cf . Cleary (2001). This decreases the computational burden as only fractured aggregates will consist of a large number of particles. However, this particle replacement method is not trivial due to mass and energy conservation considerations cf. Guo et al (2020). Instead, the stiffness reduction method will be used (Olsson and Larsson 2015) where a fractured aggregate experiences a reduction in its stiffness as (17) $E_{\mathrm{fractured}}=C_w\cdot E_{\mathrm{initial}}$(17) where $C_w$ represents the weakening coefficient and is taken as a calibration parameter. This is done to simulate local settling happening after a particle is broken. By reducing the $E_{\mathrm{initial}}$ of fractured particles, the contact stiffness between the particles is going to decrease, see Equation (8). Hence, the particles within the perimeter of those fractured particles will have a higher overlapping value in order to get into steady state, thus allowing surrounding particles to fall towards fractured particles.

In Olsson and Larsson (2015), this method provided a realistic response where a few particles fractured, and no large movements of particles occur. This is exactly the case studied presently, where a fraction of the particles is expected to fracture, and the particle movement is restricted by the confining pressure. One of the intentions with this work is to investigate, if such a computationally efficient model can capture the response of UGMs under compressive and confined loading.

Based on the particle contact and failure laws defined through Equations (3)–(17), four material parameters need to be identified: $E$, ${\sigma }_w$, $C_w$ and $m$. Among them $E$, ${\sigma }_w$ and$m$ may be measured in instrumented indentation tests, as presented by Celma Cervera et al. (2017). These parameters are reported by Olsson et al. (2019a, 2020) and also shown in Table 1.

Table 1. Parameters for the UGM models measured from instrumented indentation tests.

$E$ (GPa)	${\sigma }_w$ (MPa)	$m$ (-)
65	386.5	3.87

In this study, the values in Table 1 are used as an initial guess. Thereafter, the material parameters, representative for the aggregates used in this study are determined to gain the best fit of experimental results obtained for one of the gradation and at one load level, while the measurements at other load levels and at a different gradation are used to validate the results, as detailed in Section 5.

3. Experimental study

To investigate the influence of aggregate gradation and maximum compressive applied loading, $P_{\max }$, on aggregate crushing in UGMs, uniaxial compression tests have been conducted on two UGMs at three maximum load levels as detailed below.

3.1. Material preparation

In this study, crushed granite aggregates are used, supplied by NCC AB. The aggregates have an LA value of 25% for the fraction of 10–14 mm. The density of aggregates has been measured, following ASTM C127 (ASTM 2015), to be 2.65$\pm$0.06 g/cm³ with 95% of confidence. Aggregate crushing in UGM specimens with two different gradations is investigated: specimens composed of coarse aggregates in the 11.2–16 mm size range and specimens composed of aggregates in the 2–16 mm size range, with the gradation following the maximum density line. The two gradations investigated are shown in Figure 5 and are referred in what follows as ‘single size’ and ‘denser gradation’ respectively.

Figure 5. Gradation curves used in experiments.

The denser gradation distribution illustrated in Figure 5 is inspired by the Fuller curve used for dense packing as shown in Equation (18): (18) $G_V={\left({\displaystyle \frac{d}{d_{\max }}}\right)}^{0.5}$(18) where $G_v$ is the cumulative volume distribution, $d$ is the aggregate size and $d_{\max }$ the maximum aggregate size. The aggregates are assumed to have same density, so the Fuller curve is converted to volume passing.

3.2. Test procedure

UGMs are put into a cylindrical steel mould with 100 mm diameter, 205 mm height and 10 mm wall thickness as shown in Figure 6(b). The material weight is set between 1400–1500 g and 1580–1590 g for single size gradation and denser gradation respectively. This resulted in initial porosities between 0.44 and 0.49, and 0.41 and 0.46 for single size and denser gradations respectively.

Figure 6. (a) Uniaxial monotonic compression testing setup and (b) cross-section of the cylindrical mould and the testing equipment.

A cylindrical loading plate with 99.72 mm diameter and 25 mm height is used to transfer the vertical load to the specimen in the confining cylinder as shown in Figure 6(a, b). The aggregates are filled into the cylinder uniformly layer by layer. Experiments are conducted, at room temperature, with an MTS servo hydraulic machine, with 100 kN load capacity. Before testing, the UGMs are preloaded up to 1 kN. Later, the load is applied in displacement-controlled mode at a rate of 0.03 mm/s, and the resulting reaction force is continuously recorded during the test. The rate is chosen following earlier studies on aggregate breakage, cf. Wang et al. (2016) and Cai et al. (2020). Once the reaction force reaches $P_{\max }$, the specimen is unloaded, and the sieve analysis is conducted.

For each of the two gradations tested, two $P_{\max }$ levels are investigated. In case of single size gradation, $P_{\max }=$ 20 and 35 kN were chosen and $P_{\max }=$ 35 and 50 kN were used for the denser gradation. The force levels of 20, 35 and 50 kN correspond to vertical pressure levels of 2.55, 4.46 and 6.37 MPa respectively. These force levels were chosen, based on trial experiments, to ensure that there is a measurable amount of aggregate breakage at each force level. For each gradation and $P_{\max }$ level, three test repetitions were performed.

To capture the change in gradation due to aggregate breakage, the gradation of the specimens before and after the test is compared. The difference between the gradations is quantified in this study with the so-called Bg-Index value, following Bjarnason et al. (2000): (19) $\mathrm{Bg}=\sum _{i=1}^{n-1}\mathrm{\Delta }{S}_a^i-\mathrm{\Delta }{S}_b^i$(19) where n is the number of sieves and $\mathrm{\Delta }{S}_a^i$ and $\mathrm{\Delta }{S}_b^i$ represent the volumes of material retained at sieve size $i$ after and before the test respectively, shown in Figure 7. $<>$ are the Macaulay brackets that include only positive resulting values into the summation.

Figure 7. Definition of quantities used for the Bg-index evaluation (exemplary figure).

To account for the effect of the machine compliance $C_m$ on the measurements, the resulting compression magnitude is calculated as follows from the piston displacement: (20) ${\mathrm{u}}_{\mathrm{z}}\left(t\right)={\mathrm{u}}_{\mathrm{z},\mathrm{i}}\left(t\right)-{\mathrm{P}}_{\mathrm{z}}\left(t\right)\cdot C_m$(20) where ${\mathrm{u}}_{\mathrm{z}}\left(t\right)$ is the compression depth, ${\mathrm{P}}_{\mathrm{z}}\left(t\right)$ and ${\mathrm{u}}_{\mathrm{z},\mathrm{i}}\left(t\right)$ represent the measured load and the piston position, and $C_m$ is the machine compliance. Machine compliance was measured by pressing the loading head against the steel plate. In this way, the relation between force and load frame deformation was found to be linear, with $C_m=8.22\cdot {10}^{-6}$mm/N.

4. Computational study

To examine the capabilities of the modelling approach outlined in Section 2, DEM models of the uniaxial compression tests on the single size and denser gradation are developed. In the DEM models, the lateral confining cylinder with $100\mathrm{mm}$ diameter and with a height of $2.0\mathrm{m}$ is created to represent the cylindrical mould of the experiments. The height of the cylinder is chosen significantly higher compared to the cylindrical mould used in the experiments so that particles can be generated within this confined space without overlaps. The particles are randomly generated first in a cylindrical confinement. The number of particles is set to simulate the UGM volume identical to the corresponding experiment. As a first loading step, the gravity is applied, and the generated particles are allowed to fall to the bottom of the cylinder. Then the generated material is compressed through a rigid loading head that transfers the load to the material. The compression is conducted in a displacement-controlled mode with a velocity of 1 mm/s to reduce computation time. It should be noted that the velocity of the loading head is somewhat higher than in the experiment. To ensure quasi-static and time independent compression process in DEM, the dimensionless inertial number, $I$, is calculated as (21) $I=\dot{ϵ}d\sqrt{{\displaystyle \frac{\rho }{{\sigma }_{\mathrm{mean}}}}}$(21) where $I$ is the dimensionless inertial number, $\dot{ϵ}=0.005{\mathrm{s}}^{-1}$ is the axial strain rate, $d$ is the biggest particle size, $\rho$ is the particle density and ${\sigma }_{\mathrm{mean}}=\left({\sigma }_x+{\sigma }_y+{\sigma }_z\right)/3$ is the mean confining stress, where ${\sigma }_x$, ${\sigma }_y$ and ${\sigma }_z$ represent stress components acting on confining cylinder walls in x, y and z directions respectively. For all the cases simulated $I$ is found to be below $2.0\cdot {10}^{-5}$ at the initiation of particle breakage. As discussed in Da Cruz et al. (2005), $I<{10}^{-3}$ indicates that inertial contribution is small. To control the resulting compressive load, the reaction force on the loading head is registered. The contact force of each particle is registered every 0.1 mm of displacement increment of loading head.

Once DEM particles fall to the bottom of the cylinder, particles do not sit well to the sides of the confining boundary due to the artefact of using particles with spherical shape. Therefore, the packing density of UGM close to the walls is lower, whereas the middle of the UGM is more representative of the experiments. Hence, the volume within 10% shorter radius inside the confining cylinder is considered to calculate porosity.

Friction coefficient, $\mu$, local damping coefficient, $\zeta$, and Poisson’s ratio of particles are chosen as 0.7, 0.1 and 0.15 respectively following Celma Cervera et al. (2017) and Olsson et al. (2020). This relatively high value of $\mu$ is chosen to represent the effect of interlock due to surface irregularities. The stiffness of the confining boundaries is chosen 200 times higher than the stiffness of the particles to ensure that confining boundaries act as rigid body.

DEM models are created for both single size and denser gradations used in the experiments. Single size gradation is represented in the modelling as UGM composed of $13.0\mathrm{mm}$ spherical particles, resulting in 495 particles. The initial porosity for the single size model is 49%. For the denser gradation model, particles in the range of 2–16 mm are generated following Equation (26) and Figure 5. The resulting amount of particles is 15003 and the initial porosity is 44%.

To randomly generate a given number of spherical particles from a given volume-based size distribution, the cumulative probability distribution of the particles, $G_N$, has to be obtained from the corresponding and measured cumulative probability distribution of volume of particles, $G_V$, which is shown in Equation (18). The corresponding probability density distributions are denoted $g_N$ and $g_V$ respectively. The volume $\mathrm{dV}$, passing when the sieve size is increased from $r$ to $r+\mathrm{dr}$, can be expressed using either of those two distributions, where $r=d/2$ (22) $\mathrm{dV}=V_0\underset{r}{\overset{r+\mathrm{dr}}{\int }}{g}_v\left(r^{\mathrm{\prime }}\right)d{r}^{\mathrm{\prime }}={\displaystyle \frac{4\pi r^3}{3}{N}_0\underset{r}{\overset{r+\mathrm{dr}}{\int }}{g}_N\left(r^{\mathrm{\prime }}\right)d{r}^{\mathrm{\prime }}}$(22) where $V_0$ and $N_0$ are the total volume and number of particles respectively. The expression for the probability density function for the number distribution $g_N\left(r\right)$ can thus be expressed by (23) $g_N\left(r\right)={\displaystyle \frac{\text{β}{g}_v}{r^3}}$(23) where $\text{β}$ is a normalising factor ensuring that (24) $\underset{r_{\min }}{\overset{r_{\max }}{\int }}{g}_N\left(r\right)\mathrm{dr}=1$(24)

The cumulative number distribution $G_N\left(r\right)$ is obtained by integration (25) $G_N\left(r\right)=\underset{r_{\min }}{\overset{r}{\int }}{G}_N\left(r^{\mathrm{\prime }}\right)d{r}^{\mathrm{\prime }}=\text{β}\underset{r_{\min }}{\overset{r}{\int }}{\displaystyle \frac{g_v\left(r^{\mathrm{\prime }}\right)}{r^{\mathrm{\prime }3}}d{r}^{\mathrm{\prime }}}$(25)

Using partial integration, the following equation can be attained which uses the cumulative volume distribution $G_v\left(r\right)$ obtained experimentally directly. (26) $G_N\left(r\right)=\text{β}\left[{\displaystyle \frac{G_v\left(r\right)}{r^3}+3\underset{r_{\min }}{\overset{r}{\int }}{\displaystyle \frac{G_v\left(r\right)}{r^4}\mathrm{dr}}}\right]$(26)

To quantitatively compare experimental and numerical observations concerning aggregate crushing, Bg-Index values need to be estimated from DEM simulations. Presently, it is done, by attributing the fractured particle volume to the next smaller sieve size, i.e. it is assumed that, for example, particles retained on 8 mm sieve break into fragments smaller than 8 mm and bigger than 4 mm. The assumed volume redistribution is substantiated, at least to some degree, by the results of crushing tests on single size gradation materials. There, it was observed that approximately 70% of crushed material has been retained on the next sieve size, with the rest of the material ending up on smaller sieve sizes. It has to be emphasised nonetheless that the assumed volume redistribution is still somewhat simplistic and more thorough insight into the issue may be obtained, e.g. from combining particle crushing tests with X-Ray Computed Tomography studies. This is planned as a part of future research.

5. Results and discussion

First, the experimental results concerning the influence of gradation and maximum compressive applied loading, $P_{\max }$, on aggregate crushing in UGMs are summarised. Then, results from a brief parametric study performed on the DEM model of single size gradation experiments at $P_{\max }=$35 kN are presented and used to identify the material parameters for contact and failure laws, $E$, ${\sigma }_w$, $C_w$ and $m$ as discussed in Section 2.3. The identified material parameters are then used in DEM simulations of the uniaxial compression tests on single size gradation at $P_{\max }=$20 kN and denser gradation at $P_{\max }=$35 and 50 kN. Finally, by comparing the results, the capability of the models to capture the experimental observations is examined.

5.1. Experimental results

In Figure 8, the vertical load $,{\mathrm{P}}_{\mathrm{z}}$, recorded during the tests is shown as a function of the vertical loading head displacement, ${\mathrm{u}}_{\mathrm{z}}$, for the single size and denser gradations. The curves presented in Figure 8 correspond to the experiments at $P_{\max }=$20 kN and 50 kN for the single size and denser gradations respectively. The results for the $P_{\max }=$35 kN levels are qualitatively similar and omitted for briefness. In Figure 8, ${\mathrm{P}}_{\mathrm{z}}\left({\mathrm{u}}_{\mathrm{z}}\right)$ is presented as mean value of three tests along with the softest and the stiffest responses measured for each material.

Figure 8. Loading head’s force–displacement graph for (a) single size gradation, $P_{\max }$ = 20 kN and (b) denser gradation, $P_{\max }$ = 50 kN from experimental results.

As seen in Figure 8, the ${\mathrm{P}}_{\mathrm{z}}\left({\mathrm{u}}_{\mathrm{z}}\right)$ responses of both materials are qualitatively similar, having approximately parabolic shape and similar stiffness. This becomes clear by comparing the results presented in Figures 8(a) and (b) for ${\mathrm{P}}_{\mathrm{z}}$ < 20 kN. For ${\mathrm{P}}_{\mathrm{z}}>5$kN and ${\mathrm{P}}_{\mathrm{z}}>45$ kN for the single size and denser gradations respectively, there are jumps in the specimens’ force–displacement responses. Those jumps are associated with the aggregate rearrangements due to frictional sliding and crushing. As it is also seen in Figure 8 that the measured ${\mathrm{P}}_{\mathrm{z}}\left({\mathrm{u}}_{\mathrm{z}}\right)$ curves are accompanied by significant scatter (in particular, at higher load levels), which is furthermore significantly higher for the single size gradation compared to the denser graded one. For single size graded material, the softest and stiffest force–displacement curves deviate from the mean curve with up to 25%, with the corresponding deviation being approximately 12% for the denser graded material.

The Bg-Index (Equation 19) calculated for both gradations and all the $P_{\max }$ levels applied is presented in Figure 9. For each gradation and $P_{\max }$, the mean Bg-Index values measured in three tests are presented along with the lowest and highest values obtained. Obviously, for both gradations investigated, increasing $P_{\max }$ by 15 kN results in approximately doubling the amount of aggregate crushing. Furthermore, the amount of aggregate crushing in the case of denser gradation is significantly lower compared to the single size one. This is attributed to the higher number of aggregate-to-aggregate contact points resulting in lower magnitude of local contact forces and therefore less damage on aggregates. Similarly to the results presented in Figure 8, the scatter in Bg-Index values obtained in case of a single size gradation is somewhat higher than for the denser graded one. The results presented in Figures 8 and 9 demonstrate that the chosen gradations and load levels result in significantly different specimen responses both with respect to the ${\mathrm{P}}_{\mathrm{z}}\left({\mathrm{u}}_{\mathrm{z}}\right)$ and the amount of aggregate crushing.

Figure 9. Measured Bg-index values.

5.2. Material parameters determination

As discussed in Section 2.3, four parameters need to be identified for describing deformation and damage of aggregates due to contact loads: Young’s modulus $\left(E\right)$, Weibull stress and exponent $\left({\sigma }_w\mathrm{and}m\right)$, and weakening coefficient$\left(C_w\right)$. Presently, parameter values identified by Olsson et al. (2019a), shown in Table 1, are used as a starting point and the experimental results obtained for single size gradation up to $P_{\max }$ = 35 kN are used to adjust these parameters for the used aggregate type.

Parameter identification procedure is as follows. First, Young’s modulus ($E$) of the particles is identified to obtain the best possible fit of the initial part of the ${\mathrm{P}}_{\mathrm{z}}\left({\mathrm{u}}_{\mathrm{z}}\right)$ curve. The initial part of the curve is used assuming the absence of aggregate damage and excessive rearrangement, therefore expecting aggregate contact stiffness controlling the measured response. In these simulations, breakage of aggregates is not modelled. Then Weibull Stress (${\sigma }_w$) and weakening coefficient ($C_w$) are identified giving the optimal agreement with respect to the Bg-index values and the ${\mathrm{P}}_{\mathrm{z}}\left({\mathrm{u}}_{\mathrm{z}}\right)$ curve at high load region. Finally, the effect of the exponential constant ($m$) on the simulated aggregate crushing is examined. In Table 2, parameter values examined in this study are summarised.

Table 2. Values used in identification of contact and damage law parameters in DEM models.

	Parameters
Values	$E$ (GPa)	${\sigma }_w$ (MPa)	$C_w$ (-)	$m$ (-)
1	65	386.5	0.9	5.00
2	55	300.0	0.7	3.87
3	45	200.0	0.5	3.20
4	–	–	0.3	2.00
5	–	–	0.1	–

In Figure 10, ${\mathrm{P}}_{\mathrm{z}}\left({\mathrm{u}}_{\mathrm{z}}\right)$ relationship is presented as obtained in simulations of single size gradation experiments at $E$ = 45, 55 and 65 GPa, along with the softest and stiffest values measured in the tests. Increasing E results in an increase of uniaxial stiffness for simulated UGM behaviour. However, the results from simulation deviate significantly from the experimental ones at ${\mathrm{P}}_{\mathrm{z}}>15$ kN, indicating that at higher loads aggregate breakage starts to affect the measured force–displacement response. As seen in Figure 10, $E$ = 45 GPa results in the best agreement with the measurements. Accordingly, it is used in further modelling, being well aware, that 45 GPa is somewhat lower than the Young’s modulus values reported for granite in the literature. For instance, indentation tests results by Celma Cervera et al. (2017) show average Young’s modulus for Arlanda and Skyttorp granites as 70 GPa or higher, while average Young’s modulus of 52 GPa for Bohus granite was reported from quasi-static compression and direct tension tests by Saadati et al. (2014). However, it should be pointed out that the modulus value used in DEM simulations should not only capture the linear elastic properties of the aggregates but also capture the effect of local contact geometries on contact stiffness as well as of the non-linear deformations in the vicinity of the contact points. Both of these effects may result in decreasing the contact stiffness, as discussed by Uthus et al. (2008) and Olsson et al. (2019a).

Figure 10. $P_z\left(u_z\right)$ comparison of single size gradation simulations with different Young’s modulus.

The effect of ${\sigma }_w$ and $C_w$ on the simulated ${\mathrm{P}}_{\mathrm{z}}\left({\mathrm{u}}_{\mathrm{z}}\right)$ and Bg-Index values is examined in Figures 11 and 12 respectively. Qualitatively, it can be stated that the ${\sigma }_w$ magnitude controls the amount of fractured particles, as reflected by the change in the Bg-Index, and $C_w$ defines the effect of particle fracture on UGM stiffness, as reflected by the change in the ${\mathrm{P}}_{\mathrm{z}}\left({\mathrm{u}}_{\mathrm{z}}\right)$ curves. In the simulations where ${\sigma }_w$ is varied, $C_w$ is set to 0.5, while in the simulations with varying $C_w$, ${\sigma }_w=300$ MPa is used. In Figure 11(a), the ${\mathrm{P}}_{\mathrm{z}}\left({\mathrm{u}}_{\mathrm{z}}\right)$ relationship is presented as obtained in the simulations with ${\sigma }_w=$200, 300, 387.5 MPa along with the softest and stiffest test values. As shown, incorporating particle breakage improves the agreement between the simulated and measured UGM responses at high loads. At the same time, the effect of ${\sigma }_w$ on the specimen stiffness is not very strong and is visible only at force levels creating initiate breakage. This is expected as ${\sigma }_w$ controls the strength of the individual particles and as broken particles represent a relatively small proportion of the aggregates in UGM. At the same time, the effect of $C_w$ on the simulated ${\mathrm{P}}_{\mathrm{z}}\left({\mathrm{u}}_{\mathrm{z}}\right)$ is quite profound as illustrated in Figure 11(b). The decrease in $C_w$ results in significant decrease of specimens’ stiffness at $P_z>5\mathrm{kN}$, i.e. when aggregate crushing starts.

Figure 11. $P_z\left(u_z\right)$ comparison of UGMs with (a) different ${\sigma }_W$ with $C_w$ = 0.5 and (b) different $C_w$ with ${\sigma }_w$ = 300 MPa.

Figure 12. Bg-index comparison of UGMs with (a) different ${\sigma }_W$ with $C_w$ = 0.5 and (b) different $C_w$ with ${\sigma }_w$ = 300 MPa.

The Bg-Index values for the same simulation cases as in Figure 11 are presented in Figure 12. It is observed that while ${\sigma }_w$ affects the amount of particle breakage profoundly, the effect of $C_w$ is much smaller. Comparing the simulation results presented in Figures 11 and 12 with the experimental findings, it is concluded that ${\sigma }_w=$ 300 MPa and $C_w=$ 0.3 results in best agreement with the range of measured data.

The Weibull exponent, $m$, controls the probability distribution of particle strengths assigned in the simulations. This is illustrated in Figure 13(a), where a cumulative probability distribution of strengths values assigned to particles according to Equation (16) is presented for different m values.

Figure 13. (a) Scatter of $F_C$ to each particle within the DEM models and (b) m comparison between Bg-index differences and mean Bg-index.

As seen in Figure 13(a), decreasing $m$ results in significantly wider distributions of particle strengths. Accordingly, it may be argued that this parameter will change the distribution of critical fracture force, $F_c$, assigned between each particle, and consecutively it will lead to variations in simulated crushing performance of UGMs generated with the different random seeds, resulting in different randomisations of strength distributions. To evaluate this, for each of the m values, four simulations were performed. In the simulations, aggregates spatial distribution was kept constant and different random particle strength distributions were assigned. In Figure 13(b), the resulting mean Bg-Index from four simulations is presented as a function of $m$ along with the difference between the softest and stiffest Bg-index, which is represented as $\mathrm{\Delta }$Bg-Index in Figure 13(b). As shown in Figure 13(b), the differences between the simulations decrease rapidly with increasing $m$, the scatter between the simulated values is approximately 4 times lower for $m$ = 3.87 as compared to the one observed at $m$ = 2. Since four simulations may not fully capture the variations in the material, it is advisable to conduct a statistical evaluation based on a larger number of simulations. Following results presented in Figure 13(b) and previous results by Olsson et al. (2019a), $m$ = 3.87 is used in the simulations of crushing tests. The determined values for contact and failure law parameters are summarised in Table 3.

Table 3. Determined parameters of the UGM models.

$E$ (GPa)	${\sigma }_w$(MPa)	$C_w$ (-)	$m$ (-)
45	300	0.3	3.87

5.3. Analysis of crushing experiments

To evaluate the capability of the developed model to capture the effects of gradation and loading on the aggregate crushing, the simulations of the single size gradation at $P_{\max }=$20 kN and of the denser gradation at $P_{\max }=$35 and 50 kN were examined as follows. For all the cases simulated, the contact and failure law parameters, presented in Table 3, were used.

Force–displacement curves, ${\mathrm{P}}_{\mathrm{z}}\left({\mathrm{u}}_{\mathrm{z}}\right)$, obtained for the single size and densely graded UGMs up to $P_{\max }$ = 20 kN and $P_{\max }$ = 50 kN respectively are shown in Figure 14 along with softest and stiffest values from the experiments. For the single size material, the simulations capture the measured response very well, which is, of course somewhat expected, given that the measurements on single size gradation have been used for the parameter identification. It may also be pointed out that jumps in modelled and measured ${\mathrm{P}}_{\mathrm{z}}\left({\mathrm{u}}_{\mathrm{z}}\right)$ curves are of approximately the same magnitude of 4 kN. This indicates that the effect of aggregate breakage and rearrangement is well captured by the model. As seen in Figure 14(b) for the denser gradation, the simulated response agrees well with the measurements for the loads ${\mathrm{P}}_{\mathrm{z}}\le 20$ kN. At higher loads, the simulated response starts to deviate, overestimating the material stiffness. The maximum difference of approximately 20% is observed at ${\mathrm{P}}_{\mathrm{z}}=50$ kN. This slightly higher deviation between experimental and modelling results may be attributed to the fact that the $C_w$ parameter determining the effect of crushing on the UGM stiffness has been identified in this study from the single size gradation measurements alone. It may be expected that smaller aggregates fragment differently as compared to the large ones. Accordingly, the effect of their fragmentation on the UGM stiffness will be different. Nevertheless, the 20% accuracy of the model for the denser gradation case is an encouraging result given that no adjustment of material parameters has been made between the simulations. It is also remarkable that the magnitude of jumps in the ${\mathrm{P}}_{\mathrm{z}}\left({\mathrm{u}}_{\mathrm{z}}\right)$ curve is significantly smaller for denser gradation, which is in qualitative agreement with experimental findings.

Figure 14. Simulated and measured $P_z$ ($u_z$) curve of (a) single size gradation UGM $P_{\max }$ = 20 kN and (b) denser gradation UGM $P_{\max }$ = 50 kN.

In Figure 15, the Bg-Index, obtained from the modelling of both gradations is presented as a function of applied compressive load. For comparison, the measured value intervals are also indicated in the figure. It may first be noted that the modelling and experimental results are in good qualitative agreement giving the same ranking for all the gradation/loading combinations. Quantitative agreement is also quite encouraging with experimental and modelling values, deviating not more than 0.5%, except for the case of denser gradation and $P_{\max }=50$ kN, where a deviation of approximately 1.2% is observed. It can be concluded that the identified contact and failure law parameters capture the aggregate crushing reasonably well for all the gradations and load cases examined, thus confirming the methodology to determine these parameters in this study. The agreement may furthermore be improved by introducing a more elaborate description of particle fragmentation and this development is planned in future studies.

Figure 15. Bg-index from DEM simulations as a function of applied compressive load.

The distributions of fractured particles within the UGM specimens are shown in Figure 16. The distributions are presented for $P_{\max }=$ 35 kN for single size and denser gradations respectively. As seen in Figure 16, most of the fractured particles are located close to the top of the specimens. This is due to the fact that the compressive forces applied by the loading head are directly transferred to aggregates on top. Accordingly, the load is distributed more uniformly between the aggregates at the bottom of the specimen compared to its top. For the case of denser gradation, it may be seen that in terms of volume the majority of cracked particles belongs to the coarse fraction. Moreover, for dense gradation a smaller quantity of fractured particle appears as compared to single size gradation. This is explored further in Figure 17, where the volume of cracked particles is presented for each sieve size, normalised with the total volume of crushed aggregates. As seen in Figure 17, for both $P_{\max }$ loadings, the absolute majority of volume of crushed material is generated by the particles in the 11.2–16 mm range. This is explained by the fact that their larger volume results in higher probability of containing critical defects, it may also be argued that small particles are subjected to smaller contact forces as the coarse aggregates play a primary role in load transfer. Accordingly, ensuring adequate crushing resistance of coarse aggregates is crucial for the performance of UGM.

Figure 16. Distribution of fractured particles (pink) in (a) single size UGM at 35 kN and (b) densely graded UGM at 35 kN.

Figure 17. Broken particle volume in simulated denser gradation UGM at two load levels.

In summary, the results presented in Figures 14 and 15 indicate that the developed modelling approach captures well the effects of gradation and loading on mechanical behaviour of the UGM in terms of both load–displacement response and aggregate crushing. Accordingly, it is a starting point to establish the link between the loads experienced by unbound pavement layers in the field and their material design. Furthermore, one advantage with the DEM modelling is that it allows to gain a micromechanical insight into material behaviour which may be difficult to obtain experimentally thus reducing experimental effort and uncertainty drawbacks. As illustrated in Figures 16 and 17, where spatial and size distribution of cracked particles has been examined, the proposed DEM model may provide a useful tool for the material design optimisation.

6. Conclusion

In this study, a new discrete element method-based (DEM) model of aggregate crushing in unbound granular materials (UGM) is presented, incorporating mechanics-based granular particle contact and damage laws into a commercial software and introducing weakening phenomena to capture macro-mechanical behaviour of UGM. The model is applicable for different gradations and loading conditions. To identify the particle contact and the damage law parameters as well as to validate the model, aggregate crushing in uniaxial monotonic compression tests has been studied both experimentally and numerically for two different gradations, single size and denser gradations, and at three maximum load levels. The measurements obtained at one load level for the single size material were used to identify the contact and damage law parameters, while the rest of the experimental results were used to validate the model.

Based on the comparison of experimental and numerical findings, it is shown that the developed model captures the effect of gradation as well as of static compressive load magnitude on the UGM response reasonably well. In particular, for all cases examined, the force–displacement responses of UGMs are captured reasonably well by the model. For the majority of cases, the agreement between the simulations and experiments is within 10%, except for the denser gradation subjected to the loads above 40 kN where deviations up to 20% are observed. Furthermore, the numerical and experimental results concerning the amount of aggregate breakage are found to be in close qualitative and quantitative agreement. The remaining discrepancies between modelling and experimental results may, at least partially, be attributed to a somewhat simplified procedure used in this study to simulate the particle fragmentation after the breakage. Accordingly, the agreement may be improved further by introducing a more accurate particle fragmentation model. This development was, however, beyond the scope of this study.

The observed agreement between experimental and numerical findings indicates that the proposed particle contact and damage laws capture the interactions in the material correctly. The model thus may become a useful tool for designing UGMs and evaluating their crushing performance in unbound layers of asphalt pavements with different layer geometries and load cases. The developed approach may furthermore be useful in optimising unconventional unbound road materials such as construction and demolition waste materials against aggregate breakage. For that, however, the model accuracy should be evaluated for a wider range of UGMs composed of different aggregates, such an evaluation is planned as a part of future research.

Acknowledgements

The authors express sincere gratitude to Hassan Fadil, Gürsel Hakan Taylan and Viktor Brolund for their help in the laboratory.

Disclosure statement

No potential conflict of interest was reported by the author(s).