Investigation on the Influencing Factors of K₀ of Granular Materials Using Discrete Element Modelling

Abstract

Earth pressure coefficient at rest K₀ is commonly estimated by empirical equations, which to date has had insufficient accuracy and universality. For better prediction, the investigation on the factors influencing K₀ is required. A series of discrete element method (DEM) simulations of oedometer tests are conducted to verify the key factors influencing K₀ of granular materials. The influences of initial fabric anisotropy, particle shape, initial void ratio, inter-particle friction angle is investigated. The evolution of microstructure is monitored during the tests to reveal the relationship between the microstructure evolution and K₀ values. The results show that the effect of fabric anisotropy exists but is limited. Particle shape, initial void ratio, and inter-particle friction angle all significantly affect the K₀ values alone. According to the DEM results, an attempt is made to propose a more reasonable empirical equation in which K₀ is a function of relative density, critical state friction angle, and “shape factor”. This new empirical equation has higher accuracy and can consider the effect of particle shape, inspiring the determination of K₀ values in practical engineering.

1. Introduction

The determination of the initial stress state in soil is essential for the design and analysis of geotechnical structures. Stress state consists of both vertical and lateral stress, and the former one can be easily obtained by the product of soil unit weight and soil depth while the latter one is often much harder to calculate accurately. In practice, lateral stress is usually estimated by the product of the coefficient of earth pressure at rest K₀ [1] and the vertical stress for simplicity. As for K₀, there are two types of definitions in the literature. One is the ratio of the effective horizontal stress ${\sigma }_h^{\prime }$ and the effective vertical stress ${\sigma }_v^{\prime }$ without lateral deformation [2], and the other one is expressed in an incremental form [3]:

$$K_0=\frac{{\sigma }_h^{\prime }}{{\sigma }_v^{\prime }}$$

(1)

$$K_0^{\prime }=\frac{\Delta {\sigma }_h^{\prime }}{\Delta {\sigma }_v^{\prime }}$$

(2)

where $\Delta {\sigma }_h^{\prime }$ and $\Delta {\sigma }_v^{\prime }$ are the instantaneous increments of the effective horizontal stress and effective vertical stress with zero lateral strain, respectively.

Due to the difficulty in measuring K₀ values directly in the field, numerous approaches have been proposed to calculate K₀ from both theoretical derivation and experimental data. For instance, Jáky proposed an equation which is most widely used:

$$K_0=\frac{{\sigma }_h}{{\sigma }_v}=\left(1+\frac{2}{3}sin{\phi }^{\prime }\right)\frac{1-sin{\phi }^{\prime }}{1+sin{\phi }^{\prime }}\approx 0.95-sin{\phi }^{\prime }$$

(3)

or approximately:

$$K_0=1-sin{\phi }^{\prime }$$

(4)

for loose deposits and normally consolidated clays, where ${\phi }^{\prime }$ indicate the effective internal friction angle of the soil. This equation has been widely adopted in engineering practice for many years. Contrary to a common misconception, as Michalowski [4] had pointed out, this equation is a result of theory deduction instead of an empirical formula. It was deduced from the stress field in a wedge prism of loose granular materials. Jáky [3] pointed out that the stress on the central vertical plane is the horizontal pressure at rest. Equation (4) is a widely accepted representation of K₀ because it is easy to use and provides good enough prediction. Both experiments [5] and analysis [6] have verified the equation’s reliability. However, after years of application and research of this formula, there is considerable controversy about ${\phi }^{\prime }$. In practical engineering or scientific applications, internal effective friction angles are difficult to measure directly. Therefore, peak friction angle ${\phi }_p^{\prime }$ [7,8,9,10] and critical state friction angle ${\phi }_c^{\prime }$ [6,7,11,12,13,14] are often used as ${\phi }^{\prime }$ to calculate K₀ values. Furthermore, it is worth noting that Jáky’s originally work aimed at loose deposit soil, for which peak friction angle ${\phi }_p^{\prime }$ equals the critical state friction angle ${\phi }_c^{\prime }$. When the equation is applied to other soil, it is hard to decide which friction angle is better.

