Discrete element simulations of two wetting effects on granular materials

To explore the microscopic responses of granular materials to wetting, the inter-particle lubrication effect and particle breakage in an odometer were simulated using a two-dimensional discrete element method. The lubrication effect was modeled by reduction of the inter-particle friction coefficient and particle breakage was initiated by decreasing the particle strength. Once the strength of the particles decreased to a threshold value, the particles began to break so that new contacts could be established to transfer the external loads. Numerical simulations successfully reproduced the additional compaction of the material and the intensification of the horizontal stress in addition to the microscopic responses of the granular assemblies. The microscopic interpretation of the earth pressure coefficient at rest and the evolution of the grain number distribution during particle breaking were also investigated.

Granular materials are extensively used in civil engineering applications such as the foundations of railways and the shells of rockfill dams. These materials are exposed to changes in climate and the environment, which inevitably leads to changes in their long-term mechanical behavior [1]. One typical environmental factor that deserves special attention when designing a granular structure is the interaction between the particles and water from both underground and surface runoff. Water infiltration into the inter-particle pores and the micro-cracks in particles has at least three effects, i.e. the lubrication effect, loss of matric suction and particle breakage [2,3]. These factors often occur in combination and lead to the so-called wetting-induced collapse [3][4][5]. Alonso et al. [3] distinguished the predominant factors according to the sizes of the particles and pointed out that the wetting-induced collapse is primarily caused by particle breakage for coarse granular materials such as rockfill and gravel, while the loss of matric suction is responsible for the collapse of fine materials with grain sizes ranging from 10 to 0.1 mm, such as sands. Both particle breakage and the loss of suction overwhelm the effect of inter-particle lubrication.
The wetting-induced collapse may result in unfavorable deformation and stress deterioration in structures; therefore, this phenomenon has attracted a great deal of attention over the past three decades. Many experimental studies [6][7][8][9] have been conducted and great efforts have been devoted to constitutive modeling of this particular behavior of granular materials with the goal of understanding and explaining phenomena observed in the engineering field and numerically solving the boundary value problems in practice [10][11][12]. Among these experimental and theoretical investigations, a phenomenological approach has generally been employed and the macroscopic responses of materials were described within the framework of continuum mechanics [5,[10][11][12]. Nevertheless, the microscopic responses on the particle level have not been explored in depth due to limitations of the experimental apparatus and idealization of the granular assembly as a continuum material.
The rapid development of computational capacity and numerical techniques such as the discrete element method (DEM) [13][14][15][16][17][18][19] have provided an effective approach for investigation of microscopic responses of granular materials during wetting. As a conventional topic in unsaturated soil mechanics, the suction effect during wetting has been investigated by many authors using DEM [20][21][22][23][24], and the results obtained have shown that DEM is a rather powerful tool in reproducing the behavior of unsaturated soils. DEM simulations of particle breakage during loading have also been conducted with an emphasis on its influence on the shear and compression behavior of granular materials [25][26][27]. However, the lubrication effect and particle breakage under constant effective stress induced by water infiltration have attracted little attention to date.
The objective of the present study was to reveal the wettinginduced responses of granular materials on the particle level. To accomplish this, the lubrication effect and particle breakage under constant stress in an odometer were investigated separately using the two-dimensional discrete element method. In the numerical calculations, the lubrication effect was simulated by reducing the inter-particle friction coefficient and particle breakage was initiated by decreasing the tensile strength of the particles. Particular attention was paid to the polar distributions of the contact angle and the contact forces before and after wetting. The grain size (number) distribution after particle breakage was also studied statistically. Figure 1 shows the scheme for the explicit DEM simulations used in this study. At any time instant, t, the overlaps between particles are identified and the contact forces are calculated using the (incremental) force-displacement relationships. The total unbalanced forces and the momentum of particles can be obtained by summing the contact forces. New translational and angular accelerations are then evaluated using Newton's second law. Numerical integration of the obtained accelerations using a small time increment, t, gives the velocities at the instant t+t/2. Eq. (1) gives the iterative formulations of the translational and angular velocities [14,19]: Herein, F x , F y and M z denote the unbalanced forces (momentum); g x and g y are components of the gravitational acceleration; m and I denote the mass and inertia of particles; ,  and  are the global damping coefficients, which are introduced to dissipate the energy so that the specimen can reach a stable state. The velocities obtained from eq. (1) can be numerically integrated again to provide the displacement from t to t+t. Having obtained the new configuration of the granular assembly, another cycle of calculations is repeated to forward the simulation until the process under investigation terminates.

The techniques for DEM simulations
In most of the DEM simulations reported in previous literatures, disks or balls were generated to model the grains of a specimen owing to the convenience and high efficiency of this geometrical simplification technique in identifying the contact status. In the present study, disks with four different radiuses, 5, 4, 3 and 2 mm, were generated to prepare the specimens. The use of disks also enabled us to model the particle breakage in a way similar to that suggested by Guerrero and Vallejo [26,27]: (1) Only those particles with a coordination number N c (the number of contact points) greater than or equal to 2 (N c 2) could break. For a two-dimensional granular assembly, those particles with a coordination number less than 2 were considered unstable and often carried very small contact forces.
(2) The normal contact forces acting on a potential breakable particle were simplified to diametrical forces with a magnitude equal to the maximum normal contact force, which is denoted by F 1 in Figure 2 [26].
(3) The tensile stress induced by the diametrical forces was calculated according to 1 c , in which R and L denote the radius and width of the concerned particle; and f c reflects the fact that a higher coordi-Figure 2 Idealized model for particle breakage [26]. (a) Normal contact forces; (b) forces simplification; (c) particle fragmentation.
nation number results in a better stress condition. McDowell [28] suggested a discount factor for three-dimensional problems, f c = (N c 2) a (a > 0). Herein, for two-dimensional problems, we used in which a is a dimensionless constant and N c denotes the coordination number. For the case N c = 2, f c = 1.0.
(4) The tensile strength of the particle was assumed to be inversely proportional to its radius, i.e.
Herein, R r and T r are the radius and tensile strength of a predefined reference particle; b is another constant. The reduction of the particle strength with radius has been validated by many experiments [29,30] and can be explained by the presence of flaws, i.e. larger particles tend to contain more and larger flaws, which results in more severe stress concentration and lower strength.
(5) Once the tensile stress obtained from eq. (2) exceeded the strength evaluated according to eq. (4), the particle fragmented into eight smaller particles as shown in Figure 2. The original particle was then deleted and eight smaller particles were inserted into the particle list. The contacts related to the breaking particle were also removed from the contact list and the new contact statuses for the deriving particles and the surrounding particles were established during subsequent cycles of calculations. The particle breakage model introduced above differs slightly from that proposed by Guerrero and Vallejo [26,27] in that only particles with a coordination number of two or three can break in the latter. The introduction of the factor, f c , in eqs. (2) and (3) made the influence of the coordination number on the stress state within a particle easy to be considered. Despite the simplicity of the highly conceptual model, it is capable of capturing the progressive breakage of brittle particles and its validity for reproducing the mechanical behavior of breakable granular materials under compression and shearing has been verified by Guerrero and Vallejo [26,27].
There are a total of 11 parameters required to run the DEM simulations, namely, the normal and tangential stiffness, k n and k t , the global damping coefficients, ,  and , the inter-particle friction coefficient, , the density of particles, , the constants a and b, and the radius and strength of the reference particles, R r and T r . The choice of DEM parameters is of great importance in reproducing the mechanical properties of granular materials, especially when the results of simulation tests are compared with those of actual laboratory tests quantitatively. Nevertheless, a large amount of DEM simulations reported in the literatures were conducted using empirically evaluated constants [15][16][17][18][24][25][26][27]31]. Since the objective of the present study was to qualitatively investigate the effects of wetting on the responses of granular materials, the parameters were also tentatively chosen on an empirical basis. Table 1 shows the values of these parameters used in the following quasi-static simulations.

