The Sensitivity of Micro—Macro Mechanical Behaviour of Sand to the Inter-Particle Properties

: Sand is a particulate material but is treated as a continuum solid in some engineering analyses. This approach is proven to be acceptable when dealing with geotechnical structures, provided an adequate factor of safety is applied so that there is no risk of failure. However, the continuum approach does not account for the effect of interparticle forces on the micro–macro behaviour of sand. Sand could be modelled as a particulate material using the discrete element method (DEM), taking into account its discrete nature. This paper shows how the microscopic contact properties between the idealised sand particles inﬂuence the macro-mechanical behaviour, highlighting the development of the fabric as the soil approaches failure. Thirty DEM biaxial tests were performed to study the sensitivity of the macro–micro mechanical properties of sand to the inter-particle properties of an idealised sand particle. The conditions of these simulations were the same (e.g., particle size distribution, number of particles, porosity after radius enlargement, boundary conditions, and rate of loading). The sensitivity of the pre-peak, peak, and post-peak behaviour of these simulations to the inter-particle properties of an idealised sand particle was studied. Two extra DEM biaxial tests under different conﬁning pressures were performed to verify the cohesionless nature of the synthetic material used for this study. Since a two-dimensional DEM is used for this study, a detailed approach to interpret the results assuming either a plane strain or a plane stress situation was discussed. This study highlighted the critical inter-particle properties and the range over which these inﬂuence macro-mechanical behaviour. The results show that Young’s modulus is mainly dependent on the normal contact stiffness, and peak stress and the angle of internal friction are greatly dependent on the inter-particle coefﬁcient of friction, while Poisson’s ratio and volumetric behaviour of particulate sand are dictated mainly by shear contact stiffness. A set of relationships were established between inter-particle properties and macro-machinal parameters such as Young’s modulus, Poisson’s ratio, and angle of internal friction. The elastoplastic parameters obtained from these tests are qualitatively in agreement with the typical medium and dense sand behaviour.


Introduction
The granular soil, such as sand, is consistently used in backfilling technologies, in particular in hydraulic backfilling, for which its compressibility is determined. The compressibility characteristics of sand are generally obtained from laboratory tests. These characterisations are then used in continuum models for geotechnical problems such as analyses of shallow and deep foundations, excavations, retaining walls, and slope stability. However, these models cannot account for the effect of microscopic inter-particle properties on the macroscale response of granular materials [1][2][3][4][5]. The bulk behaviour of sands is different from usual solid materials during shearing and presents a number of micro-mechanical properties of sand to the inter-particle properties are not yet explored comprehensively in the literature. This paper aims to study this gap using DEM by conducting a general framework for DEM simulation considering the 2D disk particle shape.

Two-Dimensional Simulations in DEM
The three-dimensional solution of a problem in continuum mechanics can be approximated by a two-dimensional solution (i.e., plane strain and plane stress) depending on the geometry of that problem, meaning the components of stress and strain tensors can be reduced to three components. To calculate stress and strain tensors within the particulate soil, an averaging method is used to estimate the average stress and strain tensors. That is, to compute the stress and strain tensors for a volume of particles within a two-dimensional DEM model using PFC 2D , two in-plane force and displacement components and one in-plane moment are required. The out-of-plane force components and out-of-plane particle displacement are not taken into account in the motion equation-i.e., the out-ofplane forces and displacements which are essential to enforcing a state of plane strain or plane stress are not present. Therefore, the interpretation of 2D DEM results in terms of either plane strain or plane stress will be a controversial issue. Figure 1 displays a flowchart on two-dimensional used in DEM modellings. Section 3 provides a more detailed explanation regarding the out-of-plane force component and out-of-plane particle displacement. The two-dimensional assumption, however, has an advantage. The dynamic response of a particulate system is greatly dependent on the number of degrees of freedom of each particle within these systems. In 3D DEM simulations, each idealised particle has six degrees of freedom, while in 2D cases, there are three degrees of freedom per particle. The computational effort in 2D DEM simulations will be less than in 3D simulations and, therefore, faster. Furthermore, the number of documented 2D DEM studies published annually shows that 2D DEM simulations are Section 3 provides a more detailed explanation regarding the out-of-plane force component and out-of-plane particle displacement. The two-dimensional assumption, however, has an advantage. The dynamic response of a particulate system is greatly dependent on the number of degrees of freedom of each particle within these systems. In 3D DEM simulations, each idealised particle has six degrees of freedom, while in 2D cases, there are three degrees of freedom per particle. The computational effort in 2D DEM simulations will be less than in 3D simulations and, therefore, faster. Furthermore, the number of documented 2D DEM studies published annually shows that 2D DEM simulations are able to capture the key complex mechanical response features of soil medium [34][35][36]. Table 1 summarises the advantages and limitations of the use of 2D and 3D DEM modelling. Table 1. Advantages and limitations in the use of 2D and 3D modelling.

Two-Dimensional DEM Model
Three-Dimensional DEM Model

Advantages Limitations Advantages Limitations
Less number of particles are required to be generated to study a problem Three degrees of freedom are used in the calculation cycles (two in-plane displacements and one in-plane moment) Six degrees of freedom are used in the calculation cycles Greater numbers of particles are required to be generated to study a problem Cost of calculations is relatively less The impact of out-of-plane forces and moments are not considered in the macro-micro responses The impact of out-of-plane forces and moments are considered on the macro-micro responses The cost of calculations is relatively high. Since DEM tracks each individual particle and its interactions over time, an increase in the number of particles increases the computational time Suitable for 2D plane strain and plane stress problems The out-of-plane forces, which are essential to enforcing a state of the plane strain and plane stress, are not present.

Produces more realistic macro-micro responses
The application of 3D DEM to produce the deformable boundary conditions is complicated Fewer cycles are required for particulate systems to reach equilibrium The application of 2D DEM to produce the deformable boundary conditions is complicated More cycles are required for particulate systems to reach equilibrium Enabling clear visualisation as particle motion is restricted to one plane Less clear visualisation as particle motion is restricted to one plane More cost-effective for the models subjected to dynamic loading due to computational capabilities Less cost-effective for the models subjected to dynamic loading due to limited computational capabilities