Although Jáky’s equation is widely accepted, his derivation has been criticized by many researchers. Handy [15] affirmed that the soil state in the wedge prism was different from that in the commonly referred K₀ state. Michalowski [4] also discovered that the stress path in a wedge-shaped sand prism could not represent the one under typical K₀ state. In the opinion of Federico [16], failure stress conditions in the middle of a wedge prism appeared quite different from soil stress conditions under K₀ state, and he thought that a new formula using the angle of internal friction mobilized ${\phi }_{mob}^{\prime }$ would be more reasonable to calculate K₀ values.

Based on the results of a series of experiments [5], researchers realized that to achieve high precision, K₀ values certainly cannot be represented by a function that contains only internal friction angle as suggested by Jáky. Enormous studies have been carried out to explore the factors affecting K₀-value for predicting it better. It is found that K₀-value is affected by void ratio, confining pressure, stress history, sample preparation method, and the mineral composition of particles, etc.; many factors, among which the first three factors are recognized as the most important. However, there has been controversy over how these factors specifically affect K₀. For instance, there is no consensus on the effect of void ratio on K₀-value. Guo [6] and Feda [8] obtained similar results, that is, that the K₀ value of granular soil increases with the increase of void ratio. Contrary to this conclusion, it was found by Chu and Gan [9] and Wanatowski and Chu [7] that the K₀ value of the sample obtained by the water sedimentation method can be considered constant with the change of void ratio. It is worth mentioning that the correlation between void ratio and K₀ values has not been well revealed so far. Wu et al. [17] reported that K₀ values are affected by the applied stress, but they argued that the small variation of K₀ values is caused by the change of porosity ratio with stress level. Bauer’s [18] laboratory test results support Wu’s conclusion. Over-consolidation ratio (OCR) is often used to quantify the effect of stress history on K₀ values. Mayne and Kulhawy [5] collected and analyzed the data of K₀ values of more than 170 soil samples in the previous literature and put forward a formula to calculate K₀ values, which can take into account the influence of OCR.

Inherently anisotropic granular soils are ubiquitous in geotechnical engineering practice. Previous studies have confirmed that the fabric of granular materials greatly affects their mechanical properties, such as strength, deformation, and failure [19,20,21], while research is limited as for its effect on one-dimensional compression of granular material [7,11,22]. Northcutt and Wijewickreme [11] studied the different responses of Fraser River sand samples under one-dimensional compression by three kinds of sample preparation methods: air pluviation, tamping, and vibration. They found that specimens prepared by different reconstitution methods exhibit different responses and produce different K₀ values, and thus proposed an empirical K₀ equation involving fabric factor. Guo and Stolle [23] experimentally studied the effect of initial fabric on K₀ values of granular materials. According to the experimental results, they pointed out that the fabric of granular materials does have a significant effect on the K₀ values of granular materials. Limited by laboratory test conditions, there is no good method in existing research to present the fabric of granular material and its evolution during loading.

In laboratory, K₀ values are usually measured by triaxial tests and oedometer tests. The radial strain variation in triaxial tests [22,24,25], as well as the loading system compliance and the side friction in oedometer tests [26], bring a great difficulty to accurately measure K₀-value. This inaccuracy creates trouble for determining the influencing factors of K₀ and their influencing mode.

In recent studies about granular materials, discrete element method (DEM) [27] appears increasingly [28,29,30], which provides a new perspective from micro-mechanism to study issues in granular materials. Additionally, the measuring error in laboratory tests can be avoided in DEM and the evolution of fabric can be easily monitored. Therefore, DEM can be a powerful tool for researchers to study the influencing factors of K₀ precisely. For instance, Gu [10] studied the effects of confining pressure, soil density, and over-consolidation ratio on K₀ values using DEM method. At the same time, the development of microstructure during K₀ loading is monitored by using the characteristics of DEM method, to explore the microscopic factors affecting K₀ values. Khalili [31] generated samples with different void ratios, different initial fabric, and different ways of packing process in DEM simulation. The development of a series of macro and micro parameters, such as fabric, void ratio, coordination number, anisotropy parameter of contact force, and K₀ values in oedometric compression process are recorded in detail. Concerning the K₀ values, Khalili established a simple formula for K₀ with the parameters of anisotropy parameters of contact force and fabric. However, the influence of initial fabric anisotropy on K₀ values is missing in these DEM modellings.