The lubrication effect during wetting
To evaluate the lubrication-related response of materials, a specimen composed of disks with radiuses of 4 mm, 3 mm and 2 mm was first compressed up to a vertical stress of 0.4 MPa (point a in Figure 3), 0.8 MPa (point c) and 1.2 MPa (point e) with an inter-particle friction coefficient of 1.482. The vertical stress was then kept constant and the inter-particle friction coefficient was shifted to a value of 0.286. The reduction of the inter-particle friction coefficient triggered a relative sliding and rearrangement of the particles, as indicated by the reduction of the void ratio shown in Figure 3.
To study the influence of the sequence of loading and lubrication, the compression curves of specimens with an inter-particle friction coefficient of 1.482 and 0.286 were also depicted in Figure 3 using solid and dashed curves, respectively. The magnitude of the reduction of the void ratio increased as the vertical stress increased. However, the final states after lubrication did not lay on the compression curve of the lubricated specimen ( = 0.286). These findings seem to indicate that the sequence of loading and lubrication should be taken into consideration.
As shown in Figure 4, the reduction of the inter-particle friction coefficient also resulted in an increase in the horizontal  The strain induced by inter-particle lubrication.

Figure 4
The horizontal stress intensification induced by the inter-particle lubrication.
stress. However, the final stress states did not coincide with the stress path of the lubricated specimen either. Using the empirical relationship between the friction angle  and the earth pressure coefficient at rest (K o ) suggested by Jaky [2], i.e. K o  1sin, we can infer that the friction angle  decreased when the inter-particle friction coefficient was reduced from 1.482 to 0.286. Similar conclusions were also drawn by Liu and Matsuoka [32] based on simulations of simple shear experiments. Figure 5 shows a plot of the probability of the contacts as a function of the contact angle  at state e and state f (Figures 3 and 4). Before the inter-particle lubrication (state e), the polar distribution diagram (Figure 5(a)) showed a slight preference for orientation to the vertical direction. However,  Figure 6 shows the distributions of the normal and the tangential forces before and after lubrication (state e and state f, respectively). Although the reduction of the interparticle friction was considerable in the current DEM simulations, only a slight change in the polar distributions of the average normal contact forces was observed, as shown in Figure 6(a) and (c). Conversely, the average tangential contact forces underwent an abrupt decrease due to the lubrication effect as shown in Figure 6(b) and (d). It is also interesting to note that the average normal contact force followed an elliptical distribution, while the average tangential contact force followed a rose distribution as follows:    (Figure 6), we can express the stress components using these statistical quantities. Taking the vertical stress T 11 acting on the horizontal plane as an example (Figure 7), it can be calculated from these statistical quantities as follows: Herein, f (  ) is the probability of contacts with a contact angle of , c N denotes the number of contacts per unit length and the following assumption is implicated, namely, the specimen is homogenous and the distributions of contact angle and average contact forces along the horizontal plane are the same as those obtained for the entire specimen. This assumption is somewhat similar to the assumption that the surface porosity is identical to the volume porosity, as used for porous media [33]. Similarly, the horizontal stress T 22 acting on a given vertical plane can be expressed as follows: Therefore, the earth pressure coefficient at rest reads Figure 7 The calculation of the vertical stress using contact forces. Figure 5 indicates that the slight fabric anisotropy can be neglected and the contact angle is almost evenly distributed within [0°, 180°] or [90°, 90°]. Adopting this assumption and substituting eqs. (5) and (6) into eq. (9), one can obtain the analytical form of eq. (9) as follows:

