Geometrical network of granular materials under isochoric cyclic shearing

We use three-dimensional particle dynamics simulations to investigate the microstructure evolution of granular material subjected to isochoric cyclic shearing, driving the system to a liquefaction state. The cyclically sheared assembly presents a realistic macroscopic response as observed in physical experiments. By analyzing the contact network evolution in the post-liquefaction period, we show that the onset of liquefaction state is characterized by a sudden drop of coordination number and a fragile particle connectivity network. The simulation suggests a critical coordination number for exiting the liquefaction state. Evolution of fabric anisotropy combined with coordination number implies the isotropic and anisotropic gain or loss of contacts at certain durations of a post-liquefaction loading cycle.


Introduction
The macroscopic behavior of liquid-saturated cohesionless granular materials under dynamic loading has been explored and modeled in the realm of soil mechanics in recent decades [1][2][3][4]. The gradual evolution of granular materials under cyclic shearing towards liquefaction, characterized by the vanishing of the mean effective stress in the solid phase, is termed "cyclic liquefaction". While soil liquefaction arises from load transfer between soil and water, this phenomenon does not need a saturated granular material or dynamic conditions, but can occur in a dry granular material under isochoric (constant volume) and quasistatic loading conditions. For this reason, isochoric cyclic shearing, without explicitly incorporating pore fluid, has often been used as a proxy for the undrained cyclic deformation of liquid-saturated granular materials subjected to isobaric (constant total stress) conditions [5][6][7][8][9].
In most related discrete element modeling (DEM) studies, the exploration of microscopic features has been limited to a short-term liquefaction process, not covering enough post-liquefaction period. The physical mechanisms governing the long-time evolution of the microstructure that underlies transition to the liquefaction state remain to a large extent unexplored. In particular, the contact network anisotropy may increase dramatically in liquefaction state, reflecting the lower values of the coordination number and unstable particle arrangements. While some studies show that the normal force anisotropy prevails in the liquefaction state as compared to the fabric and tangential force anisotropies [10,11], these features have never been analyzed during shear cycles before and after liquefaction and discussed in connection with the coordination number.
In this paper, we analyze the long-time evolution of the contact network using an extensive long-time DEM simulation of the 3D packing composed of spherical particles under isochoric cyclic simple shearing. The simulation was carried out with many time steps for adequately covering the whole process before and after initial liquefaction. We first describe the numerical procedures. Then, we present the macroscopic results, followed by the microscopic investigation of the geometrical network.

Numerical procedure
A three-dimensional (3D) DEM program, named GR-Flow3D [12], was used in this work. The granular assembly was simulated using spheres interacting via softparticle laws [13,14], including normal collision, tangential sliding, rolling, and torsion. The essential quantity is the elastic deflection between particles, from which the corresponding force can be evaluated using a linear springdashpot model. The simulation involves two steps: preparing a particle assembly via isotropic compression condition and applying cyclic simple shear mode to the assembly under isochoric conditions.
The sample consists of spheres with low polydispersity, i.e., d max /d min = 2 where d min = 1 mm and d max refer to the minimum and maximum particle diameters, respectively. Between d min and d max , the particle size follows a uniform distribution of particle volumes, so that the number of particles belonging to a class of diameter d is proportional to d −3 [15]. Once the particles are generated, they are placed randomly on a 3D sparse lattice to avoid overlap. This 3D lattice is contained in a rectangular biperiodic cell whose top and bottom sides are rigid walls, and four lateral sides are periodic boundaries. The samples are compressed isotropically by translationally moving the six sides of the cell and setting gravity to zero. The tangential friction coefficient µ t is tuned to achieve a certain void ratio e, defined as the ratio of the total pore volume to the solid volume. We adopted a simple computational procedure modified from [16] to prepare samples comparable with the laboratory ones. One can refer to [17] for the details. Fig. 1(a) displays an isotropically compressed sample of 8000 spheres with void ratio e = 0.647 (medium dense) under a mean pressure of 100 kPa.
In cyclic shearing, the sample volume is kept constant by fixing four lateral sides and the bottom wall and keeping the sample height constant. The shearing is undertaken by moving the top wall horizontally at a constant velocity of v x . To reduce possible slippage between the walls and the sample, a layer of particles is glued to the top and bottom walls, respectively, as indicated by gray spheres in Fig. 1(b). The shear direction is reversed each time the shear stress τ extracted from the calculated stress tensor, as explained below, reaches the target amplitude τ amp . We adopted the inertial number I =γd ρ/p to maintain a quasistatic loading, whereγ = |v x |/h is the shear rate with h the sample height, ρ the density of particles, andd the mean particle diameter. The shear is nearly quasistatic if I < 10 −3 [18]. We tried v x between 0.005 m/s and 0.01 m/s and did not see noticeable change in the macroscopic response even in the liquefaction state. In this study, v x = 0.01 m/s is used, corresponding to a shear ratė γ ≈ 0.38 s −1 and consistent with [9]. This choice makes a faster simulation and also shearing is nearly quasistatic in the non-liquefaction states. Even in the liquefaction state, I does not increase beyond 0.025.
The simulation parameters are given in Table 1. The rolling and torsion stiffnesses and friction coefficients were set to a small nonzero value in order to make rotations slightly dissipative as a simple way to account for the effects due to aspherical particle shape [14].

