DEM Investigation of Particle-Scale Mechanical Properties of Frozen Soil Based on the Nonlinear Microcontact Model Incorporating Rolling Resistance

Although frozen soil is in nature the discrete material, it is generally treated as the continuummaterial. The mechanical properties of frozen soil are so complex to describe adequately by conventional continuum mechanics method. In this study, the nonlinear microcontact model incorporating rolling resistance is proposed to investigate the particle-scale mechanical properties of frozen soil. The failure mechanism of frozen soil is explicated based on the evolution of contact force chains and propagation of microcracks. In addition, the effects of contact stiffness ratio and friction coefficient on stress-strain curve and energy evolution are evaluated. The results show that the nonlinear microcontact model incorporating rolling resistance can better describe the experimental data. At a higher axial strain, the contact force chains near shear band which can give rise to the soil arch effect rotate away from the shear band inclination but not so much as to become perpendicular to it. The propagation of microcracks can be divided into two phases.The stress-strain curve is strongly influenced by contact stiffness ratio. In addition, friction coefficient does not significantly affect the initial tangential modulus. Compared with frictional coefficient, the effect of contact stiffness ratio on stress-strain curve and energy evolution is greater.


Introduction
In the past 50 years, a large number of engineering constructions in cold regions and artificial ground freezing projects have been built throughout Europe and East Asia [1].There is already growing evidence, on the basic of ongoing researches, that the mechanical properties of frozen soil are crucial to the stability of infrastructures constructed in cold regions [2].It is well known that frozen soil is composed of solid mineral particles, air, unfrozen water, and ice.Therefore, frozen soil is in nature the discrete material.However, it is generally treated as the continuum material [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19].
Although the continuum mechanics method plays an important role and is widely used, the mechanical properties of frozen soil are so complex that the discrete features of frozen soil should not be neglected.The Discrete Element Method (DEM) is regarded as a powerful tool to investigate the particle-scale mechanical properties of frozen soil due to the fact that DEM can control the complex responses of an assembly of discrete materials by very simple contact laws [20][21][22][23][24][25][26][27].The microcontact model is the key to investigating the particle-scale mechanical properties of granular materials accurately.Three kinds of microcontact models were widely used to study the particle-scale mechanical properties of granular materials: the linear microcontact model [28], the nonlinear microcontact model [29], and the simple bond microcontact model [30].It is widely known that the particles can roll after the bonds are broken.However, only a few investigations have been conducted on the effect of rolling resistance on particle-scale mechanical properties of granular materials [31][32][33].Hence, a better understanding with respect to the particle-scale mechanical properties of frozen soil based on the microcontact model incorporating rolling resistance is a subject of prime importance.

Mathematical Problems in Engineering
It should be pointed out that although few comprehensive studies on DEM investigation of particle-scale mechanical properties of frozen soil have been conducted, a lot of attempts on particle-scale mechanical properties of unfrozen soil have been done.Zhao et al. [34] carried out a serious of numerical triaxial test on the dense sand using the threedimensional particle flow code with a linear contact model and the influence of normal stiffness on stress-strain curve of sand was analyzed.Ng and Dobry [35] explored the use of discrete element simulations to model granular soil response and the results indicated that the stress-strain curve observed in laboratory tests on sands was well reproduced by the numerical simulations.Hosseininia [36] pointed out that the evolution of normal contact force chains revealed the development of new contacts.Jiang et al. [37] investigated the propagation of microcracks of sands by two-dimensional distinct element method.The results show that the number of microcracks reached its maximum value during strain softening.Zhou et al. [38] found that the energy allocation manner change extraordinarily within the small-strain range while the relatively steady energy allocation can be achieved at large strain.
The main objectives of this paper are to (1) present the nonlinear microcontact model incorporating rolling resistance for frozen soil; (2) explicate the failure mechanism of frozen soil based on the evolution of contact force chains and propagation of microcracks; (3) evaluate the effects of contact stiffness ratio and friction coefficient on stress-strain curve and energy evolution.