In this paper, we first investigate the effect of initial fabric anisotropy on the K₀ values of granular materials. The inherent anisotropic specimens are prepared in DEM simulation and then an oedometer test is performed to determine K₀ values. K₀ values and evolution of microstructure of specimens with different initial fabric are analyzed and compared. Then, a series of DEM simulations are carried out to analyse the effects of particle shape, initial void ratio, and angle of internal friction between particles on the value of K₀ values of granular materials. The results of DEM tests show that the effect of initial fabric anisotropy on the K₀ of granular materials is minor. Particle shape, initial void ratio, and internal friction angle all significantly affect K₀-value. Based on these results, the rationality of the existing empirical formulas of K₀ is discussed, and a new formula of relatively high accuracy is proposed in this paper.

2. DEM Modelling

The DEM simulation in this paper is carried out using the well-recognized Particle Flow Code(PFC^3D) [32]. Approximately 27,000 particles in a cubic space of 110 mm³ constitute the specimen (Figure 1), which is confined by three pairs of frictionless walls. The particle size distribution of the specimen mimics Toyoura sand (Figure 2) with a mean particle diameter of 0.214 mm and the coefficient of uniformity of 1.425. The nonlinear Hertz-Mindlin contact law [33] is used and the shear modulus and Poisson’s ratio is 5.8 GPa and 0.15, respectively. The friction between the walls and particles is set to be zero. The parameters of the simulation in this paper are listed in

Table 1. Parameters in the DEM simulation.

Specimen Dimensions (mm)	4.8 × 4.8 × 4.8
Particle density (kg/m3)	2650
No. of particles	27,000
Mean particle diameter (mm)	0.214
Inter-Particle friction coefficient	0.5
Wall-particle friction coefficient	0
Wall stiffness (N/m)	1 × 1010
Contact law	Hertz-Mindlin
Shear modulus of particles (GPa)	5.8
Poison’s ratio of particle	0.15

.

The clump of aspect ratio 1.7 (Figure 3) is adopted to produce inherent anisotropy. Initially, the anisotropic specimen is prepared with reference to the method adopted by Zhao and Guo [34]. The particles are generated with a coincident orientation and undergo two-stage isotropic consolidation with 10 kPa consolidation pressure. The orientations of the particles are fixed in the first stage and are then set to be free in the second stage. As illustrated in Figure 4, the angle δ between the long axis orientation of clumps and horizontal plane represents the inherent anisotropy principal direction. Five specimens of different anisotropy principal direction (δ = 90°, 60°, 45°, 30°, and 0°) but the same anisotropy degree is prepared. The contact normal distributions of these five specimens are plotted by 3D rose graphs in Figure 5. For specimens with angle δ equaling 90°, 60°, 45°, 30°, and 0°, the scalar anisotropy parameter a_r are 1.291, 1.336, 1.354, 1.361, and 1.311, respectively, representing the nearly identical anisotropy degree. From the rose graph, specimens with different initial anisotropic principal directions can be produced effectively by this specimen preparation method. As a reference object, an isotropic specimen is made, with random orientation in the generation stage and no restraint of orientation at any stage. In the isotropic compression stage, the inter-particle friction angle is 0.15 while in the K₀ loading stage the angle is 0.5. After consolidation, the void ratios are 0.623, 0.622, 0.626, 0.626, 0.627, and 0.624 respectively, for the isotropic specimen and the other five ones of δ = 90°, 60°, 45°, 30°, and 0°. The effect of initial void ratio has been eliminated as much as possible since their values are roughly equal.

