The relationship between contact network and energy dissipation in granular materials

A micromechanical study of the relationship between contact force networks and energy dissipation is presented. A series of drained triaxial compression tests with different stress paths have been simulated using the discrete element method. Two existing contact force network partitioning methods have been used for analysing the energy dissipation, one based on the average contact force magnitude and the other based on the contribution of contact forces to the global deviator stress. For both methods, energy dissipation in neither the strong nor weak contact networks is negligible. When the average contact force partitioning method is used, over 70% of the energy dissipation occurs in the weak contact network, but the dissipation per sliding contact is higher in the strong contact network because the tangential contact force is higher. When the contact network is partitioned based on the contribution of forces to global deviator stress, more than 60% of the energy dissipation occurs in the strong contact network. A new normal contact force threshold for splitting energy dissipation is identified. Specifically, over 93% of energy dissipation occurs at contacts with a normal contact force below 2 times the average normal contact force. There is very small energy dissipation in contacts with higher normal contact force because there is little particle sliding.


Introduction
Granular materials are widely seen in nature. They are the second-most manipulated material in the industry after water. Sand is a typical material that is frequently encountered in geotechnical engineering. The mechanical behaviour of sand (e.g., shear strength and dilatancy) at the macroscale is dependent on the interaction between sand particles at the micro-scale [32,4,17,29]. There has been extensive research on bridging the gap between the micro and macro responses of sand [1,[17][18][19]28]. In particular, it is found that the energy dissipation in granular materials which occur at the sand particle level is closely related to the macroscale response, such as shear strength, critical state and dilatancy [4,6,16,21].
Though some theoretical research has studied the relationship between energy dissipation and sand behaviour, there is insufficient understanding of the mechanism of energy dissipation at the microscale [9,11,30,31,33]. Discrete element simulations have shown that interparticle sliding is more likely to occur in contacts with smaller contact forces [29]. This is in agreement with the hypothesis made by Radjai et al. [27] that most energy dissipation occurs in weak contact networks. These hypotheses have been made based on the average force network partitioning method that considers the magnitude of contact forces [27]. The average force network partitioning method divides the contact force network into two complementary subnetworks: a strong network that carries an above-average normal contact force and a weak network carrying a below-average normal contact force. However, there is evidence to suggest that separating contact force in this way may not be appropriate. For example, Peters et al. [25] demonstrated that only half of all the strong contacts contribute to force transmission. Huang et al. [15] used DEM simulations to examine the average force network partitioning method. They concluded that using average normal force to separate the networks is not a robust method because some of the weak contacts carrying below-average normal contact forces also participate in stress transmission and contribute to structural anisotropy.
Later studies by Collins and Kelly [5] and Collins and Hilder [7] suggested that energy dissipation during normal compaction is predominantly caused by plastic deformation of the strong contact networks. However, locked-in elastic energy is mostly located at weak contacts, where stress is insufficient to cause plastic deformation. Conversely, energy is dissipated mostly in weak contacts during shear. Even though the relationship between energy dissipation and contacts networks has been discussed in many studies, there is no adequate understanding of how much energy is dissipated by the sliding of the contact networks because none of the previous hypotheses has been examined. In addition, the previous hypotheses were developed based on the average force network partitioning method. If another partitioning method is used, it is unclear whether this will affect the relationship between energy dissipation and contact networks.
The objective of this work is to quantitatively investigate the relationship between energy dissipation and contact networks. Two methods for partitioning the contact force network will be used. The first one is based on the average contact force [27] and the second one considers the contribution of the contact network to the global deviator stress [15]. A new method for partitioning the contact force networks based on energy dissipation will be presented.