Plane Strain and Plane Stress Behaviour in DEM Simulations
For most geotechnical systems such as tunnels, retaining walls, earth dams, strip foundations, deep excavation, and slopes, soil behaviour is assumed to be dependent on the vector displacement field in two dimensions, and the effect of the vector displacement field in the third dimension is not very evident on the behaviour. This situation in continuum mechanics is called the plane strain, where a three-dimensional problem is analysed as a two-dimensional problem where values of the strain tensor components in the out-of-plane dimension are set to zero (i.e., ε 22 = τ 21 = τ 23 = 0). In some cases, in continuum mechanics, the stress tensor field is also two-dimensional. In this case, the stress in the out-of-plane dimension is the intermediate stress. This situation is called plane stress. This does not apply in DEM analyses because no stress and strain tensor exist in such models. Instead, in two dimensions, only two in-plane force components and one in-plane moment are present (see Figure 1d). Section 4.1 presents the set of exemplary calculation formulas to calculate the two in-plane force components and one in-plane moment. In the averaging method, the average stress and strain tensors are computed. However, only the in-plane forces and displacements are used to calculate the average stress and strain tensors. The out-of-plane forces and displacements are not applied in calculating the average stress and strain tensors. Therefore, the out-of-plane constraint, which is essential to enforce a state of plane strain or plane stress, cannot be present. The formulations that describe plane stress and plane strain conditions will be explained, and the macro parameters (i.e., Elastic modulus, Poisson's ratio, and angle of internal friction) obtained from these two situations are compared in order to assess the differences between both interpretations.

DEM Implementation in PFC 2D
DEM is an advanced numerical tool originally proposed by Cundall and Strack [37]. It dynamically simulates and tracks the micro-macro-scale behaviour of granular materials subjected to quasi-static and dynamic loads. The central explicit finite difference scheme is used in this method to solve the dynamic equilibrium of each particle at each time step. In this paper, the DEM simulations are implemented using the software PFC 2D 4.1 [38]. The rolling resistance is not adopted in this study. It has been experimentally observed by O'Sullivan et al. [39] that if the size of a particle is bigger than 0.1 (mm), the size of particle asperity will be negligible in comparison to the particle inertia. Therefore, the surface roughness will have a minor effect on the material behaviour in comparison to the particle inertia.

Contact Law
The contact force (F i ) applied on a disk particle in PFC 2D code is decomposed into normal force (F n i ) and tangential force (F s i ). The former is directed toward the particle centre, and the latter is directed along the tangent to the particle.
To preserve the geometry of particles during loading, it is assumed in PFC 2D that particles are rigid with soft contact, which means that a contact overlap between two particles (e.g., a and b) is applied rather than a contact deformation. The magnitude of this overlap is computed by a contact law. In the present study, a linear elastic contact law is applied to calculate the components of contact forces. The normal and tangential displacements at time step ∆ t are calculated as follows: where . are scalar translational and rotational particle velocities of particles a and b, respectively. R a and R b are the particle radius. n and l are the normal and tangential unit vectors.
The magnitude of the normal and tangential contact forces is calculated via the following: The total contact shear force is compared to the Coulomb sliding friction or sliding capacity criterion (i.e., µ F n i t ) to check whether sliding has occurred. When the resultant force and torque in the z-direction (calculated by multiplying tangential contact force by the distance from the particle centre to the contact location) are computed for each particle, the local damping force F d i will be added to them: where α, |F| and .
x are damping constant, resultant force on the particle, and particle velocity, respectively. The computed resultant force and torque acting on the particle are used to determine the change in particle velocity via Newton's second law for the next time step. Once the resultant force is calculated for a particle, a damping force will be added to this at the end of each time step to reduce the kinetic energy of the particle. The damping force is calculated by multiplying a damping constant by the resultant force. The direction of the damping force is opposite to the particle velocity. The source of Equations (1)- (7) is [37,38].

Soil Fabric
In granular mechanics, soil fabric refers to the size, shape, and arrangement of soil particles. Fabric quantities include either particle orientation if 2D disk particles are used, or contact orientation and branch vector orientation for non-circular particles. These fabric quantities are often presented as an average or graphically (e.g., polar diagram of contacts). Rothenburg et al. [6] proposed a closed-form solution to estimate the polar diagram of contacts.
where a represents "fabric anisotropy" in a granular system, depending on the number and density of unit normal vectors in principles axes. E(θ) shows the deviation between the geometry of contact force distribution and the isotropic geometrical contact force distribution. For example, if a = 0, E(θ) will be a circle such that the state of the system being considered is in an isotropic state. θ a is the direction of anisotropy. Parameters, a and θ a are obtained by the following equations: θ a = 1 2 tan −1 ∑ n g g=1 N g sin((2g − 1)∆θ ∑ n g g=1 N g cos((2g − 1)∆θ (10) where N is the total number of contact, ∆θ = 360 n g , n g the number of segments and N g is the number of contacts within the gth segment. In fact, the fabric anisotropy parameter shows the ability of granular systems to create the anisotropy state in normal contact distribution.
One of the key microscopic parameters, which are defined at the particle level, is the average coordination number which increases with densification [40,41]. This parameter is the average number of contacts per particle within a specific volume of a particulate assembly, and consequently, it provides a measure of packing density or packing intensity of fabric at the particle level. For a volume of particulate assembly with N p particles and total number of contacts, N c , the definition of average coordination number C n is given by the following: Since each contact is shared between two particles, the actual number of contacts is multiplied by 2. Rothenburg and Kruyt [10] and Maeda [42] have shown that the average coordination number should be at least three for idolised disk particles when a granular system is in quasi-static equilibrium.

The Typical Behaviour of Dry Sand
As schematically shown in Figure 2, the general macro-mechanical behaviour of sand subjected to the static deviatoric loading in a standard triaxial test is characterised by the following: The slope of the volumetric strain vs. axial strain curve at a strain corresponding to half of the peak stress is used to compute Poisson's ratio, ν, for both plane strain and plane stress.
In the case of dense and medium sand, the characteristic point, , corresponds to = 0, where the dilation of the sample starts. The shear strain corresponding to this point varies between sands, as presented by Atkinson [14]. Small strains: Within this range, the stress-strain curve is assumed to be linear. In the case of dense and medium sand, the peak stress ratio or deviatoric stress occurs within this range.
Large strains: these correspond to strains exhibited when the soil is approaching failure, and the shear stiffness becomes small. Figure 2. The typical behaviour of medium and dense sand (after Momeni [15]). Figure 2. The typical behaviour of medium and dense sand (after Momeni [15]).
The initial Young's modulus, E 0 (the value of Young's modulus at very small shear strains), is an important characteristic of soil deformability and plays an important role in the dynamic and static response of the soil. This is usually estimated by using laboratory or field tests that are related to seismic wave propagation and stiffness degradation curves based on cyclic tests conducted in the laboratory, which are costly, as shown by Okur and Ansal [43]. Where E 0 is not available, E 50 , the secant modulus at 50% of peak stress, is often used to predict ground movements, as presented by Holtz et al. [44].
The slope of the volumetric strain vs. axial strain curve at a strain corresponding to half of the peak stress is used to compute Poisson's ratio, ν , for both plane strain and plane stress.
In the case of dense and medium sand, the characteristic point, G, corresponds to ∂ε v ∂ε 11 = 0, where the dilation of the sample starts. The shear strain corresponding to this point varies between sands, as presented by Atkinson [14]. Small strains: Within this range, the stress-strain curve is assumed to be linear. In the case of dense and medium sand, the peak stress ratio or deviatoric stress occurs within this range.
Large strains: these correspond to strains exhibited when the soil is approaching failure, and the shear stiffness becomes small.

Elastic Parameters
In soil mechanics, it is assumed that at the start of a triaxial test, the material is linear elastic. In the principal stress space, the behaviour is as follows: In terms of plane strain, the stress-strain relation is as follows: In terms of plane stress, the stress-strain relation is as follows: The elastic modulus and Poisson's ratio are important characteristics when predicting ground movements, as it is often assumed that soil is isotropic and homogenous and behaves elastically. There are numerous methods available to determine these properties, but the most common is the triaxial test. In this test, using local strain measurements, it is possible to measure the stress-strain response of a sample of soil subjected to a variety of load paths (e.g., Head [45]). The values of stiffness obtained are stress path dependent and vary with the strain range over which they are measured. Stiffness is also a function of the confining stress, the particle geometry, and the density of packing. In this study, only monotonic compressive loading is monitored. Furthermore, the triaxial test is three-dimensional, whereas the DEM analysis used in this study is two-dimensional. However, since the soil is assumed to be isotropic and homogenous, the material properties obtained from the triaxial test can be applied to the two-dimensional analysis. In twodimensional studies, it is necessary to consider plane stress or plane strain conditions, which lead to small differences in values of stiffness.
In plane strain situations, the secant Elastic modulus is as follows: In a plane stress situation, the Elastic modulus is as follows: When Poisson's ratio is zero, the elastic modulus is equal for the plane stress and plane strain situations. Additionally, when Poisson's ratio is 0.5, the elastic modulus in the plane strain is 75% of the elastic modulus in the plane stress [46].
Poisson's ratio is obtained from the slope of the horizontal strain vs. axial strain curve. In two-dimensional analysis, the volumetric strain is ε v = ε 11 + ε 33 . In the plane stress condition, Poisson's ratio is obtained from the following equation: where ε 33 is the horizontal strain and ε 11 is the axial strain. In a plane strain situation, it is as follows:

Plastic Parameters
Angles of friction, θ, are related to the plastic and failure behaviour of non-cohesive material. For non-cohesive sand, the following relationship can be used to compute the angle of friction: t s = Tan 2 45 + θ 2 (19) where t is the deviatoric stress and s is the isotropic stress.

The Methodology of Sensitivity Analysis
Three inter-particle properties, including normal and shear stiffness and inter-particle coefficient of friction, were considered in this research. Shear contact stiffness is presented as a ratio to normal stiffness. Values of the inter-particle properties of sand, including friction, normal, and shear contact stiffness, are listed in Figure 3 and Table 2.
where t is the deviatoric stress and s is the isotropic stress.

The Methodology of Sensitivity Analysis
Three inter-particle properties, including normal and shear stiffness and inte cle coefficient of friction, were considered in this research. Shear contact stiffness sented as a ratio to normal stiffness. Values of the inter-particle properties of sand, ing friction, normal, and shear contact stiffness, are listed in Figure 3 and Table 2. . The values of inter-particle friction for quartz sand (redrawn after Lambe and W [29]). It was necessary to establish a methodology to measure the effect of interproperties on the micro-macro mechanical behaviour of particulate idealised sa pressed in terms of the angle of internal friction and stiffness of the sample. Each inter-particle properties was varied, keeping the others constant to determine the on the macro mechanical properties. The inter-particle properties for each biaxial obtained from Table 2. Table 3 summarises the input data used for this study. F presents the particle size distribution profile used for this study. . The values of inter-particle friction for quartz sand (redrawn after Lambe and Whitman [29]). Table 2. The micro-mechanical parameters used for the sensitivity analysis (after Momeni [15]).  It was necessary to establish a methodology to measure the effect of inter-particle properties on the micro-macro mechanical behaviour of particulate idealised sand, expressed in terms of the angle of internal friction and stiffness of the sample. Each of the inter-particle properties was varied, keeping the others constant to determine the impact on the macro mechanical properties. The inter-particle properties for each biaxial test are obtained from Table 2. Table 3 summarises the input data used for this study. Figure 4 presents the particle size distribution profile used for this study. Range of Particle Size Distribution PSD (mm) 0.5-3 (refer to Figure 4) k n (N/m) (Normal contact stiffness) Refer to Table 2 k s (N/m) (Shear contact stiffness) Refer to Table 2 µ (Inter-particle coefficient of friction) Refer to Table 2 α (Damping constant) 0.8   Porosity after radius enlargement 0.12 Range of Particle Size Distribution PSD (mm) 0.5-3 (refer to Figure 4) (N/m) (Normal contact stiffness) Refer to Thirty biaxial tests in total were conducted to investigate the sensitivity of inter-particle properties on the micro-macro mechanical behaviour. The initial condition, such as the initial geometry of the biaxial chamber, particle size distribution, particle shape, contact model, porosity and isotropic stress state condition, and the lateral boundary condition for all these tests, were similar. Initially, non-overlapping disk particles in the range of 0.5 and 3 (mm) corresponding to the well grade of sand were randomly placed at half their size within four enclosed rigid walls. Their radii were then expanded to reach a porosity of 0.12. Next, the rigid boundaries of the biaxial cell were moved based on the applied strain control to approach the stress at the boundaries at 100 (kPa). Once the sample was isotopically consolidated to 100 (kPa), further cycles were needed to reach system Thirty biaxial tests in total were conducted to investigate the sensitivity of inter-particle properties on the micro-macro mechanical behaviour. The initial condition, such as the initial geometry of the biaxial chamber, particle size distribution, particle shape, contact model, porosity and isotropic stress state condition, and the lateral boundary condition for all these tests, were similar. Initially, non-overlapping disk particles in the range of 0.5 and 3 (mm) corresponding to the well grade of sand were randomly placed at half their size within four enclosed rigid walls. Their radii were then expanded to reach a porosity of 0.12. Next, the rigid boundaries of the biaxial cell were moved based on the applied strain control to approach the stress at the boundaries at 100 (kPa). Once the sample was isotopically consolidated to 100 (kPa), further cycles were needed to reach system equilibrium. It is to be noted that the porosity 0.12 mentioned in Table 2 is the porosity after radius enlargement. Next, the confining pressure on the vertical rigid boundaries was kept constant while the top and bottom rigid boundaries moved towards each other to apply deviatoric stress. The strain rate (i.e., . ε) applied for this test was 2% min , such that the incremental acceleration of each particle at each time step is small. All the imposed energy generated during the simulation was dissipated through both frictional sliding between particles and loss of contacts. In the sensitivity analysis, the critical parameters and the range over which the parameters impact the macro-mechanical behaviour are investigated. From each combination of normal stiffness, shear stiffness, and inter-particle coefficient of friction (i.e., k n , k s , and µ), the macro-mechanical parameters comprising Young's modulus, Poisson's ratio, and angle of internal friction were computed. These values will then be compared with typical values of elastic modulus, Poisson's ratio, and angle of internal friction of sand obtained from the literature (see Table 4). To investigate the sensitivity of the micro-macro-mechanical behaviour of particulate sand to the inter-particle coefficient friction, the ratio between normal contact stiffness and shear contact stiffness is assumed to be unity. By doing this, the number of biaxial tests required to be implemented was reduced to fifteen tests. The inter-particle properties used for these tests are listed in Table 5. Figure 5a illustrates the variations of deviatoric stress with axial strain for different inter-particle coefficients of friction where normal and shear contact stiffness is 8.45 × 10 7 (N/m). As seen, for a fixed value of contact stiffness, the peak stress is profoundly controlled by inter-particle frictions, while the axial strain corresponding to the peak stress together with the slope of the stress-strain is less influenced by inter-particle frictions. Figure 5b shows that the contraction behaviour of particulate sand is less impacted by inter-particle frictions, though the model with bigger inter-particle coefficient friction (i.e., 1.2) shows slightly more dilative behaviour. Table 5. The input data for sensitivity of sand system to the inter-particle coefficient of friction.  Figure 6a summarises the sensitivity of peak stress to the various inter-particle coefficient of friction for all these fifteen tests. As seen for each fixed value of contact stiffness, peak stress is significantly governed by inter-particle frictions. For example, for a fixed contact stiffness value of 8.45 × 10 7 (N/m), changing the inter-particle coefficient of friction from 0.5 to 0.9 and 1.2 results in an increase in peak stress by 37% and 53%, respectively, resulting in an increase in the angle of internal friction from 20 • to 26 • and 29 • , for the inter-particle coefficient of friction of 0.5 to 0.9 and 1.2, respectively (see Figure 6b). The angle of friction between 28 • and 37 • is typical for medium and dense sand. Based on the data provided for this study, a relationship between the inter-particle coefficient of friction and the angle of friction can be developed as follows: where θ and µ are the angle of friction and inter-particle coefficient of friction. Note that this relationship is developed for inter-particle of coefficient between 0.5 and 1.2. It can be concluded that inter-particle friction mainly controls the peak stress and internal angle of friction.
shear contact stiffness is 8.45 × 10 7 (N/m). As seen, for a fixed value of contact stiffness, the peak stress is profoundly controlled by inter-particle frictions, while the axial strain corresponding to the peak stress together with the slope of the stress-strain is less influenced by inter-particle frictions. Figure 5b shows that the contraction behaviour of particulate sand is less impacted by inter-particle frictions, though the model with bigger inter-particle coefficient friction (i.e., 1.2) shows slightly more dilative behaviour.
(a) (b) Figure 5. The sensitivity of macro-mechanical behaviour of sand to inter-particle coefficient of frictions: (a) deviatoric stress with axial strain and (b) volumetric strain with axial strain where normal and shear contact stiffness is 8.45 × 10 7 (N/m). Table 5. The input data for sensitivity of sand system to the inter-particle coefficient of friction. Geotechnics 2023, 3 427 Figure 6a summarises the sensitivity of peak stress to the various inter-particle coefficient of friction for all these fifteen tests. As seen for each fixed value of contact stiffness, peak stress is significantly governed by inter-particle frictions. For example, for a fixed contact stiffness value of 8.45 × 10 7 (N/m), changing the inter-particle coefficient of friction from 0.5 to 0.9 and 1.2 results in an increase in peak stress by 37% and 53%, respectively, resulting in an increase in the angle of internal friction from 20° to 26° and 29°, for the inter-particle coefficient of friction of 0.5 to 0.9 and 1.2, respectively (see Figure 6b). The angle of friction between 28° and 37° is typical for medium and dense sand. Based on the data provided for this study, a relationship between the inter-particle coefficient of friction and the angle of friction can be developed as follows: (20) where and are the angle of friction and inter-particle coefficient of friction. Note that this relationship is developed for inter-particle of coefficient between 0.5 and 1.2. It can be concluded that inter-particle friction mainly controls the peak stress and internal angle of friction.

kn = ks (N/m) Inter-Particle Coefficient of Friction
(a) (b) Figure 6. The sensitivity analysis: (a) inter-particle coefficient of friction vs. peak stress and (b) interparticle coefficient of friction vs. peak stress. Table 6 summarises E50 with an inter-particle coefficient of friction for both plane stress and plane stress conditions, indicating that Young's modulus is significantly governed by contact stiffness rather than the inter-particle coefficient of friction. For instance, changing the inter-particle coefficient of friction from 0.5 to 0.9 and 1.2 for a fixed contact stiffness value of 8.45 × 10 7 (N/m), the plane strain E50 goes up by 8% and 0%, respectively. In the case of plane stress, E50 goes up by 5% and 0%, where the inter-particle coefficient friction changed from 0.5 to 0.9 and 1.2. Young's modulus 25 (MPa) to 80 (MPa) is typical for medium and dense sand. Therefore, the contact stiffness above 17 × 10 7 (N/m) produces unrealistic Young's modulus values. Figure 6. The sensitivity analysis: (a) inter-particle coefficient of friction vs. peak stress and (b) inter-particle coefficient of friction vs. peak stress. Table 6 summarises E 50 with an inter-particle coefficient of friction for both plane stress and plane stress conditions, indicating that Young's modulus is significantly governed by contact stiffness rather than the inter-particle coefficient of friction. For instance, changing the inter-particle coefficient of friction from 0.5 to 0.9 and 1.2 for a fixed contact stiffness value of 8.45 × 10 7 (N/m), the plane strain E 50 goes up by 8% and 0%, respectively. In the case of plane stress, E 50 goes up by 5% and 0%, where the inter-particle coefficient friction changed from 0.5 to 0.9 and 1.2. Young's modulus 25 (MPa) to 80 (MPa) is typical for medium and dense sand. Therefore, the contact stiffness above 17 × 10 7 (N/m) produces unrealistic Young's modulus values.  Table 7 shows that increasing the inter-particle coefficient of friction leads to an increase in Poisson's ratio for both plane strain and plane stress conditions. However, the calculated Poisson's ratio for plane stress is slightly greater than that calculated for plane strain. For instance, for inter-particle coefficient friction of 1.2 where the contact stiffness is 150 × 10 7 (N/m), the calculated Poisson's ratio for plane strain and plane stress is 0.27 and 0.37, respectively. An increase in inter-particle coefficient friction leads to growth in both inter-particle sliding capacity and inter-particle forces, meaning the contact deformations and particle displacements rise. This leads to a rise in the lateral deformation of the system. The Poisson's ratio between 0.2 and 0.3 is a typical value for medium and dense sand. Table 7. The sensitivity analysis: inter-particle coefficient friction vs. Poisson's ratio for plane strain and plane stress. The trend of stress-strain behaviour presented in Figure 5a is in good agreement with average anisotropy with axial strain shown in Figure 7a. Figure 7a shows the average fabric anisotropy with axial strain to the various inter-particle coefficient of friction when a fixed value of contact stiffness is applied. The average fabric anisotropy increases until a maximum at the peak stress and then reduces with further strain for all inter-particle coefficients of friction. This clearly shows that the inter-particle coefficients of friction have a significant effect on the results. The maximum average fabric anisotropy, which is approximately 0.37, shows how much the contact arrangement drifts from the isotropic state (i.e., a = 0). This term shows how much the system being loaded can develop anisotropy in contact networks. It is also a variance term that statistically shows how well the contact networks are changing during loading. The more average fabric anisotropy there is, the more shear strength capacity can be attained [6,8]. Figure 7b shows that an increase in the inter-particle coefficient of friction has little effect on the average coordination number to peak deviatoric stress up to the axial strain of 1.5%, which is representative of the peak deviatoric stress (refer to Figure 5a). The coordination number is then slightly increased though the value depends on the inter-particle coefficient of friction.
Geotechnics 2023, 3 429 in contact networks. It is also a variance term that statistically shows how well the contact networks are changing during loading. The more average fabric anisotropy there is, the more shear strength capacity can be attained [6,8]. Figure 7b shows that an increase in the inter-particle coefficient of friction has little effect on the average coordination number to peak deviatoric stress up to the axial strain of 1.5%, which is representative of the peak deviatoric stress (refer to Figure 5a). The coordination number is then slightly increased though the value depends on the inter-particle coefficient of friction.

The Sensitivity of Sand System to Normal Contact Stiffness
A series of 15 biaxial tests were performed to determine the effect of the normal contact stiffness on the micro-macro properties of the particulate sand samples. The input data for these 15 tests are listed in Table 5. The ratio between normal contact stiffness and shear contact stiffness is assumed to be unity, while the inter-particle coefficient of friction varies between 0.5, 0.9, and 1.2. The sensitivity of combined deviatoric stress and volumetric strain with axial strain for a wide range of contact stiffness values is presented in Figure 8. It can be seen that the axial strain corresponding to the peak stress notably becomes smaller by increasing the contact stiffness, showing the particulate sand system behaves denser. However, the particulate systems with the lower contact stiffness demonstrate more softening strain behaviour at post-peak, while the particulate systems with higher contact stiffness show more hardening strain behaviour up the peak stress. Additionally, the peak stress for higher values of contact stiffness becomes sharper. Figure 8 also summarises the sensitivity of peak stress to the various normal contact stiffnesses and inter-particle coefficients of friction. Increasing the contact stiffness has less impact on the peak deviatoric stress. For instance, an increase in the normal contact stiffness from 8.45 × 10 7 (N/m) to 150 x10 7 (N/m), where a fixed value of the inter-particle coefficient of friction 0.9 is applied, results in a slight increase to the peak stress by about 2%. Figure 8 shows that the particulate sand system with the lower value of contact stiffness (e.g., 8.45 × 10 7 (N/m) and 17 × 10 7 (N/m)) demonstrates more contraction behaviour due to the shearing followed by dilation behaviour, which is similar to the normal volumetric bulk behaviour of medium and dense sand. In contrast, the model with a higher value of contact stiffness (e.g., 47 × 10 7 to 150 × 10 7 (N/m)) demonstrates notable dilation behaviour with small initial contraction.

The Sensitivity of Sand System to Normal Contact Stiffness
A series of 15 biaxial tests were performed to determine the effect of the normal contact stiffness on the micro-macro properties of the particulate sand samples. The input data for these 15 tests are listed in Table 5. The ratio between normal contact stiffness and shear contact stiffness is assumed to be unity, while the inter-particle coefficient of friction varies between 0.5, 0.9, and 1.2. The sensitivity of combined deviatoric stress and volumetric strain with axial strain for a wide range of contact stiffness values is presented in Figure 8. It can be seen that the axial strain corresponding to the peak stress notably becomes smaller by increasing the contact stiffness, showing the particulate sand system behaves denser. However, the particulate systems with the lower contact stiffness demonstrate more softening strain behaviour at post-peak, while the particulate systems with higher contact stiffness show more hardening strain behaviour up the peak stress. Additionally, the peak stress for higher values of contact stiffness becomes sharper. Figure 8 also summarises the sensitivity of peak stress to the various normal contact stiffnesses and inter-particle coefficients of friction. Increasing the contact stiffness has less impact on the peak deviatoric stress. For instance, an increase in the normal contact stiffness from 8.45 × 10 7 (N/m) to 150 × 10 7 (N/m), where a fixed value of the inter-particle coefficient of friction 0.9 is applied, results in a slight increase to the peak stress by about 2%. Figure 8 shows that the particulate sand system with the lower value of contact stiffness (e.g., 8.45 × 10 7 (N/m) and 17 × 10 7 (N/m)) demonstrates more contraction behaviour due to the shearing followed by dilation behaviour, which is similar to the normal volumetric bulk behaviour of medium and dense sand. In contrast, the model with a higher value of contact stiffness (e.g., 47 × 10 7 to 150 × 10 7 (N/m)) demonstrates notable dilation behaviour with small initial contraction.  As highlighted in Figure 8, an increase in contact stiffness results in the slope of deviatoric stress with axial strain becoming steeper, meaning the bulk stiffness of particulate sand rises. It is because the increase in contact stiffness value and inter-particle friction leads to an increase in the sliding capacity of the contacts. That is, the risk of losing the contacts per particle reduces, and the number of contacts contributing to taking the load is high. Figure 9 and Table 6 show the plane stress and plane strain secant elastic modulus of the particulate sand systems for a wide range of contact stiffness and inter-particle coefficient of friction. A linear relationship can be established for each inter-particle coefficient of friction. In the case of plane strain, E = 3 × 10 −7 k n + 13.406 for inter-particle coefficient of friction 0.5 E = 3 × 10 −7 k n + 12.089 for inter-particle coefficient of friction 0.9 (21) E = 4 × 10 −7 k n + 9.2826 for inter-particle coefficient of friction 1.2 leads to an increase in the sliding capacity of the contacts. That is, the risk of losing the contacts per particle reduces, and the number of contacts contributing to taking the load is high. Figure 9 and Table 6 show the plane stress and plane strain secant elastic modulus of the particulate sand systems for a wide range of contact stiffness and inter-particle coefficient of friction. A linear relationship can be established for each inter-particle coefficient of friction. In the case of plane strain, = 3 × 10 + 13.406 for inter-particle coefficient of friction 0.5 = 3 × 10 + 12.089 for inter-particle coefficient of friction 0.9 (21) = 4 × 10 + 9.2826 for inter-particle coefficient of friction 1.2 In the case of plane stress, another set of linear relationships can be established for each inter-particle coefficient of friction: = 3 × 10 + 12.776 for inter-particle coefficient of friction 0.5 = 4 × 10 + 11.28 for inter-particle coefficient of friction 0.9 (22) = 4 × 10 + 8.0763 for inter-particle coefficient of friction 1.2 The values of normal contact stiffness between 8.45 × 10 7 (N/m) and 17 × 10 7 (N/m) lead to values of E50, which fall within the range typical for medium and dense sand, i.e., between 25 (MPa) and 50 (MPa) and 50 (MPa) and 80 (MPa), respectively.
(a) (b) Figure 9. The sensitivity of E50 to the various particle normal contact stiffnesses: (a) plane strain and (b) plane stress. Figure 9. The sensitivity of E 50 to the various particle normal contact stiffnesses: (a) plane strain and (b) plane stress.
In the case of plane stress, another set of linear relationships can be established for each inter-particle coefficient of friction: E = 3 × 10 −7 k n + 12.776 for inter-particle coefficient of friction 0.5 E = 4 × 10 −7 k n + 11.28 for inter-particle coefficient of friction 0.9 (22) E = 4 × 10 −7 k n + 8.0763 for inter-particle coefficient of friction 1.2 The values of normal contact stiffness between 8.45 × 10 7 (N/m) and 17 × 10 7 (N/m) lead to values of E 50 , which fall within the range typical for medium and dense sand, i.e., between 25 (MPa) and 50 (MPa) and 50 (MPa) and 80 (MPa), respectively. Figure 10 and Table 7 show that an increase in the normal contact stiffness leads to an increase in Poisson's ratio with a ranged value between 0.2 and 0.4. The contact stiffness of particles between 8.45 × 10 7 (N/m) and 17 × 10 7 (N/m) produces typically ranged values of Poisson's ratio for medium sand if the inter-particle coefficient of friction is between 0.9 and 0.5. For the higher value of the inter-particle coefficient of friction (i.e., 1.2), the interpreted value of Poisson's ratio is less than typical values. A linear relationship can also be established between Poisson's ratio and the inter-particle coefficient of friction (in the case of plane strain): v = 3 × 10 −11 k n + 0.2098 for inter-particle coefficient of friction 0.5 v = 3 × 10 −11 k n + 0.1870 for inter-particle coefficient of friction 0.9 (23) v = 2 × 10 −11 k n + 0.1809 for inter-particle coefficient of friction 1.2 and 0.5. For the higher value of the inter-particle coefficient of friction (i.e., 1.2), the interpreted value of Poisson's ratio is less than typical values. A linear relationship can also be established between Poisson's ratio and the inter-particle coefficient of friction (in the case of plane strain): = 3 × 10 + 0.2098 for inter-particle coefficient of friction 0.5 = 3 × 10 + 0.1870 for inter-particle coefficient of friction 0.9 (23) = 2 × 10 + 0.1809 for inter-particle coefficient of friction 1.2 A linear relationship can also be established between Poisson's ratio and the interparticle coefficient of friction (in the case of plane stress): = 5 × 10 + 0.2650 for inter-particle coefficient of friction 0.5 = 4 × 10 + 0.2299 for inter-particle coefficient of friction 0.9 (24) = 4 × 10 + 0.2205 for inter-particle coefficient of friction 1.2 (a) (b) Figure 10. The sensitivity of Poisson's ratio to the various particle normal stiffness: (a) plane strain and (b) plane stress. Figure 11 demonstrates the evolution of contact distribution during the shearing of particulate and where a fixed value is applied to the inter-particle properties. At this isotropic state, the average fabric anisotropy is 0.0007, and a circle presents the polar diagram distribution of contacts. At this stage, the number of contacts is distributed almost equally within each segment. If the polar diagram is entirely circled, it shows that the distribution of contacts is in an isotropic state. At the peak bulk deviatoric stress, the orientation of the contacts will be toward the direction of major stress, indicating the contacts within the particulate sand are fully mobilised to take the load. The average fabric anisotropy rises to 0.35. At the post-peak state where = 10%, the orientation of contact points will be mainly toward the confining stress, and average fabric anisotropy will be reduced to 0.27, showing the particulate system has collapsed. A linear relationship can also be established between Poisson's ratio and the interparticle coefficient of friction (in the case of plane stress): v = 5 × 10 −11 k n + 0.2650 for inter-particle coefficient of friction 0.5 v = 4 × 10 −11 k n + 0.2299 for inter-particle coefficient of friction 0.9 (24) v = 4 × 10 −11 k n + 0.2205 for inter-particle coefficient of friction 1.2 Figure 11 demonstrates the evolution of contact distribution during the shearing of particulate and where a fixed value is applied to the inter-particle properties. At this isotropic state, the average fabric anisotropy is 0.0007, and a circle presents the polar diagram distribution of contacts. At this stage, the number of contacts is distributed almost equally within each segment. If the polar diagram is entirely circled, it shows that the distribution of contacts is in an isotropic state. At the peak bulk deviatoric stress, the orientation of the contacts will be toward the direction of major stress, indicating the contacts within the particulate sand are fully mobilised to take the load. The average fabric anisotropy rises to 0.35. At the post-peak state where ε 11 = 10%, the orientation of contact points will be mainly toward the confining stress, and average fabric anisotropy will be reduced to 0.27, showing the particulate system has collapsed. An increase in contact stiffness from 8.45 × 10 7 (N/m) to 150 × 10 7 (N/m) leads to an increase in average fabric anisotropy (see Figure 12a). Raising the ability of particulate systems to develop the fabric anisotropies results in an increase in their shear capacity.
Comparing Figures 8 and 12a shows that both have a similar trend, and the maximum fabric anisotropy takes place around peak deviatoric stress. An increase in contact stiffness from 8.45 × 10 7 (N/m) to 150 × 10 7 (N/m) leads to an increase in average fabric anisotropy (see Figure 12a). Raising the ability of particulate systems to develop the fabric anisotropies results in an increase in their shear capacity.
Comparing Figures 8 and 12a shows that both have a similar trend, and the maximum fabric anisotropy takes place around peak deviatoric stress. Figure 12b indicates the sensitivity of the average coordination number to the normal particle stiffness when the inter-particle coefficient of friction is 0.9. The average coordination number after peak deviatoric stress (see Figure 8) for the models using k n = 17 × 10 7 (N/m) to k n = 150 × 10 7 (N/m) falls below three contacts, which is less than the value required for static equilibrium. It is because an increase in contact stiffness or inter-particle friction leads to an increase in the inter-particle contact forces and sliding capacity between particles. The developed chain contact forces at each contact during shearing rise to maintain the particulate sand in a static equilibrium and carry the load. Therefore, a lower average coordination number is expected to develop for higher contact stiffness to resist shearing. The average coordination number increases for the lower values of contact stiffness due to the dilatant behaviour of the particulate sand. For instance, for a model using k n = k s = 8.45 × 10 7 (N/m), the system contracts till axial strain of around 0.5% (see Figure 8), demonstrating that the particulate sand becomes more compacted, and the tendency of particles to lose their contacts decreases. Therefore, the average coordination number is still high up to that point. By commencing the dilation behaviour, the tendency of particles to lose their contacts rises. That means the average coordination number drops significantly once dilation behaviour begins. By Reducing the dilation behaviour, the ability of the particulate sand to develop higher average fabric anisotropy in comparison to that for higher contact stiffness decreases.  Figure 12b indicates the sensitivity of the average coordination number to the normal particle stiffness when the inter-particle coefficient of friction is 0.9. The average coordination number after peak deviatoric stress (see Figure 8) for the models using kn = 17 × 10 7 (N/m) to kn = 150 × 10 7 (N/m) falls below three contacts, which is less than the value required for static equilibrium. It is because an increase in contact stiffness or inter-particle friction leads to an increase in the inter-particle contact forces and sliding capacity between particles. The developed chain contact forces at each contact during shearing rise to maintain the particulate sand in a static equilibrium and carry the load. Therefore, a lower average coordination number is expected to develop for higher contact stiffness to resist shearing. The average coordination number increases for the lower values of contact stiffness due to the dilatant behaviour of the particulate sand. For instance, for a model using kn = ks = 8.45 × 10 7 (N/m), the system contracts till axial strain of around 0.5% (see Figure 8), demonstrating that the particulate sand becomes more compacted, and the tendency of particles to lose their contacts decreases. Therefore, the average coordination number is still high up to that point. By commencing the dilation behaviour, the tendency of particles to lose their contacts rises. That means the average coordination number drops significantly once dilation behaviour begins. By Reducing the dilation behaviour, the ability of the particulate sand to develop higher average fabric anisotropy in comparison to that for higher contact stiffness decreases.

The Sensitivity to the Shear Particle Stiffness
To study the sensitivity of the micro-macro mechanical behaviour of particulate sand to the shear contact stiffness, an additional fifteen biaxial tests were carried out. The ratio between normal contact stiffness and shear contact stiffness for these additional tests is

The Sensitivity to the Shear Particle Stiffness
To study the sensitivity of the micro-macro mechanical behaviour of particulate sand to the shear contact stiffness, an additional fifteen biaxial tests were carried out. The ratio between normal contact stiffness and shear contact stiffness for these additional tests is assumed to be 0.5, while the inter-particle coefficient of friction varies between 0.5, 0.9, and 1.2. The input data for these additional fifteen biaxial tests are listed in Table 8. Table 8. The inter-particle parameters for additional fifteen biaxial tests to investigate the sensitivity of system to the shear contact stiffness.  Figure 13 compares the sensitivity of combined deviatoric stress and volumetric strain with axial strain to the shear contact stiffness, for a wide range of contact stiffness and a fixed value of the inter-particle coefficient of friction of 0.9. The ratio of shear contact stiffness to the normal contact stiffness varies between 1 and 0.5. Changing the shear contact stiffness to half of the normal contact stiffness has a minor impact on the peak stress. The peak stress is about 260 (kPa), though the model with a contact stiffness of 150 × 10 7 (N/m) shows a slightly higher peak stress. For cohesionless particulate systems that are under the shearing, the number of contacts contributing to taking the load plays a major role in the magnitude of the peak stress that these systems can achieve. The sliding capacity of contact is a function of normal contact force and inter-particle coefficient of friction. The normal contact force is a function of contact stiffness and contact deformation. The latter can be dictated by the strain rate of loading applied to the sample. The higher the strain rate, the higher the peak stress that can possibly be gained, as presented by Yamamuro et al. [49]. As the input data and initial conditions of these tests are similar (e.g., they are isotopically consolidated to 100 (kPa) and subjected to the same strain rate), the models with higher combined contact stiffness are stiffer and have a lower number of losing contacts during loading. That is, the stiffer models demonstrate higher hardening strain behaviour compared to those using k s k n = 0.5. It also can be observed from Figure 12 that shear contact stiffness significantly controls the pre-peak. A change in shear contact stiffness from k s k n = 1 to k s k n = 0.5 plays an important role in the slope of the backbone stress-strain curve and hardening strain behaviour. Chi et al. [50], in their sensitivity study, showed that the slope of stress-strain behaviour of sand is influenced by shear contact stiffness. This alteration is more evident for the samples with contact stiffness 8.45 × 10 7 and 17 × 10 7 (N/m), while for models with higher contact stiffness, this change has less impact on the slope of the backbone stress-strain curve. This is because the models with k s k n = 0.5 are less stiff, and the number of lost contacts can be larger during shearing for them in compression with the k s k n = 1 models. Therefore, the models with k s k n = 0.5 experience more axial displacement than k s k n = 1 models under a similar stress level. As seen, the axial strain corresponding to the peak stress rises for the models with k s k n = 0.5. Changing the shear contact stiffness to half of the normal contact stiffness also plays a major role in the post-peak and softening strain behaviour of particulate sand, such that those models using k s k n = 0.5 show more softening strain behaviour post-peak. Reducing the shear contact stiffness to half of the normal contact stiffness value results in significant changes in the contraction and dilation volumetric behaviour of sand, in particular for those two models using lower normal contact stiffness, 8.45 × 10 7 (N/m) and 17 × 10 7 (N/m). Table 9 summarises the values of Young's modulus for the wide range of inter-particle coefficients of friction and contact stiffnesses where the shear contact stiffness is half of the normal contact stiffness. As explained above and from comparing Tables 6 and 9, a change in the shear contact stiffness from k s k n = 1 to k s k n = 0.5 results in an alteration in the value of Young's modulus for both plane strain and plane stress. For instance, the value of Young's modulus at plane strain conditions for the model with k s k n = 1 and k s k n = 0.5, where the normal contact stiffness is 8.45 × 10 7 (N/m), and the inter-particle coefficient of friction is 1.2, calculated as 35 (MPa) and 29 (MPa), respectively. This is about a 17% reduction in Young's modulus. This reduction in Young's modulus at plane strain conditions for contact stiffness 17 × 10 7 (N/m), 46 × 10 7 (N/m), 133 × 10 7 (N/m), and 150 × 10 7 (N/m) when shear contact stiffness changed from k s k n = 1 to k s k n = 0.5 is 12%, 21%, 13%, and 17%, respectively. A similar exercise was carried out to investigate the sensitivity of Young's modulus shear contact stiffness where a fixed inter-particle coefficient friction of 0.9 is applied. The plane strain Young's modulus value drops by 20%, 12%, 18%, 12%, and 17% for contact stiffness 8.45 × 10 7 (N/m), 17 × 10 7 (N/m), 46 x10 7 (N/m), 133 × 10 7 (N/m), and 150 × 10 7 (N/m). Changing the shear contact stiffness from k s k n = 1 to k s k n = 0.5 leads to an average reduction of 16% in plane strain Young's modulus for the particulate sand.    Tables 7 and 9 for k s k n = 1 and k s k n = 0.5 models were compared. Comparison between Tables 7 and 10 shows that a reduction in shear contact stiffness to half of the normal contact stiffness leads to a notable increase in the value of Poisson's ratio for both plane strain and plane stress. Coetzee and Els [51] and Belheine et al. [31] showed in their parametrical study that a change in the k s k n ratio has an influence on the materials' Poisson's ratio. This reduction in Poisson's ratio at plane strain conditions and fixed inter-particle coefficient friction of 0.9 for contact stiffness 8.45 × 10 7 (N/m), 17 × 10 7 (N/m), 46 × 10 7 (N/m), 133 × 10 7 (N/m), and 150 × 10 7 (N/m) when shear contact stiffness changed from k s k n = 1 to k s k n = 0.5 is 28%, 26%, 25%, 17%, and 12%, respectively. The impact of this reduction in shear contact stiffness value to the plane strain Poisson's ratio was studied where a fixed inter-particle coefficient friction of 1.2 was used. The plane strain Poisson's ratio value increases by 14%, 14%, 8%, 8%, and 7% for contact stiffness 8.45 × 10 7 (N/m), 17 × 10 7 (N/m), 46 × 10 7 (N/m), 133 × 10 7 (N/m), and 150 × 10 7 (N/m). Comparing the numbers shows that a reduction in shear contact stiffness to half of the normal contact stiffness has a significant impact on the models with lower contact stiffness (e.g., 8.75 × 10 7 (N/m) and 17 × 10 7 (N/m)).  As highlighted above, an increase in k s k n from 0.5 to 1 has a small impact on the peak deviatoric stress. Using Equation (19), the value of the angle of internal friction can be calculated from peak stress and the applied confining pressure. The peak stress for the models using k s k n = 0.5 and k s k n = 1 is taken from Figures 8 and 13, respectively. Table 11 compares the angle of internal friction for k s k n = 0.5 and k s k n = 1. As seen, an increase in k s k n from 0.5 to 1 plays a minor role in the angle of internal friction. For instance, the angle of friction for the model with contact stiffness 17 × 10 7 (N/m) and the inter-particle coefficient of friction of 0.9, where shear contact stiffness is reduced to half of the normal contact stiffness, changed from 27 • to 26 • . Table 11. The sensitivity of angle of friction to the various particle shear stiffness. The polar diagram of contact distribution during the shearing for the model with a contact stiffness of 8.45 × 10 7 (N/m), where k s k n = 0.5 and a fixed value of 0.9 is applied to the inter-particle coefficient friction, is presented in Figure 14. Comparing Figures 11 and 14 show that at peak stress, the diameter of the polar diagram in the major direction, representing the number of contacts per segment, reduced from around 3400 to about 3200. This reduction in the number of contacts results in the particulate system becoming less stiff and taking less load. Additionally, the average fabric anisotropy drops from 0.35 to 0.33, indicating the ability of the system with k s k n = 0.5 to generate the number of contacts to take more load reduces. The orientation of contacts at the post-peak is also toward the confining stress, with average fabric anisotropy reduced to 0. 26. This reduction in the number of contacts is observed in Figure 15, where the sensitivity of the average coordination number to a reduction in shear contact stiffness from k s k n = 1 and k s k n = 0.5 for two sets of contact stiffness (i.e., 8.75 × 10 7 (N/m) and 17 × 10 7 (N/m)) is compared. At ε 11 = 0.0, the average coordination number for models with the lower shear contact stiffness value is approximately similar to that obtained from the model with higher shear contact stiffness. From ε 11 = 0.0 to the axial strain corresponds to the peak stress (see Figure 14), and the average coordination number for all models reduces. However, the models with lower shear contact stiffness experience a higher reduction in average coordination number. At post-peak, the average coordination number significantly decreases for the model with lower shear contact stiffness. resenting the number of contacts per segment, reduced from around 3400 to about 3200. This reduction in the number of contacts results in the particulate system becoming less stiff and taking less load. Additionally, the average fabric anisotropy drops from 0.35 to 0.33, indicating the ability of the system with = 0.5 to generate the number of contacts to take more load reduces. The orientation of contacts at the post-peak is also toward the confining stress, with average fabric anisotropy reduced to 0.26.  stress (see Figure 14), and the average coordination number for all m ever, the models with lower shear contact stiffness experience a hig age coordination number. At post-peak, the average coordination decreases for the model with lower shear contact stiffness. To verify the cohesionless nature of the synthetic material used biaxial tests with different confining pressure values were implemen the inter-particle properties along with the size of the sample, initia ber of particles, and confining pressures used for this study. These s cally consolidated under different confining pressures before being under a strain-control approach, as explained in Section 4.4.  To verify the cohesionless nature of the synthetic material used for this study, three biaxial tests with different confining pressure values were implemented. Table 12 presents the inter-particle properties along with the size of the sample, initial porosity, PSD, number of particles, and confining pressures used for this study. These samples were isotopically consolidated under different confining pressures before being subjected to shearing under a strain-control approach, as explained in Section 4.4.  Figure 16a shows the variation of deviatoric stress with axial strain. As expected, an increase in confining pressure leads to an increase in the peak stress. The ratio of peak stress for confining pressure of 300 (kPa) to that of 100 (kPa) is about 3, the ratio of peak stress for confining pressure 300 (kPa) to that of 200 (kPa) is almost 1.5, and the ratio of peak stress for confining pressure of 200 (kPa) to that or 100 (kPa) is approximately 1.0. This shows that these ratios are in good agreement with the ratio 300 (kPa)/100 (kPa), 300 (kPa)/200 (kPa), and 200 (kPa)/100 (kPa), respectively. The Mohr circle of each test, together with the Mohr-Coulomb failure envelope of these three tests, are presented in Figure 16b. The angle of internal friction extracted from the Mohr-Coulomb failure envelope (i.e., 27 • ) is similar to the angle of internal friction reported in Table 10 for the model using the contact stiffness of 8.45 × 10 7 (N/m) and the inter-particle coefficient of friction of 0.9. It can be seen that the nature of the synthetic material used for this study is cohesionless. stress for confining pressure of 300 (kPa) to that of 100 (kPa) is about 3, the ratio of peak stress for confining pressure 300 (kPa) to that of 200 (kPa) is almost 1.5, and the ratio of peak stress for confining pressure of 200 (kPa) to that or 100 (kPa) is approximately 1.0. This shows that these ratios are in good agreement with the ratio 300 (kPa)/100 (kPa), 300 (kPa)/200 (kPa), and 200 (kPa)/100 (kPa), respectively. The Mohr circle of each test, together with the Mohr-Coulomb failure envelope of these three tests, are presented in Figure 16b. The angle of internal friction extracted from the Mohr-Coulomb failure envelope (i.e., 27°) is similar to the angle of internal friction reported in Table 10 for the model using the contact stiffness of 8.45 × 10 7 (N/m) and the inter-particle coefficient of friction of 0.9. It can be seen that the nature of the synthetic material used for this study is cohesionless.

Conclusions
In this paper, thirty-three biaxial DEM tests with rigid walls were implemented to study the critical inter-particle properties and the range over which these properties impact the micro-macro mechanical behaviour. The modelling results show that the elastic parameter, Young's modulus, with particle diameters varying between 0.5 (mm) and 3.0 (mm), is greatly dependent on the normal contact stiffness, while the angle of friction is strongly dependent on the inter-particle coefficient of friction. It was observed that the inter-particle coefficient of friction between 0.9 and 1.2 produced an angle of internal friction between 27° and 32° such that the relationship between them seems to be linear. These values are compatible with typical ranges of the angle of friction of medium sand. That is, to produce the angle of friction for dense sand, a higher inter-particle coefficient friction (i.e., more than 1.2) should be applied. The angle of internal friction obtained from interparticle coefficient of friction 0.5 is between 18° and 20°, which is not within the angle of friction of sand. The values of normal contact stiffness between 8.45 × 10 7 and 17 × 10 7 with ratios = 1 and = 0.5 results in producing a range of Young's modulus values of medium sand. However, Young's modulus produced based on = 0.5 are about 15% less than that estimated from the models applied = 1. The higher normal contact stiffness (e.g., 46 × 10 7 (N/m) and 150 × 10 7 (N/m)) produces an unrealistic Young's modulus for sand. It was observed that Young's modulus computed from plane stress conditions is about 5% bigger than that calculated from plane strain conditions. Reduction in the value of shear contact stiffness to half of the value of normal contact stiffness led to a decrease in the total number of contacts in load bearing as the particulate system sheared. This indicates the particulate system became less stiff. This change in shear contact stiffness

Conclusions
In this paper, thirty-three biaxial DEM tests with rigid walls were implemented to study the critical inter-particle properties and the range over which these properties impact the micro-macro mechanical behaviour. The modelling results show that the elastic parameter, Young's modulus, with particle diameters varying between 0.5 (mm) and 3.0 (mm), is greatly dependent on the normal contact stiffness, while the angle of friction is strongly dependent on the inter-particle coefficient of friction. It was observed that the inter-particle coefficient of friction between 0.9 and 1.2 produced an angle of internal friction between 27 • and 32 • such that the relationship between them seems to be linear. These values are compatible with typical ranges of the angle of friction of medium sand. That is, to produce the angle of friction for dense sand, a higher inter-particle coefficient friction (i.e., more than 1.2) should be applied. The angle of internal friction obtained from inter-particle coefficient of friction 0.5 is between 18 • and 20 • , which is not within the angle of friction of sand. The values of normal contact stiffness between 8.45 × 10 7 and 17 × 10 7 with ratios k s k n = 1 and k s k n = 0.5 results in producing a range of Young's modulus values of medium sand. However, Young's modulus produced based on k s k n = 0.5 are about 15% less than that estimated from the models applied k s k n = 1. The higher normal contact stiffness (e.g., 46 × 10 7 (N/m) and 150 × 10 7 (N/m)) produces an unrealistic Young's modulus for sand. It was observed that Young's modulus computed from plane stress conditions is about 5% bigger than that calculated from plane strain conditions. Reduction in the value of shear contact stiffness to half of the value of normal contact stiffness led to a decrease in the total number of contacts in load bearing as the particulate system sheared. This indicates the particulate system became less stiff. This change in shear contact stiffness plays a major role in the volumetric behaviour of these systems, including contraction and dilation, meaning Poisson's ratio is mainly governed by shear contact stiffness. The results show reducing shear contact stiffness resulting in an increase in Poisson's ratio values for both plane strain and plane stress conditions. The values of Poisson's ratio for plane strain conditions obtained for those particulate systems applying the normal contact stiffness between 8.45 × 10 7 and 17 × 10 7 and k s k n = 1 and k s k n = 0.5 are within the typical values of Poisson's ratio of medium sand. That is, to produce Poisson's ratio values within the typical range of dense sand, the ratio k s k n should be less than 0.5. The macro-mechanical behaviour of particulate sand is dictated by inter-particle properties such that a change in one of these parameters results in a significant change in the macro-mechanical behaviour, such as stress-strain and volumetric behaviour. A set of relationships was established between inter-particle properties (including normal contact stiffness and inter-particle coefficient of friction) and macro-machinal parameters such as Young's modulus, Poisson's ratio, and angle of internal friction. Reviewing the polar diagram of contact distribution shows that the peak stress and hardening strain behaviour for a particulate system are a function of the number of contacts developed in the major stress direction. The ability of a system to develop the number of contacts is tied to its fabric anisotropy. It was seen that the models use the higher inter-particle coefficient of friction value to develop higher fabric anisotropy and subsequently produce higher peaks and angles of internal friction. This paper only covers three components of ideal disk inter-particle properties, e.g., normal contact stiffness, shear contact stiffness, and particle friction. The impact of other particle-scale parameters, such as particle shape, rolling resistance, and lateral boundary conditions on the macro-mechanical properties, are not covered in this study. This study presents a set of relationships between a wide range of inter-particle properties and elastoplastic parameters for medium and dense sands. These relationships can be used to estimate the inter-particle properties values from lab test data of sands, including elastoplastic parameters for DEM modelling of geotechnical problems. For instance, the first author of this paper applied these relationships to estimate the input data of inter-particle properties of seabed sand from lab test data to simulate the lateral behaviour of the Gravity Base Foundation of an offshore windfarm turbine using PFC 2D for an industrial project. Nomenclature . x Translational particle velocity .