Next, for evaluating K₀-value, one-dimensional compression is performed in the DEM simulation with the stress path, similar to that in laboratory triaxial tests. Due to the difficulty of measuring lateral stress in oedometer tests, triaxial tests are usually performed and an initial isotropic stress is applied to stand the sample before the K₀ tests. Therefore, all the samples are initially isotropic consolidated, after which K₀ loading is conducted. Meanwhile, numerical drained triaxial tests are also conducted to determine the strength and critical state parameters. The existing DEM simulation results [10] verified that the K₀ values determined in one-dimensional compression with high vertical pressure are reliable and that the K₀ values obtained from triaxial tests are consistent with those obtained from oedometer tests.

3. Influence of Inherent Anisotropy

3.1. Macroscopic Results

Figure 6 shows the stress path of each specimen during K₀-loading. The number of the legend represents the bedding angle δ of specimen, and for the isotropic specimen, it is labeled by “iso”. Consistent with existing simulation results [6,25] and discussion [14], all stress path lines are slightly bent and then gradually changed into straight lines. Meanwhile, for the specimens with δ = 60°, 45°, and 30°, their stress path lines are close with each other and higher than the line representing the isotropic specimen when the vertical stress is relatively large. In contrast, for the specimens with δ = 90° and 0°, the lines are similar but lower than the “iso” line. This clearly shows the effect of initial anisotropy on the evolution of horizontal stress with vertical stress.

The evolution of two definitions of K₀ values along with stress variations are shown in Figure 7. Consistent with previous expectations, K₀ decreases with the increase of vertical stress, while $K_0^{\prime }$ rise to a peak and then decrease [10]. When the vertical stress is high enough, both $K_0^{\prime }$ and K₀ tend to be stable. However, the value of $K_0^{\prime }$ fluctuates while the value of K₀ is stable. Therefore, the stable values of K₀ at high pressure are treated as the representative K₀ values of the specimen in the following paragraphs. To assess the effect of initial anisotropy on K₀, the change of the representative K₀ values with δ is plotted in Figure 8. K₀ increases first and then decreases with the change of δ. The K₀ value reaches a maximum of 0.343 when δ = 45° and a minimum of 0.308 when δ = 90°. From the above result, the initial anisotropy of the specimen does have a certain effect on the K₀ value, consistent with some previous studies [23]. It’s worth noting that the difference between the maximum and the minimum is only about 10%. Additionally, due to the particularity of the sample preparation method in this simulation, the anisotropy degree of the specimens in this simulation is much higher than that of the naturally generated soil anisotropy (mainly due to gravity deposition) [35]. Therefore, it can be inferred that the difference of K₀ caused by anisotropic principal direction difference in actual engineering or laboratory test is much smaller than that of 10%. The initial fabric anisotropy is not a main influencing factor on K₀ in engineering practice.

3.2. Microscopic Results

Granular materials bear and transmit load through contact between adjacent particles for their particulate nature. Therefore, the microstructure of the granular material, such as the number of contacts, contact force, and distribution, plays a crucial role in the properties of granular soils. The investigation of the microstructure evolution of granular soils under K₀ loading is significant for studying the behaviour of granular soils under oedometer test at the particle level.

3.2.1. Fabric Anisotropy Description

To quantify the initial fabric anisotropy of the specimens and observe the microstructure evolution in the loading process, a second-order fabric tensor [36,37] which is from the statistics of spatial distribution of contact normal is introduced here. The fabric tensor is defined as follows:

$$R_{ij}={\displaystyle {\int }_{\mathsf{\Omega }}E\left(\mathsf{\Omega }\right)n_i}{n}_{j}d\mathsf{\Omega }=\frac{1}{N}{\displaystyle \sum _{c\in N}{n}_i}{n}_j$$

(5)