Macroscopic response
The stress tensor σ of the granular assembly is determined via the microscopic interactions between particles over a selected volume V: where l c is the branch vector connecting the centres of two particles for interior contact or connecting the particle centre and the contact point for exterior contacts, f c is the contact force, ⊗ denotes the dyatic tensor product, and the summation runs over all the contacts N c in the selected volume V. The selected volume shares the same center and occupies 80% of the whole sample, excluding the boundary layer effect from the top and bottom. The shear stress τ and mean stress p is obtained from σ, i.e., τ = σ zx and p = (σ xx + σ yy + σ zz )/3. The variation of p compared with the initial mean stress p 0 , is denoted by "excess pore pressure" ∆u as ∆u = p 0 − p. It can be normalized by p 0 , thus introducing excess pore pressure ratio r u , i.e., r u = ∆u/p 0 = 1 − p/p 0 . The shear strain γ is defined by x w /h, where x w is the cumulative horizontal displacement of the top wall. Fig. 2 presents the typical macroscopic behavior of the constant volume cyclic simple shear test with τ amp = 25 kPa, described in terms of stress path and stress-strain curve. The stress path in Fig. 2(a) moves cyclically leftward from p = 100 kPa and τ = 0 kPa, with a decreasing p (increasing r u ) due to the system contraction  tendency. The first time that r u reaches 0.99 is termed "initial liquefaction", and its corresponding number of cycles is denoted as N IL . the cyclic shearing process is divided into two periods, before and after initial liquefaction, namely pre-and post-liquefaction periods, as colored in gray and red in Fig. 2, respectively. The shear strain in pre-liquefaction period is negligibly small but develops significantly in post-liquefaction period, especially when τ vanishes and r u approaches 1.0. Hereafter, we assume that the system gets into liquefaction state when r u > 0.99 and it exits liquefaction state when r u < 0.99.
A post-liquefaction cycle C is highlighted in Fig. 2, with six characteristic states. C 0 is the first time τ ≥ 0 distinguishing loading from unloading, C 0 is the exit from liquefaction state, and C 1 is τ reaching τ amp . C 2 , C 2 , and C 3 are similar to C 0 , C 0 , and C 1 , respectively, in terms of shearing along the negative x direction (see Fig. 1(b)).

Granular microstructure
In this section, the granular microstructure is investigated in terms of particle connectivity and fabric anisotropy. The coordination number z g , considering non-floaters (particles with contacts) and force-bearing contacts, is adopted: where N p is the total number of particles, N 0 p is the number of floaters, and N c is the number of contacts. Fig. 3 displays the evolution of z g with the fractional number of cycles N (another way representing time), where the time history is colored by the value of r u . The number of cycles reaching initial liquefaction N IL is marked by the vertical dashed line. We see that z g decreases from its initial value z g 4.76 with small oscillations in pre-liquefaction period and drops below 4.0 when the system approaches initial liquefaction. In the post-liquefaction period, z g stays below 4.0 and fluctuates significantly down to values as low as 1.5 with a negative static redundancy, implying insufficient constraints to hold the system stable. A zoomed-up window is added to present details of z g evolution in cycle C. The horizontal dashed line for z g = 3.6 corresponds to the inflection point of z g as a function of N, coinciding with the exit of liquefaction state (C 0 and C 2 ). It may be considered the percolation threshold of the particles allowing for force transmission across the system through the contact network.
The connectivity of particles P c is defined as the proportion of particles with exactly c contacts. Its distribution at the characteristic states of cycle C shown in Fig. 4 provides more detailed information about the microscopic state than z g . C 0 and C 2 exhibit a high proportion of particles with c < 4, implying a fragile contact network. This fragile network disappears only when the system exits the liquefaction states as shown by {P c } at C 0 and C 2 . Given large shear strain development between C 0 and C 0 or C 2 and C 2 , one can infer that sample deformation rebuilds the fragile network resulting from unloading (compare C 1 and C 2 ) although p does not increase markedly. The system stays stable at C 1 and C 3 while C 0 and C 2 represent the states with the weakest contact network in the postliquefaction cycle.
We use a scalar a c to quantify the contact network anisotropy, defined as where S c = a c : s/( √ a c : a c √ s : s) is the normalized first joint invariant between the fabric anisotropy tensor a c , as   Figure 4: The connectivity P c of the contact network in the post-liquefaction cycle C.
defined in [19], and the deviatoric stress tensor s. The sign of scalar a c determines the coaxiality between the deviatoric stress tensor and the fabric tensor. Fig. 5 presents a c versus z g during the cyclic shearing process, a picture of the reorganizations of the contact network in response to external loading [20], with cycle C being highlighted. The diagram starts from the rightmost part with a high value of z g and a c 0, and follows a path towards to the left with oscillations between a c 0 and a maximum value of a c that increases gradually with N. The non-liquefaction states belong to the region of z g ≥ 3.6. In liquefaction states, both z g and a c vary significantly and follow long paths exemplified by that from C 0 to C 0 , corresponding to anisotropic gain of contacts. When the system exits liquefaction state and reaches τ amp (between C 0 and C 1 ), or unloading back to liquefaction (between C 1 and C 2 ) while z g either increases or declines, a c is nearly constant, implying isotropic gain or loss of contacts.

Conclusions
This paper investigated the evolution of geometrical network during isochoric cyclic simple shearing using discrete-element numerical simulations. The macroscopic behavior is characterized by gradual reduction of mean stress until the system enters a liquefaction state. In the transition to liquefaction state, the coordination number drops significantly due to isotropic loss of contacts, and the force-bearing network collapses. Large deformation in the liquefied states leads to rebuilding the contact network and exiting the liquefaction state at the coordination number of 3.6 corresponding to the geometrical percolation threshold of the contact network.