A comparative assessment and unification of bond models in DEM simulations

Bonded contact models have been increasingly used in the discrete element method (DEM) to study cemented and sintered particulate materials in recent years. Several popular DEM bond models have been proposed in the literature; thus it is beneficial to assess the similarities and differences between the different bond models before they are used in simulations. This paper identifies and discusses two fundamental types of bond models: the Spring Bond Model where two bonded particles are joined by a set of uniform elastic springs on the bond’s cross-section, and the Beam Bond Model in which a beam is used to connect the centres of two particles. A series of cantilever beam bending simulation cases were carried out to verify the findings and assess the strength and weakness of the bond models. Despite the numerous bond models described in the literature, they can all be considered as a variation of these two fundamental model types. The comparative evaluation in this paper also shows that all the bond models investigated can be unified to a general form given at a predefined contact point location.


Introduction
During the last two decades, there has been an increase in the number of attempts to apply the Discrete Element Method (DEM) to simulate the behaviour of cementitious materials like concrete, rocks and fibers [1][2][3][4][5][6][7]. A key reason for this growth is, unlike finite element method simulations, the ability of DEM simulations to naturally take into account the discontinuity or microstructure in the material behaviour allowing detailed analysis of the strength and failure mechanism of a cementitious material [8][9][10]. In general, the contact models in DEM can be classified into one of three groups: cohesionless, cohesive and bonded contacts. DEM was originally developed for cohesion-less systems [11], in which the contact model can be either a linear model like a Hookean spring-dashpot contact model [11] or a non-linear model such as the Hertz-Mindlin contact model [12]. In recent years, there have been increased efforts in applying DEM to simulate more complex systems such as using a cohesive contact model to simulate cohesive systems [13][14][15][16][17][18][19][20] or using a bonded contact model to simulate cementitious or bonded systems [1][2][3][4]21]. A cohesive contact model is often referred to particle-particle interaction caused by an adhesive force such as van der Waals type force or liquid bridges. While these models often have some finite contact area, the contribution of the contact area to the twisting or bending resistance of the contact is typically ignored. A cohesive contact exists at any time two particles are within a close proximity of each other for the cohesive force to come into play. Unlike the cohesionless family of contact models, particles do not necessarily need to be in physical contact for a cohesive contact to develop. These contacts can break and re-form. New cohesive contacts can be formed with any nearby particles. The final family of inter-particle contact models are bonded contact models. Bonded contacts, on the other hand, are formed at a bond initialisation timestep in a simulation and the breakage of a bond is irreversible. A bonded contact will be removed permanently upon failure, which is in contrast to the cohesive contacts which can re-form. Bonded contact models also typically provide bending and twisting resistance in the contact, although some cohesionless [22] and cohesive [13,23] also provide some bending and twisting resistance. The cross-sectional area of the bond is used to determine the bending and twisting resistance of the bonded contact and the use of a full rotation matrix allows the bending history from the initial free state to be considered. For example, a bonded contact model can easily capture the bending and twisting forces present in a static, deformed beam. A bonded-particle model is a common approach used to simulate a cementitious material or rock-like material [24][25][26][27][28][29]. Bonded-particle models consist of an assembly of particles that are bonded together to create a virtual material that approximates the bulk behaviour. Cracking is represented by the failure and breakage of bonds. The fracture process is captured by the joining of multiple bond breakages. The micro-mechanisms of the failure of materials can be investigated since the formation and growth of micro-cracks on a particle scale are captured, which makes bonded-particle models an attractive tool. The bonded contact model, from hereafter referred to simply as the bond model, is a hugely important component of the bonded-particle model methodology.
A number of different bonded contact models have been developed to study the failure of cementitious materials and these have been categorised into two main types of bond model -spring and beam models. The simplest type of bond model is the single spring bond model or 'pin' bond model, where an infinitesimal interface ("point glue") is used to bond particles [2]. The physical behaviour of the point bond can be envisioned as a spring with constant normal and shear stiffness acting at a point, which is capable of resisting the relative displacement of compression and tension actions but cannot resist bending or twisting actions [30]. Note that the single spring bond model is also called the "linear contact bond model" in the commercial software Particle Flow Code (PFC) [31]. Another bonded contact model available in PFC is called linear parallel bond model which is developed by Potyondy and Cundall [1]. The linear parallel bond model (referred to as parallel bond model or PBM for short hereafter) is envisioned as a "reinforced" pin bond model by installing a finite-size "glue" at the contact so that it can resist twisting and bending moments. Potyondy and Cundall [1] summarised the motivation for the development of the parallel bond model and described a parallel bond as a set of elastic springs, uniformly distributed over a circular crosssection, centred at the contact point. Therefore, the parallel bond model is classified as a multi-spring bond model. The advantages and limitations of using the parallel bond model for predicting the damage of rocks have been discussed elsewhere [3,32,33]. Besides springs, beam elements (as another fundamental structure element) are also used in bonded-particle models to mimic the behaviour of the cemented bond [9,[34][35][36][37][38]. In these beam bond models, either a Euler-Bernoulli beam or a Timoshenko beam is assumed to link the centres of bonded particles. Carmona [35] studied the detailed development of the fragmentation processes of brittle agglomerates using a beam bond model that is based on an extension of 2D Euler-Bernoulli beam theory. André [9] used a Euler-Bernoulli beam bond model (referred to as EBBM for short thereafter) to study the micro-macro laws for homogeneous and isotropic materials. Obermayr et al. [39] proposed a beam bond model based on Timoshenko beam theory and demonstrated that the behaviour of the bond was equivalent to a linear finite-element Timoshenko beam element with reduced integration. Brown [40] developed a Timoshenko beam bond model (referred to as TBBM for short thereafter) that was able to produce the dynamics response of various structural elements such as simply supported beams, cantilever beams, multi-storey plane frames and thin plates. Further successful applications of the TBBM to more complicated processes such as impact loading of cementitious materials and loading of fibre reinforced polymers bonded to concrete have also been demonstrated [40,41].
Although bonded-particle models are increasingly used in the simulations of cementitious or sintered materials, there is limited effort on comparative studies of the advantages and disadvantages of different bonded contact models. The similarities and differences between the bonded contact models are unclear. Due to the differences in the mathematical expressions and model implementations, it can be confusing to the end-user whether a proposed model is a completely new model or just a slight variation of an existing model. In this paper, we identify and discuss two fundamental types of bonded contact models, namely, the multi-spring bond model and the beam bond model as shown schematically in Fig. 1. A series of verification cases were proposed to evaluate the efficacy of these commonly used bonded contact models and establish any discrepancy in the model predictions.

Description of DEM bonded contact model
The ability for the bonded particle models to produce a realistic representation of the mechanical behaviour of a cemented granular assembly mainly depends on the implementation of an inter-particle bonded contact model. This paper focuses on clarifying the micro-mechanics of the bonded contact behaviour before bond failure. The names of different bond models originate from the differences in their elastic components. The damping component and plastic deformation in the bond contact model will also be introduced. A bonded contact model typically has some defined criteria for bond failure which will not be discussed in this paper.

Spring bond model
The Spring bond model considers that two bonded particles are joined by a set of uniform elastic springs on the bond's cross-section. An example of the spring bond model is the parallel bond model in PFC3D. The parallel bond model in PFC3D incorporates a bonded contact component and a non-bonded contact component. When the two bonded particles are in physical contact, both bonded and non-bonded contact components are active in parallel. If two particles are bonded together but are not in physical The bonded contact component can be considered as a set of elastic springs uniformly distributed over the circular cross-section of the particle-particle contact plane [1,31]. These springs are centred at the contact point and have constant normal and shear stiffnesses. Since the bonded contact is activated over a finite area, this contact can resist both forces and moments.
The calculation of the forces and moments in a bond model is usually carried out in a local coordinate system [9,39,42]. Figure 2 shows a schematic of the contact description of the parallel bond model. The ̂ x axis of this local coordinate system is defined by the central axis of the bond, which joins the centres of the bonded particles. The other two local axes ̂ y ,̂ z lie normal to each other, as well as normal to the central axis. A transformation matrix is used to transform the contact vectors in a global coordinate system to the local coordinate system [31,38]. The contact point of the parallel bond model, c , is located at the centre of the interaction volume (gap or overlap) as shown in Fig. 2, which is calculated as, where g c is the length of contact gap or overlap, x is the centroid of particle , R is the radius of particle , L b is the length of the bond, between the centre of particles and . In order to generate a bonded contact at the bond initialisation time, a reference gap is specified by the user. A bonded contact is formed when the contact gap g c is less than or equal to the reference gap at the bond initialisation Note that a similar strategy of using the contact radius concept is implemented in EDEM [43] and YADE [44] to bond neighbouring particles that are not in physical contact, which will be discussed in more detail in the beam bond model section below. The translational velocity of particle at the contact point is calculated as, where U is the centroid translational velocity of particle , is the rotational velocity of particle and r c is the contact vector. The relative translational velocity ( v cr ) (Eq. 6) and rotational velocity (Eq. 7) at the contact point are given by where the subscript of cr represents the relative velocity at the contact point. The bond contact forces and moments in the parallel bond model, that can be calculated incrementally, are given in the Table 1 where k n , k s , k tor , k ben are the normal, shear, twisting and bending stiffness, respectively. It can be seen that the bond contact stiffness depends on the bond material properties E b , b , and bond geometry L b , r b , which are the Young's modulus E b , Poisson's ratio b , length L b and radius r b of the bond respectively, while is the ratio of normal stiffness to shear stiffness. The bond radius is calculated as the minimum radius of the particle contact pair times a constant multiplier. Note that the default input value of normal stiffness in PFC3D is an average normal stiffness over the cross-section of the bond, i.e. k n = E b ∕L b .
Each contact stores a force and moment that acts at the contact location in an equal and opposite direction on the two particle centres, while the shear forces at the contact point cause an additional torque at the particle centres that are not necessarily equal. Table 2 summarises the parallel  bond contact force and moment calculations at particle centres.
While the parallel bond model by Itasca (described above) might be the most common implementation, various changes to its formulation have been proposed by researchers to address some of its shortcomings or limitations. Brendel et al. [45] presented a visco-elastic bond model for caked spheres. In this model the elastic part of the contact is derived from the contact elasticity, which has the same formulation as the parallel bond model, as shown in Table 1 except that the ratio of normal to shear stiffness given by in which G b is the shear modulus of the bond. In addition to the elastic part, various authors have added a viscous damping component to the bond forces and moments with coefficients for critical damping. The formula of bond damping is analogous to contact damping and can be calculated as follows [45,46] and I is the reduced moment of inertia given by I = The viscous bond damping coefficient b determines the energy dissipation rate as the bond deforms. The bond damping coefficient is given the same value for all the forces and moments here for simplicity. In general, the damping coefficients in axial, shear, twisting and bending are not necessarily equal. Note that while the particle-particle contact component of the parallel bond model contains a damping component, this is only active when there is physical contact between the two 2k tor I cr particles. It does not apply if the particles are bonded but with a physical gap between the two particles. For static and quasi-static problems, a so-called "local damping" [31], "global damping" [40] or "background friction" [13] can be introduced to accelerate the system convergence to a steady-state solution. To avoid confusion with damping occurring at the contact, the term "particle damping" will be used hereafter instead of global, local or background friction as it is a clearer description of the type of damping. The opposing damping forces F ′′ and moments M ′′ are applied to each of the six degrees of freedom for each particle given as follows, where F i and M i are the sum of the force and moment acting on particle i , U and are the translational and angular velocity, and x , y and z represent each degree of freedom. b is a dimensionless coefficient that defines the magnitude of the particle damping.
Rojek et al. [47] proposed a modification of the stiffness to take into account the non-uniform cross-sectional area consisting of two segments. The segment area and length are proportional to the particle size in their method. Therefore, the calculated stiffness of a bond in the polydisperse case will be different from the parallel bond model in PFC3D. Shen et al. [48] also proposed a modification of the bond strength and stiffness to capture the concave end geometry of the cylindrical bond.
Twisting moment Bending moment

Beam bond model
Another class of DEM bonded particle models is the beam bond model [9,35,39,40,42,49] which proposes a beam element to connect the centres of two bonded particles, as shown in Fig. 3. In contrast to the parallel bond model, the beam bond model typically directly calculates the forces and moments at the two ends of the bond (particle centres) instead of the bond centre.
In Euler-Bernoulli beam theory it is assumed that under deflection of a beam, plane sections remain plane and normal to the neutral axis after deformation. This means that in long beams with very small deflections there will only be very minor (almost negligible) strains in the outermost (compressive and tensile) fibres of the beam which generate very small transverse shear forces along the longitudinal axis. However, in shorter beams or with larger deflections where the bending is more pronounced, the strains in the outer most fibres are significant and no longer negligible, and the plane section can no longer be considered to be perpendicular to the axis. It is this rotation of the plane and the resulting transverse shear forces that is considered by the more comprehensive Timoshenko beam theory.
Timoshenko beam theory is adopted in the formulation of Timoshenko beam bond model (TBBM) [38,40] because it accounts for the effects of transverse shear deformations in beams and is better suited than Euler-Bernoulli beam theory to study beams that are very short or undergo large deformations. In most cemented materials, where particle packings are dense, bonds tend to be short and stocky. Timoshenko beam theory will also correctly capture the behaviour of longer, more slender bonds. The TBBM model assumes the mechanical behaviour of a bonded contact under axial, shear, twisting and bending will follow the Timoshenko beam theory.
In the TBBM, a bond is formed when the contact radii of two particles touch or overlap at the bond formation time as shown in Fig. 4. The contact radius of a particle is defined as the product of its physical radius and a radius multiplier r . A contact radius multiplier of unity specifies that a contact only exists if there is physical contact between two particles. By increasing the value of r above 1, the contact radius becomes larger than the physical radius and a contact is formed between particles that are not physically touching. This will allow the bonding of particles to occur between particles not in direct contact, resulting in a denser matrix of bonds being formed at the bond formation time. The contact radius multiplier thus allows different degrees of cementation to give a better representation of a cementitious matrix by increasing the bond density when needed.
The forces and moments at particle centres are calculated in an incremental manner that is similar to the parallel bond model.
The incremental displacements and rotations at the bond ends (i.e. particle centres) are related to the incremental forces and moments using a stiffness matrix given as follows, where ΔF and Δ are the incremental forces (moments) and displacement (rotations) at the two ends of the beam.
The term ΔF contains forces (F) and moments (M) at the two ends of the bond, Δ contains translational displacements (d) and rotations ( ) that are calculated based on the translational velocity (U) and rotational velocity ( ) of particles and , respectively. The elements of the stiffness matrix, K (Eq. (14)), are given as follows [50]: The Timoshenko shear coefficient is a dimensionless quantity that is introduced to account for the shear stresses and strains no longer being uniformly distributed over the cross-section of the beam. The shear coefficient is dependent on the ratio of bond's radius to length and Poisson's ratio. If the Timoshenko shear coefficient Φ s becomes zero, the behaviour of the Timoshenko beam is reduced to an Euler-Bernoulli beam and the corresponding bonded particle model becomes the Euler-Bernoulli beam bond model (EBBM). When Φ s becomes zero, the shear stiffness of EBBM equals k 2 which is given by The EBBM was used as the cohesive beam model [9] and the cohesive discrete element method [37,51,52]. However, it should be noted that the word "cohesive" usually represents contacts that can be separated and reformed in DEM simulations as explained in the introduction. By contrast, a bonded contact can form once at an initialisation time step and its breakage is irreversible and cannot be reformed in DEM simulations.

Test cases to evaluate common DEM bond models
The bending behaviour of a cantilever beam under a point load is used as a reference case for a comparative evaluation of the different bond models described above as the preditions can be verified with the well founded solutions of the Euler-Bernoulli and Timoshenko beam theories. The cantilever beam, which is anchored at one end and free to deflect and rotate at the other, is modelled by bonding a number of spherical particles in a straight line to form the beam. An example of a beam with four particles is shown in Fig. 5.
Since the two ends of a single bond are fixed at the centres of two neighbouring particles, the relative movements of the particles result in the deformation of the bonds. Fig. 5 Schematic of the bending of a cantilever beam and bonded particle model. a physical model of cantilever beam b bonded particle model For the reference problem, a cantilever beam is modelled by bonding together a row of 11 particles. Each particle is just in physical contact with its neighbour particles; thus, the gap of the contact is zero. This case is similar to the verification example of parallel bond model (referred to as PBM for short) provided by Itasca for PFC3D® 5.0 [31] and similar cases have also been studied in the literature [39,49,53].
To model the bending of a cantilever beam, the extreme left particle is fixed against all translational and rotational displacement while the extreme right particle is subjected to a point load perpendicular to the axis of the beam. Gravity is not considered in the test cases. A quasi-static loading process is achieved by applying a small particle damping coefficient ( b = 0.001 ) to dissipate the kinetic energy so that the system comes to rest. The particle damping coefficient is used to make the system reach the steady state faster. A preliminary test has been carried out to confirm that it has a negligible effect on the steady state result, which is also verified by the excellent agreement between numerical simulation and analytical solutions of the cantilever beam test. The beam Young's modulus E b is 200 GPa and the Poisson's ratio b is 0.3. The applied external loading force F is 100 kN. The quasi-static loading is processed by increasing the load gradually at a rate of 100 kN/s for the first 1 s and then keep the load constant for the rest of simulation time. The system is assumed to reach the quasi-static state when the maximum particle velocity is below 10 -9 m/s.
In order to study the responses of the different bond models, a series of tests (cases 2-5) are proposed where the length and configuration of the cantilever are varied by changing the bond radius, the bond length and the number of particles/bonds representing the beam. As this study only focuses on the behaviour of the bond models, our analysis intentionally exclude the non-bonded particle contact model that kicks in after a bond breakage. It should be noted that the PBM model also includes the non-bonded contact before the bond breakage which is also excluded in this evaluation. The configuration parameters for the reference problem and four other test cases are listed in Table 3. A schematic of these simulation cases is shown in Fig. 6.
In Case 2, the effect of a beam's slenderness on its shear deformation is studied by varying the cantilever beam length while keeping the beam cross-section constant. In Case 3, the effect of bond length is investigated by varying the number of bonds used to form the cantilever beam whilst keeping the other parameters the same as the reference Case 1. Note that in Case 3, unlike the reference case, the particles are in not in physical contact-the gap between the particles is non-zero. This case is designed to study the sensitivity of the predictions to the beam resolution i.e. the number of particles representing the beam. In Case 4, both the number of bonds and the particle radius are varying. The particle radius is adjusted to fill in the beam so that each particle is still in physical contact with its neighbour particle. This case is to study if the gap or the number of bonds affects the predictions. Finally, for Case 5, the influence of poly-dispersity is studied using particles of different sizes to form the beam.
All the test cases are simulated using three bond models, i.e. the PBM, the EBBM and the TBBM as described in the section above. The simulation cases of the PBM were carried out using the commercial software PFC3D® 5.0 with the built-in parallel bond contact model [31], while EBBM and TBBM codes were implemented in the commercial software EDEM® using the API [43]. Meanwhile, we also implemented the PBM through the API of EDEM to double check our simulation results and the correctness of the formulation discussed in this work. The normal and shear strengths of the respective bond models are set to an extremely large value so that for all the cases the bonds cannot break. The normal, shear stiffness and frictional coefficient of non-bonded contacts is set to 0 in the PBM so that only stiffness in the bonded contact is active. The shear stiffness in the PBM is calculated using Eq. (8), which is based on the verification example of PFC3D® [31] and also adopted elsewhere [45,46,54] except in case 5 where the shear stiffness is intentionally adjusted as will be clarified later. Note that if the shear stiffness of the PBM is calculated using Eq. (8), then the twisting stiffness will be the same as the TBBM and the EBBM. Additionally, the default input value of stiffness in PFC3D is an average stiffness over the cross-section of the bond, i.e. k = k∕A b . The input stiffness will be multiplied by the cross-section area of the bond in the background [31]. The governing equations for a Timoshenko cantilever beam with a circular cross-section are as follows [55,56], where x is the distance from the fixed end, f s is the form factor for shear and equals 10/9 for a circular cross-section. By integrating the equations and substituting the boundary conditions, one can get the theoretical deformation given as follows, By substituting x = L 0 into the solution, the deflection at the beam tip can be determined as where the first term on the right-hand size is the bending contribution while the second term is the shear contribution which can be significant for short beams. If the second term is ignored, the solution reduces to be the theoretical solution of Euler-Bernoulli beam theory. Compared with Euler-Bernoulli beam theory, Timoshenko beam theory is known to be superior in predicting the response of beams, especially for thick beams. Therefore, the analytical solution of Timoshenko beam theory is used as the benchmark solution to compare with the numerical predictions in this paper. The relative error of the models to the Timoshenko beam theory analytical solution in Eq. (29) is used to evaluate the respective performance.

Case 1: reference problem
The schematic of the reference cantilever beam used in this case is shown in Fig. 7. The numerical predictions of the tip deflection for the PBM, EBBM and TBBM are shown in Table 4 along with the relative error from Timoshenko beam theory. All the three bond models produce satisfactory result in this case where the relative errors are all less than 1%. In particular, the EBBM prediction is very close to the Euler-Bernoulli theory analytical solution and TBBM prediction is close to the Timoshenko theory analytical solution. The final positions in the DEM simulations are used to determine the deflection at points along the beam and shown in Fig. 8a where the predictions by all the three models are also found to match the analytical values well. In general, the relative errors of the TBBM predictions are less than 1% while relative errors of PBM and EBBM predictions Fig. 7 Schematic of the simulated cantilever beam in the reference case

Case 2: slender and deep beams
According to the Eurocode 2 on the designing of concrete structure [57], any beam with a slenderness ratio less than 3 should be considered a deep beam. The slenderness ratio of a beam ( 0 ) is calculated as When a beam is classed as a deep beam, the Timoshenko beam theory is more applicable than Euler-Bernoulli theory. Since shear deformation is important in the case of short, deep beams, case 2 is designed to study the effect of this changing slenderness ratio in the various models. A series of simulations were conducted where the length of the beam is gradually decreasing as shown in Fig. 9. As the cantilever beam shortens, the slenderness ratio of the cantilever beam increases. Each beam is formed by bonding a row of particles with equal radius. The bond radius is set to be the same as the particle radius and each particle is in physical contact with its neighbour.
Since the bond radius does not change in these cases, shortening the cantilever beam length will lead to a decrease of the deformation under same loading force according to Eq. (29) The simulation results are compared in Fig. 10 where the relative displacement is calculated as the displacement at the tip normalized by the beam length. The reference case results appear in Fig. 10 as the smallest beam slenderness ratio, with λ 0 equal to 0.1. It is observed that all the  shortening the cantilever beam length three models follow the correct trend with increasing beam slenderness ratio. It is not surprising that TBBM predictions have close agreement with the analytical solutions of Timoshenko beam theory and that the EBBM predictions closely match the Euler-Bernoulli theorical solution. However, for deep beams, the Euler-Bernoulli theory predicts a lower relative displacement than the Timoshenko theory due to neglecting the shear deformation. At this point it is worth considering that the EBBM DEM models used in the literature such as André [9] may not be producing the correct behaviour of beam bending especially when the bond lengths are short, which is common in many DEM applications where densely bonded configurations are involved. The predictions of the PBM are observed to be in between Euler-Bernoulli theory and Timoshenko theory analytical solutions for this series of cases.

Case 3: effect of the discretization of the beam
Bonded-particle models are usually employed in the study of a cementitious material such as rock or concrete. For these cases, the particles in the bonded-particle models are typically considered as coarse-grained particles that are larger than the realistic constituent particle in the material. A critical issue of using Bonded-particle models to study the mechanical behaviour of a cementitious material is the need to choose an appropriate number of particles to represent the macroscopic fabric. This choice on the level of fidelity of the model will determine the level of detail and type of feature that can be studied with that model. As the number of particles in a simulation increases, the resolution of the DEM simulation and its ability to capture phenomena occurring at the microscale will increase. However, due to the computational cost associated with DEM simulations there is a compromise between using a sufficient number of particles to capture the bulk properties of the material and using large enough number of particles to study micro-scale phenomena such as cracking.
Simulation resolution can be described as the amount of detail that a simulation holds with a higher resolution meaning more detail is captured. Resolution of a DEM simulation can be defined in a number of ways, but is usually described in terms of the ratio of the particle size to the system feature of interest. In the case of beam bending, the resolution may be taken as the ratio of the length of an individual bond to the overall beam length. Therefore, in this case 3, the effect of the resolution of the beam model is studied by varying the number of particles (and bonds) used to represent the beam. The schematic of this series of cases is shown in Fig. 11.    Figure 12 compares the numerical predictions of the different bond models with the analytical solution of maximum displacement at the tip. An excellent agreement is achieved between TBBM predictions and the Timoshenko theoretical solutions, irrespective of the number of bonds used to form the cantilever beam. A similar trend is observed for the EBBM predictions whereas the PBM prediction is found to be significantly affected by the beam resolution in this simple test case. However, the relative error observed with the PBM prediction with respect to the Timoshenko theory reduces with increasing the number of bonds. Note that this dependency was also reported by Guo [53] who found that in order to reduce the error it was necessary to use 10 particles to form a fibre. Because refining the beam resolution will inevitably increase the computational cost, this series of cases indicates that the TBBM implementation is more attractive for simulating structural elements such as beams and flexible materials such as fibres or sheets.

Case 4: effect of gap length
As case 3 indicates that the PBM prediction is sensitive to the beam resolution, it is not clear if it is because of the decrease of the number of bonds in the beam or the increase of the gap between a bond pair of two neighbour particles. In case 4, we increase the particle radius when the number of particles is decreasing so that a bonded pair of particles remain in physical contact. Note that the bond radius is kept constant so that both the beam radius and beam length do not change. The schematic of this series of cases is shown in Fig. 13. Figure 14 shows the simulation results of this series of cases. The same trends as the previous case are observed for all the three models. Therefore, it is confirmed that the deterioration of PBM prediction is because of the lack of resolution-a certain resolution (number of bonds) is required to produce the theoretical result. The fundamental reason is the difference of shear stiffness calculation between the spring bond model and beam bond model, which will be further illuminated in the later discussions section. Note that the bond radius is intentionally set to be half of the reference case, which is to demonstrate that the conclusion will still hold when the bond radius is not the same as the particle radius.

Case 5: effect of polydispersity of bonded particles
This case is designed not only to test if one could improve the prediction of PBM with a low beam resolution, but also to investigate the different bond models' prediction for polydisperse cases. In this case, the bond shear stiffness is calculated using the Euler-Bernoulli beam theory in advance and is then used as the input for the parallel bond model. Because the cases studies here are pure bending loading, the twisting moments are zero. Therefore, one does not need to consider the difference of twisting stiffness caused by changing the shear stiffness in this test case. The schematic of the series of tests involved in this test cases is shown in Fig. 15. A beam is formed by three particles with the size of the middle particle varied. The polydispersity index of the simulation system is defined as the dispersion factor here, which can be calculated as follows,  where R max , R min and R avg are the maximum, minimum and average discrete particle radius values. Figure 16 shows the relative errors of numerical predictions of PBM, EBBM and TBBM for these three simulations. Note that the shear stiffness used in PBM is now calculated using the Euler-Bernoulli beam theory, i.e. Equation (26), and due to this, the predictions from the PBM and EBBM are almost identical for the monodisperse configuration. It can be seen that for the monodisperse case, the relative error of PBM is 0.62%, which is over 10 times smaller than the same two-bond layout in case 3 (see Figs. 12 and 14) because the shear stiffness has been specified correctly. However, with an increase in the size of the middle particle, the PBM prediction starts to deteriorate, even though the shear stiffness is already specified according to beam theory and remains unchanged. In contrast, the relative errors of the TBBM predictions are less than 0.1% for all the three simulations as shown in Fig. 16. This case highlights a crucial difference between the PBM and TBBM (or EBBM). Recalling Figures 2, 3 we can see that there is a difference in how the contact point location is calculated in the respective models. In the poly-disperse case, the contact point location is adjusted to be closer to the smaller particle in the PBM while it remains exactly in the centre of the beam connecting the two-bonded particles in the TBBM and EBBM. In the parallel bond model, the contact point shifts further away from the beam centroid as the size ratio of the particles increases, leading to increased errors. This case highlights a key limitation of the PBM where bond bending and twisting are important as it will only provide accurate results for monodisperse particle pairs.

Bond slenderness ratio effect
It has been shown above in the reference problem that all three bond models (i.e. PBM, EBBM and TBBM) can predict the bending behaviour of a thin cantilever beam consisting of a number of contacting particles. In the investigation of the effects of beam resolution (case 3) and gap length (case 4) the simulation results indicate that the predictions of the EBBM and TBBM are insensitive to the beam resolutions whereas the PBM can give rise to a large relative error if there is not a large enough number of constituent bonds in the beam. Therefore, the EBBM and TBBM can be considered to be independent of the number of bonds used. Another difference between the models is the shear stiffness as mentioned in Sect. 2 above. The ratio of the axial to shear stiffness adopted in the PBM, the EBBM and the TBBM models can be summarised as follows, where the bond slenderness ratio is defined as It can be observed that the axial to shear stiffness ratio of EBBM is a function of the bond's geometry alone whereas that of PBM is only a function of the Poisson's ratio of the bond but not the bond's geometry. The stiffness ratio of the TBBM, however, includes both the bond's geometry and the Poisson's ratio. André [9] reported that the bond Poisson's ratio has a negligible effect on the bulk scale behaviour of the material in EBBM simulation. This comparative study has shown that it is because of the lack of inclusion of the bond Poisson's ratio in the bond shear stiffness in EBBM. The bond Poisson's ratio only plays a role in the torsional stiffness in EBBM. 5 presents the values of axial stiffness and shear stiffness of each bond calculated in case 3 where different numbers of particles are used to discretise the beam (Fig. 12). As the length of bond decreases with increasing number of bonds in the beam, both the axial and shear stiffnesses increase for all three models. In the case of a single bond, the shear stiffness of the PBM is the largest while that of the TBBM is the smallest among these three models, which means that the deformation of TBBM under same point load would be largest. In particular, the shear stiffness of the PBM is around 52 times larger than that of TBBM for the lowest beam resolution simulation. With such a high shear stiffness, the deflection is significantly under-estimated in the PBM. As the number of the constituent bonds in the beam increases, the relative difference of the shear stiffness between PBM and TBBM decreases. When the number of bonds increases to 10, the shear stiffness of the parallel bond is 1.6 times larger than that of TBBM. The shear stiffness of the EBBM is 3.1 times larger than that of TBBM at this beam resolution. It is also interesting to note that the relative difference of the shear stiffness between EBBM and TBBM increases when the discretisation of the beam is increasing.
The ratio of axial to shear stiffness of the three bond models with increasing number of constituent particles representing the beam is shown in Fig. 17. As the number of constituent particles in the beam increases, the bond length decreases which leads to a fall in the axial to shear stiffness ratio for the EBBM and TBBM (Eqs. (33), (34) respectively). However, this stiffness ratio remains a constant in the PBM calculation (Eq. (32)). The difference between the PBM and TBBM (or EBBM) decreases with increasing the number of constituent particles in the beam, which explains the trend observed in case 3 and case 4. This constant stiffness ratio, due to a failure to include the bond geometry information in the bond shear stiffness, leads to the large errors in deflection seen at low beam resolutions (Fig. 12). In case 5a, we further showed that if the shear stiffness of PBM was adjusted using beam theory, the PBM can also correctly predict the bending deformation when only monodisperse particles are used. However, the predictions of the PBM for polydisperse particles still lead to a significant error as shown in case 5b and 5c. This is attributed to the difference of the intrinsic characteristic between a spring bond model and beam bond model. To be more specific, the implementations of the spring bond model and beam bond model are also different apart from the difference in shear stiffness calculations. It can be shown that each model uses a different method of calculating the contact point location for the forces and the moments. The PBM uses the contact point from the contact radii while the TBBM (or EBBM) always uses the centre point of the beam connecting the two particles. This difference was confirmed in case 5b and case 5c, where polydisperse particles were used with beams of the same length connecting them.

Beam slenderness ratio effect
It is also worth addressing the differences between the Timoshenko theory and Euler-Bernoulli theory since they are both used in beam bond models. In Case 2, where the effect of beam slenderness is investigated, a Fig. 17 The ratio of axial stiffness to shear stiffness for different bond models in case 3 difference between the TBBM and EBBM predictions has been shown for varying beam slenderness ratios (Fig. 10) -the Euler-Bernoulli and Timoshenko solutions diverge as the depth of the beam increases. The key advantage of the Timoshenko theory is that it takes into account the shear deformation, making it suitable for describing the behaviour of both thick and slender beams and vibration of beams under high frequencies. As shown in Fig. 18, the Timoshenko theory allows a rotation between the cross-section and the bending axis. This rotation comes from a shear deformation that is more noticeable for stocky beams. The Euler-Bernoulli theory assumes that the cross-section perpendicular to the neutral axis of the beam will remain both plane and perpendicular after a deflection. This assumption is known to be applicable for slender beams when the deflections are small compared to the depth [56,58]. According to Eq. (29), the Euler-Bernoulli theory neglects the shear contribution term which results in a stiffer beam. This behaviour can be observed even in the reference case (case 1). As shown in Fig. 8b, the displacement prediction of EBBM along the beam deteriorates as it approaches the fixed end. The load force is equal along the beam while the "effective beam length" becomes shorter as the position gets closer to the fixed end. This indicates that the point near the fixed end has a larger effective beam slenderness ratio, which makes EBBM prediction worse compared with the Timoshenko theory. It can be also confirmed through Eq. (28), as the distance from the fixed end ( x ) becomes smaller, the relative contribution of the second term (shear contribution) becomes more important. The relative error of Euler-Bernoulli theory compared with Timoshenko theory in the bending case studied here is related to the slenderness ratio and Poisson's ratio of the beam. The theoretical relative error of Euler-Bernoulli theory can be calculated as follows, Figure 19 presents the analytical difference between Euler-Bernoulli and Timoshenko theory for different slenderness ratios. It can be seen that the maximum relative error can be as high as 38% when the beam diameter equals the beam length. Recall that the shear bond stiffness in EBBM is always larger than TBBM (Table 5), which results in an under prediction of the beam deformation. Additionally, the damage caused by shear is an important failure mechanism for cementitious materials [59,60]. Therefore, for stocky beams such as those concerned in the majority of DEM of bonded materials, the TBBM should be superior to EBBM for the failure simulations of the cementitious materials.

General formulation of bond models
In order to better understand the differences between the parallel bond model and the beam bond models, it is first required to convert the beam bond models into a similar   14.88 format to that of the parallel bond model where force and moments at the contact point and at the particle centroids are presented in Tables1 and 2 respectively. Here we start with the Timoshenko beam bond model. Substituting the stiffness matrix into Eqs. (13) and (14), one can get the following relations after some re-arrangements: As shown in Fig. 3, we define that the contact point is located at the centroid of the beam, i.e., Combined with Eqs. (6) and (7), both the force and moment relations at the contact point can now be rearranged to give the concise forms outlined in Table 6. Furthermore, the force and moment relations calculated at the centroid of each particle in the bonded pair can also be rewritten and are given in Table 7 [61].
Comparing Tables 1 with 7, it becomes evident that whilst the formulations of the PBM and TBBM bond models appear quite different in the literature, they are actually of a similar formulation, with the key exception being, the definition of the contact point location for each model. The parallel bond model defines that the contact point is centred at the gap or overlap of the bonded particle pair while the contact point of the TBBM can be considered as located at the centroid of the beam. By substituting Eqs. (2) into (1), it can be shown that the contact points for the parallel bond model and TBBM are coincident when the two bonded particles have the same physical and contact radii (Test cases 1-4). However, if the sizes of bonded particles are different, the contact point locations for these two models are different (see Test case 5). The contact point location is of major importance because it will affect the calculation of relative velocity, displacement and lever arms for the bending moments. The difference lies in the original assumptions of the model. Because the TBBM focuses on the mechanical behaviour of the beam that connects the bonded particles, the contact point is considered to be located at the centroid of beam. In contrast for the PBM, the contact point location will adjust to be closer to the smaller particle if the sizes of bonded particles are different. In this sense, the TBBM can be more suitable for modelling continuous materials since it has less dependence on the constituent particle size chosen in the discretisation of the domain.
The other significant difference between spring bond model and beam bond model is the calculation method for Twisting moment ΔM cx = k tor cr,x Δt k tor = 0.25k n r 2 b ∕ 1 + b Bending moment ΔM cy = k ben cr,y Δt k ben = 0.25k n r 2 b ΔM cz = k ben cr,z Δt the shear stiffness, which in turn will also affect the calculation of torsional stiffness in the model. The shear stiffness in a spring bond model is typically calculated as a ratio of the normal to shear stiffness, as shown in Table 1. This ratio is often defined with respect to the bond material properties and according to Equation (8) [45,46,54]. The shear and torsional stiffnesses of TBBM are derived from the Timoshenko beam theory which takes into account the shear deformation of a beam under loading. Obermayr [39] developed another beam bond model and showed that it was equivalent to a linear finite-element Timoshenko beam element with reduced integration of the shear deformation for small deformations. Absolute displacements and rotations using quaternion algebra framework were adopted in the force and moment calculations in their proposed model. However, we find that the elastic part of the force and moment calculations can also be written in a similar form to Tables 6 and 7  where s is a simplified shear coefficient. At this point it should be noted the Timoshenko beam model proposed by Obermayr [39] is not resolution independent as there is a significant difference in the deflections calculated between a beam consisting of 2 particles and 10 particles. While they compared a 10-particle solution to the solution using ANSYS BEAM188 finite element method and found a relatively close agreement, no comparison was made with the actual Timoshenko theory.
In summary, all the aforementioned bond models can be unified in a generic form as shown in Table 8, provided that the relative translational velocity and rotational velocity at each contact point location are calculated using Eqs. (6) and (7) respectively. The differences between the bond models lie in the definition of contact point location and contact stiffness calculations.
Following the unification of the different bond models into the generic formulation, the key aspects (contact point location and bond stiffness calculations) of these bond models can be summarised as shown in Table 9. There are several important conclusions to be drawn from this. The first is that while the parallel bond model has a different definition of the contact point, it will collapse to the same definition as both beam bond models if monodisperse particles are being used. With this in mind, it is possible to then compare the various model stiffnesses to see where the difference lies.
All three models have been found to have the same normal and bending stiffness while the main difference between all the models is the shear stiffness. The parallel bond model has the simplest shear stiffness model which is simply related to the normal stiffness, a method that is typically employed in the cohesionless and cohesive contact model families. The Euler-Bernoulli and Timoshenko models use a definition that includes both the bond geometry (radius and length) and the material properties (Young's Modulus). The Timoshenko shear stiffness is modified by the Timoshenko shear coefficient, which was described in Eq. (25), to account for the shear deformations in deep beams. The twisting stiffness is the same for both beam models and is again dependent on geometric definition (bond radius) and material property (Poisson's ratio), whereas for the parallel bond model, it is again simply related to the normal stiffness. This difference is the modification proposed by Brendel et al. (Eq. (8)).

Conclusions
In this paper, a detailed assessment of the common bond models used in DEM bonded particle simulations is presented. Two fundamental types of bonded contact models Twisting moment ΔM cx = k tor cr,x Δt ΔM x = −ΔM x = −ΔM cx Bending moment ΔM cy = k ben cr,y Δt ΔM y = −ΔM cy − r c ΔF z ΔM cz = k ben cr,z Δt ΔM y = ΔM cy + r c ΔF z ΔM z = −ΔM cz + r c ΔF y ΔM z = ΔM cz − r c ΔF y Contact point Twisting stiffness k tor = 0.5k s r 2 b k tor = 0.25k n r 2 b ∕ 1 + b Bending stiffness k ben = 0.25k n r 2 b are identified, namely the spring bond model and the beam bond model. The spring bond model calculates the bond forces and moments at the contact point location and then transfer them to the centres of the bonded particles whereas the beam bond model directly calculates the bond forces and moments acting at the centres of the bonded particles based on beam theory. A series of cantilever beam bending cases have been carried out to evaluate the performance of the three typical bond models, namely, the parallel bond model (PBM), the Euler-Bernoullli beam bond model (EBBM) and the Timoshenko beam bond model (TBBM). It is found that all the three models are capable of quantitatively predicting the bending behaviour of a slender cantilever beam. However, the PBM is very sensitive to the number of constituent particles used to represent a beam and could only predict the deflection if the number of constituent particles is large enough-at least 10 bonds are required to predict the beam response to the same results as the beam bond models. As expected, the EBBM correctly predicts the deformation for slender beams but underpredicts the deformation for deep beams due to neglecting the shear deformation contribution. The TBBM is found to have excellent predictions of beam deformations for all the studied cases due to the rigorous nature of the Timoshenko Theory upon which it is based.
A generic formulation is presented for the different bond models and used to assess the key differences between the models. It is found that the contact point definition varies between the parallel bond model and the beam bond models, however, in the case of monodisperse particles the contact point could be shown to be the same and a direct comparison between the models could be made for this configuration. Several similarities and differences were found in the various stiffness components of the models. All three models share the same normal and bending stiffness while the torsional stiffness and shear stiffness differ. The Euler-Bernoulli and Timoshenko modes share the same definition of torsional stiffness. Finally, it is worth mentioning that all the bond models can provide reasonably good results with right sets of input parameters, i.e., by providing enough resolution and considering the right beam slenderness ratios that are likely to be encountered in the study phenomena. Therefore, appropriate calibration, verification and validation are vital before applying DEM in solving real-world problems [62].
Future study on the use of these models with bonded particle models to study complicated cases like failure of concrete or rock needs to be carried out to understand the effects of the fundamental difference in these models to their ability to predict real-world cases.