where n_i is the unit contact normal in the i-direction, N is the total number of contacts, and $E\left(\mathsf{\Omega }\right)$ is the distribution function on the unit sphere $\mathsf{\Omega }$.

Similarly, the distribution of normal contact force and tangential contact force in space can also be expressed by second-order tensor as follows [38,39]:

$$F_{ij}^n=\frac{1}{4\pi }{\displaystyle {\int }_{\mathsf{\Omega }}{\overline{f}}^n\left(\mathsf{\Omega }\right)n_i{n}_{j}d\mathsf{\Omega }}={\displaystyle \sum _{c\in N}\frac{f^n{n}_i{n}_j}{N\left(1+a_{kl}^r{n}_k{n}_l\right)}}$$

(6)

$$F_{ij}^t=\frac{1}{4\pi }{\displaystyle {\int }_{\mathsf{\Omega }}{\overline{f}}^t\left(\mathsf{\Omega }\right)t_i{n}_{j}d\mathsf{\Omega }}={\displaystyle \sum _{c\in N}\frac{f^t{t}_i{n}_j}{N\left(1+a_{kl}^r{n}_k{n}_l\right)}}$$

(7)

where ${\overline{f}}^n\left(\mathsf{\Omega }\right)$ and ${\overline{f}}^t\left(\mathsf{\Omega }\right)$ are the spatial distributions of normal contact force and shear contact force, respectively, and $a_{ij}^r$ is the second-order anisotropy tensor of contact normal, which can be deduced from the deviation of fabric tensor ($R_{ij}^{\prime }$):

$$a_{ij}^r=\frac{15}{2}{R^{\prime }}_{ij}$$

(8)

Then, scalar anisotropy parameters are used to quantify the degree of contact normal anisotropy [36,38,39]:

$$a_r=\sqrt{\frac{3}{2}{a}_{ij}^r{a}_{ij}^r}$$

(9)

Similarly, we use two other parameters a_n and a_t to quantify normal contact force and shear contact force anisotropy:

$$a_{ij}^n=\frac{15}{2}{F}_{ij}^{n^{\prime }}/{\overline{f}}_0^n$$

(10)

$$a_{ij}^t=5F_{ij}^{t^{\prime }}/{\overline{f}}_0^n$$

(11)

$$a_n=\sqrt{\frac{3}{2}{a}_{ij}^n{a}_{ij}^n}$$

(12)

$$a_t=\sqrt{\frac{3}{2}{a}_{ij}^t{a}_{ij}^t}$$

(13)

3.2.2. Evolution of Coordination Number

Contact number is one of the most important characteristics of granular materials and is often represented by the coordination number CN, which is defined as the average number of contacts per particle. In the definition of the mechanical coordination number, only particles with no less than two contacts are used to calculate CN value [40]. Previous studies have shown that void ratio is a macro quantity closely related to CN. Rothenburg and Kruyt [41] discussed the relation between CN and void ratio in the biaxial test simulated by DEM. They pointed out that anisotropy must be considered when establishing a relationship between void ratios and CN.

The evolution of CN under growing stress is shown in Figure 9. Before loading, the CN values of all specimens were very close to each other. We found that the CN value of all samples decreased slightly at the beginning and then increased steadily. This is the same as Khalili’s simulation results [31]. Based on Rothenburg and Kruyt’s [41] theory, the decrease of CN value in the initial stage is induced by the rearrangement of particle contact, producing fabric anisotropy to bear the applied stress. After the completion of rearrangement, CN increases steadily with the densification of the specimen. For specimens with large δ (δ = 60° and 90°), we can find that the decrease of CN is greater than that of other specimens.

3.2.3. Evolution of Anisotropy

The evolutions of the scalar anisotropy parameter of contact normal a_r of different initial anisotropy specimens are illustrated in Figure 10. When δ = 0° and 30°, a_r initially showed an upward trend and then stabilized. However, when δ = 45°, 60°, and 90°, a_r firstly decreased and then tended to be stable. This phenomenon can be explained in this way: under the action of increasing vertical stress, the contact is redistributed, and more contacts are formed in the vertical direction to withstand greater stress. From the previous initial contact normal distribution rose graph, vertical contacts will increase the anisotropy of contact normal for smaller δ while decreasing it for larger δ. The stability of a_r in the subsequent stage indicates less contact redistribution and more frequent variations of contact force under the condition of larger vertical stress during K₀ loading.