Nonlinear Microcontact Model for Frozen Soil
The microcontact model in this study is similar to that proposed by Jiang et al. [32] and consists of a normal contact model, a tangential contact model, and a rolling contact model, as shown in Figure 1.
The basic mechanical elements are spring, bond, slider, divider, dashpot, and roller.The elastic relationship between contact force and relative displacement is reflected by spring.The cementation at contact is reflected by bond.The frictional sliding is reflected by slider.The divider reflects that the tensile force between particles is zero once the bonds are broken.The energy dissipation is reflected by dashpot.The roller reflects the rolling resistance when a particle rotates against another particle through a small angle.
The normal contact force, tangential contact force, and rolling contact moment are calculated as follows.
2.1.Normal Contact Force. Figure 2 presents the mechanical response of normal contact model.The normal contact force under compression can be described by where   is the normal stiffness;   is the normal displacement.
Before the bond breaks, the normal contact force under tension can be also described by (1).Once the bond breaks,  the particles are separated and the normal contact force under tension abruptly reduces from normal bonding strength  nb to zero.

Tangential Contact Force.
The mechanical response of tangential contact model is presented in Figure 3.
Before the bond breaks, the tangential contact force can be expressed by where   is the tangential stiffness;   is the tangential displacement.
Once the bond breaks, the two contacting particles may slip relative to one another and the tangential contact force  under tension abruptly reduces from tangential bonding strength  tb to friction associated with normal contact force and friction coefficient .

Rolling Contact Moment.
The rolling contact moment can be divided into two parts: the rolling contact moment before the bond breaks and that after the bond breaks.