Reinspection of
in which x = b/a, y = c/a. At state f in the current case, a, b and c were estimated to be 4.5 kN, 9.0 kN and 0.96 kN, respectively. Consequently, the earth pressure coefficient at rest was approximately 0.67, which is of the same order as that determined from Figure 4 (K o ≈ 0.60). Figure 8 plots the contours of K o predicted by eq. (10) in the two-dimensional space (x, y). For a given x (the ratio of b to a), the decrease of y, which can be attributed to the lubrication effect here, will result in an obvious increase in the earth pressure coefficient at rest. The increase of x for a given y will also result in a change of K o . However, the sensitivity of K o to x is not as evident as that to y, especially when y is close to 0.5. Figure 8 also indicates that the necessary conditions for K o = 1.0 for granular materials under confined compression are frictionless contact between particles (y→0) and a circular distribution of normal contact force (x→1). A particular case in point is water, which can be interpreted as an aggregate of uncountable tiny particles without inter-particle friction.

Particle breakage induced by wetting
The infiltration of water into granular materials not only results in degradation of the inter-particle frictional strength, but also degradation of the particle strength due to the complex interaction between water and minerals. Once the strength of the particles decreases to the magnitude of stress induced by external loadings, particles start to break and the materials undergo additional volume compaction. To reproduce this phenomenon, a specimen composed of disks with radiuses of 5 mm, 4 mm and 3 mm and a referential strength of 10 MPa was first compressed to a vertical stress of 0.4 MPa (point a in Figure 9), 0.8 MPa (point c) and 1.2 MPa (point e) with an inter-particle friction coefficient of 0.286. The vertical stress was then kept constant and the tensile strength of the reference particle was reduced to a value of 1.0 MPa.
The breakage-related collapse was simulated using a breaking-rearranging logic. First, the breaking criterion was checked for all particles in a stable step, and those particles meeting the criterion were then replaced by eight new particles. Second, conventional cyclic calculations were performed with the new granular assembly composed of the original (unbreakable) particles and the deriving particles until a new static configuration was achieved. The entire breaking process was modeled by alternately implementing the above two operations until a final stable state without particle breakage was reached. Since the number of particles will increase considerably and very tiny particles will emerge due to particle breakage, the number of particles was initially small and the mass scale technique [15,31] was adopted so that the numerical simulations could be finished within an acceptable timescale.
As shown in Figure 9, the reduction in particle strength triggered grain breakage and resulted in an abrupt decease of the void ratio (e.g. c→d and e→f). However, a threshold stress under which no particle breakage would occur existed (point a). Comparison of Figure 9 and Figure 4 also indicates Figure 9 The strain induced by particle breakage.