The evolution of contact normal force anisotropy parameters a_n and contact tangential force anisotropy parameter a_t are shown in Figure 11 and Figure 12, respectively. In Figure 11, it is observed that the contact force anisotropy parameter a_n increases sharply at first and then stabilizes with the increase of vertical stress. It is worth noting that a_n value of the specimen with δ = 90° is significantly higher than that of other specimens, and the development of a_n value of other specimens is very similar. As for a_t value illustrated in Figure 12, a_t value is very small and fluctuates violently during K₀ loading, but similar to a_n, it also reaches a stable value at the end of the test. Interestingly, a_t values of the specimens with δ = 30° and δ = 45° in the stable phase are the same, while those of the specimens with δ = 60° and δ = 90° are the same.

4. Effects of Particle Shape, Initial Void Ratio, and Inter-Particle Friction Angle

In addition to the above-mentioned simulations concerning the influence of inherent anisotropy, the effects of particle shape, initial void ratio, and inter-particle friction coefficient on K₀-value are studied by a further series of simulations. Meanwhile, variations of these factors can induce changes of strength parameters (peak strength and critical state strength) of the specimen, thus the relationship between K₀ values and strength parameter can be better investigated.

Particle number and particle size distribution of all specimens in this section are identical to those described in Section 3. Three different kinds of shape (shown in Figure 3) were used to investigate the effect of both particle shape and potential particle rotation. Aspect ratio is represented by AR below for conciseness. Specimen formation and consolidation mode are the same as the isotropic specimens in Section 4. The initial void ratio of these specimens is controlled by adjusting inter-particle friction coefficient during the consolidation stage. Summary information about all specimens in these simulations is listed in

Table 2. Summary of the test program.

Aspect Ratio of Particles	Friction Coefficient in Specimen Preparation	Initial Void Ratio
1	0.001	0.588
	0.01	0.598
	0.15	0.670
	0.4	0.701
	0.5	0.744
1.3	0.001	0.493
	0.01	0.530
	0.15	0.591
	0.4	0.646
	0.5	0.727
1.7	0.001	0.499
	0.01	0.560
	0.15	0.600
	0.4	0.709
	0.5	0.775

. After completion of specimen formation, five different inter-particle friction coefficients (0.3, 0.35, 0.45, 0.5, and 0.7) are adopted during loading stage for each specimen. Besides the one-dimensional compression experiment described in Section 4, triaxial tests are also conducted to obtain corresponding macroscopic strength parameters.

Figure 13 and Figure 14 show the variations of K₀-values with void ratio from two aspects. Figure 13 focuses on the impact of inter-particle friction coefficient whereas Figure 14 focuses on particle shapes. According to Figure 7, K₀-value tends to reach a stable value when vertical stress is large enough. Therefore, K₀-values were determined at the end of the loading process (at vertical stress of 900 kPa). The void ratio was determined after consolidation completed before one-dimensional loading.

The results presented in Figure 13 and Figure 14 show an increasing trend of K₀-value, with increasing specimen initial void ratio for specimens with any inter-particle friction coefficient and particle shape in the tested range.