DEM simulations
A series of drained triaxial tests have been carried out using a modified version of LAMMPS [26]. DEM samples were created containing 8262 spherical particles with a minimum diameter of 0.1 mm and a maximum diameter of 0.14 mm. The number of particles was varied to ensure a representative volume element (RVE) was obtained. Figure 1 shows the particle size distribution (PSD). A simplified Hertz-Mindlin contact model was used and the model parameters are shown in Table 1. The samples were enclosed within periodic boundaries to ensure uniform deformation [8]. It was found that these boundaries eliminate inhomogeneities at the sample's periphery [14]. Two samples were produced in loose and dense states, respectively. The loose samples were first prepared by placing the particles randomly as a nontouching cloud in the periodic cell and then the cell was isotopically compressed to the target mean normal stress p 0 = 250, 500 and 1000 kPa with a friction coefficient μ of 0.3. Once the target stress was reached the samples were allowed to equilibrate. The dense sample was prepared using the twostep process [12]. The particles were placed randomly in a periodic cell and the sample was equilibrated with μ = 0 to reach the target p 0 . The μ was then set to 0.3 and the sample was left to equilibrate again. Drained triaxial compression tests were then carried out, following two different stress paths, one with constant radial stress σ r and the other with constant mean effective stress p . Loading was controlled by moving the upper boundary while the lower boundary was fixed. The shear rate was constant with the inertial number of I = 1 × 10 −4 , which indicates quasi-static loading. Local damping of 0.3 was used during the preparation but set to 0 during shearing, which ensures that all dissipation is due to interparticle friction during triaxial compression. Note that energy dissipation due to rolling friction or resistance is neglected because there is still dispute in this issue. Figures 2 and 3 show the stress-strain relationship of dense and loose sand. The dense sand shows strain softening and volume expansion while the loose sand shows strain hardening and volume contraction during the tests. Note that the true strain is used in this study. Furthermore, Figure 4 shows the stress ratio against axial strain, it can be noticed that regardless of the stress path applied, there is a unique relationship between a and q∕p when the initial void ratio is the same. Figure 5 shows the critical state line (CSL) in the e − log(p) plot and q − p plot. The critical state stress

Stress-strain relationship
The CSL is nonlinear in the e − log(p) plot, which is consistent with previous results [11,13,23,24]. It can be seen that the critical state void ratio is almost constant when p < 1000 kPa, which is different from that for real sand. There are several plausible reasons for this. First, spherical particles cannot represent the shape and movement of real sand particles in the 3D space [22]. Secondly, the contact model used here may not be able to produce the real contact behaviour between sand particles.

Energy dissipation
In this study, the energies were traced using the method implemented in Granular LAMMPS by [11]. The error in the energy balance during the simulation is examined by using the equation below as described in [11].
where W sn and W st are the normal and tangential strain energies respectively; W kt and W kr are translational and rotational kinetic energies respectively; W and W f are the boundary work and energy dissipation respectively; indicate the timestep. The error in the energy balance during triaxial compression is illustrated in Fig. 6a and b. The maximum error is observed at the beginning of shear, which is about 0.06%. This is caused by the sudden movement of the particles when the loading starts. Similar results have been reported by [11]. Furthermore, Bernhardt and Biscontin [3] have indicated that an increase in μ during shearing could (1)  produce unrealistic results at low strains because the contact was prevented from sliding at the start of loading.
It is well known that interparticle sliding in granular materials is the main source of energy dissipation [10]. The other sources of energy such as strain energy and kinetic energy are very small in the system due to quasistatic loading conditions [9][10][11]. In this study, the main focus is on the energy dissipation associated with sliding within the contact networks. During loading, the energy dissipation is continuously monitored. When sliding occurs, the increment of energy dissipation W f is calculated using the equation below where denotes a time step, Δs ← −1 is the incremental tangential displacement at the sliding contact between timestep ← − 1 . F c s is the tangential component of the contact force of sliding contacts; V is the volume of the sample and N c is the number of sliding contact in the system. Note that the energy increment has been normalized by the total volume because the response of a representative soil element is considered. The total frictional energy dissipation as accumulation is The incremental and cumulative energy dissipation is shown in Fig. 7a and b respectively. The energy dissipation with different soil densities and stress paths is also shown in the figures. It should be noted that the results presented in Fig. 7 are based on the dissipation of all contact sliding, which does not provide any information about the contribution of weak and strong contacts to energy dissipation. A method for determining the contribution of strong and weak contacts to energy dissipation will be discussed later. Figure 7a shows that the energy dissipation reaches a constant rate at the critical state. The initial void ratio and stress path do not affect the mechanism of energy dissipation except that more energy is dissipated when samples are sheard with constant r . Figure 7b shows the cumulative energy dissipation compared with boundary work, which are almost equal, especially when the samples reach the critical state. This is in agreement with previous studies [10,11]. (2)