Figure 10
The horizontal stress intensification induced by particle breakage. that the collapse induced by particle breakage was much more evident than that induced by lubrication, both of which increased with applied vertical stress. Moreover, as shown in Figure 10, a reduction of the particle strength will also lead to an intensification of the horizontal stress T 22 if particle breakage presents in the specimen. Figure 11 shows the distributions of the contact angle before and after particle breakage. An obvious preference of contact orientations toward the vertical and horizontal directions is clearly shown in Figure 11(a), which indicates the existence of initial fabric anisotropy. However, this fabric anisotropy was alleviated as a result of particle breakage, which tends to result in an isotropic distribution of the contact angle, as shown in Figure 11(b). Particle breakage also leads to substantive decreases in the average normal and tangential contact forces in all directions, as shown by comparing the polar distributions of the average contact forces before and after the particle breakage in Figure 12. Figure 13 shows the grain size distributions before and after the particle breakage under two vertical stresses. The Figure 11 The polar distributions of the contact angles before and after particle breakage. (a) State e; (b) state f.

Figure 12
The polar distributions of the normal and tangential contact force before and after particle breakage. (a) Normal contact force (state e); (b) tangential contact force (state e); (c) normal contact force (state f); (d) tangential contact force (state f).

Figure 13
The grain size distributions of the specimen before and after particle breakage. results revealed that the grain size distribution broadened considerably as a result of particle breakage. Furthermore, the grain size distribution of the specimen wetted under a vertical stress of 1.2 MPa (point f in Figure 9) was broader than that of the specimen wetted under a vertical stress of 0.8 MPa (point d). These findings indicate that particles break more sufficiently under a higher vertical stress, which leads to a higher content of fine particles. In fact, when the wetting-induced breakage was finished, the original 90 grains fragmented into approximately 1500 smaller particles and 2500 smaller particles under the vertical stress of 0.8 MPa and 1.2 MPa, respectively. Figure 14 shows the configurations of the specimen at different states. As shown in Figure  14(b) and (c), large particles were surrounded and protected Figure 14 The grain size distributions of the specimen before and after particle breakage. (a) Initial state; (b) state d; (c) state f. by many tiny particles emerging from particle breakage and were less vulnerable to crushing. Conversely, tiny particles generally had less contacts with neighboring particles and were therefore more prone to suffer progressive breakage. These results are in good agreement with the experimental findings of Nakata et al. [30] and the theoretical analysis conducted by McDowell [34].
McDowell [28,34] suggested that the geometry of the broken specimens tends to be self-similar or fractal, which can be mathematically expressed as Herein, N (d > D) is the number of particles with a diameter larger than D; A is a proportional constant and B is the fractal dimension. B is close to 2.5 for most materials under pure crushing in three-dimensional problems [34]. For twodimensional problems, McDowell proposed the following empirical relationship [28]: in which m is the so-called Weibull modulus, which indicates the variability of the particle strength. For most materials, m ranges from 5 to 10. Thus, the fractal dimension B is expected to range from 1.2 to 1.4 in two-dimensional cases. Figure 15 shows the final distributions of the particle numbers when the wetting-induced breakage was completed. The distribution of the numbers of particles with a diameter greater than 0.5 mm could be approximated well by a fractal distribution, and the fractal dimension was calibrated to be 1.36, which lies in the expected range. The results also show that, as grains broke more sufficiently, the actual particle number distribution became closer to the theoretical fractal distribution. Therefore, it is reasonable to postulate that there is an ultimate fractal distribution of the grain numbers that can be approached when the particles break sufficiently. This assumption would be very useful when establishing a constitutive model for breakable granular materials, and it has already been used by Einav [35,36] in Figure 15 The grain number distributions of the specimen after particle breakage.
construction of the theory of breakage mechanics.