Furthermore, according to Figure 13, K₀-value decreases with inter-particle friction coefficient increasing. This trend matches with existing laboratory tests [23,42]. For real granular material specimens, the difference of inter-particle friction coefficient between particles reflects the difference in mineral composition and particle grinding roundness. The effect of particle shape on K₀ can be observed from Figure 14. Although there is a slight difference in the initial void ratio range corresponding to each particle shape, we can still find that the K₀-value decreases with the increase of the aspect ratio. This trend is well observed in previous literature [23,42]. Guo [23] believed that this trend stemmed from the fact that changes in the aspect ratio of particles caused changes in the branch vector (see its literature for a specific explanation), thus changing the mechanical properties of the material. Lee [42] thought that the change of self-locking effect caused by the change of particle shape was the main reason. The stronger the self-locking effect [43] is, the weaker the relative sliding between particles is, thus less lateral stress is required under the same vertical stress. However, the difference of K₀-value between spherical particle and clump particle with AR = 1.3 is much larger than that between particle with AR = 1.3 and particle with AR = 1.7. Thus, it can be inferred that the change of K₀-value caused by the aspect ratio is more drastic when the aspect ratio is relatively small.

For granular materials, particle shape, inter-particle friction coefficient, and initial void ratio control their strength parameters (peak strength and critical state strength). By the combination of the variations of these three factors, enough wide range of strength parameters can be obtained to study the correlation between strength parameters and K₀-value. Figure 15 shows the variation of K₀-values with the critical state friction angle and peak friction angle for all specimens in this section. The prediction line of the Jáky’s [44] simplified formula is also drawn in the figure. From the figure, it is more general to adopt the peak friction strength parameter in this simplified formula. The critical friction strength is constant as the initial void ratio changes, contradicting the phenomenon that K₀-value varies with initial void ratio. However, there is still a certain gap between the measured values and the prediction results using peak strength parameters. Northcutt and Wijewickreme [11], Michalowski [4], and Wanatowski and Chu [7] have also indicated that Jáky’s formula does not provide a concise estimation of K₀-value.

Therefore, we do not think K₀ can be predicted precisely enough by just one strength parameter. Based on the results presented in this section, we try to give the following formula for estimating the K₀-value of granular materials as a function of critical state friction angle (φ_cs), particle shape, and the relative density of initial state (Dr). In the calculation of Dr, $e_{max}$ and $e_{min}$ are obtained by the method proposed by Wood and Maeda [45]:

$$K_0=S(15-Dr-15\times sin{\phi }_{cs})-1$$

(14)

where S is the particle shape factor. For particles with AR = 1.0, 1.3, and 1.7, the shape factor is 0.154, 0.179, and 0.183, respectively.

In the formula above, we adopt critical state friction angle instead of peak friction angle, because the later one is affected by too many factors such as initial void ratio, confining pressure, mineral composition, etc., to establish a more accurate formula of multiple factors. The critical state friction angle is generally accepted to be controlled only by mineral composition [46]. Additionally, relative density was used to describe the initial state characteristics of the material [47], and particle shape parameter S was used to characterize the effect of particle shape on K₀ value. Figure 16 illustrates the fitting effect of the formula. As we can see from the figure, the proposed formula is accurate enough to predict the DEM simulation results in this section.

5. Conclusions

DEM simulations were carried out to investigate the factors influencing the K₀ value of granular materials, including the initial fabric anisotropy, initial void ratio, inter-particle friction coefficient, and particle shape. The main findings of this study can be summarized as follows:

(a) When the initial anisotropy degree is high, the initial anisotropy does have a certain effect on the value of K₀, but the effect is limited. We consider that the effect of initial anisotropy on K₀ can be neglected under laboratory or engineering conditions as the initial anisotropy degree is smaller in these conditions.

(b) The evolution of microstructure (CN, a_r) in specimens with different initial anisotropy shows different trends during the K₀ loading. The larger the angle between the principal direction of the initial anisotropy and the principal stress is, the larger the fabric redistribution will occur in the initial stage of K₀ loading.

(c) The shape of particles, the initial void ratio, and the coefficient of inter-particle friction all have significant effects on K₀ values.

(d) It is more appropriate to adopt the peak friction angle in Jáky’s simplified formula for calculating K₀ values. However, this equation still cannot provide accurate predictions. Based on the simulation results in this paper, we establish a prediction formula of K₀ values which relate to initial relative density, shape parameter, and critical state friction angle. The new proposed formula gives better results of K₀.