Partitioning of contact networks
This paper utilized two methods for partitioning contact networks. The first method is based on average normal force magnitude [27], and the second is based on the contribution of contact networks to deviator stress [15]. Though the average contact force partitioning method is widely used, it was not considered a robust method when the contribution to global deviator stress is considered for 3D modelling [15]. To illustrate this, this work examines the average contact force partitioning method based on how the contact network contributes to deviator stress.
The deviator stress is calculated from the contact force, where the stresses are obtained based on the stress tensor equation in Bagi [2] where V is the volume of the domain considered f c i is the component of the contact force acting on the particle and l c j is the component of the branch vector from the centre of the contacting neighbour to the centre of the particle and N c is the number of contacts. To determine the contribution of weak and strong networks to total deviator stress, the total average stress tensor can be decomposed into contributions from the strong  the branch vector of strong and weak contacts, and N S c and N w c are the number of strong and weak contacts respectively. Figure 8 depicts the cumulative contribution of the contact force to deviator stress as a function of the normalised contact force, where f n is the normal contact force and ⟨f n ⟩ is the average contact force. The figure evaluates the contribution in several positions during the test: near peak (2%, 3% and 5% axial strain) and the critical state (30% axial strain and 50% axial strain). It is evident that the weak contacts (contact force below average) also contribute to deviator stress, with approximately half of the weak contacts contributing to the overall deviator stress which is in agreement with the findings by [15]. In addition, it can be noticed from Fig. 8 that approximately half the mean normal contact force ( f n ∕⟨f n ⟩ = 0.5) gives a reasonable separation for the contact network contribution to the overall deviator stress at all strain levels. Therefore, f n ∕⟨f n ⟩ = 0.5 is used as the second partitioning method in this study. This partitioning definition is based on the assumption that weak contacts do not contribute to deviator stress and strong contacts are responsible for carrying out all deviator stress. This method will thus be called a deviator stress partitioning method. As shown in Fig. 9, the value of f n ∕⟨f n ⟩ = 0.5 used for the deviator stress partitioning is valid for all the tests with different initial densities and stress paths. Therefore, the value of f n ∕⟨f n ⟩ = 0.5 is used for partitioning the force networks based on contribution to the deviator stress in a 3D system. Figure 10 illustrates the contribution of the strong and weak force networks to the total deviator stress and compares both partitioning methods [29]. It is apparent from the figure that all deviator stress is carried by strong contacts when the deviator partition method is used, which confirms that the load will be transferred via the strong contact only.
However, when the average force partition method is used it is apparent that some of the weak contacts also contribute to deviator stress.

Relationship between contacts network and energy dissipation
To understand the relationship between energy dissipation and contact networks, results at the microscale are needed. In this study, the objective is to identify the amount of energy dissipated by each type of contact network (strong and weak) separately. Note that recording data for each contact throughout the simulation is unfeasible due to the large amount of data produced. In addition, incremental  dissipation at consecutive timesteps must be calculated to assess the contribution of the strong and weak network, as contacts form and disappear and change in their normal force magnitude. Identifying the strong and weak contacts costs excessive computation time. Therefore, the energy dissipation was examined for a subset of every 1% axial strain. As shown in Fig. 11, restart files were created for every 1% strain. To evaluate the energy dissipation for a given strain level, a probe test was performed, which involves applying a small axial strain increment Δ a = 1 × 10 −5 from each restart file. Each probe was then analysed and the energy dissipation was calculated using the equations below.
where E c,fd is friction dissipation per single contact at timestep β, E −1 c,fd is friction dissipation per single contact at timestep β-1, E tot is total increment energy dissipation across all contact in timestep β, and ∅E tot is total cumulative energy dissipation at timestep β. These equations were implemented in MATLAB.
Since all the energy is dissipated by the sliding of contact networks, it is initially convenient to assess the percentage of contacts that are sliding during the shearing process. In this study, the particles considered to be sliding if F t = f n × where f n is the normal contact force and F t tangential force. The contribution of strong and weak networks to the proportion of total contact sliding is illustrated in Fig. 12. The figure also compares both partitioning methods used in this study. The results show that the highest proportion of sliding occurs in weak contacts, but some of the strong contacts are sliding as well. The same behaviour is observed when both partitioning methods are used. However, it is noticeable that there is an increase in the proportion of sliding in the strong contact when the deviator stress partitioning method is used, from around 18% to around 29%. Furthermore, there is also an increase in the proportion of sliding in the weak contact, from 55% to around 65%. Previous studies showed that all sliding would occur in weak contacts and all non-sliding contacts were strong contact [27]. In contrast, the present study proved otherwise, showing that the amount of sliding at strong contacts is not negligible. Energy dissipation was traced for strong and weak contact networks separately. Figure 13a and b illustrates the contribution of the strong and weak contacts to the total energy dissipation. Using the average force partitioning method, Fig. 13a shows that approximately 70% of dissipation occurs in weak contacts, while 30% occurs in strong contacts. In the deviator stress partitioning method, however, Fig. 13b shows that the energy dissipation occurs more in the strong contact with approximately 65%, while the weak contact is approximately 35%. Table 2 shows the split of energy dissipation at different strain levels for a dense sample. Note that the portion of sliding contacts is calculated with reference to a specific contact network (e.g., strong or weak). Contribution to energy dissipation is defined as the ratio of energy dissipation in a contact network and the total dissipation in the sample. The examination is specifically located before and after the critical state at 2%,3%,5%,30% and 50% axial strain. The   shows that although strong contact sliding was around 15% when the average force partitioning method was used, this caused around 30% of total energy dissipation. Around 53% of weak contacts slide when applying the average force partitioning method, causing around 70% of the total energy dissipation. Using the deviator stress partitioning method leads to an increase in the proportions of strong contacts sliding to around 25%, which results in their contribution to energy dissipation being around 65% of the total dissipation.
Additionally, using deviator stress partitioning increased the proportion of sliding in weak contacts significantly. However, the contribution of weak contacts to energy dissipation decreased to around 35%. The results show that the relationship between energy dissipation and contact networks is not influenced only by the number of sliding contacts but also by how much force each contact network carries. When the normal and tangential components of contact force in a strong contact network are large, the dissipation during sliding becomes greater. A similar relationship was identified between contact networks and energy dissipation with other soil densities or stress paths.
It is evident that there is significant energy dissipation in both the strong and weak networks when both partitioning methods are used. This indicates the importance of investigating energy dissipation in contact networks based on their normal force magnitude. To examine this, it is convenient to begin by examining the distribution of the contacts network based on f n ∕⟨f n ⟩ . Figure 14 shows the distribution of contact forces at the critical state using the probability density function PDF. Though there is an exponential decay after the peak, the number of contacts with f n ∕⟨f n ⟩ between 1 and 3 is not negligible. Figure 15 shows the rate of energy dissipation for contacts with different ranges of f n ∕⟨f n ⟩ near peak (5% axial strain) and critical state (30% axial strain). It is evident that the energy dissipation rate is much higher when f n ∕⟨f n ⟩ ≤ 2 . Figure 16 depicts the cumulative energy dissipation for several contact groups throughout the simulation. Results similar to Fig. 15 are found, with the majority of energy dissipation occurring at the contact f n ∕⟨f n ⟩ ≤ 2 . Furthermore, Table 3 shows the portion of sliding contacts and energy   2  19  62  29  71  29  73  64  36  3  16  54  26  74  25  64  63  37  5  15  52  27  73  24  62  62  38  30  15  50  29  71  24  60  66  34  50  15  51  28  72  25  62  64  36 dissipation for contacts with f n ∕⟨f n ⟩ ≤ 2 and f n ∕⟨f n ⟩ > 2 . It can be seen that the majority of energy dissipation (> 93%) occurs in contact with f n ∕⟨f n ⟩ ≤ 2 . The primary reason is that the portion of sliding contacts is much higher when f n ∕⟨f n ⟩ ≤ 2 (Table 3).