Conclusions
Two wetting effects on granular materials in an odometer, the inter-particle lubrication and the particle breakage, were investigated using the discrete element method. Although the vertical stress remained constant during wetting, additional strain and rearrangement of the particles were observed and the microscopic responses of the materials were investigated in both cases. The conclusions are summarized as follows.
(1) Both the lubrication effect and the particle breakage initiated by wetting resulted in additional compaction of the material and intensification of the horizontal stress in the odometer. However, the magnitude of the strain induced by particle breakage exceeded that induced by the lubrication. In both cases, the strain increased as the vertical stress increased.
(2) Both the lubrication effect and the particle breakage eliminated the preference for the contact orientation or the fabric anisotropy, and tended to yield an isotropic distribution of the contact angle. However, the disturbance induced by particle breakage was more evident than that induced by inter-particle lubrication.
(3) The wetting-induced inter-particle lubrication resulted in an abrupt decrease of the average tangential contact forces, while an evident decrease in the average normal contact forces was also observed as a result of the particle breakage.
(4) Microscopic investigation of the earth pressure coefficient at rest (K o ) indicated that the necessary conditions for K o = 1.0 for granular materials under confined compression are frictionless contact between particles and an isotropic distribution of the normal contact force.
(5) Statistical investigations of the number of particles after breakage confirmed the existence of a theoretical frac-tal distribution for the grain numbers. However, the fractal distribution could only be achieved when the particles broke sufficiently.