Rolling Contact Moment before Bond Breaks.
For the conventional microcontact model which is widely used to investigate the particle-scale mechanical properties of granular materials, two adjacent particles are in contact at a discrete point.The ESEM micrograph of bonded granular materials shows that two adjacent particles are in contact over a width [39].For the microcontact model in this paper, two adjacent particles are in contact over a width.We assume that a series of bonds are distributed continuously over the bond contact width B, as shown in Figure 4. Particle  rotates counterclockwise against particle  through a small angle  before bond breaks.Since  is extremely small, the contact force at the lefthand side   and right-hand side   can be calculated by (3), as shown in Figure 5: ( The rolling contact moment before bond breaks is determined by

Rolling Contact
Moment after Bond Breaks. Figure 6 shows that particle  rotates counterclockwise against particle  through a small angle  after bond breaks.
The rolling contact moment after bond breaks can be divided into two parts: the rolling contact moment purely due to bond   and that purely due to normal force   .
(1) Rolling Contact Moment Purely due to Bond.When the rolling contact moment is purely caused by bond, the maximum contact force and minimum contact force can be calculated by ( 5), as shown in Figure 7: where  is the contact width that is still distributed with the remaining bonds.Since  is extremely small, the rolling contact moment is determined by Let  = , the rolling contact moment can be written as follows: (2) Rolling Contact Moment Purely due to Normal Force.
The rolling contact moment purely due to normal force is analogous to that purely due to bond and can be described by The total rolling contact moment is a linear summation of the rolling contact moment before bond breaks and that after bond breaks for simplicity.So the total rolling contact moment can be calculated as follows: Equation ( 9) shows that the rolling contact moment changes nonlinearly with  and the nonlinearity can result in poor computational efficiency.To solve the problem, it is necessary to develop the mechanical response of rolling contact model in a simple way.We assume that the rolling contact moment does not move along the curve ACD but moves along the polyline ABCEF after the bond breaks, as shown in Figure 8.It is reasonable to assume that the energy dissipation in the simplified model is equal to that in the original model.Hence, the area of ABC is equal to that of CDFE based on the energy equivalence principle.
Hence, the total rolling contact moment can be expressed by Figure 9: DEM sample.

Model Calibration.
For the calibration of microcontact model for frozen soil, the DEM sample is generated according to the test conditions of our previous laboratory experiment on the mechanical characteristics of frozen soil in [40] and then the simulated results of microcontact model are compared with the test data in [40].

DEM Sample Generation and Loading
Procedure.The microcontact model for frozen soil is implemented in the two-dimensional particle flow code (PFC2D) provided by Itasca Consulting Group.The processes of DEM sample generation are performed as follows.(1) Soil particles are randomly generated inside the rectangular frictionless boundaries which are simulated by four rigid walls.The system then arrives at the equilibrium state by using an equilibrium ratio limit.Note that the grain size distribution curve of soil in laboratory tests is adopted as that in numerical tests and the assembly of soil particles has a width of 39.1 mm and height of 80 mm.(2) The DEM sample is subsequently subjected to an isotropic constant confining pressure.In this stage, a slight change in sample height which can be neglected may occur.(3) A few "floating" soil particles are eliminated to obtain a denser network of contact, as shown in Figure 9. (4) The bonds between all soil particles are then installed throughout the assembly with the bond parameters.
(5) Finally, the DEM sample is vertically compressed while the confining pressure acting on the left and right walls is kept unchanged.During this stage, the energy evolution and the number of microcracks are monitored.It should be noted that the loading rate is 1.072 mm/min and the confining pressure is 0.2 MPa.

Model Calibration.
It is well known that not only the internal structure of frozen soil changes constantly but  also the physical and mechanical properties of ice are not constant during the loading process.The strain is considered as an important index reflecting the dynamic process.It is reasonable to assume that the normal bonding strength  nb changes with strain  in the loading process and other microparameters are invariable for simplicity.Hence, the nonlinear microcontact models for frozen soil are adopted.
A series of triaxial compression tests on the static mechanical properties of frozen soil were performed to provide the experimental data for model calibration.The laboratory tests are introduced briefly here in order to explain the numerical results further.The soil used in laboratory tests is Fujian standard sand.The grain size distribution curve of Fujian standard sand is presented in Figure 10.Other physical indices such as maximum dry density, coefficient of curvature, and maximum void ratio are listed in Table 1.
The test device is an electric-hydraulic servo-controlled triaxial material testing machine which is equipped with an automatic numerical control system and a data collection system, as shown in Figure 11.The triaxial compression test conditions of specimen DTJwd2 are summarized in Table 2 in detail.The stress-strain curve of specimen DTJwd2 is shown in Figure 12 and more details about our previous laboratory tests can be found in [40].
Due to the restrictions of experiment equipment, the microparameters such as bonding strength cannot be obtained easily even by the advanced technologies such as X-ray, the stereophotogrammetric technique and particle image velocimetry.The "trial and error" is considered as the best way to determine the microparameters at present.The microparameters such as normal bonding strength in this paper are determined by "trial and error."In this study, the normal bonding strength in linear model is 10 kN and the  relationship between  nb and  in nonlinear model can be defined by Other microparameters of nonlinear and linear models are listed in Table 3.The comparison of our previous test data of specimen DTJwd2 in [40] with the simulated results of nonlinear and linear models is shown in Figure 12.
Before the deviatoric stress reaches its peak value, for the nonlinear and linear models, the specimens exhibit strainsoftening behavior which can be found from the laboratory test results and the stress values are close to the experimental data for a certain strain.Hence, the nonlinear and linear models can better describe the laboratory data before the deviatoric stress reaches its peak value.However, it is obvious that the stress value of laboratory test for a certain strain is greater than that of linear microcontact model and is close to that of nonlinear microcontact model after the deviatoric stress reaches its peak value.The conclusion that the simulated results of nonlinear microcontact model show a good agreement with the test data is achieved.The reason is that the normal bonding strength does not vary greatly before the deviatoric stress reaches its peak value, as shown in Figure 13.When the deviatoric stress reaches its peak value, the normal bonding strength varies dramatically since the change of internal structure of frozen soil is major.Then, the normal bonding strength gradually increases as the internal structure of frozen soil gradually reaches the steady state.

DEM Investigation of Particle-Scale Mechanical Properties of Frozen Soil
3.1.Failure Mechanism of Frozen Soil.The evolution of contact force chains and propagation of microcracks can well account for the failure mechanism of frozen soil since the evolution of contact force chains can illustrate how the applied loads are transmitted within frozen soil and the formation of shear bands in frozen soil is associated with the propagation of microcracks.

Evolution of Contact Force
Chains.The contact force chains can be divided into two parts: the strong and weak contact force chains.Almost all the applied loads transmit along the strong contact force chains.Figures 14(a)-14(d) present the contact force chains at various axial strains (1%, 4%, 9%, and 15%).The contact force chains are plotted as lines with thickness proportional to contact force magnitude.At the initial state ( = 1%), the contact force chains are distributed uniformly throughout the frozen soil due to the relatively small applied loads, as shown in Figure 14(a).Furthermore, an obvious webbed pattern of weak contact force chains can be observed.Hence, all soil particles can jointly bear the applied loads.When the axial strain reaches 4%, the thickness of contact force chains increases due to the increase of applied loads and the columnar pattern of strong contact force chains takes the place of the webbed pattern of weak contact force chains, as shown in Figure 14(b).In addition, the contact force chains are distributed mainly in the vertical direction.At a higher axial strain (9%), although the columnar pattern is still the main pattern of contact force chains, the contact force chains near the upper and lower parts of the specimen are not oriented along the vertical direction and rotate away from the shear band inclination, as shown in Figure 14(c).The contact force chains near shear band can give rise to the soil arch effect and the increase of porosity.At the final state ( = 15%), several shear bands can be observed and the contact force chains near the shear bands are almost perpendicular to them.Among these shear bands, only one becomes extraordinarily thick which is in good agreement with the experimental results, as shown in Figures 14(d) and 14(e).It is worth noting that although other shear bands are not observed from experimental results, the existence of these shear bands is undeniable.

Propagation of Microcracks.
It is common knowledge that there are two kinds of microcracks in frozen soil: original cracks and microcracks caused by the breakage of bond (ice).In this paper, the original cracks are not taken into account for simplicity.Figures 15(a)-15(d) presents the propagation of microcracks at various axial strains (1%, 4%, 9%, and 15%).
At the initial state ( = 1%), it is obvious that only a few isolated microcracks can be observed throughout the frozen soil due to the relatively small applied loads, as shown in Figure 15(a).When the axial strain reaches 4%, although the number of microcracks increases, the amount of microcracks is still comparatively small, as shown in Figure 15(b).At a much higher axial strain (9%), the number of microcracks increases dramatically.Moreover, the microcracks are concentrated within some narrow regions, tending to coincide with the shear bands, as shown in Figure 15(c).The reason is that a huge number of bonds are broken due to the increase of contact force chains near the shear bands.Furthermore, the microcracks grow and coalesce to form some larger cracks.At the final state ( = 15%), the larger cracks develop to an array of continuous fracture planes.Among these fracture planes, only one becomes a macroscopic crack which is in good agreement with the experimental results, as shown in Figures 15(d) and 15(e).It is worth noting that although other fracture planes are not observed from experimental results, the existence of these fracture planes is undeniable.
Figure 16 shows the number of microcracks at various axial strains.It becomes clear that the propagation of microcracks can be divided into two phases.In the first phase (OA), point O indicates the initial state before loads are applied; point A can be defined as the "yielding point" of propagation of microcracks since this is where a large number of microcracks begin to be observed.In the second phase (AB), a roughly parabolic increase of the number of microcracks is observed.The number of microcracks increases from 232 to 1405 which means a 505 percent growth as the strain increases from 4% to 9%.Moreover, the final number of microcracks is 1940 at the end of test.

Influence of Contact Stiffness
Ratio.The contact stiffness ratio  is the ratio of normal contact stiffness to tangential contact stiffness.The influence of contact stiffness ratio on stress-strain curve is studied by three numerical tests with different contact stiffness ratios (1, 1.3, and 1.5), as shown in Figure 17.It is evident that the stress-strain curve is strongly influenced by contact stiffness ratio.The stiffness of sample increases with the increase of contact stiffness ratio.Figure 17 also demonstrates that the peak stress decreases from 8.5 MPa to 6.7 MPa which means a 21 percent decline as the contact stiffness ratio decreases from 1.5 to 1.One reason may contribute to the occurrence of this phenomenon: the higher contact stiffness ratio can lead to the greater capability of resisting normal deformation which is presented as low compressibility of sample at macro level.Furthermore, little change of residual stress can be found at the end of test.

Influence of Friction
Coefficient.Three numerical tests with various friction coefficients  (0.1, 0.3, and 0.5) are performed to study the impact of friction coefficient on stressstrain curve, as shown in Figure 18.The friction coefficient does not significantly affect the initial tangential modulus.It can be seen that all the samples present obvious strain softening and the peak stress increases from 6.2 MPa to 6.9 MPa as the friction coefficient increases from 0.1 to 0.5.The results may be due to the fact that more energy is needed to overcome the slip between particles when the friction coefficient becomes larger and therefore greater stress is needed in order to assure a certain loading rate.

Energy Evolution of Frozen Soil.
The energies considered in this study are boundary, bonding, frictional, rolling, strain, and kinetic energies.The boundary energy is the total accumulated work done by all boundaries on the assembly and can be defined by where   is the boundary energy,   is the number of boundaries,   and  3 are the resultant force and moment of boundaries, and Δ  and Δ 3 are the displacement and rotation of boundaries.The strain energy is the total energy stored at all contacts and can be defined by where   is the strain energy,   is the number of contacts, |   | and |   | are the magnitudes of the normal and tangential components of the contact force,   and   are the normal and tangential contact stiffness.
The frictional energy is the total energy dissipated by sliding and can be defined by   where   is the frictional energy, ⟨   ⟩ is the average shear force, and (Δ   ) slip is the increment of slip displacement.The rolling energy is the total energy dissipated by rolling and can be defined by where   is the rolling energy, ⟨  ⟩ is the average moment, and (Δ  ) roll is the increment of rotation.
The kinetic energy is total kinetic energy of all particles, accounting for both translational and rotational motion and can be defined by where   is the kinetic energy,   is the number of particles,   is the mass of particle,   is the inertia moment of particle, u   is the translational velocity of particle, and θ   is the rotational velocity of particle.
The bonding energy   is the total energy dissipated by breakage of bonds and can be defined by (18) based on the law of conservation of energy: 3.3.1.Influence of Contact Stiffness Ratio.The impact of contact stiffness ratio on energy evolution is studied by three numerical tests (1, 1.3, and 1.5), as shown in Figure 19.At the initial state, the boundary energy increases slightly with the increasing contact stiffness ratio and then the increasing rate of boundary energy becomes quite large, as shown in Figure 19(a).The bonding energy is hardly influenced by contact stiffness ratio, as shown in Figure 19(b).At the end of test, the bonding energy decreases from 6311 J to 4522 J which means a 28 percent decline as the contact stiffness ratio decreases from 1.5 to 1. Figure 19(c) demonstrates that the relationship between contact stiffness ratio and strain energy approximately can be divided into two phases.In the first phase, when the strain range is from 0 to 3%, and the strain energy is slightly influenced by contact stiffness ratio.
In the second phase, when the strain range is from 3% to 15%, the strain energy decreases with the increasing contact stiffness ratio.As the contact stiffness ratio increases, more soil particles can slip and roll due to the increasing bonding energy.As a result, the frictional and rolling energies increase with the increase of contact stiffness ratio, as shown in Figures 19(d) and 19(e).Figure 19(f) indicates that the kinetic energy is extraordinarily low and can be ignored.

Influence of Friction
Coefficient.The influence of friction coefficient on energy evolution is studied by three numerical tests (0.1, 0.3, and 0.5), as shown in Figure 20.The boundary energies of samples with different friction coefficients tend to increase with the increasing strain, as shown in Figure 20(a).Moreover, the boundary energy increases with the increase of friction coefficient.As shown in Figure 20(b), the bonding energy increases with the increasing friction coefficient.At the end of test, the bonding energy increases from 2333 J to 2839 J which means a 21 percent growth as the friction coefficient increases from 0.1 to 0.5.The shape of strain energy in Figure 20(c) is similar to that of stressstrain curve in Figure 16.Compared with frictional energy, rolling energy is relatively low, as shown in Figures 20(d) and 20(e).Moreover, the frictional energy does not obviously change when friction coefficient is greater than 0.3.The kinetic energy is extraordinarily low due to the small value of translational and rotational velocities of soil particle, as shown in Figure 20(f).

Conclusions
The nonlinear microcontact model incorporating rolling resistance for frozen soil is proposed.Then, the failure mechanism of frozen soil is explicated based on the evolution of contact force chains and propagation of microcracks.Finally, the effects of contact stiffness ratio and friction coefficient on stress-strain curve and energy evolution are evaluated.The salient findings are summarized as follows: (1) Although the nonlinear and linear models can better describe the laboratory data before the deviatoric stress reaches its peak value, only the simulated results of nonlinear microcontact model show a good agreement with the test data after the deviatoric stress reaches its peak value.(2) The evolution of contact force chains and propagation of microcracks can well account for the failure mechanism of frozen soil since the evolution of contact force chains can illustrate how the applied loads are transmitted within frozen soil and the formation of shear bands in frozen soil is associated with the propagation of microcracks.(3) The stiffness of sample increases with the increase of contact stiffness ratio.The friction coefficient does not significantly affect the initial tangential modulus.
In addition, almost all energies increase with the increasing contact stiffness ratio and friction coefficient.Compared with frictional coefficient, the effect of contact stiffness ratio on stress-strain curve and energy evolution is greater.

Figure 2 :
Figure 2: Mechanical response of normal contact model.

Figure 5 :Figure 6 :
Figure 5: Contact forces at the left-hand side and right-hand side before bond breaks.

Figure 10 :
Figure 10: Grain size distribution curve of Fujian standard sand.

Figure 11 :
Figure 11: The static triaxial material testing machine.

Figure 13 :
Figure 13: Relationship between normal bonding strength and strain.

Figure 17 :
Figure 17: Relationship between stress-strain curve and contact stiffness ratio.

Figure 18 :
Figure 18: Relationship between stress-strain curve and friction coefficient.

Table 1 :
Other physical indices of Fujian standard sand.

Table 2 :
The triaxial compression test conditions of specimen DTJwd2.

Table 3 :
Other microparameters of nonlinear and linear models.