Conclusion
In this study, a DEM-based investigation is conducted to investigate the mechanism of energy dissipation in granular materials. The relationship between contact networks and energy dissipation is investigated for dense and loose sand in drained triaxial compression with different stress paths. Two partitioning methods for the contact force networks were used. The first partition method is based on the contact force magnitude and the second is based on the contribution to global deviator stress. The main conclusions are: • When the average contact force partitioning method is used, more interparticle sliding and energy dissipation occur in the weak networks. However, the strong contact network has a five times higher energy dissipation per sliding contact due to higher contact forces. • When the deviator stress partitioning method is used the sliding still occurs more in the weak contact. Nonetheless, there has been a slight increase in the proportion of strong contact sliding compared to the average contact force partitioning results. This increase of sliding in strong contacts has led to more energy dissipation in the strong contacts. • It is observed that the energy dissipation that occurs in both strong and weak contacts is not negligible when both partitioning methods are used. Therefore, A new threshold of f n ∕⟨f n ⟩ for partitioning the contact network has been identified. Specifically, f n ∕⟨f n ⟩ = 2 can be used to determine if the contact contributes to energy dissipation. Over 93% of the total energy dissipation occurs in contact with f n ∕⟨f n ⟩ ≤ 2.
The results of this study confirm that both the proportion of sliding contacts and the magnitude of tangential contact force between particles play a significant role in the relationship between contact networks and energy  dissipation. For simplicity, this work has used only a uniform particle size distribution. Mukwiri et al. [20] demonstrated that particle size distribution influences the total amount of energy dissipation. It is expected that the particle size distribution does not have a significant influence on the conclusions of this study. But more research is needed in this regard.  Table 3 Sliding and energy dissipation in contacts with different normal contact forces (dense sample with constant r =250 kPa) Axial strain % Portion of sliding contacts (%) Contribution to total energy dissipation (%) f n ∕⟨f n ⟩ > 2 f n ∕⟨f n ⟩ ≤ 2 f n ∕⟨f n ⟩ > 2 f n ∕⟨f n ⟩ ≤ 2   5  5  42  4  96  30  5  40  8  92  50  5  40  7  93