Applied Science, Faculty of
Civil Engineering, Department of
DSpace
UBCV
Banerjee, Sounik Kumar
Doctor of Philosophy - PhD
University of British Columbia
Cyclic liquefaction of granular soils is influenced by particle-level parameters such as particle size distribution and particle shape. These parameters are difficult to isolate in real soils, whereas use of Discrete element method (DEM) can be an important step to resolve this problem. It can also be used to model these individual effects on the undrained cyclic response. Samples with increasing particle size distributions, denoted by an increase in coefficient of uniformity are prepared under isotropic compression at two different relative densities. Constant volume, cyclic simple shear simulations at a low relative density under various loading intensities show that increasing the coefficient of uniformity initially increases, then decreases, the cyclic liquefaction resistance. This observation is largely reversed when the relative density is increased, in which case the liquefaction resistance decreases initially and then increases as the coefficient of uniformity increases. Microscopic examinations of the samples at their initial state are conducted to provide insight into their liquefaction resistance. Microparameters based on inter-particle contacts and normal forces reveal reasonable correlations with the macroscopic response. Furthermore, the initial state-based state parameter is found to significantly correlate with the cyclic resistance ratio irrespective of coefficient of uniformity and relative density. Finally, the determinate structure of particle assemblies (isostaticity) at the advent of cyclic liquefaction is studied in terms of average contacts. A similar approach is used to study the effect of particle shape under constant volume cyclic loading. Samples with a variety of simplified particle shapes denoted by aspect ratio and blockiness are isotropically compressed under two relative densities. Samples with non-spherical particles have lower liquefaction resistance than those with spherical particles at low relative density. In contrast, at a medium relative density, the non-spherical samples show comparatively higher liquefaction resistance. It has been discovered that these macrolevel observations under the influence of particle shape can be adequately explained by microparameters previously used to study the particle size distribution. The determinate structure is also identified here and its variation is studied for the range of particle shape and relative densities.
EFFECTS OF PARTICLE SIZE DISTRIBUTION AND PARTICLE SHAPE ON CYCLICLIQUEFACTION RESPONSE OF GRANULAR MATERIALSbySounik Kumar BanerjeeB.E., Indian Institute of Engineering Science and Technology Shibpur, 2013M.Tech., Indian Institute of Technology Kharagpur, 2015A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinTHE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES(Civil Engineering)The University of British Columbia(Vancouver)April 2022© Sounik Kumar Banerjee, 2022The following individuals certify that they have read, and recommend to the Faculty of Graduateand Postdoctoral Studies for acceptance, the thesis entitled:Effects of particle size distribution and particle shape on cyclic liquefaction response of granular materialssubmitted by Sounik Kumar Banerjee in partial fulfillment of the requirements forthe degree of Doctor of Philosophyin Civil EngineeringExamining Committee:Dr. Mahdi Taiebat, Civil Engineering, The University of British ColumbiaSupervisorDr. Anthony Wachs, Mathematics and Chemical and Biological Engineering, The University ofBritish ColumbiaSupervisory Committee MemberDr. John Howie, Civil Engineering, The University of British ColumbiaSupervisory Committee MemberDr. Arcesio Lizcano, Industry Expert, SRK ConsultingSupervisory Committee MemberDr. Davide Elmo, Mining Engineering, The University of British ColumbiaUniversity ExaminerDr. Carlos Ventura, Civil Engineering, The University of British ColumbiaUniversity ExaminerDr. Jidong Zhao, Civil and Environmental Engineering, Hong Kong University of Science andTechnologyExternal ExamineriiAbstractCyclic liquefaction of granular soils is influenced by particle-level parameters such as particle sizedistribution and particle shape. These parameters are difficult to isolate in real soils, whereas useof Discrete element method (DEM) can be an important step to resolve this problem. It can alsobe used to model these individual effects on the undrained cyclic response.Samples with increasing particle size distributions, denoted by an increase in coefficient ofuniformity are prepared under isotropic compression at two different relative densities. Constantvolume, cyclic simple shear simulations at a low relative density under various loading intensitiesshow that increasing the coefficient of uniformity initially increases, then decreases, the cyclicliquefaction resistance. This observation is largely reversed when the relative density is increased,in which case the liquefaction resistance decreases initially and then increases as the coefficient ofuniformity increases. Microscopic examinations of the samples at their initial state are conductedto provide insight into their liquefaction resistance. Microparameters based on inter-particlecontacts and normal forces reveal reasonable correlations with the macroscopic response.Furthermore, the initial state based state parameter is found to significantly correlate with thecyclic resistance ratio irrespective of coefficient of uniformity and relative density. Finally, thedeterminate structure of particle assemblies (isostaticity) at the advent of cyclic liquefaction isstudied in terms of average contacts.A similar approach is used to study the effect of particle shape under constant volume cyclicloading. Samples with a variety of simplified particle shapes denoted by aspect ratio andblockiness are isotropically compressed under two relative densities. At low relative densities,samples containing non-spherical particles exhibit less liquefaction resistance than samplesiiicontaining spherical particles. In contrast, at a medium relative density, the non-spherical samplesshow comparatively higher liquefaction resistance. It has been discovered that these macrolevelobservations under the influence of particle shape can be adequately explained bymicroparameters previously used to study the particle size distribution. The determinate structureis also identified here and its variation is studied for the range of particle shape and relativedensities.ivLay SummaryCyclic loading from earthquakes, waves, and wind can liquefy saturated granular soils, causingmassive human and infrastructure damage. Although the particle size distribution and shape ofgranular soils can influence liquefaction, their effects are difficult to study individually due todifferent mineral composition and surface features. Discrete element modeling can solve theproblem of overlapping effects often present in laboratory experiments. Particle size distributionand particle shape effects are thus isolated, and their individual effects on the cyclic liquefactionare analyzed in the idealized granular assemblies. This detailed particle-level study providesinsight for mechanism governing the cyclic liquefaction response in sandy soils.vPrefaceI, Sounik Kumar Banerjee, am the primary and sole author of all six chapters in the dissertation.Professor Mahdi Taiebat, my research supervisor have offered guidance and direction on thedissertation.The research is conducted using discrete element simulations and the main contributions aredivided into two main parts. The first part contains the results and discussions from Chapter 4, inwhich the effect of particle size gradation on the cyclic response of isotropic granular materialsis investigated. The second, which is described in detail in Chapter 5, investigates the effect ofsynthetic particle shapes on the cyclic response. Both chapters are undergoing revisions as journalmanuscripts.viTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiList of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .xviiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .xxiii. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.3 Layout of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.4 Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10vii2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.2 Effects of particle size distribution on the monotonic response . . . . . . . . . . . 102.2.1 Experimental testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.2.2 DEM testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.3 Effects of particle size distribution on the cyclic response . . . . . . . . . . . . . . 142.3.1 Experimental testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.3.2 DEM testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.4 Effects of particle shape on the monotonic response . . . . . . . . . . . . . . . . . 172.4.1 Experimental testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.4.2 DEM testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.5 Effects of particle shape on the cyclic response . . . . . . . . . . . . . . . . . . . . 202.5.1 Experimental testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.5.2 DEM testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253.2 DEM program . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253.2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.2.2 Contact laws . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.2.3 Input file structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 283.2.4 Use of high performance computing . . . . . . . . . . . . . . . . . . . . . 283.3 Sample preparation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 293.4 Shearing process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 323.5 Description of particle size distribution . . . . . . . . . . . . . . . . . . . . . . . . 333.6 Scalability or choice of optimum particle numbers . . . . . . . . . . . . . . . . . . 343.7 Choice of optimum particle classes . . . . . . . . . . . . . . . . . . . . . . . . . . 38viii3.8 Description of particle shape . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 413.8.1 Shape parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 413.8.2 Relation with roundness and sphericity . . . . . . . . . . . . . . . . . . . . 433.9 Validation with laboratory tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . 453.10 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 484 Effect of particle size distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 494.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 494.2 Extreme void ratios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 494.3 Test design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 514.4 Macroscopic response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 524.5 Cyclic response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 554.5.1 Evolution of excess pore pressure . . . . . . . . . . . . . . . . . . . . . . 564.5.2 Cyclic liquefaction resistance . . . . . . . . . . . . . . . . . . . . . . . . . 574.6 Relation between CRR15 and initial state . . . . . . . . . . . . . . . . . . . . . . . 604.6.1 Relation with z0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 604.6.2 Relations with inter-particle force and contact distributions . . . . . . . . . 624.7 Relation between CRR15 and state parameters . . . . . . . . . . . . . . . . . . . . 694.8 Evolution of the non-load bearing particles . . . . . . . . . . . . . . . . . . . . . . 724.9 Examination of the isostatic stage . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.10 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 765 Effect of particle shape . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 785.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 785.2 Extreme void ratios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 785.3 Test design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 805.4 Macroscopic response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 815.4.1 Evolution of excess pore pressure . . . . . . . . . . . . . . . . . . . . . . 81ix5.4.2 Cyclic liquefaction resistance . . . . . . . . . . . . . . . . . . . . . . . . . 835.5 Relation between CRR15 and initial state . . . . . . . . . . . . . . . . . . . . . . . 885.5.1 Relation with z0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 885.5.2 Relation with contact normal force distributions . . . . . . . . . . . . . . . 905.6 Relation between CRR15 and state parameters . . . . . . . . . . . . . . . . . . . . 955.7 Evolution of zm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 975.8 Analysis of the isostatic stage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 985.9 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1006 Summary, Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . . . . 1026.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1026.2 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1046.2.1 Effect of particle size distribution . . . . . . . . . . . . . . . . . . . . . . 1046.2.2 Effect of particle shape . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1056.3 Recommendations for future work . . . . . . . . . . . . . . . . . . . . . . . . . . 106A LIGGGHTS CODE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133xList of TablesTable 3.1 DEM parameters involving interparticle and particle-wall interactions. . . . . . 29Table 3.2 The effects of total particles Np and random positioning seed Sn on the void ratioe, mechanical coordination number zm, and liquefaction resistance indicator NIL. 38Table 4.1 Simulated constant volume cyclic simple shear tests. . . . . . . . . . . . . . . . 52Table 5.1 Simulated constant volume cyclic simple shear tests on samples withp0 = 100 kPa. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81Table 5.2 Roundness, sphericity and overall-regularity of the samples. . . . . . . . . . . . 87xiList of FiguresFigure 1.1 Cyclic shear load results in contact loss, which eventually results in cyclicliquefaction [1]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1Figure 1.2 Effects of cyclic liquefaction: (a) landslide [2], (b) building collapse [3], (c)sinking of transportation infrastructure [4], and (d) lateral spreading [5]. . . . . 2Figure 1.3 Microphotographs of different types of particle size distributions ranging from(a) uniformly graded to (b) well graded samples [6]. . . . . . . . . . . . . . . . 4Figure 1.4 Microphotographs of various particle shapes found in real sands: (a) and (b)ASTM 20/30 sand, (c) glass beads, (d) and (e) crushed granite, and (f) ticinosand [7]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5Figure 2.1 Effect of particle size distribution on liquefaction resistance (modified from [8]). 14Figure 2.2 Effect of particle shape on liquefaction resistance (a) at σv = 100 kPa and, (b)at a significantly higher σv = 1600 kPa (modified from [9]). . . . . . . . . . . . 21Figure 2.3 Effect of simulated particle shape of real sands on liquefaction resistance for a(a) loose sample, and (b) dense sample (modified from [10]). . . . . . . . . . . 22Figure 3.1 Particle arrangements and boundary conditions for an assembly with 12600particles and sample size ≈ 20×D50: (a) at the end of isotropic compressionwith bi-periodic lateral walls, and (b) during volume cyclic shearing. The greyparticles at the top of the sample form the shearing plate, while the particles onthe bottom form the rigid plate. . . . . . . . . . . . . . . . . . . . . . . . . . . 30xiiFigure 3.2 Snapshots of the samples representing relevant particle size distributionswhere, (a) Cu = 1, (b) Cu = 1.5, (c) Cu = 2, and (d) Cu = 3. . . . . . . . . . . . 31Figure 3.3 Snapshots of samples with different particle shapes where (a) AR = 0.5, B = 2,(b) AR= 1, B = 2, (c) AR= 1.5, B= 2, (d) AR= 1, B= 4, and (e) AR= 1, B= 8. 32Figure 3.4 The effect of total number of particles (Np) in the assembly on the initial valuesof (a) e, (b) zm, and (c) NIL. The effects from random particle distributions(change in seed number) on packing properties and liquefaction resistance aredemonstrated by the error bars. . . . . . . . . . . . . . . . . . . . . . . . . . . 37Figure 3.5 Effect of (a) total particle classes (Nc), and (b) total particles in each class(Ncp), on the void ratio e of samples with Cu = 3, p0 =100 kPa. . . . . . . . . . 39Figure 3.6 Design of the log-linear particle size distribution with ten particle classes (C1-C10) for a typical case of Cu = 2 and D50 = 5 mm. . . . . . . . . . . . . . . . . 40Figure 3.7 A typical superquadratic ellipsoid with marked half lengths along major axis cand minor axes a and b. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42Figure 3.8 Superquadratic particle shapes with different aspect ratio AR and blockiness B.Highlighted yellow and green portions represent the ranges of AR and B studied. 42Figure 3.9 Relation between (a) aspect ratio AR, and (b) blockiness B, of superquadraticswith roundness and sphericity calculated from the current analysis. . . . . . . . 45Figure 3.10 The typical (a) stress path, and (b) stress-strain curve, for a sample subjectedto cyclic loading in a simple shear test (current study). . . . . . . . . . . . . . 46Figure 3.11 The typical (a) stress path, and (b) stress-strain curve, for a sample subjectedto cyclic loading in a hollow-cyclinder laboratory test [11]. . . . . . . . . . . . 46Figure 3.12 Qualitative comparisons between a typical CSR curve for Cu = 1.5 and σv =100 kPa from (a) current study using DEM, and (b) an undrained cyclic hollowcylinder test by Wichtmann [12]. . . . . . . . . . . . . . . . . . . . . . . . . . 47Figure 4.1 Particle size gradation curves with five levels of Cu used in the study. . . . . . . 50xiiiFigure 4.2 Maximum and minimum simulated void ratios from current study andcomparisons with literature. . . . . . . . . . . . . . . . . . . . . . . . . . . . 51Figure 4.3 Macroscopic shear response of a sample with Cu = 2, Dr = 50% and CSR =0.35 where: (a) stress path, (b) stress-strain curve, (c) derived pore pressureevolution, and (d) shear strain development. . . . . . . . . . . . . . . . . . . . 54Figure 4.4 Evolution of excess pore pressure ratio ru given at the beginning each cycle(τ ≃ 0) for samples with different levels of Cu at: (a) Dr = 30%, at τ ≃ 0,and (b) Dr = 50%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56Figure 4.5 Cyclic liquefaction resistance curves for samples with different values of Cuand Dr. The hollow markers represent results for Dr = 30%, and the solidmarkers for Dr = 50%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58Figure 4.6 Variation of cyclic liquefaction resistance with Cu, where (a) liquefactionresistance is represented as CRR15; and (b) as NIL. The hollow markersrepresent results for Dr = 30%, and the solid markers for Dr = 50%. . . . . . . 59Figure 4.7 Variation of z0 for samples with different Cu. Those packed at Dr = 30% areshown with hollow markers, and for Dr = 50% with solid markers. . . . . . . . 60Figure 4.8 Correlation between z0 with CRR15 for samples with different Cu and Dr. Theexponential best fit is shown as the blue dotted line. . . . . . . . . . . . . . . . 61Figure 4.9 Probability distribution functions of normalized contact normal forces fn/⟨ fn⟩for samples with different Cu at (a) low Dr with wide force distribution, and (b)medium dense Dr with narrow force distribution. . . . . . . . . . . . . . . . . 64Figure 4.10 The relation between the rate of decay β of the strong force networkfn/⟨ fn⟩ > 1 with CRR15. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64Figure 4.11 Snapshots of normal force chains of the samples at the initial state at Dr = 30%for: (a) Cu = 1, (b) Cu = 2, and (c) Cu = 3. The zoomed insets are provided forthe sake of clarity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66xivFigure 4.12 Snapshots of normal force chains of the samples at the initial state at Dr = 50%for: (a) Cu = 1, (b) Cu = 2, and (c) Cu = 3. The zoomed-in insets are includedfor clarity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67Figure 4.13 The fraction Pc of particles with c contacts for different values of Cu at (a)low Dr with significant presence of non-load bearing contacts in higher Cu (b)medium-dense Dr with suppressed non-load bearing contacts. . . . . . . . . . . 68Figure 4.14 Drained monotonic simple shear test on samples with different Cu at Dr =30% and 50%: (a) stress-strain response, (b) evolution of void ratio, and (c)evolution of mechanical coordination number. . . . . . . . . . . . . . . . . . . 70Figure 4.15 Cyclic liquefaction resistance against state parameters for samples withdifferent particle size gradations and relative density derived from: (a) voidratio, and (b) mechanical coordination number. The hollow markers representDr = 30%, and solid markers represent Dr = 50%. . . . . . . . . . . . . . . . . 71Figure 4.16 Evolution of ffr with N/NIL at the end of each loading cycles (τ ≃ 0) atdifferent levels of Dr. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73Figure 4.17 The evaluation of ziso at initial liquefaction based on the loss of σv: (a) Thedegradation of zm with loss of σv leading to the vertical liquefied zone at σv ≃0, and (b) comparisons of N at liquefaction corresponding to ziso and ru = 0.95. 74Figure 4.18 The effects of particle size distribution and relative density on ziso. . . . . . . . 75Figure 5.1 The maximum and minimum simulated void ratios for samples with (a)different aspect ratios AR and same blockiness B, and (b) different B and sameAR. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79Figure 5.2 Log-linear particle size distribution of samples with Cu = 2 for all values of ARand B considered in this study. . . . . . . . . . . . . . . . . . . . . . . . . . . 80xvFigure 5.3 The effect of particle shape on the evolution of accumulated ru at CSR = 0.25and τ ≈ 0 for: (a) AR = 0.5 to 1.5 and Dr = 20% (b) AR = 0.5 to 1.5 and Dr =50%, (c) B = 2 to 8 and Dr = 20%, and (d) B = 2 to 8 and Dr = 50%. . . . . . . 82Figure 5.4 Cyclic liquefaction strength curves of loose and medium-dense samples withdifferent values of AR, when B is kept constant at 2. . . . . . . . . . . . . . . . 84Figure 5.5 Cyclic liquefaction strength curves of loose and medium-dense samples withdifferent values of B when AR is kept constant at 1. . . . . . . . . . . . . . . . 85Figure 5.6 Variation of liquefaction resistance in terms of (a) interpolated values of CRR15with AR, (b) NIL with AR, (c) interpolated values of CRR15 with B, and (d) NILwith B. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86Figure 5.7 Variation of cyclic liquefaction resistance indicators (a) CRR15, and (b) NILwith overall regularity (OR). . . . . . . . . . . . . . . . . . . . . . . . . . . . 88Figure 5.8 The relation between z0 with CRR15 for samples with different B and AR. . . . 89Figure 5.9 Probability distribution functions of normalized contact normal forces fn/⟨ fn⟩for samples with different B and AR at (a) AR = 0.5, 1 and 1.5, Dr = 20%, (b)AR = 0.5, 1 and 1.5, Dr = 50%, (c) B = 2, 4 and 8, Dr = 20%, and (d) B = 2, 4and 8, Dr = 50%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90Figure 5.10 The relation between the rate of decay β of the normalized strong contactnormal force network fn/⟨ fn⟩ > 1 for different AR and B with CRR15. . . . . 91Figure 5.11 Snapshots of normal force chains of the samples with different aspect ratios atthe initial state for: (a) AR = 0.5, Dr = 20%, (b) AR = 1, Dr = 20%, (c) AR =1.5, Dr = 20%, (a) AR = 0.5, Dr = 50%, (e) AR = 1, Dr = 50%, and (f) AR =1.5, Dr = 50%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92Figure 5.12 Snapshots of normal force chains of the samples with different blockiness atthe initial state for: (a) B = 2, Dr = 20%, (b) B = 4, Dr = 20%, (c) B = 8, Dr =20%, (d) B = 2, Dr = 50%, (e) B = 4, Dr = 50%, and (f) B = 8, Dr = 50%. . . . . 93xviFigure 5.13 Cyclic liquefaction resistance indicator CRR15 against state parameters forsamples with different particle shapes and relative density derived from: (a)void ratio, and (b) mechanical coordination number. . . . . . . . . . . . . . . . 96Figure 5.14 The effect of particle shape on the evolution of zm for (a) AR = 0.5, 1, 1.5 at Dr= 20%, (b) AR = 0.5, 1, 1.5 at Dr = 50%, (c) B = 2, 4, 8 at Dr = 20%, and (d) B= 2, 4, 8 at Dr = 50%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98Figure 5.15 The value of ziso is described with (a) the corner in the zm curve with an arrow,and (b) position of ziso curve with the corresponding value of N. . . . . . . . . 99Figure 5.16 Variation of ziso at different relative densities with (a) different AR at B = 2,and (b) different B at AR = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 99xviiList of Symbols⊗ Dyadic product⟨⟩ Absolute value or norm of a vectortr(·) Trace operatorai Radius of particle iAR Aspect ratioB BlockinessCc Coefficient of curvatureCu Coefficient of uniformityCSR Cyclic stress ratioCRR15 Cyclic resistance ratio corresponding to 15 loading cyclesd Particle diameterdmax Maximum particle diameterdmin Minimum particle diameterD10 Particles finer than this size are 10% of the total assembly volumeD50 Mean particle diameter of an assemblyxviiiD60 Particles finer than this size are 60% of the total assembly volumeDr Relative densitye Void ratioe0 Void ratio at the initial stateec Void ratio at the critical stateemax Maximum void ratioemin Minimum void ratiofn,Fn Normal contact forcefr Fraction of floaters and rattlers in an assemblyFV Percentage finer by particles volumeFnel Elastic part of the contact normal forceFndis Dissipative part of the contact normal forceF t Total contact tangential forceF tel Elastic part of the contact tangential forceF tdis Dissipative part of the contact tangential forceI Inertial numberkn Spring normal stiffnesskt Spring tangential stiffnessK0 Earth pressure coefficient at restxixME Margin of errorM0el Mechanical spring component of the inter-particle torqueM0dis Mechanical dissipative component of the inter-particle torquen Sample size based on seedsN Number of loading cyclesNIL Number of loading cycles to initial liquefactionNc Total number of contactsNp Total number of particlesN0p Number of particles without contactsN1p Number of particles with one contactOR Overall regularityp Mean effective stressp0, σvo Initial mean effective stressPc Fraction of particles with c contactsP( fn) Probability distribution of contact normal forcesru Derived excess pore pressure ratiorc Position of contact cr p Center position of particle prc,p Vector connecting rc and r pxxR RoundnessReff Effective radius of the particles in contactS Sphericityzm Mechanical coordination numberz0 Mechanical coordination number at the initial statezc Mechanical coordination number at the critical statezcr Mechanical coordination number at the isostatic statezα Quantile corresponding to a confidence level of 95%β Slope of the strong contact normal force distributionγ Shear strain along y axisγamp Shear strain along y axis at τamp in a loading cycleδ[ ] Overlap between two particles∆u Pore pressureε Strain tensorµt Tangential friction coefficientµr Rolling friction coefficientρ Particle densityσ Standard deviationσxx x-component of stress tensorxxiσyy y-component of stress tensorσzz z-component of stress tensorσyz yz-component of stress tensorσxz xz-component of stress tensorσxy xy-component of stress tensorσv Effective vertical compressive stressσ Stress tensorτ Shear stressτamp Maximum value of τ in a loading cycleτx Shear stress along x axisτy Shear stress along y axisτu,lim Shear stress corresponding to ru,limψ State parameterψe State parameter based on void ratioψz State parameter based on mechanical coordination numberxxiiAcknowledgementsI, Sounik Kumar Banerjee, would like to express my heartfelt appreciation to my researchsupervisor, Professor Mahdi Taiebat, for his consistent guidance and encouragement throughoutmy doctoral research journey. In short, I have convincingly felt that his approach to supervision isto concentrate on the research fundamentals and to be rightfully suspicious of any data wegenerate. Once honed, the style of remaining neutral to both sides of the coin can direct anindividual to concentrate on any given problem with a moderate time investment. Thisunderstanding was strengthened in my mind and actions as a result of the collaborative assistanceby Dr. Ming Yang during the brief but important final period of my Ph.D. I sincerely thank bothof them and hope that this training will enable me to conduct deep and important research withhigh-quality publications in a variety of fields in geomechanics.Special thank to my supervisory committee members Dr. Arcesio Lizcano, Professor AnthonyWachs and Professor John Howie for their support during the several stages of my journey. Iwould also like to express my appreciation to the members of my research group, namely Andre´sReyes Pererra, Dr. Andres Barrero, Eric Lino, Keith Perez, and Sheng Zheng, for their constructivecomments during individual discussions and/or group meetings.Additionally, I would like to appreciate my mother Runa Banerjee and my grandmotherKalyani Mukherjee for their unwavering support in my life. Besides, I am indebted to my wifeAsmita Chakraborty for being an extremely reassuring partner and true companion.Finally, I would like to express my gratitude for the primary financial assistance provided bythe Four Year Doctoral Fellowship from University of British Columbia (UBC) and the associatedtop up grant by Natural Sciences and Engineering Research Council of Canada (NSERC).xxiiiDedicated to my motherChapter 1: Introduction1.1 BackgroundGranular soils tend to change volume in response to any type of shearing. Cyclic shearing arisingfrom repetitive seismic loadings on saturated granular soils may not allow adequate time for porefluid drainage. This may lead to the development of excess pore water pressure for both loose anddense soils. The pressure in the pore fluid leads to reduction and, in some cases, complete lossof the inter-particle force transmission as shown in Figure 1.1, thus leading to partial or completeloss of shear resistance of the material. This phenomenon of loss in shearing resistance is largelyknown as cyclic liquefaction.Figure 1.1: Cyclic shear load results in contact loss, which eventually results in cyclic liquefaction [1].Cyclic liquefaction is noted in a variety of earthquakes across the world. The 1964 Niigataearthquake, the 1964 Alaskan earthquake, the 1989 Loma Prieta earthquake, the 1994 Northridge1earthquake, 2010 El-mayor Cucapeh earthquake, and the 2011 Christchurch earthquake, to namea few, are among the most devastating earthquakes that caused catastrophic damage to people andinfrastructure due to soil liquefaction. The total of 10,000 landslides, one of which is shown inFigure 1.2(a), triggered by the Northridge earthquake, caused a total $20 billion in propertydamage. The liquefaction may result in foundation failures in tall buildings without anysuperstructure damage as evidenced in Figure 1.2(b) or sinking of transport vehicles as observedafter the Christchurch earthquake (refer to Figure 1.2(c)). Liquefaction in the inner layers mayalso result in lateral spreading of the top soil as evidenced in the El-mayor Cucapeh earthquakeand shown in Figure 1.2(d).(a) (b)(c) (d)Figure 1.2: Effects of cyclic liquefaction: (a) landslide [2], (b) building collapse [3], (c) sinking oftransportation infrastructure [4], and (d) lateral spreading [5].Cyclic liquefaction can be influenced by several key factors: (a) the loading type, (b) statevariables [13, 14], and (c) intrinsic factors [15, 16]. Undrained cyclic biaxial test (CBT) [17],Undrained cyclic triaxial test (TXT) [18], Undrained cyclic direct simple shear test (DSS) [19], andUndrained cyclic hollow cylinder test (HCT) [20] are common laboratory based element tests used2to investigate cyclic liquefaction of cohesionless soils. In these standardized tests, a representativevolume element (RVE) of soil is prepared and then sheared under truly undrained or constant-volume conditions often under a simple harmonic shearing to approximate the irregular periodicloading caused by an earthquake. In general, each of the laboratory tests mentioned above involvestwo distinct types of cyclic shearing, depending on the load increment condition: stress control andstrain control tests. The terminology used to briefly discuss them here is derived from a DSS test,as its loading mode and the rotation of principal stresses are very similar to that of an earthquake[21] and will be the current mode of operation. In a stress control test, usually the shear stressincrement is uniform over time and for a cyclic load application, the loading direction is reversedwhen the shear stress reaches a predetermined level.Whereas in the case of strain-limit test, the shear strain is incrementally applied to a soil sample,and reversed when either a predefined shear stress or strain limit is reached. The stress controlledtests generally allow for higher instability in the soil sample during cyclic degradation since thereis no control over the rate of shearing. The strain control tests with strain limit [22] though allowstable increment by controlling the shear velocity, however, it prevents shear strain from exceedinga predefined limit even after shear stress degradation. Whereas strain controlled tests with a shearstress limit are more interesting because they not only allow for strain control even in the liquefiedphase [23], but also permit modeling of large deformations.Regardless of loading technique, a granular soil sample will initially contract, followed bydilation-contraction loops until it approaches initial liquefaction. Initial liquefaction may be mainlydefined when the mean effective stress is near zero or when shear strain development is significant,apart from various other definitions as discussed by many researchers [24, 25] .The state variables of a soil sample mainly include, sample density, confining pressure,isotropic or anisotropic (fabric) condition, and over-consolidation ratio. These state variables canhave an effect on the cyclic response of the system. For example, the increase in sample densityalone, mainly measured in terms of decrease in void ratio (e) or increase in relative density (Dr)results in an increase in liquefaction resistance and vice-versa. Meanwhile, other state variables3like confining pressure can also significantly impact the liquefaction resistance [26, 27].Similarly, changes in fabric properties alone can largely affect the cyclic liquefaction [28].Therefore, in order to investigate a specific state variable, others should be comparable acrosssamples in order to offset their effects.It is also necessary to determine the unique, intrinsic (particle-level) factors in sands [14]. Thefactors may include particle size distribution, particle shape, mineralogy, inter-particle friction,among others. The particle size distribution of soil is generally classified based on the meanparticle size D50, the content of fine particles FC, and the coefficient of uniformity Cu, amongother indicators. The fine content refers to particles finer than 0.075 mm [29]. The value of Cu isdefined as D60/D10, where the D60 and D10 refer to the maximum diameter of the particles that is60 % and 10% of particles by mass of the total sample, respectively. There are mainly three typesof size distributions found in soils, as shown in Figure 1.3. A well graded soil contains particleswith a smooth distribution of all representative particle sizes, and can be denoted by Cu greaterthan one. A uniformly graded soil has soil grains of similar sizes with a Cu equal to unity. Gapgraded soils on the other hand can have an excess or a deficiency of particle sizes, or they can bedevoid of specific grain sizes.(a) (b)Figure 1.3: Microphotographs of different types of particle size distributions ranging from (a) uniformlygraded to (b) well graded samples [6].Vaid et al. [30] showed that particle size distribution can significantly influence cyclicliquefaction of granular soils. Further studies in this regard have been conducted by someresearchers [8, 31]. The main challenge has been to isolate the effect of Cu from other factors4(a) (b) (c)(d) (e) (f)Figure 1.4: Microphotographs of various particle shapes found in real sands: (a) and (b) ASTM 20/30sand, (c) glass beads, (d) and (e) crushed granite, and (f) ticino sand [7].such as particle shape, particle size, and mineralogy. To avoid the overlapping effects of the statevariables on particle shape and particle size distribution, samples must be prepared at the samerelative density and using the same sample preparation method. Adhering to such a level ofprecision across multiple parameters is challenging in experimental studies [32].Particle shape is an intrinsic property of sand that can also influence its mechanical behavior,including cyclic liquefaction [33]. Although the particle shapes are random and difficult to quantify[34–36], there are two widely used indicators in the literature to describe particle shapes, namelyroundness and sphericity [35, 37, 38]. They will be described qualitatively here in the form of achart, and in greater detail in Chapter 3, Methodology 3.8.1. Visual inspection in Figures 1.4(a)to 1.4(a) infer that the sand particles do not have sharp edges or corners on the surface, whichindicate that the particle are round. However, the particle in 1.4(a) in terms of its geometry possesselongated shape compared to the perfect spherical shape in 1.4(a), which indicate that the earlierhas less sphericity. Similar visual examination for the set of particles in the lower row i.e. Figures1.4(d) to 1.4(f) have distinct sharp edges. Among the several billions of particles present even5in an element-level soil sample, each has its unique size, and shape that is not perfectly roundedor spherical. Apart from the fact that quantifying the size of individual grains complicates sizedistribution analysis, it is also extremely difficult to evaluate the shape of enumerable grains usingvisual charts or even with statistical interventions with a relatively small sample size. As a result ofthis problem, separate effects from these intrinsic factors are often not taken into account in mostgeotechnical practice. However, there are several particle image or scanning based techniques[39], complex fourier analysis [40] and others to scan the image of particle shape. Nevertheless,the measurements remain tedious and are viable for a small number of particles. Apart frommeasuring particle shape, focusing on this parameter in laboratory experiments by excluding otherintrinsic and state variables is a significant challenge [32].The discrete element method (DEM) [41] is capable of numerically assembling particles witha desired grain size distribution and/or shape while isolating other variables as necessary [32, 42].Though laboratory tests deal with samples in lab which are mostly reconstituted, someexperimentalists may argue that DEM simplifies real-world complexity because it does not usereal soil particles. However, eliminating randomness via simulation enables a researcher to studyeach of these parameters independently, which may be a significant advantage. Furthermore,DEM enables the study of micromechanics by tracking individual particles and their interactionsin a simulation, which is not possible in a laboratory test. DEM has been used to capture theeffect of particle shape on the monotonic shear strength [43–45] and particle interactions [46–48]of granular assemblies. It uses various techniques to model particle shape in three dimensions(3D), such as spherical clusters, polyhedrons, and superquadratic particles, among others. It canalso be used to model undrained response in granular materials as a constant-volume cyclicresponse that leads to cyclic liquefaction [49, 50]. Although, to the best of knowledge, noprevious research has examined the cyclic liquefaction resistance considering the effects fromparticle size distribution and non-spherical particles using 3D-DEM.Therefore it is interesting and novel to study the effects of both particle size distribution andparticle shape isolated from other factors as is possible using DEM discussed before in a constant-6volume, cyclic direct simple shear framework. The framework will look at the cyclic liquefactionresistance and its possible connection with certain macro level and micro level state variablesmeasured both at the pres-shearing stage and at the critical state [14].1.2 ObjectivesThe purpose of this dissertation is to investigate the effects of particle size distribution andparticle shape on the cyclic liquefaction response of granular materials. These problems areaddressed using numerical modeling of a cyclic constant volume direct simple shear test based on3D-DEM. The objectives are analyzed at the macroscopic level, which mainly includedetermining the stress-strain response and the number of cycles required for each sample to reachliquefaction. Additionally, state-based variables at the macroscopic level, based on void ratio, andat the microscopic level, based on inter-particle contacts and contact forces, are used to establishcorrelations with the cyclic liquefaction resistance of all samples. Additionally, this procedure istested at various relative densities in the hope of achieving some generalization.1.3 Layout of the thesisApart from the Introduction in Chapter 1, i.e., the current chapter, the dissertation consists of fiveadditional chapters which are organized in the following order:Chapter 2 focuses on the literature review, identifying the gaps in literature. The chapter willdiscuss the state of knowledge regarding the effect of particle size distribution and particle shapeon the response of granular materials under both monotonic and cyclic loading, using examplesfrom laboratory experiments and DEM tests.The methodology is discussed in Chapter 3. It covers the major highlights of the simulationsoftware used for the research. Additionally, it discusses how to design samples with a varietyof particle size distributions and particle shapes while keeping other interfering variables constant.Hereafter, the design of the sample preparation procedure is included followed by the cyclic simple7shear procedure. Numerous sensitivity analyses are performed to optimize the design parametersof the sample components, such as random particle orientations, particle numbers, and particleclasses each containing a particular particle size, in order to optimize accuracy and computationalefficiency.Chapter 4 examines the effect of particle size distributions on the cyclic response of granularmaterials. At two different relative densities, five different size distributions are prepared, asmeasured by the uniformity coefficient (Cu). The cyclic liquefaction resistance curves are thendetermined for 58 constant-volume cyclic simple shear simulations. At a low relative density, ithas been discovered that an increase in size distribution may increase the liquefaction resistanceuntil a certain limit is reached, followed by a decrement with further increases in sizedistributions. However, the observations are different when the same effect of gradation onliquefaction resistance is assessed at a higher relative density. Two state parameters, based on (a)void ratio (macro-based) and (b) coordination number (micro-based), are used to describe thetrend in the cyclic liquefaction resistance irrespective of size distribution and relative density.Micro-parameters are further used to describe the samples before and during the constant-volumecyclic loading.In Chapter 5, the effect of particle shape on the cyclic liquefaction response of granularmaterials is explored. Five distinct particle shapes, quantified by their aspect ratio (AR) andblockiness (B), are prepared at two different relative densities, encompassing both loose anddense assemblies. The cyclic liquefaction strength curves are then determined using 41constant-volume cyclic simple shear simulations. It is discovered that particle shape increases ordecreases the liquefaction resistance depending on the relative densities of the samples. The stateparameters derived using the critical state theory for initial void ratio and coordination number,exhibit a strong correlation with liquefaction resistance. Additionally, the evolution of the internalstructure is monitored to characterize the solid-fluid transition as a granular assembly approachesliquefaction after subjected to undrained loading cycles.In Chapter 6, the dissertation is concluded with summary, conclusions and recommendations8for future work.1.4 NotationsScalar variables are denoted in this dissertation by characters with normal letters (e.g., p for meaneffective stress), while second-order tensors are denoted by boldface characters (e.g.σ for stresstensor). The dyadic product between two tensors is denoted by a ⊗ b, alternatively expressedin index notation as ai jblm. The symbol : between two second order tensors A and B meanssummations over the consecutive pairs of indices in reverse order and is given as tr(AB) = A : B =Ai jB ji.9Chapter 2: Literature Review2.1 IntroductionThe literature review consists of summarizing the existing knowledge on the effects of a particlelevel parameter like particle size distribution and particle shape on the shearing behavior ofgranular materials in terms of both laboratory tests and DEM models. The extract from the reviewpoints out that the most fertile of these studies, which deal with a thorough investigation of aparticle level parameter, are primarily concerned with simple monotonic shearing in the case ofexperiments. In the case of DEM simulations, the majority of researchers have concentrated onconfining pressure, density, and sample preparation effects under both monotonic and cyclicshearing. However, the effects from size distribution and shape on cyclic liquefaction has beenexplored on a limited scale. In that small pool of studies, fundamental limitations of theexperimental outlook have been pointed out and adequate rationale is provided to justify the useof DEM in studying the effects of particle size distribution and particle shape individually underconstant volume cyclic shearing, which includes the isolation of these particle level parametersand detailed microscopic investigation otherwise improbable using real soils.2.2 Effects of particle size distribution on the monotonic responseThere are numerous methods for quantifying particle size distribution [51]. The difference betweenthe largest and smallest particles in an assembly can be used to denote size distribution [52]. Also,a particle diameter corresponding to a specific percentage by mass of a sample can be used as the10effective size of the distribution following the sieve test. For instance, D50 denotes the grain sizeat which 50% of the sample passes by weight [53]. Size distribution can also be measured by thefraction of a fine material (less than 0.075 mm) in a coarser material, and is often expressed asfines content (FC) [54, 55]. Particle disparity ratio (D/d), where D and d are the effective particlesizes of the base sand and the fines respectively [56], can be used to indicate the size distribution.Another well known size distribution parameter used in classifying soil is called the coefficient ofcurvature defined as Cc = (D30)2/(D10)(D60), where D10, D30, and D60 are the particle diametercorresponding to 10%, 30%, and 60% passing by weight, respectively. Finally, the coefficient ofuniformity (Cu) is a simple and well-known parameter used to measure size distribution of wellgraded soils defined as D60/D10 [30, 51, 57, 58]. Higher values of Cu correspond to samples withwider distribution of particle sizes.2.2.1 Experimental testingParticle size distribution may alter the frictional and dilatant characteristics of the monotonicbehavior of sands. This can be made evident by standardized triaxial, direct shear and simpleshear tests as followed by many researchers. Zelasko [59] has used dense Franklin and Ottawasands to assess the effect of Cu on monotonic behavior using triaxial tests. It is concluded thatdistribution has little effect on peak friction angle but a considerable effect on dilatancy orvolumetric behavior. Well graded samples used in this study were found to be more compressiblethan uniformly graded samples.However, Simoni and Houlsby [60] using Leighton Buzzard sand in drained direct shear tests,found that an increase in Cu can increase the peak friction angle, the critical friction angle, and thedilative behavior of the material.Carrera et al. [61] utilizing a series of triaxial tests on stava tailings, underlined that the effect ofparticle size distribution on drained resistance may come from overlapping parameters like particleshape. It is pointed out that if such a factor is excluded by design, distribution does not influence11the monotonic shear resistance at the critical state.Yang and Luo [62] conducted both drained and undrained triaxial tests using Fujian sand, andconcluded that size distribution has no effect on critical state friction angle. A contractiveresponse, on the other hand, is observed in samples with comparatively larger particle sizedistribution. This is reflected in the undrained response where well graded samples showed higherliquefaction susceptibility under undrained monotonic loading.Thus, experimental studies are divided about whether distribution has effects on strengthparameters in both peak and critical states, but they all concur that wider size distribution resultsin higher contraction or tendency of contraction than uniformly graded samples. Moreover, thestudy of size distribution should be conducted with caution [62], to avoid the risk ofunderestimating the influence of other parameters like particle shape and relative density.2.2.2 DEM testingNumerous researchers have also re-examined the results of monotonic tests conducted inlaboratories using DEM to examine the effect of particle size distribution on the aspects of shearstrength and volume change.Azema et al. [63] used bi-axial tests and Cantor et al. [64] used 3D-simple shear tests toinvestigate the monotonic response of dense disks and spherical assemblies, respectively. Thecritical state friction angle was found to be unaffected by size distribution in both studies. However,it has been reported that expanding the particle size distribution results in increased volumetriccontraction or a decrease in dilation during shearing. It should be noted that the samples hadvarying relative densities at the time of shearing, which may have interfered with the effect fromdifferent size distribution. These findings were later corroborated by Wood and Maeda [65] usingdrained triaxial simulations.Jiang et al. [42] prepared isotropic specimens at the same void ratio but with varyingdistribution values. The samples are assumed to be medium dense from the drained triaxial shear12response. The drained tests revealed that samples with a higher Cu have a lower peak shearstrength but a comparable shear strength at critical state. Additionally, as the particle sizedistribution was increased, the volumetric strain was found to gradually change from dilative tocontractive. The undrained tests concluded that increasing distribution may increase thesusceptibility to static liquefaction.Certain limitations may have impacted the observations of these studies. For example, all ofthe studies above either have not discussed or not considered samples with varying particle sizedistributions at the same relative density at the initial state. Given the effect of size distributionon extreme void ratios, comparisons at even the same void ratio may be impractical [13]. Thereare also other concerns that parameters, such as particle shape and mineralogy, could overlapthe effect of size distribution on shear response of granular materials. Finally, the density of thesamples examined under monotonic loading varies from moderate to dense, and it is unknown howsamples of varying sizes behave under shear in loose packings.Nonetheless, the literature indicates that particle size distribution has an effect on the peakshear strength of granular material but not on the shear strength at the critical state. In addition, anincrease in size distribution has consistently been shown to increase the volumetric contraction ofthe material in drained monotonic tests. As a result, higher size distribution has been associatedwith a higher static liquefaction potential of granular materials in undrained monotonic tests. Therole of size distribution in influencing shear resistance is clear from these observations usingsimple monotonic shearing. While the scope of work is restricted to monotonic shearing, theirsignificance cannot be overlooked. This has prompted some researchers including the currentresearch to investigate the effect of particle size distribution under more complex loadingconditions, such as cyclic undrained loading.132.3 Effects of particle size distribution on the cyclic response2.3.1 Experimental testingThe effect of particle size distribution on cyclic liquefaction resistance has been studied in someexperimental tests. Vaid et al. [30] explored the effect of Cu on the cyclic liquefaction resistance ofgranular assemblies using cyclic undrained triaxial tests on Earls Creek sand. To prepare samplesat different Cu, three different distributions of the same material were prepared with similar D50,Dr, and mineral composition. For the specific type of sand considered, samples that had higherCu were found to demonstrate higher cyclic liquefaction resistance in medium Dr (40-45%). Incontrast, in high Dr (greater than 60%) increase in Cu resulted in decrease in the liquefactionresistance. Meanwhile, the reason/s behind such contrasting cyclic response with change in Drwas/were beyond the scope of the study. Apart from that, there were no measurements of particleshape reported, and the practise of using the slurry deposition method should have introduceddissimilar inherent anisotropy, biases in the direction of particle deposition [8, 66], and possiblesegregation [67]. But these are possibilities limited by the data provided.Kokusho et al. [8] studied the effect of Cu on cyclic liquefaction resistance of river sands anddecomposed granite soils using triaxial tests, prepared at similar Dr, but different D50 and mineral100 101 102Cu0.00.10.20.30.4CSRDr = 20%Dr = 35%Dr = 50%Dr = 60%Figure 2.1: Effect of particle size distribution on liquefaction resistance (modified from [8]).14composition . The effect of Cu on cyclic liquefaction resistance was found to be negligible both inlow and high Dr values. The specific responses for river sands are shown in Figure 2.1. However,these results should be viewed in light of the following: The river sand itself is reported to becomposed of five different minerals, each of which can interfere with the effect of Cu. The highlyangular decomposed granite may have been subjected to particle breakage, resulting in changesin Cu during cyclic shearing. Moreover, the observations were based on the wet tamping methodof sample preparation, which may have introduced dissimilar inherent anisotropy into the samples[68]. Additionally, the differences in mean size and particle shape among different batches of riversand and granite soil could influence the aforementioned findings.Yilmaz and Mollamahmutoglu [31] investigated the effect of Cu on cyclic liquefactionresistance using triaxial tests on dense reconstituted Pasabahce sands, which share the same D50and mineralogy. The sample preparation method based on the ASTM D5311 standard [69]resulted in isotropic samples at a Dr = 60%. The liquefaction resistance was found to decreasewith an increase in Cu. However, the samples were gap-graded, and hence the D50 in thedistribution curve could not represent the mean size of the sample. In addition, the particle shapeof the fines was dissimilar to the coarse material and was not quantified, which can cast doubt onthe effective similarity of particle shape between samples.Wei et al. [56] investigated the effect of Cu on cyclic liquefaction resistance in silica-basedOttawa, Fujian, and Toyoura sands with non-plastic silty fines using triaxial tests. This study foundthat there is no direct link between Cu and cyclic liquefaction resistance. Apart from possessingdifferent Cu values, the sands also had different D50, particle shape and Dr, which might haveinfluenced the liquefaction response. The moist tamping method used in the study can result innon-uniform specimens [68].Capellaro et al. [70] conducted a series of simple shear tests to investigate the effect of Cuon cyclic liquefaction resistance using Christchurch sand. For samples prepared at the same Dr,increase in Cu was found to reduce the cyclic liquefaction resistance, and this reduction becomesmore pronounced as the Dr is increased.15To summarize, the researchers investigated the effect of size distribution on cyclic liquefactionresistance using a variety of materials, sample preparation methods, and cyclic tests. Apart fromthe overlapping parameters mentioned, the inability to quantify particle shape in this context, anabsence of unified or clear understanding of the effects of size distribution on cyclic liquefactionmay be due to the differences in the sample preparation and mode of shearing.The review of the literature reveals considerable differences in the way Cu affects cyclicliquefaction resistance, which might result from each or a combination of the following factors:(a) increasing the fines content in a base material to increase Cu can alter D50 which cansignificantly alter the cyclic liquefaction response [31, 71, 72] without any effect from systemcompliance [73], (b) samples with different size distributions can have different particle shapes,and (c) geological uncertainty. The overlapping effects can also arise from the state parametersamong compared samples like: (a) differences in relative density at the initial state i.e., sampleswith significantly different void ratios may show similar liquefaction resistance as Cus differ [70],clearly indicating void ratio is not suited for comparison of density, (b) sample preparationmethods can result in inherent anisotropy and segregation, and (c) type of loading like triaxial andsimple shear tests can make a significant difference [74, 75] in the shearing response. Thus, it isimportant to recognize that the overlapping effects of multiple parameters complicate theseparation of the effect of particle size distribution on cyclic liquefaction in laboratoryexperiments. Finally, due to the difficulty of measuring and tracking actual particle-particleinteractions, physical experiments can only hypothesize their significance in explaining theeffects of particle size distribution on cyclic liquefaction. A more holistic investigation shouldinclude the ability to decouple the effect of size distribution and the ability to record macro-level(stress-strain) and micro-level (particle-particle) responses in order to establish causality.162.3.2 DEM testingAt the time of investigation, it is unknown if any published study exist that has used either 2Dor 3D-DEM to investigate the effect of Cu on cyclic liquefaction resistance. Hence, this topic isinvestigated as the first major topic of the dissertation and presented in Chapter 4.2.4 Effects of particle shape on the monotonic responseThe evaluations of shape include ocular representation, which classifies it into four main types:rounded, sub-rounded, sub-angular, and angular. There are also manual measurements with thehelp of charts [76] and sieving techniques to measure the elongation index, also known as theaspect ratio. Other parameters that are fundamental in the domain of particle shape measurementsare roundness and sphericity [35, 37]. Since manual measurements of the 2D-representation ofparticles are both overwhelming and oversimplified, sophisticated 3D image analysis [77, 78] on acomputer base has been used in the classification of particle shapes, which allows dealing with alarger number of particles. The effect of particle shape is hereby studied under monotonic shearing.2.4.1 Experimental testingParticle shape has been found to have a considerable effect on the shear resistance of granularmaterials. Different researchers have used experimental methods to investigate the effect ofparticle shape on shear response parameters such as peak shear strength, friction angle, dilatancy,and critical state, using monotonic drained tests [76, 79, 80]. As demonstrated in these studies,assemblies composed of angular and elongated particles exhibit greater drained peak and criticalfriction angles than rounded sands.Shinohara et al. [80] conducted drained triaxial tests on angular and rounded granularassemblies of steel powders and concluded that angular assemblies exhibit a greater peak frictionangle than rounded assemblies, and that this enhanced effect of angularity becomes lesspronounced as the void ratio decreases. However, in order to study the shape effect, the17researchers have not specified whether the initial sample densities are comparable across samples.While it is critical that samples have comparable density conditions to avoid shape interference,there are disagreements regarding the appropriate density parameter.Numerous researchers have noted [13, 76, 81] that samples should be compared using the samerelative density rather than the same void ratio, as the former takes into account the extreme voidratios of the samples, affected by particle shape. However, these extreme void ratios are not uniquebecause they are dependent on the testing method and reproducibility of the test [13].Apart from the issue of reproducibility with relative density, sample preparation has emergedas a significant factor that may complicate the effects of particle shape [39]. In general, thestudies demonstrate that particle shape does have an effect on the monotonic shear response underdrained loading. To isolate the effect of shape, it is critical to eliminate certain particle-level andstate parameters. By understanding these design philosophies, other investigators may be able toadvance the study of the effects of particle shape under monotonic loading without drainage.Georgiannou [82] has conducted undrained monotonic triaxial tests to investigate the stabilityof sand matrix mixed with angular fines. The results from samples prepared at the same “quasi-natural” relative density indicate that angular and flaky mica fines can cause instability in themonotonic response, whereas rotund silty fines demonstrated a stable response.Ghadr and Langroudi [83] have looked at the effect of shape on the monotonic response ofloose samples prepared with the same initial void ratio. They discovered that sub-angular F-131sands with a similar D50 value to sub-rounded Urmia lake sands have comparable shear strengthsat the point of phase transformation.Tsomokos and Georgiannou [74] have conducted similar undrained triaxial tests onsub-angular Ham river sand and rounded M-31 sand at the same relative density. They observedthat for a Dr ≃ 41%, angular sands exhibited a stable response (constant increase in strength),whereas rounded sands exhibited an unstable response (decrease in strength after a transientpeak). Increased relative density at 63% increased the dilative tendency of sub-angular samples.In summary, observations indicate that as particles become more rounded, shear strength18decreases, but the effect is more pronounced at higher relative densities. Undrained experimentalsamples have studied relative densities greater than 40%. More importantly, a systematic study ofparticle shape which covers both the roundness and sphericity aspects of particle shape is notfound.2.4.2 DEM testingDEM studies provide an alternative and often simplified method for simulating particle shape andits effect on the shear behavior of assemblies. It can also help to separate the effect of particleshape from other parameters such as particle size distribution, relative density, and changes in soilstructure (fabric), to name a few. Furthermore, DEM can also adequately expose micro-mechanics,which is rare in experiments.Nouguier-lehon et al. [43] discovered that polygonal and elongated shaped discs indicatedsuperior shear resistance to their circular counterparts when subjected to monotonic loading atmoderate densities. However, they noted that polygonal discs were significantly more compressibleat low densities, despite their higher shear strengths compared to the circular disks. One mightargue that increased compressibility or volume reduction for polygonal disks may result in a higherpore pressure development in constant volume shear tests. Therefore, increased compressibilityindicates that such non-circular particles may be more susceptible to liquefaction under undrainedmonotonic and cyclic loading. The macro observations of peak shear strength were also discussedbased on the differences in initial anisotropy and the particle rotations occuring from the change inparticle shape.Anthony and Kuhn [84] demonstrated that ellipsoidal particles exhibit greater monotonic shearresistance than spherical particles in medium dense samples. Additionally, the study emphasizesthe effect of contact anisotropy in such differences in shear resistance. It implies the formationof more contacts across the principal shear plane in the case of ellipsoidal particles, mobilizinga greater number of them along the plane, contributing to the higher shear resistance. However,19little emphasis is placed on the overlapping effects of initial conditions interfering with the effectof particle shape on the monotonic shear resistance.Jiang et al. [42] conducted monotonic constant-volume tests using clumped spheres wheresamples posses similar void ratio (not at the same relative density), and identical values of Cu andD50. It has been discovered that elongated clumped particles have a greater undrained monotonicshear resistance. Nonetheless, the parameters defining the shape of the particles are not variedsystematically, and the loading criterion is monotonic.To summarize, it can be concluded that researchers have shown an adequate appetite forseparating the effect of particle shape from other factors such as particle size distribution, surfacetexture, mean particle size, and following an identical sample preparation, among others. Theclear intention was also found to investigate and link the micro-structure with the macro-response,a lesson which can be extrapolated in undrained studies of particle shape as well. Finally, fromthe conclusions drawn from one of the studies, it is evident that angular particles are morecompressible in loose densities and, hence, if a potential undrained test is conducted, increasedvolume reduction may result in a reduction in shear strength. But, it is an informed hypothesisthat has been tested in the thesis under cyclic shearing.2.5 Effects of particle shape on the cyclic response2.5.1 Experimental testingVaid et al. [9] determined the liquefaction resistance of rounded Ottawa and angular, flaky Tailingssands using constant-volume cyclic direct simple shear (DSS) tests. The samples had similar D50,but different Cu and mineral composition. In Figure 2.2(a), at low confining pressures and a relativedensity greater than 40%, the angular sand is more resistant to liquefaction than the rounded sand.However, at high confining pressure as shown in Figure 2.2(b), angular sands may be less resistantthan rounded sands. The particle shape of the two materials were not quantified and the gradationcurves before and after liquefaction for both ottawa and tailings sand were not compared to check2035 45 55 65 75Dr [%]0.080.100.120.140.160.18CRR10σv = 200 kPa(a)45 55 65 75 85Dr [%]0.080.100.120.140.160.18CRR10σv = 1600 kPaOttawaTailings(b)Figure 2.2: Effect of particle shape on liquefaction resistance (a) at σv = 100 kPa and, (b) at asignificantly higher σv = 1600 kPa (modified from [9]).for change in size distribution due to breakage.Hubler et al. [85] investigated the cyclic response of three uniform gravels using DSS testsunder K0 condition. The samples were characterized by similar mineralogy and Dr but differentD50. It was concluded that increasing particle angularity do not have a significant effect on theliquefaction resistance of medium dense (Dr = 47%) and dense (Dr = 86%) samples for small cyclicstress ratio (CSR). The particle shape was quantified using roundness, but other aspects of shapelike sphericity and aspect ratio were not studied. Moreover, the K0 method of sample preparationusing dry tamping might have caused significant differences in contact anisotropy. These factorscombined could have a remarkable interference with the effect of particle shape on cyclic response.Rui et al. [86] used cyclic triaxial tests to compare the liquefaction resistance of calcareoussand, quartz sand, and steel balls at a high relative density of 70%. It was found that non-sphericalsand particles accumulate excess pore pressure more slowly and require more cycles to liquefythan spherical particles. However, they do not completely exclude the effects of mean particle size,particle mineralogy, surface roughness, and particle size distribution on cyclic liquefaction. Thismay cast doubt on their conclusions.These studies also only looked at medium or dense samples. It is unclear how particle shape2110−1 100 101 102 103NIL0.00.10.20.30.4CSRGlass beadsOttawa angular(a)10−1 100 101 102 103NIL0.00.10.20.30.4CSRGlass beadsOttawa angular(b)Figure 2.3: Effect of simulated particle shape of real sands on liquefaction resistance for a (a) loosesample, and (b) dense sample (modified from [10]).affects cyclic liquefaction of loose samples. Moreover, in light of the fact that the laboratory testsonly obtained macroscopic responses, what occurs at the particle level and how microscopic datacan be linked to macroscopic data merits further investigation.2.5.2 DEM testingThese gaps can be addressed numerically using the three-dimensional discrete element method(3D-DEM). It consists of (a) separating particle shape from other parameters, (b) preparingsamples covering a range of packing densities, and (c) exploring microscopic descriptors linked tothe macroscopic finding. Generally, there are two ways to consider particle shape in DEM. Thefirst one is to introduce a rolling resistance model [87] to spherical particles, which is not thefocus of this dissertation, and the other explicitly adopts non-spherical geometries. To date,particle shape geometries consist of superquadratics [88], polyhedrons [89], agglomerates [46],non-uniform rational B-Splines based particles [90], and real particle shape mapping via level-setimaging [91].Using 2D-DEM, Ashmawy et al. [10] clumped a series of overlapping circular elements toapproximate real particle shapes and found that the sample possessing the highest level of particle22angularity (surface curvature) exhibits the highest resistance to liquefaction, whereas for loosesamples particle angularity showed little influence on liquefaction resistance as shown in Figure2.3. However, the shape parameters are varied arbitrarily and overlapping effects are introducedby the dissimilar particle size distribution among samples.Kuhn et al. [92] and Morimoto et al. [93] have investigated the liquefaction resistance ofclumped spheres, but neither has varied the shape parameter to investigate its effect on cyclicliquefaction resistance. Whereas, the current research objectives include studying the effect ofnon-spherical particles on cyclic liquefaction resistance. Particle shape can be defined based oncertain features viz., elongation, angularity, and convexity. In this research, only elongation andangularity are captured, using a non-spherical 3D descriptor known as superquadratics [94].Particles composed of the superquadratics family allow us to maintain computational efficiencywhile studying the different aspects of particle shape [95] like roundness and sphericity. Thisstudy investigates superquadratic particles with two shape parameters, aspect ratio andblockiness, which correspond to the geometrical properties and surface curvature of the particles,respectively. Additionally, the current shape parameters allow for the proper separation of particleshape from other parameters such as particle size distribution and surface roughness, as well asensuring that each assembly contains particles with distinct shapes. The granular assemblies arethen subjected to constant-volume cyclic DSS loading to investigate the effect of aspect ratio andblockiness primarily on liquefaction resistance.2.6 SummaryThe extensive review of experimental data on particle size distribution indicates that otherparameters such as particle shape, particle size, and particle mineralogy, among others, canobstruct its effects on shearing, including cyclic shearing. Due to the fact that these effectsoverlap, it is difficult to attribute differences in monotonic and cyclic response to a singleparameter. More such effects can come from non-identical sample preparation, unequal density23among samples. This is also true when it comes to isolating particle shape effects. Whilelaboratory experiments demonstrate the importance of considering particle size distribution andshape, they do not provide the level of control over these parameters that DEM provides. DEMalso allows in significant monitoring of the sample to the particle level, which can raise thepossibility of connecting even inter-particle contacts to the macroscopic response of cyclicliquefaction. According to the review of the literature on 3D-DEM studies, it is clear that particlesize distribution and shape have been studied primarily under monotonic loading conditions.Considering cyclic liquefaction modeling may signal to significant increase in computation timeaggravated with the consideration of particle shape and particle size distribution. This may haveengaged the research community from abstaining from direct research into these particle levelparameters.The current study uses advanced hi-performance computational facilities to study the effectsof particle size distribution and particle shape on cyclic liquefaction of granular materials. It uses3D-DEM and aims to bridge this gap in literature both from the viewpoint of cyclic liquefactionand connecting the trends with particle level observations.24Chapter 3: Methodology3.1 IntroductionThis chapter starts with the introduction to the DEM tool for analysis, the basic structure of its inputscript, and the contact laws governing the inter-particle interactions. The tool is then utilized todescribe the details of sample preparation and shearing protocol followed for all simulations. Theparticle size distribution and particle shape models used in the DEM tool are described, and they areoptimized based on initial sets of simulations. The optimization schemes mainly involve the choiceof adequate numbers of particles in the system, discrete classes of particles each with a specificparticle size in order to design the particle size distribution, and particles in each class. Finally,a small set of constant volume cyclic tests are modeled and their characteristics are adjudgedqualitatively based on experimental results from literature.3.2 DEM programThe numerical platform is chosen based on open source systems that can be easily accessed, havean active community and documentation, provide capacity for modeling particle size distributionand particle shape in particular, and capable of conducting massively parallel molecular dynamicssimulations.253.2.1 IntroductionThe modeling is conducted using the open source version (public) of LIGGGHTS [96], whichstands for LAMMPS Improved for General Granular and Granular Heat Transfer Simulations,where LAMMPS is abbreviated from Large-scale Atomic/Molecular Massively ParallelSimulator. The versions of open-source LIGGGHTS prior to 2017 can be used to study cyclicliquefaction with spherical particles, but the incorporation of superquadratic particle function inthe latest version 3.8.0 [97] made the study of particle shape possible. LIGGGHTS usesNewtonian laws of motion to simulate the kinetics of the particles, while the kinematics originatesfrom the inter-penetration between colliding particles. The contact detection algorithm based onnon-spherical particles are provided in Podlozhnyuk et al. [97].3.2.2 Contact lawsThe contact interactions are based on normal collision, tangential sliding, and rotations, whereeach interaction has an elastic and an energy dissipative component. The contact laws used for thecurrent study are briefly described below.For total normal force Fn, the Hertz model damping is adopted from [98] which consists of anelastic and a dissipative component described as follows:Fnel =23f1(E,ν)√Reffξ 3/2n , (3.1)where Fnel is the elastic component of the normal force, f1 = 2E/3(1−ν2)is a function of Young’sModulus E and Poisson’s ratio ν , Reff is the equivalent radius of the particles in contact, ξn is thenormal component of the particle overlap. The dissipative component of the normal force Fn(dis) isgiven as:Fn(dis) = f2(E,ν ,ε)√ReffA√ξnξ˙n, (3.2)where, f2 = E/(1−ν2)depends on E, ν and the coefficient of restitution ε , A is a function of E26and ν , and ξ˙n is the tangential displacement rate between particles in contact.For tangential forces, the Mindlin-Deresiewicz formulation [99] is adopted where the totaltangential force F t is divided as:F t(el) = f3(E,ν)√Reffξ 3/2t , (3.3)where F t(el) is the elastic component of the tangential force, f3 = 8Geff and1/Geff=[2(2− v1)(1+ v1)+2(2− v2)(1+ v2)]/E ,and the dissipative component of thetangential contact force F t(dis) is given asF t(dis) = f4(E,ν ,ε)√ReffA√ξnξ˙t . (3.4)The actual expression of f4 is referred to Di Renzo et al. [100] for the interest of the reader. Thevalue of F t has a maximum value governed by the Mohr-Coulomb sliding criterion F tmax = µFn.The modified elastic-plastic spring dashpot model (EPSD2) is used for including rollingmoment for spherical particles in the current study in Chapter 4, whereas in Chapter 5 rollingfriction model is not considered since it may interfere with the rolling resistance offered by thenon-spherical particles. The total rolling torque Mr in this model has two components: amechanical spring torque (Mr(el)) and a viscous dissipative torque (Mr(dis)) as shown below:Mr(el) = 2JnReffFn, (3.5)where Jn is a dimensionless coefficient which varies theoretically from 0.25 to 0.5, andMr(dis) = 0. (3.6)The maximum value of the rolling torque is bounded by Mrmax = µrReffFn. Because the currentstudy in both Chapters 4 and 5 do not deal with dynamic flow, the viscous dissipative torque is27disabled in the EPSD2 model [45, 101]. This also avoids the complication of calibrating a largenumber of unknown parameters, such as the rolling viscous parameter. For the same reason, addingto its limited effect on the shear resistance [87], the twisting torque is also disabled as followed inrecent studies using LIGGGHTS [45, 97].3.2.3 Input file structureThe structure of a typical LIGGGHTS input script is mainly composed of four portions:1. Define: This is the initial stage of the script, where the particle size and geometry aredefined followed by random insertion of particles within a bounding box with nointer-particle contact. Finally, the contact model parameters for the inter-particle andparticle-wall interactions are described.2. Direct: This stage provides commands to register particle motion, contact detection, andcontact forces in relation to applied compression as well as cyclic loading and other boundaryconditions.3. Derive: At this stage, the recorded particle interactions are used to calculate macroscopicparameters such as shear stress and shear strain, as well as microscopic variables such ascoordination number and anisotropy, at specified time intervals.4. Dump: At this stage, at predefined time-steps, snapshot of both macro and micro parametersare printed to one or more data files.3.2.4 Use of high performance computingLIGGGHTS has been used in a number of HPCs to generate data relevant to this thesis. The mostimportant systems used are: (a) DesignSafe [102] provided by the Texas Advanced ComputingCentre at the University of Texas at Austin and (b) UBC-Sockeye provided by Advanced ResearchComputing at the University of British Columbia.283.3 Sample preparationThe objective of sample preparation is to compress an assembly of particles isotropically withperiodic lateral boundaries and rigid walls enclosing the top and bottom sides. This setting isidentical to that of a bi-periodic simulation box, which has been used by a number of researchers[103–105]. The material parameters for the simulation are presented in Table 3.1. The values of E,ν and ε are adopted from Sufian et al. [106]. A low εpw is used to dissipate the kinetic energy of theparticles and achieve a fast and steady compression. During compression, the value of µ is variedto obtain different sample densities. During shearing, however, µ is held constant at 0.5, a fairlycommon value used by modelers since higher values neither affect the the maximum void ratio norincrease the sheraing resistance [104, 107]. The sample preparation involves two major steps:Table 3.1: DEM parameters involving interparticle and particle-wall interactions.Description ValueParticle density, ρ 2500 kg/m3Young’s modulus, E 70 GPaPoisson’s ratio, ν 0.25Coefficient of tangential friction, µ 0.5∗Coefficient of rolling friction, µr 0.1†, 0.0‡Coefficient of restitution, ε 0.2∗ used during shear; † 0.1 is used to study particle size distribution; and ‡ 0.0 for particle shape.(a) particle insertion and (b) compression. Particles with mean particle size D50 of 0.0025 m arerandomly inserted into a 0.15 m cubical box in a randomly dispersed state [108] with no particle-particle contact. The scalability of the system is tested using different total number of particlesinserted to the system as discussed later in Sec. 3.6. The size of the insertion box is chosen toavoid any inter particle forces thus ensuring there is no overlap. To ensure that all particles fitinto a small cubical box, the particles are inserted in order of decreasing particle size. Two rigidplates are also inserted inside the box, at the top and bottom of the dispersed sample along the zdirection. Subsequently, the rigid plates and lateral periodic boundaries are used as servo walls to29(a)(b)Figure 3.1: Particle arrangements and boundary conditions for an assembly with 12600 particles andsample size ≈ 20×D50: (a) at the end of isotropic compression with bi-periodic lateral walls, and (b)during volume cyclic shearing. The grey particles at the top of the sample form the shearing plate, whilethe particles on the bottom form the rigid plate.compress the assembly to 10% of the final target value of 100 kPa. As a result, the sample densityis increased and a stable contact network is formed. The value of µ is varied between 0.0 and 0.5in each of the trials to obtain different target densities. Further increases in µ had no effect on voidratios, and thus the highest value of 0.5 is used during the final stage of compression and maintainedthroughout subsequent shearing. The second stage of compression applies compression p0 equal to100 kPa, which completes the sample preparation process. Figure 3.1(a) depicts a typical samplewith spherical particles after sample preparation for Cu = 2. The snapshots of the samples aftersample preparation for different particle size distributions are presented in Figure 3.2, while for30(a) (b)(c) (d)Figure 3.2: Snapshots of the samples representing relevant particle size distributions where, (a) Cu = 1,(b) Cu = 1.5, (c) Cu = 2, and (d) Cu = 3.different particle shapes are shown in Figure 3.3.Because the sample preparation procedure compresses the samples by moving all sixboundaries with identical servo control, the contacts formed at the end of compression (contactdistribution) are assumed to be unbiased in a particular direction. It is assumed that thishomogeneous, or frequently referred to as isotropic, distribution of contacts exists in all threedimensions. Several researchers [105, 109, 110] have used DEM to follow this procedure of sixwall servo control for preparing isotropic samples with spherical particles. However, becausenon-spherical particles were used in a part of this study, it is unknown whether the contact forcesand distribution of contacts are truly isotropic, which could be the subject of future research.31(a) (b) (c)(d) (e)Figure 3.3: Snapshots of samples with different particle shapes where (a) AR = 0.5, B = 2, (b) AR = 1, B =2, (c) AR = 1.5, B = 2, (d) AR = 1, B = 4, and (e) AR = 1, B = 8.3.4 Shearing processFollowing compression, the sample is prepared for shearing under simple shear conditions. Asillustrated in Figure 3.1(b), a top layer of (black) particles is used as the shearing plate extendinginto the sample to a diameter of D50, while a similar layer at the bottom acts as a fixed base. Thespecified thickness of the particle plate is adequate to represent a single particle layer, as pointedout by Zhang and Kuhn [111]. The convex and concave surfaces of the particle plates for shearresemble the shearing process found in natural soils and is a good alternative to rigid flat plateswhere slippage may occur [105]. The movement is restricted along the z axis to avoid changingthe volume of the sample and thus simulating a constant-volume condition.Cyclic simple shearing is undertaken by moving the top plate horizontally at a constant velocity(vy) of 0.01 m/s in order to shear the sample. The shearing direction is reversed each time the shearstress τ reaches its target amplitude of τamp. This is referred to as the uniform amplitude cyclicshear test with strain control [92]. This is consistent with laboratory tests using cyclic simple shear32with symmetric and cyclic loading [112]. The ratio of τamp to p0 is referred to as the “cyclic stressratio” (CSR)CSR =τampp0. (3.7)To maintain a quasi-static shear regime, the inertial number I = γ˙D50√ρ/p is monitored,where γ˙ = |vx|/h is the shear strain rate, h is the sample height, ρ is the density of sphericalparticles, and D50 is the mean diameter of the particle assembly. A value of I < 10−3 indicatesthat the sample is in quasi-static regime during cyclic shearing [113]. The movement of the topplate at a velocity of 0.01 m/s ensures a steady γ˙ ≈ 0.125s−1 throughout the shear resistant state,complying with the quasi-static criterion suggested by Martin [114]. The higher values of I atthe liquefied state arise from unstable deformation, which is unavoidable while modeling cyclicliquefaction and does not depend on the rate of shearing [105].Notably, the current study makes use of the widely used constant volume condition tosimulate the fully undrained state of saturated granular materials during cyclic shearing. Whileshear behavior of completely saturated granular materials under undrained loading can bereplicated using constant-volume conditions [49, 50, 115], fluid-particle interactions in the porespaces of saturated granular materials, direct measurement of the pore pressure and total stresscan be calculated only when the pore fluid is directly modeled using DEM and computationalfluid dynamics (CFD) coupling [116, 117]. It was recently noted [118] that realisticconsiderations such as pore air compressibility, surface tension, and air solubility in the pore fluidcan be modeled using 3D-DEM and CFD coupling while being deficient in constant volumeconditions to simulate undrained loading.3.5 Description of particle size distributionThere are several methods for quantifying particle size distribution [51]. The difference betweenthe largest and the smallest particle size in a sample, referred to as the range, can be used for thispurpose. However, it is problematic to find the largest and smallest particles. Other techniques33involve finding the particle diameter on the grain size distribution curve that corresponds to aspecific percent of soil mass passing through a sieve. For instance, D50 denotes the grain diameterat which 50% of the sample by weight passes through a sieve. The percentage of fines in a coarsermaterial, often expressed as fine content (FC) can also be used to describe particle size distribution.Particle size disparity, which reflects the size difference between the coarser fraction and the finefraction, can also be characterized by the disparity ratio (D/d), where D and d are the effectiveparticle sizes of the base sand and the fines, respectively [56]. The coefficient of curvature (Cc) isanother well-known quantifier for particle size distribution. The coefficient of uniformity (Cu) is asimple and well-known parameter used to measure size distribution of well graded soils defined as[30, 51, 58] :Cu =D60D10(3.8)where, D60 and D10 are defined in the current study as the sieve sizes at which 60% and 10% ofthe particles are passing. Samples with higher values of Cu indicate wider distribution of particlesizes.In the current research, Cu is used to quantify particle size distribution. Due to the assumptionthat all particles in the current DEM model have the same density, a volume-based particle sizedistribution curve is used in place of the mass-based sieve procedures typically used inexperiments. The particle size distributions are modeled as a log-linear distribution curve,replicating the experimental procedure of sieve analysis and the cumulative frequency diagram.The subsequent subsections lay out the procedure in detail.3.6 Scalability or choice of optimum particle numbersA soil element in a laboratory contains overwhelming number of particles to be considered in a3D-DEM platform even with the advent of modern computing, and further, performing basic testslike analyzing packing properties, and monotonic tests [119]. It is even more challenging tosimulate cyclic liquefaction using a large number of virtual particles since it involves a lower34strain rate than monotonic tests and often large shearing cycles to reach liquefaction. Nonetheless,there are techniques for reducing the particle count in order to facilitate the modeling of cyclicliquefaction. By using lateral periodic boundaries, an “optimum number” of a few thousandparticles can emulate a shearing test in DEM just as well as an assembly of millions of particlesby minimizing the edge effects [103]. The argument could be based on the impact of totalparticles in an assembly on key metrics like void ratio and coordination number. Additionally, thecriteria for choosing an optimum number may include comparing the mechanical response toshearing of samples with varying total number of particles. On that same basis, a small number ofparticles less than the optimum number may be insufficient and thus susceptible to theseimportant metrics. Thus, in the current DEM simulations approximate optimal range of the totalnumber of particles in a sample are determined.Numerous tests are conducted on assemblies with varying particle counts to determine theirpre-shearing void ratio (e), and the mechanical coordination number (zm) defined as [120]zm =2Ncv−N1pNp−N0p −N1p, (3.9)where Ncv denotes the total number of contacts in the assembly of volume V ; N0p and N1p denotethe number of particles with no contact and particles with only one contact, respectively. Theparticles with zero contact is defined here as floaters [121] and those with one contact as rattlers[114]. They are expected to not contribute to the force network and are thus excluded from theaverage contact count [122, 123]. Additionally, because the current research is concerned withcyclic liquefaction, the number of cycles required to achieve initial liquefaction i.e. ru = 0.95, istaken into account when determining the optimal particles. Figure 3.4 illustrates the effect of thetotal particle number (Np) on the e, zm of the compressed sample, and the number of cycles requiredto reach initial liquefaction NIL. All samples for this objective are prepared with spherical particlesat Cu = 3, coefficient of tangential friction µ = 0.3 , where other material properties consistent withTable 3.1. The mean effective stress at the pre-shearing stage (p0) defined as tr(σ )/3 is kept equal35to 100 kPa. Over a given volume (V ), the stress tensor (σ ) [124] is calculated as:σ =1V ∑i∈Ncvl i⊗ f i, (3.10)where l i is the branch vector joining the centroids of two neighbouring particles for interior contactor connecting the particle centroid and the contact point for exterior contacts, f i is the contact force,⊗ represents dyadic product, and the summation ∑ covers Ncv. The samples are subsequentlysheared under a CSR = 0.25.It is observed that as Np is increased from a small value of 2000 to 8000, e increases rapidly(refer to Figure 3.4(a)) which results in decrease in zm as shown in Figure 3.4(b). Subsequently,when Np values higher than 8000 i.e., 16000 and 24000 are considered, the values of both e andzm appears to stabilize. Similarly, as shown in Figure 3.4(c), NIL is also found to stabilize asNp exceeds a value of 8000. Therefore, it can be concluded that changes in the key parametersused for the analysis like e, zm and NIL have insignificant effect on assemblies with a minimum of8000 particles for Cu of 3. It should be noted that lower values of Cu requires even lower numberof particles for this optimization [125]. As discussed later in 3.7, further increase in number ofparticles with may be required to ensure that all discrete classes and particles within each class areadequately represented for the intended particle size distribution.The effect of the random positioning of particles during the sample preparation stage iscaptured using the Park Miller pseudo randomization technique initialized by a seed number Sn asfollowed in LIGGGHTS (refer manual). Please note that the seed number is typically a primenumber above 10000, but for simplicity the three seeds considered here are referred as 1, 2, and 3.The margin of error for the three combinations of seed are evaluated for e, zm, and Np and shownin Table. 3.2 in columns 4, 5 and 6 for a total of 15 samples. The sensitivity of Sn on the results isshown in terms of margin of error (ME) in % [126] ME = zα ×√σ2n where, α is the confidencelevel equal to 95%, z is the quantile, σ is the standard deviation, and n is the sample size, which isequal to 3 in this case. ME is found to decrease sharply as Np increases, implying that random362 4 8 16 24× 103NP0.520.540.560.580.60e(a)0 5 10 15 20 25Np5.055.105.155.205.255.305.35z m(b)2 4 8 16 24× 103Np070140210280350NIL(c)Figure 3.4: The effect of total number of particles (Np) in the assembly on the initial values of (a) e, (b)zm, and (c) NIL. The effects from random particle distributions (change in seed number) on packingproperties and liquefaction resistance are demonstrated by the error bars.positioning has little effect on assemblies because they are made up of a larger number ofparticles. When the effect of ME is compared among e, zm and NIL, it is observed that zm is theleast affected by change in Sn followed by e. In both cases, the ME value is well below a toleranceof 5%. While NIL is found to be most sensitive to Sn at Np = 4000. However, as Np increasesbeyond 4000, ME decreases to less than 5%. In general, this sensitivity analysis indicates that avalue of Np ≥ 8000 is also suitable for removing the effect of random particle positioning.37Table 3.2: The effects of total particles Np and random positioning seed Sn on the void ratio e, mechanicalcoordination number zm, and liquefaction resistance indicator NIL.Test Sn NpME (%)e zm NIL1-312000 2.9 0.75 7.30234-614000 0.8 0.33 15.10237-918000 0.6 0.19 3.102310-12116000 0.38 0.13 4.352313-15124000 0.5 0.11 3.2523Please note random positioning changes with seed number (Sn).Higher ME signifies increased effect from change in Sn.3.7 Choice of optimum particle classesA number of particles with identical particle size can constitute a class of particles. A smallelement of granular material may contain different particle sizes and thus can possess severalparticle classes (Nc). In laboratory-based sieve tests, only a limited number of discrete sieve sizesare used to retain specific classes. Such discrete classes can also be designed to represent aspecific size distribution using DEM modeling. These discrete classes should address twoprimary design considerations: first, determining the optimal particle classes to represent a givensize distribution; second, each discrete class should have a sufficient number of particles to ensurethat the class is representative. Viovret et al. [127] demonstrated using 2D-DEM tests that thestatistically representativity of particle size classes for a well calibrated polydispersed samplescan be determined on the basis of void ratio of a packing. It is found that increasing the number ofclasses in a particle assembly beyond an optimum value has no effect on e at the initial state. It382 6 10 14 18Nc0.440.460.480.500.52e(a)0 10 20 30 40 50 60Ncp0.540.580.620.660.70e(b)Figure 3.5: Effect of (a) total particle classes (Nc), and (b) total particles in each class (Ncp), on the voidratio e of samples with Cu = 3, p0 =100 kPa.was also identified that a broad size distribution with well-represented populations of differentsizes requires many more particles than a narrow distribution. Using 3D-DEM tests, Mutaburakaet al. [128] specified that an optimum number of Nc = 10 is required to correctly represent aparticle size distribution. In the current study the optimum number of Nc is investigated for a3D-DEM sample with a Cu = 3 at µ= 0.0 with 12000 particles, where the other material propertiesare depicted in Table3.1. It can be observed in Figure 3.5(a), that the effect of Nc on the initialvoid ratio e of the sample saturates after Nc ≃ 9. To be on the safe side, a Nc value of 10 is usedto represent any sample designed in the dissertation. A typical distribution of particle sizesconsidered in the current dissertation is shown in Figure 3.6. There are ten classes in thelog-linear distribution considered in this study, as marked from C1 to C10 with the arrowsshowing the possible value of D10, D50, and D60 for a value of Cu = 2. Following the work ofViovert [127], it is assumed that the particle volume in individual classes are identical to oneanother. With that assumption it is obvious that C10, with the largest sized particles will containthe minimum particles. As shown in Figure 3.5(b), Ncp = 30 is sufficient to represent the classC10, as any further increase in the value of Ncp will not effect the initial value of e. This is lessthan the findings of Mutaburaka et al. [128], who recommended an optimum value of Ncp = 50.392 4 8 12D [mm]020406080100FV [%]C1C105010Figure 3.6: Design of the log-linear particle size distribution with ten particle classes (C1-C10) for atypical case of Cu = 2 and D50 = 5 mm.Expectedly, the other classes, C1 to C9, will have more particles than 30 and will thus beadequately represented by definition in the size gradation. Overall, regardless of the sizegradation and shape considered, the optimum criteria established herein for Nc and Ncp are alwaysfollowed when developing or shearing a sample. A few points must be noted here after the designof particle size distribution is discussed:• It should be noted that the process of sample preparation with various particle sizedistributions can be replicated for other PSD curve shapes than the log-linear used in thisthesis as long as the design element of the optimum number of classes and particles in eachclass is respected.• Additionally, all samples in this study are designed to have the same D50, hence in eachsample the total volume of coarse and fine particles is the same. This results in a significantincrease of the number of fine particles while exploring a larger size distribution, whena minimum number of particles is to be present in each class size. This design featuredistinguishes this study from some others in which a fine granular material is mixed witha coarse material to alter the size distributions. A more detailed and targeted study in the40future may address the micromechanical aspects of the explanations in studies such as thatof Yamamoto et al. [129].3.8 Description of particle shapeTo come close to the complexity of real particle shapes, shape mapping via level set imaging,multisphere models, and other methods may be required. However, a much simplified family ofshapes known as superquadratics can capture the main particle shape aspects like roundness andsphericity with significantly advantage in computational costs [130]. The current study modelsshape with two indicators related to the superquadratic family, namely aspect ratio (AR) andblockiness (B) as discussed in detail in the following.3.8.1 Shape parametersThe class of superquadratic particles can be expressed as [94]:F(x,y,z)≡(∣∣∣∣xa∣∣∣∣n2 + ∣∣∣∣yb∣∣∣∣n2) n1n2+∣∣∣∣ zc∣∣∣∣n1 −1 = 0, (3.11)where a, b, c are the half-lengths of the particles along the principal axes as illustrated in Figure3.7, and, n1 and n2 are blockiness parameters.By setting a = b and n1 = n2 = B, we can reduce the five shape parameters (a, b, c, n1, n2) totwo independent parameters viz., Aspect ratio (AR = c/a) representing the particle geometry andBlockiness (B) representing the surface curvature. All values of AR and B discussed hereafter aretaken into account for analysis in Chapter 5, and are highlighted in yellow and green, respectively,in Figure 3.8. The values of AR considered in this study are 0.5, 1 and 1.5. To exclusively studythe effect of AR, the value of B is set constant at 2. A AR value of 0.5 denotes an oblate spheroid,while a value of 1 represents a perfect sphere, and finally AR of 1.5 depicts a prolate ellipsoid. As anext step, in order to study the individual effect of the variable B, AR is kept constant at 1. When B41equals 2, the particle resembles a sphere; as B increases, the particle resembles a cube. To study theeffect of B values of 2, 4, and 8 are chosen. It should be noted that while it is possible to considerB values greater than 8 using the adopted implementation of superquadratics, doing so may resultin simulation convergence problems [97]. The figure 3.8 also illustrates additional convex shapesFigure 3.7: A typical superquadratic ellipsoid with marked half lengths along major axis c and minoraxes a and b.2 4B0.51.01.5AR8Figure 3.8: Superquadratic particle shapes with different aspect ratio AR and blockiness B. Highlightedyellow and green portions represent the ranges of AR and B studied.42which are not highlighted, that can be created using Equation 3.11. These parameters, namely ARand B, are related to well-known particle shape indicators such as roundness (R) and sphericity (S)in the following.3.8.2 Relation with roundness and sphericityIn general, the roundness and sphericity of a particle are calculated using a two-dimensional (2D)projection of the particle, proposed by Sloss [38]. This is a simplification of an actualthree-dimensional (3D) particle. Recently, Zheng et al. [131] has devised a computationalprocedure to calculate 3D roudness which can automatically determine the inscribed spherediameters. However, the process is more suited for complex particle shapes and hence a fairlysimple numerical algorithm proposed by Goldman [132] can be used to calculate the inscribedradius as follows:• Calculate the Gaussian curvature (KG) and mean curvature (KM) around the enclosed surfaceof a particle and choose the unique maximum value for a certain shape.KG =∇F×H(F)×∇FT|∇F |4 , M(F) =∂F∂x∂F∂y∂F∂ z , ∇F =∂ 2F∂x2∂ 2F∂x∂y∂ 2F∂x∂ z∂F∂y∂x∂ 2F∂y2∂ 2F∂y∂ z∂ 2F∂ z∂x∂ 2F∂ z∂y∂ 2F∂ z2 (3.12)KM =∇F×H(F)×∇FT −|∇F |2 Trace(H(F))2|∇F |2 (3.13)where, F(x,y,z) here is the implicit superquadratic function, and H is called the Hessianmatrix.• Principal curvatures(ki and k j)are thereby calculated from the values of KG and KM as:ki,k j = KM±√K2M−KG.43• Equivalent radius which is the reciprocal of curvature is calculated as req = 1/KM, whereaverage mean curvature for the two particles (i, j) is KM = (k1+ k2)/2.• With this process of calculating req, Wadell’s definition [37] of roundness can then beemployed to calculate R as:R =∑Ni=1∣∣reqi∣∣Nreqic(3.14)where the sum is taken over all the N corners of the projection, reqi is the radius of curvatureof the ith corner, and reqic is the radius of the largest circle that can be inscribed in the section.Particle sphericity is calculated in this study as the ratio of the radius of the maximum inscribedsphere (rins) and the radius of the minimum circumscribed sphere (rcir), a standard definition ingeotechnical engineering [76] as:S =rinsrcir. (3.15)The values of R and S evaluated from the above procedures are presented in Figure 3.9 toestablish links with AR and B. In Figure 3.9(a) it can be observed that at AR = 1, both the values ofR and S are = 1. As AR is varied to either greater than or less than 1, the values of R and S decrease.The decrease in S is found to be significantly higher as compared to R. The link between B, R, andS is illustrated in Figure 3.9(b). As B increases, both R and S are found to decrease. In comparisonto AR, which has a greater influence on S, B is found to have a stronger influence on R.The current research does not seek to replicate actual soil grains, but rather to capture shapecharacteristics such as angularity and elongation. These two characteristics, which can bequantified in terms of roundness and sphericity, are investigated here using particles withsimplified geometry but identical shape features to those found in real soil grains.440.5 1.0 1.5AR0.50.60.70.80.91.0RSB/2 R S(a)2 4 6 8B0.20.40.60.81.0AR R S(b)Figure 3.9: Relation between (a) aspect ratio AR, and (b) blockiness B, of superquadratics withroundness and sphericity calculated from the current analysis.3.9 Validation with laboratory testsCharacteristics of cyclic liquefaction response were successfully modeled using 3D-DEM as earlyas 1994 by Ng [50]. Thereafter, several DEM studies have been conducted that have modelledthe cyclic liquefaction response [10, 110, 133]. A constant-volume cyclic response from a simpleshear test conducted in the present analysis using an isotropic compressed sample with bi-periodiclateral boundaries is shown in Figure 3.10. This response is compared with a Figure 3.11 whichrepresents a typical mechanical response of an Ottawa sand sample subjected to cyclic undrained45hollow-cylinder loading [11]. Due to the fact that external stress is not applied on the principalplanes in both simple-shear and hollow cylinder mode, along with the control the stress in all sixdirections, at least under the current design of simple shear, these two test types can be argued tosimilar. It can be observed from the two cyclic responses that the DEM response is qualitativelycomparable with the experimental result, however, it appears to be comparatively stiffer. This maylikely occur due to the inability to model the uncertainty in sample preparation.Apart from assessing the qualitative agreement between the DEM and laboratory-based0 20 40 60 80 100σv [kPa]−30−20−100102030τ[kPa](a)−4 −2 0 2 4γ [%]−30−20−100102030τ[kPa](b)Figure 3.10: The typical (a) stress path, and (b) stress-strain curve, for a sample subjected to cyclicloading in a simple shear test (current study).0 20 40 60 80 100σv [kPa]−30−20−100102030τ[kPa](a)−6 −4 −2 0 2 4 6γ [%]−30−20−100102030τ[kPa](b)Figure 3.11: The typical (a) stress path, and (b) stress-strain curve, for a sample subjected to cyclicloading in a hollow-cyclinder laboratory test [11].4620 40 60 80 100NIL0.10.20.30.40.5CSRDr = 30%Dr = 50%(a)20 40 60 80 100NIL0.050.100.150.20CSRDr ≈ 28%Dr ≈ 58%(b)Figure 3.12: Qualitative comparisons between a typical CSR curve for Cu = 1.5 and σv = 100 kPa from(a) current study using DEM, and (b) an undrained cyclic hollow cylinder test by Wichtmann [12].stress-strain responses, the effects of two significant parameters, namely Dr and CSR, arepresented for the current study and qualitatively compared with the experimental study byWichtmann [12]. Figure 3.12(a) depicts the response of samples with Cu = 1.5 to two differentlevels of Dr, namely 30% and 50% compressed at p0 of 100 kPa using DEM. It is clear thatincreasing Dr and CSR results in an increase and decrease in the liquefaction resistance indicatorNIL, respectively. Similar observations can be observed from the undrained triaxial laboratorytests, in Figure 3.12(b), for the Karlsruhe sand with Cu = 1.5 at Dr of ≈ 28% and 58%. Thequalitative observations also match well with various other cyclic studies on different types of47sands [19, 33, 134]. Since quantitative comparison is extremely rare in DEM modeling due todifficulties in recreating the initial conditions, such characteristic observations have been criticalparameters to adhere to in DEM modeling over the last few decades. As a result, it can beconcluded that the analysis carried out in terms of relative density and cyclic stress ratio usingDEM in the present analysis is qualitatively consistent with the experimental analysis.3.10 SummaryThe methodology describes the detailed design of particle size distribution Cu and particle shape(AR, B). The value of Cu ranges from 1 to 3 when particles are limited to only spherical, whereasthe particle shape indicators AR is varied between 0.5 and 1.5, and B from 2 to 8 where Cu is keptfixed at 2. Isotropic samples are prepared compressed at a p0 of 100 kPa at Dr of 30% and 50%for studying the effect of particle size distribution and 20% and 50% for particle shape. Finally, asmall set of current cyclic test results are compared qualitatively to existing literature to show thatDEM follows the qualitative trends of experimental analysis.48Chapter 4: Effect of particle size distribution4.1 IntroductionThe preceding chapter had mainly provided the details of isotropic sample preparation in DEMfor different particle characteristics. This chapter discusses nuances in preparing the samples atdifferent relative densities, followed by the details of test design. As the samples with variousparticle size distributions are investigated under cyclic loading, their liquefaction resistance isexamined under different relative densities and loading intensities. The macroscopic findings arethen discussed in detail and their cause and effects are explored in terms of macro and microparameters based on the initial state and the state parameter that leads to establishingmacro-micro links. The evaluation of the micro structure is also used to arrive at the solid-fluidtransition that occurs during initial liquefaction (isostatic system) and is evaluated primarily basedon the changes of particle size distribution and relative density.4.2 Extreme void ratiosIn order to prepare samples with different values of coefficient of uniformity Cu at the same relativedensity, the minimum (emin) and maximum void ratios (emax) are determined. This is achieved byvarying the inter-particle tangential friction coefficient µ , at the first stage of isotropic compressionto 0.0 for the emin and 0.5 for the emax. The µ is varied in between and largely with the help oftrial and error; the intermediate void ratios and thereby relative densities are achieved for sampleswith five levels of Cu including 1, 1.2, 1.5, 2, and 3. This technique for generating dense and492 4 8 16D [mm]020406080100FV[%]Cu = 1.0Cu = 1.2Cu = 1.5Cu = 2.0Cu = 3.0Figure 4.1: Particle size gradation curves with five levels of Cu used in the study.sparse samples by varying µ has been used by a number of researchers in the field of DEM [103,104, 135]. The typical variation in the particle size distribution for all the samples with the finerpercentage by volume (FV) of particles passing given diameters in mm is presented in Figure 4.1.All of these samples are prepared with the optimized number of at least 30 particles per class witha total of ten classes. The samples all have the same mean particle size D50 of 5 mm, resemblingcoarse sand, and are compressed under a common mean effective stress p0 of 100 kPa.Since the sample preparation uses all the six boundaries as moving with identical servocontrol, to compress the samples, the contacts forming during at the end of compression (contactdistribution) is assumed to be unbiased in a specific direction. It is assumed that thishomogeneous or often termed isotropic distribution of contacts is prevalent along all thedirections in three dimensions irrespective of the particle shape. Finally, this isotropiccompression may not guarantee isotropic distribution of contact forces and can be a subject ofinvestigation for future DEM works.Figure 4.2 presents the maximum and minimum simulated void ratios obtained in this study.The values of emax are compared to the void ratios of random loosest packing (RLP) published inthe literature both in experimental database [76] and in 3D-DEM simulations [15, 136]. Similarly,the values of emin are compared to the void ratios of random densest packings (RDP) cited in501.0 1.5 2.0 2.5 3.0Cu0.40.50.60.70.8eemaxCho et al. (2006)Zheng and Hryciw (2016)Salerno et al. (2018)Present studyeminCho et al. (2006)Zamponi (2008)Salerno et al. (2018)Mutaburaka et al. (2019)Taylor et al. (2019)Present studyFigure 4.2: Maximum and minimum simulated void ratios from current study and comparisons withliterature.relevant experimental [76], and several 3D-DEM tests [76, 125, 128, 136, 137]. This could bebecause several studies used different sample preparation methods than the current study, owingprimarily to a general lack of standardized sample preparation in DEM. A decreasing overall trendin both emax and emin is observed with the increase of Cu. This largely demonstrates that higherCu can result in higher result in more densely packed specimens [128, 138]. The samples thusprepared at extreme void ratios demonstrate a close resemblance with similar data observed in pastresearch which have followed varied protocols of sample preparation, thus providing a strong basisfor further exploration using 3D-DEM.4.3 Test designTable 4.1 presents the 58 types of constant volume cyclic direct simple shear tests conducted forthe current study. Each simulation is assigned an individual ID, shown in the first column ofTable 4.1. To elaborate, Cu2 Dr50 CSR35 is a sample with a Cu value of 2, a Dr of 50%, and issubjected to a CSR value of 0.35. The second column displays the different Cu , the third columnpresents the total number of particles (Np) used in each sample, which varies between 8070 and24000, depending on the Cu level. It should be noted that in order to achieve the optimal classes51of particle sizes and particles per class, the total number of particles must significantly increaseas the particle size distribution increases; this is why the current study is limited to Cu of 3.0.For instance, increasing Cu from 3.0 to 4.0 increases the number of particles in the sample from24,000 to 110, 000, significantly increasing the computational time even with parallel processing.Higher values of Cu may well be targeted in the future with improved computational resources.The criteria for selecting a larger number of particles with a higher uniformity coefficient werediscussed in Chapter 3 and are therefore not repeated here. The samples are prepared in eitherof the two relative densities, Dr = 30% and 50% as shown in the fourth column, under a uniformtarget stress of p0 = 100 kPa. Following that, the samples are subjected to shearing to a rangeof CSRs as shown in the final column. For Dr = 30%, the CSRs are varied from 0.05 to 0.30,Table 4.1: Simulated constant volume cyclic simple shear tests.Test ID† Cu Np Dr (%) CSRCu10 Dr** CSR** 1 8070 30 0.10, 0.15, 0.20, 0.25, 0.3050 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60Cu12 Dr** CSR** 1.2 7960 30 0.15, 0.20, 0.25, 0.3050 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60Cu15 Dr** CSR** 1.5 9030 30 0.15, 0.20, 0.25, 0.3050 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60Cu20 Dr** CSR** 2 12600 30 0.15, 0.20, 0.25, 0.3050 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60Cu30 Dr** CSR** 3 24000 30 0.15, 0.20, 0.25, 0.3050 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60, 0.80, 1.00† Test ID with Cu12 Dr30 CSR20 refers to a sample with Cu = 1.2, Dr = 30%and CSR = 0.20.while for Dr = 50%, the range of CSRs are higher and varied between 0.25 and 0.60, for most ofthe samples. These simulations will be used to explore the effect of uniformity coefficient on thecyclic liquefaction response in the sequel, from both macro- and microscopic views, respectively.4.4 Macroscopic responseThe average stress tensor can be used to characterize the overall mechanical response of anassembly under constant volume cyclic shearing. Over a given volume V , the stress tensor σ is52calculated as an average of inter-particle forces and branch vectors:σ =1V ∑i∈Ncvl i⊗ f i, (4.1)where the branch vector l i connects the centers of two particles, f i is the contact force, ⊗ refer tothe tensor dyadic product, and the summation includes all the contacts Ncv in the selected volumeV . The shear stress τ , the vertical stress σv, and the mean effective stress p can also be calculatedfrom the stress tensor in the simple shear test, as τ = σzy, σv = σzz, and p = tr(σ )/3, respectively.The shear strain γ is also an important part of the mechanical response of the assembly whichcan be evaluated as:γ =ywh, (4.2)where, yw refers to the horizontal displacement of the top plate along y direction.While the pore water is not explicitly simulated in this study, the derived excess pore pressurein the constant volume test can still be calculated as the reduction in vertical stress relative to itsinitial value σvi as:∆u = σv−σvi. (4.3)The derived pore pressure ratio (ru) can be calculated by normalising ∆u with σvi:ru =∆uσvi= 1− σvσvi. (4.4)To present the evolution of a quantity such as ru, a “fractional cycle number” N′ is used to replace53the current running time t by interpolation between two successive cycles:N′ = N+t− tNT (N)(4.5)where N is the cycle number, tN is the starting time for the Nth cycle, and T (N) is the duration ofcycle N. The value of N′ coincides with N at t = tN . To avoid confusion in the dissertation, N isused to represent fractional cycle number defined by N′.0 20 40 60 80 100[kPa]−40−2002040[kPa]v(a)−2 −1 0 1 2[%]−40−2002040[kPa](b)0 10 20 30 40 50N0.00.20.40.60.81.0r uPre-liquefactionPost-liquefactionNIL(c)0 10 20 30 40 50N−2−1012γ[%](d)Figure 4.3: Macroscopic shear response of a sample with Cu = 2, Dr = 50% and CSR = 0.35 where: (a)stress path, (b) stress-strain curve, (c) derived pore pressure evolution, and (d) shear strain development.544.5 Cyclic responseTo illustrate the typical macroscopic behavior of granular materials under constant-volume cyclicshearing, the simulation Cu20 Dr50 CSR35 is presented in Figure 4.3 in a four-way plot,including stress path, evolution of derived pore pressure ratio, stress-strain loops, and shear straindevelopment. The progression of the stress path in terms of shear stress τ and effective verticalstress σv is shown in Figure 4.3(a). The beginning of the stress path occurs at σv = 100 kPa andτ = 0 kPa i.e., without any shear bias or asymmetry in shearing. As repeated cycles of shearingare applied, the stress path oscillates and leads to a decrease in the σv as evident from the leftwardmovement of the stress path. This irreversible loss of σv coincides with the increase in the level ofshear strain γ generation as shown in Figure 4.3(b). Eventually, the decrement in σv reaches avery low value corresponding to the derived pore-pressure ratio ru reaches a value equal to 0.95 asevidenced from 4.3(c) when increment in γ as referred in Figure 4.3(d) increases steadily into anunbounded nature with number of shearing cycles. In this study, the attainment of ru = 0.95 isdefined as the condition of initial liquefaction of a granular sample. It is noticeable that at initialliquefaction both τ and σv are very close to zero with the number of complete cycles of shearingto reach liquefaction (NIL) is counted as one of the prime indicators of cyclic liquefactionresistance. With further shearing, the sample demonstrates shearing resistance and follows a fixedstress path in the shape of a butterfly [139]. This phenomenon of restoring shearing resistanceafter initial liquefaction is known as cyclic mobility. Overall, the findings of the 4-way plot matchthose of the experimental undrained cyclic simple shear test on sands [23, 140, 141]. This pointdivides the loading process into pre- and post-liquefaction stages, shown in blue and red colors inFigure 4.3, respectively. The corresponding number of cycles is denoted as NIL, a measure forcyclic strength or cyclic liquefaction resistance. As shown in Figure 4.3(b), cyclic shearing causesa slow accumulation of shear strain amplitude in each cycle. As the loading progresses, the shearstrain amplitude begins to accumulate at a faster rate. It increases unboundedly until the test isterminated in the post-liquefaction stage, as shown in Figure 4.3(d). Overall, the findings of the550.0 0.2 0.4 0.6 0.8 1.0N/NIL0.00.20.40.60.81.0r uDr = 30%CSR = 0.15Cu = 1.0Cu = 1.2Cu = 1.5Cu = 2.0Cu = 3.0(a)0.0 0.2 0.4 0.6 0.8 1.0N/NIL0.00.20.40.60.81.0r uDr = 50%CSR = 0.35(b)Figure 4.4: Evolution of excess pore pressure ratio ru given at the beginning each cycle (τ ≃ 0) forsamples with different levels of Cu at: (a) Dr = 30%, at τ ≃ 0, and (b) Dr = 50%.4-way plot match those of the experimental undrained cyclic simple shear test on sands[23, 140, 141].4.5.1 Evolution of excess pore pressureGiven the rather qualitative description of the cyclic behavior of granular materials, a naturalquestion which arises is how Cu affects the evolution pattern of ru. Figure 4.4 presents theaccumulation of ru from ten simulations chosen in Table 4.1, covering both relative densities,against the normalized number of cycles by the corresponding NIL. To avoid cycle fluctuations,56only the values of ru accumulated at the start of each new cycle of loading when τ ≃ 0 for the firsttime are presented. Wichtmann [12] also presented similar excess pore-pressure evolutionfollowing each cycle of loading. The simulations with NIL ≥ 10 are chosen to cover sufficientcycles in the pre-liquefaction stage allowing for a more smooth continuous trend rather than thedrastic jump caused by considering high CSRs that result in fewer cycles to liquefaction.Surprisingly, for both relative densities of 30% and 50%, Cu does not have noticeable effects onthe evolution trend of ru.4.5.2 Cyclic liquefaction resistanceLiquefaction resistance is expressed in terms of the number of complete shearing cycles (NIL)required to bring the ru for the first time to a value of 0.95. The liquefaction resistance is assessedbased on 58 cyclic tests (refer to Table 4.1), that differ by Cu, Dr, or CSR, and are presented inthe form of typical CSR vs NIL curves as shown in Figure 4.5. For samples prepared at a Dr =30%, given a certain CSR the tests where Cu = 1.5 show the the maximum value of NIL or inother words possess the highest liquefaction resistance. Samples with Cu values greater or lessthan 1.5 that were considered in this study demonstrated a decrease in its value. The liquefactionresistance curves does not intersect at this relative density implying that the above observation ismade regardless of the level of CSR. However, for the medium dense case with Dr = 50%, the effectof Cu on the liquefaction resistance curves becomes complex, as they are found to be dependenton the CSR level. In the range of CSR considered for medium dense samples, i.e., 0.30-0.60, it isfound that among all samples sheared at the same CSR between 0.40-0.60, one with Cu = 3 showsthe maximum NIL and this trend is amplified as the magnitude increases from 0.4 to 0.6. However,between CSRs of 0.30 and 0.35, the higher NIL for Cu = 3 is replaced gradually by curves for lowerCu values of 1.0, 1.2, and 1.5. The effect of Cu of 1, 1.2, 1.5, and 2 on the NIL is found to benoticeable at low CSRs while their cyclic resistance curves tend to converge beyond CSR of 0.4.This last observation, together with the figure at Dr = 50%, may imply that for dense samples,57100 101 102NIL0.10.20.30.40.50.60.7CSR30 50= 1.0= 1.2= 1.5= 2.0CuCuCuCuCu = 3.0Dr (%)Figure 4.5: Cyclic liquefaction resistance curves for samples with different values of Cu and Dr. Thehollow markers represent results for Dr = 30%, and the solid markers for Dr = 50%.the vast majority of the Cu values considered in this study have a marginal effect on liquefactionresistance at CSR values greater than 0.5.The effect of Cu on liquefaction resistance is further illustrated based on certain macro andmicro-parameters after a detailed analysis of the CSR-NIL curve. During the liquefactionassessment procedure, Seed et al. [142] define the cyclic resistance ratio CRR as a simplifiedmacro-level quantity that can be used in lieu of the CSR-NIL plot. As stated by Robertson andWride [143], the CRR is usually evaluated at 15 cycles of uniform loading to represent anequivalent earthquake magnitude of 7.5, denoted by CRR15. Such extractions from the CSR-NILcurve using CRR can be referred as cyclic strength [144]. The CRR15 is plotted against variouslevels of Dr and Cu in Figure 4.6(a). For simulations with Dr = 30%, CRR15 follows a “A” shapedrelation within the interval of Cu, reaching its peak at Cu = 1.5 before plateauing at Cu = 2. Forhigh Dr = 50%, the CRR15 trend is “V” shaped with the minimum occurring at Cu = 2. In Figure4.6(b), NIL corresponding to a common CSR = 0.30 is used as the liquefaction resistanceindicator, to study the effect of Cu. Similar trends in the liquefaction resistance are observed atdifferent Dr, as found using CRR15. Because the results obtained at Dr = 50% are CSR581.0 1.5 2.0 2.5 3.0Cu0.10.20.30.40.50.60.7CRR1530 50= 1.0= 1.2= 1.5= 2.0CuCuCuCuCu = 3.0Dr (%)(a)1.0 1.5 2.0 2.5 3.0Cu100101102NIL(b)Figure 4.6: Variation of cyclic liquefaction resistance with Cu, where (a) liquefaction resistance isrepresented as CRR15; and (b) as NIL. The hollow markers represent results for Dr = 30%, and the solidmarkers for Dr = 50%.dependent, the trend for Figure 4.6(b) may slightly differ if a CSR other than 0.30 is used, but thebasic “V” shape will remain. In summary, Figure 4.6 demonstrates a significant effect of Cu onliquefaction resistance, with the characteristics of the influence varying according to the level ofDr at which the samples are compared. The underlying assumption in any particulate mechanicsis that the source of any macrobehavior can be traced back to the microlevel. The macroscopic59observations will then be analyzed from a microscopic perspective of the particle system in itsinitial or pre-shearing state. For example, mechanical coordination number is used to representthe contacts, and the value β is used to quantify the strong contact normal forces among others,and compared at the initial state of the samples with their liquefaction resistance indicator CRR15.4.6 Relation between CRR15 and initial stateThe cyclic liquefaction resistance of samples with different uniformity coefficients is found to besignificantly influenced, and the analysis of the microstate based on the contacts, contact forces,among others, can provide a better understanding of the behavior of the granular materials.4.6.1 Relation with z0To gain a better understanding of the macroscopic response, attempts are made to trace theirparticle-scale origins at the initial state. The average number of contacts in an assembly devoid ofnon-load bearing particles known as the mechanical coordination number (zm) [145] is used as amicro-mechanical parameter to explain cyclic liquefaction resistance. To transmit force through a1.0 1.5 2.0 2.5 3.0Cu3.23.74.24.75.2z 030 50= 1.0= 1.2= 1.5= 2.0CuCuCuCuCu = 3.0Dr (%)Figure 4.7: Variation of z0 for samples with different Cu. Those packed at Dr = 30% are shown withhollow markers, and for Dr = 50% with solid markers.60particle (load-bearing) in a pseudostatic system, a collinear force transmission path with at leasttwo inter-particle contacts must exist. Hence, the non-load bearing particles in an assembly ofgranular particles either have zero (floaters) or just one contact (rattlers). The value of zm is givenas:zm =2Ncv−N1pNp−N0p −N1p(4.6)where, Ncv denotes the total number of contacts in the assembly of volume V ; N0p and N1p are thenumber of particles with 0 and 1 contact, respectively. Following an evaluation of all possible typesof coordination number, it is determined that only zm at the initial state (z0) adequately accountsfor cyclic liquefaction resistance. In light of Figure 4.6, the variation of z0 is plotted against therange of Cu for the two levels of Dr in Figure 4.7. The current study clearly shows that z0, in itsinitial state, exhibits a similar “A” and “V” shaped trend with variation in Cu for Dr = 30% and50%, respectively, as CRR15 and NIL ascribe for the same range of Cu. In light of the similarity3.8 4.2 4.6 5.0z00.10.20.30.40.50.60.7CRR1530 50= 1.0= 1.2= 1.5= 2.0CuCuCuCuCu = 3.0Dr (%)Figure 4.8: Correlation between z0 with CRR15 for samples with different Cu and Dr. The exponentialbest fit is shown as the blue dotted line.between the variation of z0 and CRR15, with Cu, they are plotted against each other for different Dras shown in Figure 4.8. It is observed that CRR15 and z0 show a promising monotonic correlationwell described using an exponential curve, shown here with the dotted line. A similar plot hasalso been utilized in the recent work by Morimoto et al. [93] for a fixed particle size distribution.61This suggests that the initial mechanical coordination number which envisage that the participatingnetwork” in a granular assembly is a good predictor of liquefaction resistance in isotropic granularsamples with varying uniformity coefficients and relative densities. As discussed in the followingsection, it is also worth noting the force transmission network and force chains of a sample thatmay be used to represent this participating network in order to study the liquefaction resistance.4.6.2 Relations with inter-particle force and contact distributionsThe contact force network can provide insight into the microstructure of a particulate assembly[47]. The current research evaluates the inter-particle normal contact forces at the initial state andattempts to link certain characteristics of this force distribution with cyclic liquefaction resistance.The probability distribution function [P( fn)] of contact normal forces is a crucial and frequentlyused parameter to analyze the contact force distributions [146, 147]. At the initial state, theprobability distributions of normal contact contact forces fn, normalized by the mean normalcontact forces ⟨ fn⟩, are compared between samples with varying uniformity coefficients andfurther investigated to determine a connection between the variations in the distributions andCRR15. Figures 4.9(a) and 4.9(b) illustrate the distribution of the normalized contact normalforces for samples with different Cu and Dr. In general, contact forces lower than the mean valuefn/⟨ fn⟩ < 1 are referred to as weak forces, which decay in a power-law relation, whereas theforces larger than the mean, are the strong contact forces which decay exponentially [148]. Thedistribution can be written in the following form:P( fn) = k(fn/⟨ fn⟩)−α, fn/⟨ fn⟩< 1keβ (1− fn/⟨ fn⟩), fn/⟨ fn⟩> 1(4.7)62where α is the slope of the power-law fit for weak forces, β is the slope of the exponential best fitfor strong forces, and k is known as the normalization factor expressed as1k=11−α +1β(4.8)Additionally, since the mean force ⟨ fn⟩ is chosen as the point of intersection between the two partsof the distribution, the relation between the exponents can be obtained as:β 2 = (1−α)(2−α) (4.9)Each sample may be described as comprising a heterogeneous force network since it containsboth strong and weak contact normal forces (non-uniform forces). The probability of the strongcontact normal force is found to decay with the increase ⟨ fn⟩; comparatively, the decay occursmuch slower as the samples with higher Cu are considered. Moreover, this decay is found to beenhanced with the increase of Dr from 30% to 50%. Considering that all samples are isotropicallycompressed at the same level of effective stress, the wider distributions of the contact normalforces with Cu increasing from 1.0 to 3.0 and with decrease in Dr from 50% to 30% indicatetowards larger forces present in the network at the expense of lower average and weak forces.This invariably points to the fact that there are fewer contact channels for force transmission andlarger fraction of particles existing in the network as non-force transmitting i.e., floaters or rattlers.Larger number of floaters or rattlers have also been linked with higher Cu by a recent study byZhang et al. [149]. Because forces are transmitted through particles that have two or morecontacts, the relative scarcity of force networks assemblies may also indicate a lower z0 value. Thedecay in the strong normal force ( fn/⟨ fn⟩> 1) distributions is further investigated by quantifyingits slope in the semi-log space with the exponent β from Equation 4.7. The lesser distributedthe strong force network is, the greater the value of β and vice-versa. In Figure 4.10, the βvalue is plotted along with CRR15 for different levels of Cu and Dr. Increases in β have been630 4 8 12 16fn/ 〈fn〉10−310−210−1100P(f n)Dr = 30%1Cu = 1.0Cu = 1.2Cu = 1.5Cu = 2.0Cu = 3.0(a)0 4 8 12 16fn/ 〈fn〉10−310−210−1100P(f n)Dr = 50%1Cu = 1.0Cu = 1.2Cu = 1.5Cu = 2.0Cu = 3.0(b)Figure 4.9: Probability distribution functions of normalized contact normal forces fn/⟨ fn⟩ for sampleswith different Cu at (a) low Dr with wide force distribution, and (b) medium dense Dr with narrow forcedistribution.found to correlate with increased cyclic liquefaction resistance of an assembly, as indicated byan exponential curve fitting except for the lone outlier representing less than 10% of the datapresented in this figure. According to the limit set by Nau [150], this is a good regression valuebetween two non-stationary variables (time-independent) which are based on well-designed tests.A comparatively higher β indicates faster decay of the strong contact normal force, thus leadingto a relatively homogeneous distribution of the contact forces. This type of greater homogeneity0.0 0.5 1.0 1.5 2.00.00.40.81.2CRR15Dr = 30%rD = 50%Figure 4.10: The relation between the rate of decay β of the strong force network fn/⟨ fn⟩ > 1 withCRR15.64in force distribution implies there are more contact pathways for force transfer arising from highercontacts per particle; thus providing a reasonable explanation for the positive correlation between βand CRR15. The linkages between homogeneity in the distribution of the contact normal forces andhigher contact pathways in an assembly have also been discussed by some researchers [149, 151],which led to higher shear resistance. Moreover, larger contact forces in a force network may alsoarise from the capacity of larger diameter particles to hold higher magnitude contact forces, whichis prevalent in samples with higher Cu values as pointed out by some researchers in 2D [127], and3D-DEM [128] simulations.The lone outlier could have been caused by the competing effect from the excessive loss ofload bearing particle fraction with polydispersity, and by the increased number of several smallerparticles filled in between few large particles in the system. Additional research using highlypolydisperse samples is necessary to determine whether the relationship between CRR and β holdstrue in a more general sense or is an artifact of the small data pool considered in this problem.Overall, it can be concluded that increased homogeneity of the strong contact normal forces mayincrease cyclic liquefaction resistance of granular assemblies, regardless of the limited type ofparticle size distribution and relative density used in this study.The networks of contact normal forces are subsequently presented in the visual form of forcechains. The force chains are visualized using colored tubes of varying sizes depending on themagnitude of the contact normal force. These tubes form a network by connecting the centroidsof particles in contact. 3D force chains are used to represent the differences in force distributionsthat arise from the changes in particle size distribution. Force chains depicting the contact | fn|values for Cu = 1, 2 and 3 are shown in Figure 4.11 for Dr = 30%, and for Dr = 50% in Figure4.12. The insets are used to zoom in a portion of the force network for analysis. The force networkfor Cu = 1 as shown in Figure 4.11(a) portrays substantial tubes with similar diameters that arepredominantly blue in color. This clearly demonstrates the visible homogeneity in contact forcesin assemblies with identical sized particles. In comparison, the force network for Cu = 2 shownin Figure 4.11(b) contain thicker and red force tubes along with thinner blue force tubes, and this65(N)(a) (b)(c)Figure 4.11: Snapshots of normal force chains of the samples at the initial state at Dr = 30% for: (a) Cu =1, (b) Cu = 2, and (c) Cu = 3. The zoomed insets are provided for the sake of clarity.trend is amplified further for Cu = 3 in Figure 4.11(c). These observations are in agreement withthe probability distributions of contact normal forces as shown in Figure 4.9(a), that also shows anincreased heterogeneity of the contact normal force network with the increase in the particle sizedistribution. Due to the increased heterogeneity, there are presence of very thin, blue force chains66(N)(a) (b)(c)Figure 4.12: Snapshots of normal force chains of the samples at the initial state at Dr = 50% for: (a) Cu =1, (b) Cu = 2, and (c) Cu = 3. The zoomed-in insets are included for clarity.placed abundantly in between the stronger force chains. For Dr = 50%, a similar observation canbe made (refer to Figure 4.12), except these are marked with fewer visible red chains compared tothe sample at Dr = 30%, consistent with Figure 4.9(b). The findings from the current examination670 2 4 6 8 10 12c0.00.10.20.30.4P cDr = 30%Cu = 1.0Cu = 1.2Cu = 1.5Cu = 2.0Cu = 3.0(a)0 2 4 6 8 10 12c0.00.10.20.30.4P cDr = 50%Cu = 1.0Cu = 1.2Cu = 1.5Cu = 2.0Cu = 3.0(b)Figure 4.13: The fraction Pc of particles with c contacts for different values of Cu at (a) low Dr withsignificant presence of non-load bearing contacts in higher Cu (b) medium-dense Dr with suppressednon-load bearing contacts.of the force networks illustrate the presence of the arching effect [152] in the granular structure,which is found to increase with increasing Cu, and is more prevalent in looser samples than indenser samples.Because the initial state provides a reasonable context into the cyclic shearing behavior, asevident from the microscopic analysis, both the participating and non-participating networks areexamined in terms of particle fraction with a specific number of contact numbers for furtherinsight into the microstructure. For the entire assembly, the contact network is analyzed using theproportion of particles with c contacts, denoted by Pc. Figure 4.13 illustrates the correlation68between Pc and c for various values of Cu and Dr. As illustrated in Figure 4.13(a), there is anon-zero fraction (about 5%) of particles with c = 0 when Cu = 1.0, which increases significantlyto around 35% as Cu increases to 3.0. This observation holds true for most particles with c = 1,though the fraction of such particles is relatively small. Thus, it is evident that as particle sizedistribution increases, a greater proportion of the granular assembly is excluded from the forcebearing network. These increasing numbers of floaters and rattlers also indicate the increasingpresence of arching effects, which corroborate from the previous understanding based on theforce chains, providing sufficient insights into the more heterogeneity in the microstructure. Theparticles with equal or more than four contacts, because of their stability form the arches. The Pcvalue of these particles decreases with increase in Cu as evident from the figure. For example,only 15% of particles in an assembly with Cu = 3.0 have been found with c = 4, compared toapproximately 25% for the monodispersed sample. Due to the fewer particles used to form thearches, the force carried along each particle-particle contact is increased, as are the sizes of thesearches formed around more floaters and rattlers. As shown in Figure 4.13(b), samples with Dr =50% have a lower number of floaters and rattlers in the assembly, but the increase in Cu has stillresulted in an increase of these particles. It can therefore be concluded that increased archingeffects at the initial state results in burgeoning heterogeneity of contact forces and this in turn isfound to reduce the liquefaction resistance. These observations simply convey that if a granularassembly is stable due to a few large arches rather than a dense network of smaller ones, it is moreprone to cyclic liquefaction.4.7 Relation between CRR15 and state parametersA state parameter is generally considered to be a more accurate to describe the “state” of a granularmaterial than the void ratio or relative density. This is because it compares the difference of theinitial state with the critical state of a material at the same confinement stress. The state parameterψe is equal to e0− ec, where e0 is the initial void ratio and ec is the critical state void ratio. This is690 10 20 30 40 50γ [%]010203040506070τ[kPa](a) 0 10 20 30 40 50[%]0.540.560.580.600.620.64eDr = 30%Dr = 50% e0 ec = e0 - ec(b)0 10 20 30 40 50[%]2.53.03.54.04.55.0z mzm0zmc = zm0 - zmc(c)Figure 4.14: Drained monotonic simple shear test on samples with different Cu at Dr = 30% and 50%: (a)stress-strain response, (b) evolution of void ratio, and (c) evolution of mechanical coordination number.claimed [14, 134, 153] to capture dilatancy more accurately than the majority of other parameters,which is thought to be the source of all shear resistance, including the cyclic strength. However,ψe cannot capture the effect of fabric in a unique relationship with CRR. Recently, Gu et al. [154]proposed zm as an alternative to e for its usefulness as a state parameter ψz which do not considerthe non-load bearing particles. In addition, it was claimed that ψz can capture the effect arisingfrom differences in fabric among samples into its unique relationship with CRR for samples withsimilar particle size distribution. Nevertheless, there is lack of evidence in literature about whethera unique relation between state parameter and the cyclic strength (expressed here as CRR15) canbe presented, which can capture the effect of particle size distribution.70In the current study, each sample is subjected to drained, monotonic, simple shear simulationswith the element boundaries adjusted to maintain a constant mean effective stress (p0). These−0.20 −0.15 −0.10 −0.05 0.00ψe0.00.20.40.60.81.0CRR15Dr = 30%Dr = 50%(a)0 1 2 3 40.00.20.40.60.81.0CRR15Dr = 30%Dr = 50%30 50= 1.0= 1.2= 1.5= 2.0CuCuCuCuCu = 3.0Dr (%)(b)Figure 4.15: Cyclic liquefaction resistance against state parameters for samples with different particlesize gradations and relative density derived from: (a) void ratio, and (b) mechanical coordinationnumber. The hollow markers represent Dr = 30%, and solid markers represent Dr = 50%.constant p0 tests enabled the calculation of state parameters for each sample from a single testwithout requiring a large number of simulations to obtain the typical e− p0 relationship. Drainedresponse of a sample with Cu = 2 is presented in Figures 4.14(a) to 4.14(c) at Dr = 30% and 50%.The values of τ , e and zm are found to evolve with γ until they reach a level of saturation. It canbe observed that both of the samples show little to high dilation with gain in e and subsequent loss71in zm due to shearing, and ultimately reaching an unique value at the critical state viz. ec and zmc,respectively. These responses are recorded for samples with other Cu values and not repeated herefor brevity. As a next step ψe and ψz , where ψz = ψz0-ψzc , are correlated with CRR15 in Figure4.15. It can be clearly seen from Figure 4.15(a) that CRR15 overall reduces non-linearly with ψefollowing an exponential correlation [155]. This implies that the effect of particle size distributionand relative density both on the cyclic liquefaction resistance can be captured by this exponentialrelation. This good correlation while using ψe as a state parameter is understandable in the currentstudy given that all samples are prepared under identical isotropic sample preparation. Figure4.15(b) demonstrates that a similar exponential but a positive relation exists between CRR15 andψz . It signifies that higher the contact loss for a sample during the stage of drained shearing, whichresults in increase in ψz , greater its liquefaction resistance. Though, higher dilation during drainedshearing captured by ψe may appear to be analogous with higher contact loss ψz in the micro scale,but loss of contacts, has been shown to be better correlated with CRR15. After examining the initialstate in terms of contact normal force and contact distribution, the contact evolution and a criticalparameter at liquefaction are examined in the following.4.8 Evolution of the non-load bearing particlesAt the initial state, the floaters and rattlers (non-load bearing particles) in a sample is found toincrease significantly with increase in Cu as studied from the combination of the Pc − c,P( fn)− fn/⟨ fn⟩ and the force chain plots. It was observed that increased presence of suchparticles may help indicate other fabric details such as heterogeneity in the force distribution andheightened presence of the arches which form the load bearing fabric. Here from, the fraction ofthese non-load-bearing particles is denoted as ffr and its evolution with N/NIL is studied till initialliquefaction in Figure 4.16. In Figure 4.16(a), it can be observed that ffr increases for all Cu butonly to a small extent, indicating the fabric as described at the initial state remains valid throughout the pre-liquefaction stage. As illustrated in Figure 4.16(b), this observation also holds true for720.0 0.2 0.4 0.6 0.8 1.0N/NIL0.00.20.40.60.81.01.2f frDr = 30%CSR = 0.15Cu = 1.0Cu = 1.2Cu = 1.5Cu = 2.0Cu = 3.0(a)0.0 0.2 0.4 0.6 0.8 1.0N/NIL0.00.20.40.60.81.01.2f frDr = 50%CSR = 0.35(b)Figure 4.16: Evolution of ffr with N/NIL at the end of each loading cycles (τ ≃ 0) at different levels of Dr.the dense sample with Dr = 50%. Meanwhile, the relative extreme positions of the evolutioncurves for different Cu can be observed in a comparatively wider band for Dr = 30 % as comparedto Dr = 50%. This indicates that the variation in fabric composition due to the differences in Cu isless prevalent with the increase in Dr.Instead, it is more interesting to consider the pivotal point marked by the transition from soild-like or hyperstatic state to a fluid-like or hypostatic state. Huang et al. [156] suggested that thistransition point, also known as the isostatic condition, can be described by the coordination number,and its value “may” depend on factors like particle shape, particle size gradation among others.73Previously, while some researchers [114, 139] emphasized the importance of studying the solid-to-fluid transition. While some researchers [156, 157] have attempted to calculate the value of zm atthe isostatic condition for spherical particle assemblies. Under such condition, when inter-particlefriction µ is infinite, it was determined that zm should be equal and not exceeding four, and lessthan four for finite values of µ .4.9 Examination of the isostatic stageThe state of initial liquefaction i.e., the first occurrence of isostaticity is examined hereafter bymeans of identifying the critical coordination number ziso. With the advancement of cyclic loading,σv reaches to near zero levels, where an abrupt drop in the corresponding zm is observed as shownin Figure 4.17(a). This decrease in zm occurs as a result of the transition from a solid-like to a fluid-like state, and its corresponding value is defined as the critical coordination number ziso, also shown0 20 40 60 80 100σ v [kPa]2345z mCu = 2Dr = 50%CSR = 0.35Liquefied zoneisoz = 3.20(a)0 5 10 15 20 25 30 35 40N2345z m z = 3.20ILN = 30(ru = 0.95)ru < 0.95ru > 0.95iso(b)Figure 4.17: The evaluation of ziso at initial liquefaction based on the loss of σv: (a) The degradation ofzm with loss of σv leading to the vertical liquefied zone at σv ≃ 0, and (b) comparisons of N at liquefactioncorresponding to ziso and ru = 0.95.in Figure 4.17(b). It should be noted that the coordination number in a fluid-like state during initialliquefaction and subsequent post cycles, follows a certain vertical band designated as the liquefiedzone. Similar drop in coordination is also observed by Huang et al. [156] while modeling cyclicliquefaction using triaxial tests using 3D-DEM. This observation is consistent with that observed74Figure 4.18: The effects of particle size distribution and relative density on ziso.at other Cu levels. These provide clear evidence that both methods are consistent in identifying thetransition of solid-like to the fluid like stage in the samples considered.As a next step, the values of ziso are plotted for different Cu and Dr in Figure 4.18. It can beobserved that ziso decreases with increase in Cu and decrease in Dr. The rate of decrease of zisois more prevalent at the lower Dr considered in this study. The reduction in ziso values indicatethat an assembly require lesser participatory contacts (two and above) to achieve a isostatic stage,which is also consistent with the observation at the initial state, where lesser particles are foundto provide the contact pathways, thus producing larger contact forces with increase in particlesize distribution. Increased floater and rattler fractions may also account for the lower ziso, andoccupying the otherwise void spaces as they increase with higher uniformity coefficients and lowerrelative densities. It should be noted that the effect of CSR on ziso has been determined to beinsignificant and is thus not presented in the study.754.10 ConclusionsThe primary objective of this chapter was to isolate and then systematically study the effect ofparticle size distribution on the resistance of isotropic samples to cyclic liquefaction when shearedunder constant volume cyclic simple loading. Furthermore, using micro-level parameters, thisstudy investigated the pre-shearing micro-structure, as well as state parameters in order todetermine their possible link with the macro-level response of cyclic liquefaction resistance.Samples with Cu ranging from 1 to 3 were prepared isotropically at Dr = 30% and 50% andsubjected to varying levels of CSRs, with cyclic liquefaction resistance determined using twodifferent indicators, CRR15 and NIL. At Dr = 30%, increase in particle size distribution resulted inan “A” shaped variation in liquefaction resistance i.e., an initial increase followed by decrement.Whereas, when Dr was increased to 50%, the observation exhibited an inverse relationship,closely resembling a “V” shaped trend. In this case, as opposed to the less dense case, theresponse was found to be CSR dependent. These observations were found consistent with boththe liquefaction indicators.The macroscopic observations were thereby judged by exploring the initial state and thecritical state of the samples. It was observed that the initial state variable mechanical coordinationnumber zm is a meaningful predictor of CRR15 irrespective of particle size distribution. It wasinteresting to note that zm does not monotonically increase with increase in Cu as generallypresumed, and follows a complex variation directly influencing the cyclic strength. Moreover,when zm was substituted with the state parameter ψz , it significantly improved the uniqueexponential relation with CRR15 for different Cu and Dr values.The fabric was also explored beyond zm into the contact force distributions at the initial state. Itwas observed that increase in Cu leads to increased heterogeneity in the force distributions arisingfrom the increased presence of small floaters and rattler fractions. This resulted in a smaller numberof large particles carrying the entire contact forces in the form of arches, causing reduction incyclic strength. Furthermore, these differences in force and contact networks with increase in Cu76were found to become less prevalent with increase in sample density. These observations aboutthe contact and force networks were found to be true at the initial state as well as for the entirepre-liquefaction stage.Finally, zm was further used to identify the transition of the assemblies from a solid like state(shear resistant) into a fluid like state with negligible shear resistance. At this point, the values ofzm referred to as ziso, appear to be unaffected by the CSR level, but they decrease in samples withhigher uniformity coefficients and lower relative density. This decrease in the ziso requirement forisostaticity is believed to be due to ability of polydispersed samples to form load bearing systemswith a smaller number of large particles, as indicated by the increased number of floaters as Cuincreases with smaller particles. The reduction may also arise from the intermittent supportprovided, as well as the increased space occupied by floaters and rattlers.These macroscopic observations made with respect to cyclic strength are quite non-intuitivewhen the effect of size distribution is completely isolated from other factors like shape, surfacefeatures, fabric, among others. Such non-intuitive macroscopic observations are found to be wellexplained with respect to the micromechanics and state parameter based on the initial state of theisotropic samples. Future studies should look at particle shape for different particle sizedistributions and loading types among others to investigate if such macro-micro linkages can bedeveloped in view of the state parameters of void ratio and coordination number.77Chapter 5: Effect of particle shape5.1 IntroductionThis chapter discusses the design and results of DEM tests performed to determine the effect ofparticle shape on the cyclic liquefaction resistance of granular assemblies. The methodologyfollowed in this chapter for isotropic sample preparation and subsequent shearing was describedin Chapter 3, and is similar to the one used to study particle size distribution in the previouschapter. Additionally, this chapter examines the effect of particle shape on the differences in themicrostructure at the initial state in terms of contacts and contact forces, as well as the stateparameter based on both initial and critical states, and thereby investigates a link between themicrostructure, the state parameter, and the macroscopic response. Finally, contact-basedanalyses are performed to ascertain the critical number of contacts that could define the solid-fluidstate of transition during cyclic liquefaction.5.2 Extreme void ratiosAs samples with different particle shapes are created, preparing them at the same void ratio is notcorrect because their densest and loosest state might vary significantly [76]. Hence, in the currentstudy relative density Dr is used and therefore simulated maximum emax and minimum emin voidratios are calculated for different aspect ratio AR and Blockiness B values, where the value ofcoefficient of uniformity Cu is kept constant at 2.0, and are presented in Figure 5.1. It should benoted that the procedures for calculating emax and emin, as well as the design of the particle class780.5 1.0 1.5AR0.30.40.50.60.7eemaxemin(a)2 4 6 8B0.30.40.50.60.7e(b)Figure 5.1: The maximum and minimum simulated void ratios for samples with (a) different aspectratios AR and same blockiness B, and (b) different B and same AR.for Cu = 2.0 and the definitions of AR and B, are covered in Chapter 3 and are thus not repeatedhere. In this study, the extreme void ratios are calculated for samples with AR of 0.5, 1, and 1.5when B is kept constant at 2.0, and at B = 2, 4, 6, and 8 when AR is kept constant at 1. Thesevalues are presented in Figure 5.1. As illustrated in Figure 5.1(a), the values of emax and emin are0.65 and 0.48 occurring at AR = 1. As AR is increased to 1.5, the corresponding values decreaseto 0.58 and 0.40. Similar decrease in emax and emin can also be observed for a reduced AR of 0.5.Some researchers [44, 158] have also noted a decrease in emax and emin in samples with ellipsoidalparticles, i.e., AR ̸= 1 as compared to spherical particles, i.e., AR = 1. The values of emax and eminare shown to decrease as samples with a higher B value are prepared, as shown in Figure 5.1(b).As can be seen, emin decreases to a higher level as compared to emax. Similar reductions in extremevoid ratios are observed in the literature with increasing blockiness [136, 159]. It should be notedthat no specific study has been discovered that investigates the packing of blocky and/or ellipsoidalparticles which considers the similar Cu, as considered here. As a result, the literature cited aboveonly discusses qualitative trends. The detailed quantitative comparisons are made for samples withspherical particles (AR = 1, B = 2) at different Cu, which can be found in Chapter 4. Based onextreme void ratios and varying the value of µ , samples are prepared at two different levels of Dr,792 4 8 12D [mm]020406080100FV[%]Figure 5.2: Log-linear particle size distribution of samples with Cu = 2 for all values of AR and Bconsidered in this study.namely 20% and 50% at a fixed Cu of 2.0 (refer Figure 5.2). As mentioned later in the chapter, aDr of 20% is chosen to be very close to the critical state void ratio, while a Dr of 50% is chosento show significant dilation when samples are subjected to a monotonic drained simple shear test.In accordance with the test protocol outlined below, the samples are then subjected to constantvolume cyclic loading.5.3 Test designEach sample designed in this study contains ten distinct classes of particles of varying sizes,corresponding to a log-linear particle size distribution with a Cu of 2, D50 of 2.5 mm (refer toFigure 5.2), and a Np of 10000 particles. The shearing tests used to investigate the effects of ARand B are referred in the Table 5.1. The effect of shape on cyclic liquefaction resistance isexamined using a total of 41 cyclic shearing tests at two different relative densities. Thenomenclature AR05 B2 Dr20 CSR15 describes the cyclic test for a sample with AR = 0.5, B = 2,and Dr = 20%, sheared under a CSR of 0.15. For the low Dr of 20%, the CSR varies from 0.10 to0.25 (with some exceptions). While the CSR ranges between 0.25 and 0.55 with an equal interval80Table 5.1: Simulated constant volume cyclic simple shear tests on samples with p0 = 100 kPa.Test ID† AR (−) B (−) Dr (%) CSR (−)AR05 B2 Dr** CSR** 0.5 2 20 0.10, 0.15, 0.20, 0.2550 0.25, 0.35, 0.45, 0.55AR10 B2 Dr** CSR** 1.0 2 20 0.10, 0.15, 0.20, 0.2550 0.25, 0.35, 0.45, 0.55AR15 B2 Dr** CSR** 1.5 2 20 0.10, 0.15, 0.20, 0.2550 0.25, 0.35, 0.45, 0.55AR10 B4 Dr** CSR** 1.0 4 20 0.10, 0.15, 0.20, 0.2550 0.25, 0.35, 0.45, 0.55AR10 B8 Dr** CSR** 1.0 8 20 0.05, 0.10, 0.15, 0.20, 0.2550 0.25, 0.35, 0.45, 0.55† Test ID consists of the particle aspect ratio, blockiness, sample relative density, and cyclic stressratio.of 0.10 for the medium high Dr of 50%. The next section discusses the cyclic results of a specifictest ID, followed by a comprehensive discussion of cyclic liquefaction.5.4 Macroscopic response5.4.1 Evolution of excess pore pressureThe evolution of the ru value is crucial because it indicates soil stiffness degradation, theaccumulation of shear strain, and the rate of progression to the liquefaction stage. Numerousstudies have demonstrated that the rate at which ru evolves is a good indicator of the level ofcatastrophy caused by liquefaction-related failures [160, 161]. The evolution of ru can be affectedby the particle shape of assemblies [86]. In Figure 5.3, the effects of ellipsoidal and blockyparticles on the evolution patterns of ru accumulated at the beginning of each loading cycle whenτ ≈ 0 are examined with normalized cycles to liquefaction (N/NIL) at Dr = 20% and 50%,respectively. The figure 5.3(a) illustrates the evolution of ru at Dr = 20% with each cycle of cyclicloading when τ ≃ 0. The value of ru gradually increases from 0.0 to 0.38 as N/NIL increasesfrom 0.0 to ≃ 0.8 when AR = 0.5. As N/NIL increases from 0.8 to 1.0, a significant rise in ruoccurs from 0.38 to 0.95. A similar trend is observed for the sample with AR = 1.5, with an initial810.0 0.2 0.4 0.6 0.8 1.0 1.2N/NIL0.00.20.40.60.81.0r u Dr = 20%AR = 0.5AR = 1.0AR = 1.5(a)0.0 0.2 0.4 0.6 0.8 1.0 1.2N/NIL0.00.20.40.60.81.0r u Dr = 50%AR = 0.5AR = 1.0AR = 1.5(b)0.0 0.2 0.4 0.6 0.8 1.0 1.2N/NIL0.00.20.40.60.81.0r u Dr = 20%B = 2B = 4B = 8(c)0.0 0.2 0.4 0.6 0.8 1.0 1.2N/NIL0.00.20.40.60.81.0r u Dr = 50%B = 2B = 4B = 8(d)Figure 5.3: The effect of particle shape on the evolution of accumulated ru at CSR = 0.25 and τ ≈ 0 for:(a) AR = 0.5 to 1.5 and Dr = 20% (b) AR = 0.5 to 1.5 and Dr = 50%, (c) B = 2 to 8 and Dr = 20%, and (d)B = 2 to 8 and Dr = 50%.gradual evolution followed by a relatively steep build-up, though this transition occurs earlier, atN/NIL ≃ 0.7. At N/NIL ≃ 0.6, the transition occurs more rapidly when AR = 1. It is vital to notethat samples with AR values of 0.5 and 1.5 (oblate and prolate ellipsoids) develop ru at a slowerrate than the spherical sample during the first 60% to 80% of cycles to N/NIL = 1. They do,however, generate ru at a faster rate than spherical assemblies during the final stages ofpreliquefaction cycles. As a result, the presence of such particle shapes may result in more abruptliquefaction than spherical samples. When the Dr is increased to 50%, the same set of samplesexhibits a slightly different trend, as illustrated in Figure 5.3(b). The sample with AR = 0.582demonstrates a gradual increase in ru until N/NIL = 0.8, followed by a sharp increase to 0.95when N/NIL = 1. In addition, this abrupt increase in ru is far less sudden than it is when Dr =20%. As AR is increased to values of 1.0 and 1.5, the evolution of ru demonstrates a gradualincrease from initial shearing to liquefaction at N/NIL = 1. A similar type of analysis isperformed on samples with different B values at Dr = 20% and 50%. At Dr = 20%, samples withB = 2 and 4 exhibit a similar trend in their evolution of ru, as illustrated in Figure 5.3(c). Thetransition from the gradual to the steeper or rapid stage of ru evolution occurs at N/NIL ≃ 0.6.Whereas in the sample with a value of B equal to 8, this transition occurs later, at N/NIL ≃ 0.7,and is followed by a rapid buildup to liquefaction. At Dr = 50% (refer to Figure 5.3(d)), theresponses are more gradual for samples with B = 2 and 4. Nevertheless, the buildup of ru for B =8 exhibits a two-stage characteristics leading to liquefaction, with the transition occurringapproximately at N/NIL = 0.8. In summary, ellipsoidal and blocky particle assemblies are foundto be more prone to abrupt liquefaction than round or spherical assemblies.5.4.2 Cyclic liquefaction resistanceIt is standard practice in geotechnical engineering to construct the liquefaction resistance curveusing NIL at various levels of CSR. This method is employed in the current study to describe theeffect of particle shapes (aspect ratio and blockiness) on cyclic liquefaction resistance. The cyclicliquefaction strength curves for samples with AR values of 0.5, 1, and 1.5 are shown in Figure 5.4.As the CSR rises, the value of NIL decreases, confirming the long-established fact that higherdemand reduces the resistance of a material to cyclic liquefaction [19]. In simulations with Dr =20%, sample with AR = 1 has the highest cyclic liquefaction resistance, followed by AR = 1.5, andsample AR = 0.5 has the lowest resistance, though the differences are found to be marginal. Sincethe liquefaction curves do not intersect, the above observations can be made regardless of theCSR. Observations contrast when a higher Dr of 50% is considered. Cyclic simulations indicatethat the sample with AR = 1 shows the least resistance to liquefaction, followed by AR = 1.5, and83the sample with AR = 0.5 shows the highest cyclic liquefaction resistance. At this relative density,1 10 100NIL0.00.20.40.6CSRDr = 20%Dr = 50%AR = 0.5AR = 1.0AR = 1.5Figure 5.4: Cyclic liquefaction strength curves of loose and medium-dense samples with different valuesof AR, when B is kept constant at 2.it is discovered that the difference in NIL values between AR is CSR dependent. While the samplewith AR = 1.5 has a significantly higher liquefaction resistance than the sample with AR = 1, thedifference gradually narrows as the CSR is increased during the cyclic simulations. Moreover,cyclic resistance of AR = 0.5 is higher than the sample with AR = 1.5, and the former is found todisplay substantially greater cyclic resistance as CSR is increased. At this density, the curves donot overlap within the examined range, indicating that there is no qualitative shift in response tovariation in CSR. In summary, when compared to spherical particle assemblies, ellipsoidalparticles are less resistant to cyclic liquefaction at low relative densities but more resistant andCSR dependent at medium densities.The cyclic liquefaction strength curves for samples with B values of 2, 4, and 8 are shown inFigure 5.5. In simulations with Dr = 20%, spherical particle samples with B = 2 have the highestcyclic liquefaction resistance. As B increases from 2 to 4, then to 8, there is a progressive decreasein cyclic liquefaction resistance, which is amplified by increases in CSR. However at a higher Dr841 10 100NIL0.00.20.40.6CSRDr = 20%Dr = 50%B = 2B = 4B = 8Figure 5.5: Cyclic liquefaction strength curves of loose and medium-dense samples with different valuesof B when AR is kept constant at 1.= 50%, as the B of the samples increases, for the same CSR, the cyclic liquefaction resistanceincreases monotonically. In this case, CSR variation does not contribute to the relative differencein liquefaction resistance among samples with varying blockiness. In conclusion, the CSR-NILplots show that non-spherical particles have lower liquefaction resistance at low Dr and elevatedliquefaction resistance at medium Dr than spherical particles.The data presented in Figures 5.4 and 5.5 are simplified and presented in Figure 5.6, by takingone point from each CSR curve at NIL equal to 15. In Figures 5.6(a) and 5.6(c), the respectiveCRR15 values are plotted with AR and B, respectively. It is affirmed that CRR15 decreases asparticles become non-spherical (AR ̸= 1, B ̸= 2) for samples at Dr = 20%, whereas, thesesamples show significant increase in CRR15 when compared at a higher Dr = 50%. Theobservations are corroborated when CRR15 is replaced by NIL to re-present the data for a CSR =0.3, as shown in Figures 5.6(b) and 5.6(d). While the liquefaction resistance curves indicate thatparticle shape is CSR dependent, the characteristic trends in particle shape and relative densityremain unchanged for the range of CSR considered in this study. Ashmawy et al. [10]850.5 1.0 1.5AR0.000.150.300.450.60CRR15Dr = 20%Dr = 50%(a)0.5 1.0 1.5AR01020304050NILDr = 20%Dr = 50%(b)2 4 6 8B0.000.150.300.450.60CRR15Dr = 20%Dr = 50%(c)2 4 6 8B01020304050NILDr = 20%Dr = 50%(d)Figure 5.6: Variation of liquefaction resistance in terms of (a) interpolated values of CRR15 with AR, (b)NIL with AR, (c) interpolated values of CRR15 with B, and (d) NIL with B.demonstrated density-dependent behavior of non-spherical assemblies on cyclic liquefactionresistance, but their study was limited by the fact that initial state variables like density,anisotropy among others were not isolated, and they used a simplified 2D-DEM method.Attempts are made to unify the different aspects of shape viz., AR and B under a single shapeparameter. Several researchers have used compound indicators like overall regularity (OR) torepresent particle shape [76, 86] defined asOR =R+S2, (5.1)86Table 5.2: Roundness, sphericity and overall-regularity of the samples.AR B R S OR0.5 2 0.67 0.5 0.591 2 1 1 11.5 2 0.92 0.67 0.791 4 0.50 0.79 0.641 8 0.25 0.63 0.44where, R and S are the roundness and sphericity of an assembly of particles. In this way CRR15 canbe directly linked with a single shape parameter and conclusions drawn can be further simplified.Hence, shape parameters used in the numerical analysis A and B are used to calculate R and S asper the procedure presented in section 3.8.2, and thereby used to calculate OR. The values arepresented here in Table 5.2. The maximum OR is found in spherical particles with B = 2 and AR =1, while the minimum OR is found in cubic particles with B = 8 and AR = 1. The values of CRR15and NIL presented earlier for AR and B are plotted with the compound shape indicator OR in Figure5.7. It is very succinct that CRR15 shows a increasing trend with the increase in OR for the low Dr =20% as shown in Figure 5.7(a). This is found to be valid for both AR and B, the two different typesof shape indicators used in the study. In contrast, at the higher Dr = 50%, the values of CRR15 seemto follow a negative i.e., decreasing slope with OR. Comparable trend is observed in the relationbetween NIL and OR as shown in Figure 5.7(b). The contrasting nature of the CRR15-OR curvewith change in Dr highlights a very important and novel contribution from the current researchwork, under the pre-condition that all other indicators like size distribution, surface roughness andso on, are completely isolated from shape. For dense sands, reduction in liquefaction resistancewith increasing OR was also noted by Rui et al. [86]. This experimental study found that non-spherical samples with higher liquefaction resistance had more inter-particle contacts, leading toincreased interlocking. In the following section, we extrapolate this “connection” between inter-particle contacts and cyclic liquefaction for different particle shapes.870.4 0.5 0.6 0.7 0.8 0.9 1.0OR0.000.150.300.450.60CRR15Dr = 20%Dr = 50%AR = 0.5AR = 1AR = 1.5B= 4B= 8(a)0.4 0.5 0.6 0.7 0.8 0.9 1.0OR0102030405060NILDr = 20%Dr = 50%(b)Figure 5.7: Variation of cyclic liquefaction resistance indicators (a) CRR15, and (b) NIL with overallregularity (OR).5.5 Relation between CRR15 and initial state5.5.1 Relation with z0Micromechanical analysis of the effects of particle shape on liquefaction resistance is used tobetter understand the particle-scale origins of the macroscopic response. Morimoto et al. [93]884 5 6 7z0.00.20.40.6CRR15Dr = 20%Dr = 50%AR = 0.5AR = 1AR = 1.5B= 4B= 80Figure 5.8: The relation between z0 with CRR15 for samples with different B and AR.recently explored the correlation between liquefaction resistance and the initial value of zmdenoted here as z0. Similarly, Figure 5.8 describes a connection between the liquefactionresistance indicator CRR15 and z0. Based on Dr, z0 values are divided into two groups, eachrepresenting the set of AR and B considered in the study. It is revealed from the figure that thevariation of CRR15 may be affected by z0. The correlation is specially noteworthy when thesamples are compared at the higher Dr of 50% with a good overall exponential fit. The correlationbetween the coordination number z (including floaters and rattlers) and the CRR15 is excluded inthis study due to its lower degree of correlation. It refers to the fact that, regardless of particleshape or relative density, the contacts established within the “participating network” should be areliable indicator of liquefaction resistance. Representations of the “participating network” interms of coordination number may be feasible using micro CT images [78, 162] from in-situsamples, and as suggested by the current study and recently by Morimoto et al. [93] could helpenrich liquefaction susceptibility assessments.890 4 8 12 16fn/ 〈fn〉10−310−210−1100P(f n)Dr = 20%1AR = 0.5AR = 1AR = 1.5(a)0 4 8 12 16fn/ 〈fn〉10−310−210−1100P(f n)Dr = 50%1AR = 0.5AR = 1AR = 1.5(b)0 4 8 12 16fn/ 〈fn〉10−310−210−1100P(f n)Dr = 20%1B = 2B = 4B = 8(c)0 4 8 12 16fn/ 〈fn〉10−310−210−1100P(f n)Dr = 50%1B = 2B = 4B = 8(d)Figure 5.9: Probability distribution functions of normalized contact normal forces fn/⟨ fn⟩ for sampleswith different B and AR at (a) AR = 0.5, 1 and 1.5, Dr = 20%, (b) AR = 0.5, 1 and 1.5, Dr = 50%, (c) B = 2, 4and 8, Dr = 20%, and (d) B = 2, 4 and 8, Dr = 50%.5.5.2 Relation with contact normal force distributionsApart from the contact network, assessments of the contact force network can provide insight intothe microstructure [47] the way which might help understand its influence under different particleshapes at the initial state, and how they connect with cyclic liquefaction resistance. In case ofDr = 20%, Figure 5.9(a) shows that for samples with AR not equal to 1, the distribution of strongcontact normal forces is wider than for spherical (AR equal to 1) samples. Considering that themean global stress is the same in all assemblies at the initial state, the wider force distributionin non-spherical (prolate and oblate ellipsoids) samples indicates fewer contact channels and thus900.4 0.6 0.8 1.0β0.00.20.40.60.81.0CRR15Dr = 20%Dr = 50%AR = 0.5AR = 1AR = 1.5B= 4B= 8Figure 5.10: The relation between the rate of decay β of the normalized strong contact normal forcenetwork fn/⟨ fn⟩ > 1 for different AR and B with CRR15.fewer force transmission networks. Due to the fact that forces are transmitted via particles withtwo or more contacts, the relative paucity of force networks in non-spherical particles at this Drmay also indicate a lower zm value. A recent study by Marteau and Andrade [151] also cites thelink between contact normal force heterogeneity (wider distribution) and lower zm. For Dr = 50%,the contact-normal force distributions change as shown in Figure 5.9(b). For samples with ARless than or greater than 1, the distribution of strong contact forces is more uniform. This showsevidence that non-spherical aspect ratios can reduce the heterogeneity of contact normal forces athigher relative densities, as opposed to lower densities. For this reason, non-spherical assembliesmay have more contact pathways and thus comparatively higher zm than compressed sphericalparticle assemblies. Further, in terms of particle shape, it has been discovered that a change inrelative density may transition the nature of force distribution.Figures 5.9(c) and 5.9(d) show force distributions for various particle shapes described byindicator B at Dr = 20% and 50%, respectively. Non-spherical samples obtained by increasingthe B values (2 to 8) at Dr = 20%, result in wider force distributions. This is qualitatively similarto the wider distribution observed when different AR values are used to achieve non-sphericity.Furthermore, the transition of contact normal distributions observed for various AR when Dr equals50%, can also be seen in Figures 5.9(c) and 5.9(d) for various values of B. In summary, both91(a) (b) (c)(d) (e) (f)Figure 5.11: Snapshots of normal force chains of the samples with different aspect ratios at the initialstate for: (a) AR = 0.5, Dr = 20%, (b) AR = 1, Dr = 20%, (c) AR = 1.5, Dr = 20%, (a) AR = 0.5, Dr = 50%,(e) AR = 1, Dr = 50%, and (f) AR = 1.5, Dr = 50%.particle shape indicators produce non-sphericity that has a similar qualitative effect on contactnormal force distributions.The decay in the strong normal force distributions is further investigated by quantifying itsslope with the exponent β . The strong normal force distribution becomes less heterogeneous as βincreases and vice-versa. The potential dependence of β on the liquefaction resistance indicatorCRR15 is studied for all levels of AR, B and Dr as shown in Figure 5.10. This figure demonstratesthat β has a positive correlation with CRR15, as indicated by an exponential curve fitting. As β92increases, greater homogeneity of contact normal forces indicates more contact pathways andcontacts per particle, possibly leading to this positive association. Overall, it can be stated thatregardless of the limited type of particle shape indicators and relative density used in this study,increased homogeneity of the strong contact normal forces may increase cyclic liquefactionresistance of granular assemblies irrespective of shape.Direct visualization studies of force distributions as “force chains” are carried out in limitedexperimental studies using photoelastic disks [146, 163]. To visualize a force chain in a DEM(a) (b) (c)(d) (e) (f)Figure 5.12: Snapshots of normal force chains of the samples with different blockiness at the initial statefor: (a) B = 2, Dr = 20%, (b) B = 4, Dr = 20%, (c) B = 8, Dr = 20%, (d) B = 2, Dr = 50%, (e) B = 4, Dr =50%, and (f) B = 8, Dr = 50%.93simulation, colored tubes of varying sizes (depending on the magnitude of the force) connect thecenters of two contacting particles that transmit forces. The force chains are used in this studyto aid in the qualitative visualization of the distribution of contact normal forces with effect fromAR, B, and Dr. They are presented in Figures 5.11 and 5.12, with insets illustrating a zoomed-in area of the force network. Smaller contact normal forces are represented by blue tubes withsmaller diameters, while larger contact normal forces are represented by red tubes with largerdiameters. In Figures 5.11(a) to 5.11(c), the sample with AR = 1 has marginally thinner red tubesthan the samples with AR = 0.5 and 1.5, indicating slightly lower heterogeneity than the ellipsoidalsamples. This observation also correspond to the wider distribution of fn/⟨ fn⟩ in Figure 5.9(a) forAR = 0.5 and 1.5, and to their β values lower than AR = 1 in Figure 5.10. With Dr = 50% asshown in Figures 5.11(d) to 5.11(f), the distributions of larger diameter force tubes appear to fadeaway at AR = 0.5 and 1.5 (as indicated by the insets), but are more prevalent with AR = 1. Thisindicates that the contact normal forces are more uniform (less heterogeneous) at AR = 0.5 and 1.5,as previously observed in Figure 5.9(b). Similar qualitative observations can also be made for B= 2, 4 and 8 in Figure 5.12 for both levels of Dr. As explained previously for AR, this figure canvisually support the observations made in the contact force distributions for various B in Figures5.9(c) and 5.9(d). To avoid repetition, the comparisons are assumed to be self-evident.The analysis based on the initial state is found to be a reasonable direction for unravelingthe cyclic liquefaction resistance of a sample. Meanwhile, initial state variables like void ratio,coordination number and intrinsic variables like particle shape influence the mechanical behaviorin terms of volumetric response. It should be kept in mind that liquefaction is a function of thevolumetric response of a specimen. However, information on the initial state alone is insufficient,as knowledge of the critical state, also known as the end state, is required. The state parameter,which includes both the initial and critical state variables, better represents volumetric response. Inthe following, state parameters are obtained with the aim of using them as a predictor of the cyclicliquefaction resistance of granular materials.945.6 Relation between CRR15 and state parametersA state parameter is floated as that parameter that can accurately describe the state of a granularmaterial better than void ratio or relative density. This is because it compares the difference of theinitial state with the critical state of a material at the same confinement stress. This is claimed tobetter capture dilatancy than the majority of other parameters, which is believed to be the sourceof all shear resistance, including the cyclic strength. However, there are other factors like particlesize distribution that can influence the dilatancy of materials, like differences in particle sizedistribution. Nevertheless, there is no direct evidence about whether a unique relation betweenstate parameter and the cyclic strength (expressed here as CRR15) can be presented, which cancapture the effect of particle size distribution.The state parameters ψe, and ψz are plotted against the cyclic strength parameter CRR15, andpresented in Figure 5.13. In Figure 5.13(a), the ten markers representing ψe are calculated fromsamples with different AR and B at Dr = 20% and 50%, respectively. An exponential correlationis found between these two parameters with a R2 value of 0.94. Evidently, as ψe decreases from0.05 to -0.15, the observed CRR15 values rise exponentially. This clear association indicates thatthe effect of the AR and B on liquefaction resistance can be reasonably captured using void ratio-based ψ regardless of relative density. This link confirms that monotonic drained sand responsecan reveal key cyclic response aspects in granular materials like liquefaction resistance [164, 165].Additionally, the presence of mostly negative ψe values indicates that the majority of the samplesconsidered in this study are dilative, with the exception of those with B = 4 and 8 at Dr = 20%.Since ψe signify the net dilation or contraction of a sample, a good correlation with CRR15 signifya direct link between volumetric tendency and cyclic liquefaction resistance of isotropic samples.Similarly, ψz is also used as an indicator, based on the micro-state, to correlate with CRR15. Theresults are shown in Figure 5.13(b), where an exponential fit to describe the data trend can be usedwith a higher regression coefficient than ψe. Thus, the state parameter derived from zm exhibits astrong correlation with liquefaction resistance for all samples. Specifically, micro-state parameter95−0.15 −0.10 −0.05 0.00 0.05ψe0.00.20.40.6CRR15(a)−0.5 0.0 0.5 1.0 1.5 2.0 2.5ψz0.00.20.40.6CRR15Dr = 20%Dr = 50%AR = 0.5AR = 1AR = 1.5B= 4B= 8(b)Figure 5.13: Cyclic liquefaction resistance indicator CRR15 against state parameters for samples withdifferent particle shapes and relative density derived from: (a) void ratio, and (b) mechanicalcoordination number.zm provides the opportunity for a strong macro-micro relationship to be established. This impliesthat samples with a greater contact loss due to dilation in drained constant p tests show higherresistance to cyclic liquefaction at the various particle shape and relative density levels examined.Despite the fact that the current set of tests are performed on isotropic samples, ψz can also be usedto establish a strong linkage with CRR 15 for anisotropic samples where ψe is not a good choice[154] for the reasons stated previously.965.7 Evolution of zmThe microstructure of a granular material can be characterized by coordination number and avariety of other fabric indicators like arrangement of particles and their contacts, arrangement ofcontact force distributions, void spaces among others [166]. This research focuses on theevolution of coordination number, a reliable indicator to represent fabric characteristics, underconstant volume cyclic loading on granular materials. The evolution of zm can be seen in Figure5.14 for all levels of particle shape and relative density considered. It can be observed that asN/NIL is increased, there is a constant decrease in zm, indicating there is an overall loss ofcontacts in the assembly. As a result more particles are generated with unstable and metastableconfigurations. This is similar to other DEM studies [105, 110]. For samples prepared at a Dr of20%, zm starts with a small range of initial values corresponding to varying AR and B levels, asillustrated in Figures 5.14(a) and 5.14(c). When samples with varying AR and B values arecyclically loaded, they appear to lose contact at a similar, gradual and nearly linear rate; with thefew cycles prior to initial liquefaction causing a rapid degradation of zm. At Dr = 50%, the initialvalues of zm for various AR (refer to Figure 5.14(b)) are quite widely distributed than at Dr =20%. The evolution plot for AR = 0.5 (black line) reveals a gradual decline in zm until N/NIL =0.8; increasing N/NIL further to 1 results in a sharp decline in zm. Whereas, zm for samples withAR values of 1 and 1.5 exhibits gradual degradation throughout the pre-liquefaction stage untilN/NIL = 1, with no evidence of a sharp decline prior to initial liquefaction. Figure 5.14(d) showsthe evolution of zm for samples at Dr = 50% and B = 2, 4 and 8. All samples show a gradualdecline in zm before N/NIL = 0.8. However, in the final few cycles prior to initial liquefaction, anincrease in B results in a change in the evolution from a gradual to a more rapid decline, implyingthat the rate of contact loss increases with increasing blockiness at this stage. In summary, whengranular samples are subjected to constant volume cyclic loading, increasing relative density andnon-sphericity in terms of aspect ratio and blockiness can result in a faster loss of contacts and“rapid” entry into cyclic liquefaction.970.0 0.2 0.4 0.6 0.8 1.0N/NIL0123456z mDr = 20%AR = 0.5AR = 1.0AR = 1.5(a)0.0 0.2 0.4 0.6 0.8 1.0N/NIL0123456z mDr = 50%AR = 0.5AR = 1.0AR = 1.5(b)0.0 0.2 0.4 0.6 0.8 1.0N/NIL01234567z mDr = 20%B = 2B = 4B = 8(c)0.0 0.2 0.4 0.6 0.8 1.0N/NIL01234567z mDr = 50%B = 2B = 4B = 8(d)Figure 5.14: The effect of particle shape on the evolution of zm for (a) AR = 0.5, 1, 1.5 at Dr = 20%, (b) AR= 0.5, 1, 1.5 at Dr = 50%, (c) B = 2, 4, 8 at Dr = 20%, and (d) B = 2, 4, 8 at Dr = 50%.5.8 Analysis of the isostatic stageThe figure 5.15(a) illustrates a specific response of a sample with a Dr = 50% and cyclically shearedat a CSR = 0.25. With progressive loading cycles, the σv value decreases, accompanied by agradual decay of the zm value. However, as σv approaches 0, an arrow is used to indicate an abruptdecrease in the zm along a precise corner. As the assembly rapidly loses contact, the corner isrepresentative of the critical point at which the solid-like behavior transitions to fluid-like behavior.The value of zm at which this transition occurs is defined as ziso. In Figure 5.15(b), the value of zisowith the backdrop of NIL.980 20 40 60 80 100[kPa]123456z mz = 2.8B = 2Dr = 50%CSR = 0.25Liquefied zoneiso(a)0 5 10 15 20 25 30N123456z mz = 2.8NIL = 17ru < 0.95ru > 0.95iso(b)Figure 5.15: The value of ziso is described with (a) the corner in the zm curve with an arrow, and (b)position of ziso curve with the corresponding value of N.Hereafter, the values of ziso are plotted for different Dr and AR in Figure 5.16(a). It is discoveredthat the value of ziso increases when Dr is increased from 20% to 50%. At Dr = 20%, ziso presentsa “A” shaped variation with aspect ratio, increasing slightly from AR = 0.5 to 1, followed by asmall decrement as AR equals 1.5. In contrast, for Dr = 50%, ziso is found to have a “V” shapedrelation with AR, where increase in AR from 0.5 to 1, results in substantial decrease in ziso, and as itincreases from 1 to 1.5, results in a notable increase in ziso. Additional observations on the relationbetween ziso and B are reflected in Figure 5.16(b). At a Dr of 20%, B has a weak influence on ziso,(a) (b)Figure 5.16: Variation of ziso at different relative densities with (a) different AR at B = 2, and (b) differentB at AR = 1.99but as the relative density of the samples rise to 50%, ziso grows considerably as B is increasedfrom 2 to 4, and slightly less augmented as B is raised from 4 to 8. It should be noted that sinceziso is determined to be unaffected by the CSR level, all findings are shown for a given CSR valueof 0.25.5.9 ConclusionsThe effect of particle shape on the cyclic liquefaction resistance of isotropic, loose, and medium-dense samples was investigated using 3D-DEM.The macroscopic observations based on the effect of particle shape, relative density and CSRwere expressed as CSR-NIL curves. The initial state and critical state based descriptors were used todescribe the characteristics of the cyclic liquefaction curves irrespective of their particle shape andrelative density. Finally, contact-based analyses were performed to ascertain the critical number ofparticle-particle contacts that could define the solid to fluid transition.At the same Dr of 20%, samples with cubic (blocky) particles exhibited less resistance thansamples with spherical particles. However, the observation greatly varied for the same particleassemblies when compared at a higher Dr = 50%. In this case, for all CSR levels considered, itwas discovered that ellipsoidal and cubic particles exhibited a significant increase in liquefactionresistance when compared to spherical samples.It was found that state parameter ψe demonstrates a well-defined exponential association withthe liquefaction resistance indicator CRR15. This resulted in the conclusion that DEM studies canmodel the exponential relation between a void ratio-based state parameter and the the liquefactionresistance of isotropic (or similar anisotropy) granular materials regardless of the particle shape andrelative density. However, the different coefficients of the exponential fit for particle shape-basedrelations and the same for particle size distribution points to the fact that a single relation involvingboth shape and size distribution is not possible. This might be due to the slight rolling friction usedfor particle size distribution, which was zero for particle shape analysis in this chapter.100It was found that CRR15 increases exponentially with the increase in ψz . This association ofψz was found to be the most robust among other parameters which were linked with CRR15 in thisstudy.The value of ziso was observed to reduce with decrease in Dr. It was also shown that a changein particle shape in terms of AR and B can appreciably alter its value.Finally, future research should build on the current liquefaction-based study by examining thepost- and re-liquefaction responses of samples with a variety of particle shapes. While investigatingthe effect of particle shape on the cyclic strength of granular materials, a variety of particle shapesand particle size distributions may also be considered.101Chapter 6: Summary, Conclusions and Future Work6.1 SummaryThis dissertation uses three-dimensional discrete element modeling (3D-DEM) to make aqualitative assessment on the effects of particle size distribution and particle shape on theresponse of granular materials under constant-volume cyclic loading, with the main emphasis oncyclic liquefaction.The first contribution investigates the effects of particle size distribution on cyclic response.The coefficient of uniformity (Cu) is used to quantify size distribution among virtual samples withthe same particle shape, mean particle size, and surface roughness. In addition, all the samples areprepared under the identical sample preparation protocol to achieve the same isotropic compressivestress. The samples are also prepared at the same relative density (Dr) to eliminate the differencesin density from affecting the shearing response. This level of control is possible via DEM modelingand further it cam be used to isolate all effects other than Cu from the following cyclic shearingregime at constant volume. All samples have been sheared under varied levels of cyclic stressratio (CSR) and Dr at-least until initial liquefaction. The effects of Cu on the cyclic liquefactionresistance are observed to be significant, and the nature of such effect depends on the Dr. Whileat a low Dr, the liquefaction resistance showed an overall decreasing trend with an increase in Cu,this effect is observed to be decreasing at a medium dense Dr. To further explore this variability,one of the widely used liquefaction resistance indicators, CRR15 is correlated with microscopicindicators that define the initial state (pre-shearing) of the compressed assemblies. In addition, twostate parameters based on the macrostate and the microstate, respectively: (a) based on void ratio102and (b) on the mechanical coordination number, are used to investigate correlations with CRR15.The goal of the correlations, particularly those based on the microstate, is to explain the trend incyclic liquefaction resistance independent of Cu and Dr, thereby establishing a macro-micro linkin the overall liquefaction analysis.In the second investigation, a similar structure of analysis has been followed while studying theeffects of particle shape. Aspect ratio and blockiness are variables used to independently vary theparticle shape among different samples. In spite of using simplified geometries like near cuboidsand perfect ellipsoids, these shape variables have been found to capture both the roundness andsphericity aspects of real, more complex particle shapes. In this regard, the Cu, mean particlesize, and surface roughness are kept constant. The sample preparation has been used to preparesamples with different particle shapes under isotropic compression. The cyclic simple shear testsconducted under constant volume conditions revealed that particle shape has a significant effect onthe cyclic liquefaction resistance, and the influence is density dependent. Further analysis based oninitial state based parameters and state parameters revealed that mechanical coordination numberpresents a significant advantage in predicting the liquefaction resistance in terms of CRR15 whichis typically difficult to explain from the current findings, through macro based techniques. This istrue even when the tests are performed at different Dr and CSR levels, as was discovered previouslyin the case of particle size distribution analysis.It should be mentioned at the end that qualitative validation was done based on the effectsof relative density and CSR on the liquefaction strength. Certain aspects of the effect of particleshape on liquefaction strength were also qualitatively comparable to experiments. For instance,the increased liquefaction resistance observed with non-spherical samples prepared at a mediumdensity state can also be observed in experiments.It should be noted here that no quantitative validation with experiments was done, and noaspects of laboratory-based sample preparation were objectively replicated in the simulations. Thegoal was to make sample preparation more controlled and qualitative aspects were studied, asopposed to a calibration based quantitative study. Also, the type of simple shearing followed is103typical of a simple shear type of load but does not replicate the simple shear apparatus that is usedin laboratory experiments.6.2 ConclusionThe main contributions of the research are as follows:• Three dimensional DEM modeling of particle shape and particle size distribution tosystematically study their effects on cyclic liquefaction resistance for isotropic samples.• Microscopic investigation into the contact network, and the normal force network to unraveltheir relation with the macroscopic response of liquefaction resistance.• State parameter based approach that can be applied independently to study the isolatedeffects of particle shape and particle size distribution based on void ratio (macro-level) andcoordination number (micro-level) in terms of their liquefaction resistance.6.2.1 Effect of particle size distributionSamples with varying particle size distributions are compared at different relative densities forcyclic liquefaction resistance. At a low relative density, samples with low polydispersity exhibithigher cyclic liquefaction resistance than samples with purely uniform particle sizes. However,with a further increase in polydispersity, the liquefaction resistance is found to decrease. Whensamples are compared at a medium relative density, the introduction of small polydispersity isfound to reduce the liquefaction resistance, while further increasing it shows increased resistance.It should be noted that the observations are based on simple shear loading, in which the principalstress rotates smoothly during the shearing process, influencing the undrained response.While the macroscopic characteristics from the current study are quite counterintuitive, theyare found to be reasonably explained by the mechanical coordination number, that is, the averagenumber of contacts in an assembly without floaters and rattlers. This microparameter, when used104to define the state parameter, is found to best explain the complex nature of cyclic liquefactionresistance for samples with different particle size distributions, irrespective of the relative densityand loading intensity. It has been worth noting that the samples with higher size distributionsalways have a lower void ratio than their more uniform counterparts, but the absence of void spacesdoes not consistently guarantee higher liquefaction strength as opposed to the current notion. Thishas been attributed to lower mechanical coordination numbers and increased floaters and rattlers,lesser dilative tendencies, and a high level of arching effects observed with higher size distributions,with fewer larger particle networks taking the bulk of the effective stress. Apart from predictingthe characteristics of liquefaction resistance, the coordination number is also used to define thesolid-fluid transition that occurs during initial liquefaction. This has provided information for theminimum number of contacts necessary to maintain a solid soil structure for completely saturatedgranular materials at different size distributions and relative densities.6.2.2 Effect of particle shapeAt various levels of relative density, the effects of particle shape variation on the cyclic responsehave been investigated using a similar approach to that used for particle size distribution. Whilethe particle shapes used are simplified in comparison to the complex shapes found in real soils,they allow for adequate modeling and variation of roundness and sphericity, two critical shapecharacteristics. The results of undrained cyclic simple shear tests on isotropically compressedsamples reveal that particle shape has a significant effect on the liquefaction resistance of granularmaterials depending on their relative density. When comparing non-spherical particle assembliesto purely spherical particle assemblies at a medium density, it was discovered that they aresignificantly more resistant to liquefaction, whereas they are noticeably less resistant whensamples are prepared at a low relative density.This is contrary to the prevailing understanding among researchers that a non-spherical shapeinvariably contributes to a better interlocked system that leads to higher shear resistance. The105current study highlights that for loose samples, that does not hold true, along with the fact thatinterlocking itself is mostly qualitative and cannot be successfully measured across all samplesto explain the liquefaction behavior. Instead, the mechanical coordination number calculated atthe initial state and more convincingly when used as a state parameter, can successfully explainthe liquefaction resistance across various shapes, relative densities, and loading intensities. Theminimum contact at the solid-fluid transition, also known as isostaticity, during initial liquefactionis found to be dependent on the particle shape and relative density, and independent of the loadingintensity of the samples.6.3 Recommendations for future workThe recommendations aim to challenge the limitations and scope of the current research with aimtowards studying a wider range of particle size distribution and shape than studied in the currentthesis. The five important items are listed in the following in the optimal order of importance andease of study:• Investigate the post-liquefaction shear strain generation under the influence ofparticle size distribution: The current dissertation focuses on the initial state,preliquefaction, and initial liquefaction of granular samples. Large and rapid shear straingeneration occurs during the post-liquefaction stage. It is debatable whether thisaccumulation spirals out of control or whether such rapid shear strain generation has anupper bound or saturation limit. This evidence of strain saturation plays a role in bothpractical and theoretical aspects of granular flows. DEM can be used to study factors likeparticle size distribution and particle shape that may induce or prohibit such an occurrence.• Effects of particle size distribution for non-spherical samples: Using only sphericalparticles, the dissertation examined the effect of particle size distribution on cyclicliquefaction. As a result, the natural extension of the work will be to investigate the effect106of particle size distribution on cyclic liquefaction when non-spherical particle geometriesare used.• Effects of particle size distribution and particle shape for K0 samples: One of the keyhighlights of the current dissertation is that relative density can significantly alter the effectsof particle size distribution and particle shape on the cyclic response of isotropic granularsamples. Since anistropy can play a key role in shearing behavior and is more commonamong natural soil deposits, it is important to extend the current approach to K0 compressedsamples to further generalize the current findings.• Effects of particle size distribution on cyclic liquefaction resistance using Coupled CFDand DEM : The current research has used constant volume condition in order to replicateundrained behavior. However, since saturated granular materials are two phases and thereis often presence of air pockets even inside a fully saturated media. Hence, a more realisticorder will be to include the pore fluid behavior modeled into the system.• Effect of particle size distribution on cyclic liquefaction under cyclic triaxial and multi-directional loading: The current study liquefies the samples using unidirectional cyclicsimple shear loading. However, the observed results may differ depending on the type anddirection of loading, for example, true triaxial loading in which the loading is axial and theprincipal stress axes exhibit instantaneous rotation in a manner distinct from that observedin simple shear tests. Similarly, when multidirectional loading is used in the study, similarstatements can be made.107Appendix A: LIGGGHTS CODEThe primary code written for the current research using LIGGGHTS consists of three stages:insertion, compression, and shearing. The three files are organized as following three items:• INSERTIONprocessors 2 2 2echo botha tom sty le superquadr icatom modify map ar rayh a r d p a r t i c l e s yes # NOTE{ Had to add l i n e due to Y mod > 1e9 }newton o f fcommunicate s i n g l e ve l yesu n i t s s iboundary p p m############################################################ Del ta t###########################################################v a r i a b l e d t equal .0000002t imestep ${dt}v a r i a b l e nStepsRun equal 300/${dt}v a r i a b l e wr i teOutEvery equal 0 .005/${dt}v a r i a b l e thermofreq equal 0 .01 / ${dt}v a r i a b l e res ta r t Eve r y equal 100############################################################ def ine box reg ion###########################################################v a r i a b l e x l equal 0v a r i a b l e xh equal 0.20v a r i a b l e y l equal 0v a r i a b l e yh equal 0.20v a r i a b l e z l equal 0108v a r i a b l e zh equal 0.20######################## I n s e r t i o n reg ion ########################v a r i a b l e x l1 equal 0.01v a r i a b l e xh1 equal 0.19v a r i a b l e y l1 equal 0.01v a r i a b l e yh1 equal 0.19v a r i a b l e z l1 equal 0.01v a r i a b l e zh1 equal 0.19################################################### Creat ion o f the S imu la t ion & I n s e r t i o n box###################################################reg ion reg block ${ x l} ${xh} ${ y l} ${yh} ${ z l} ${zh} u n i t s boxreg ion reg1 block ${x l1} ${xh1} ${y l1} ${yh1} ${z l1} ${zh1} u n i t s boxcreate box 1 reg # Number o f atom types used ( check l a t e r i f i t s ok )# DEACTIVATE FOR RESTART########################neighbor 0.002 b inneigh modi fy delay 0#################################################### Ma te r i a l P rope r t i es######################### c o e f f i c i e n t R e s t i t u t i o nv a r i a b l e cf ww equal 0.9v a r i a b l e cf pw equal 0.9v a r i a b l e c f pp equal 0.9################################################## Phys ica l parameters ################################################################# c o e f f i c i e n t p a r t i c l e − p a r t i c l e F r i c t i o nv a r i a b l e cof pp equal 0.35# c o e f f i c i e n t R o l l i n g F r i c t i o nv a r i a b l e corf ww equal 0.0v a r i a b l e cor f pw equal 0.0v a r i a b l e co r f pp equal 0.0109# Young ’ s Modulusv a r i a b l e Yw equal 70e+9v a r i a b l e Yp equal 70e+9######################### Block iness ############################v a r i a b l e b lock iness1 equal 2v a r i a b l e b lock iness2 equal 2v a r i a b l e ar equal 1.0############################################################### Ma te r i a l p r o p e r t i e s assignment f o r new p a i r s t y l e s ####################################################################f i x m1 a l l p roper ty / g loba l youngsModulus peratomtype ${Yw}f i x m2 a l l p roper ty / g loba l poissonsRat io peratomtype 0.25f i x m3 a l l p roper ty / g loba l c o e f f i c i e n t R e s t i t u t i o n peratomtypepai r 1 ${c f pp}f i x m4 a l l p roper ty / g loba l c o e f f i c i e n t F r i c t i o n peratomtypepai r 1 ${cof pp}f i x m5 a l l p roper ty / g loba l c o e f f i c i e n t R o l l i n g F r i c t i o n peratomtypepai r 1 ${co r f pp}#######################Contact Law######################## Her tz ian w i thou t cohesionp a i r s t y l e gran model he r t z t a n g e n t i a l h i s t o r y r o l l i n g f r i c t i o n epsd2 sur face superquadr ic l i m i t F o r c e onp a i r c o e f f * *################################################################# P a r t i c l e Morphology ( Size , Shape and D i s t r i b u t i o n ) ####################################################################v a r i a b l e n r a d i i equal 10v a r i a b l e rad ius1 equal 0.005223v a r i a b l e radius1a equal ${ rad ius1 }*${ar}v a r i a b l e rad ius2 equal 0.004547v a r i a b l e radius2a equal ${ rad ius2 }*${ar}v a r i a b l e rad ius3 equal 0.003958v a r i a b l e radius3a equal ${ rad ius3 }*${ar}v a r i a b l e rad ius4 equal 0.003446v a r i a b l e radius4a equal ${ rad ius4 }*${ar}v a r i a b l e rad ius5 equal 0.003000v a r i a b l e radius5a equal ${ rad ius5 }*${ar}v a r i a b l e rad ius6 equal 0.002612v a r i a b l e radius6a equal ${ rad ius6 }*${ar}v a r i a b l e rad ius7 equal 0.002274110v a r i a b l e radius7a equal ${ rad ius7 }*${ar}v a r i a b l e rad ius8 equal 0.001979v a r i a b l e radius8a equal ${ rad ius8 }*${ar}v a r i a b l e rad ius9 equal 0.001723v a r i a b l e radius9a equal ${ rad ius9 }*${ar}v a r i a b l e radius10 equal 0.001500v a r i a b l e radius10a equal ${ radius10 }*${ar}v a r i a b l e f rac1 equal 0.1v a r i a b l e f rac2 equal 0.1v a r i a b l e f rac3 equal 0.1v a r i a b l e f rac4 equal 0.1v a r i a b l e f rac5 equal 0.1v a r i a b l e f rac6 equal 0.1v a r i a b l e f rac7 equal 0.1v a r i a b l e f rac8 equal 0.1v a r i a b l e f rac9 equal 0.1v a r i a b l e f rac10 equal 0.1##################################group pts11 i n i t i a l i z egroup pts12 i n i t i a l i z egroup pts13 i n i t i a l i z egroup pts14 i n i t i a l i z egroup pts15 i n i t i a l i z egroup pts16 i n i t i a l i z egroup pts17 i n i t i a l i z egroup pts18 i n i t i a l i z egroup pts19 i n i t i a l i z egroup pts20 i n i t i a l i z ef i x pts1 pts11 p a r t i c l e t e m p l a t e / superquadr ic 10487 atom type 1 dens i t y constant 2500shape constant ${ rad ius1} ${ rad ius1}${ radius1a} b lock iness constant ${block iness1} ${block iness2}f i x pts2 pts12 p a r t i c l e t e m p l a t e / superquadr ic 11887 atom type 1 dens i t y constant 2500shape constant ${ rad ius2} ${ rad ius2}${ radius2a} b lock iness constant ${block iness1} ${block iness2}f i x pts3 pts13 p a r t i c l e t e m p l a t e / superquadr ic 11923 atom type 1 dens i t y constant 2500shape constant ${ rad ius3} ${ rad ius3}${ radius3a} b lock iness constant ${block iness1} ${block iness2}111f i x pts4 pts14 p a r t i c l e t e m p l a t e / superquadr ic 12589 atom type 1 dens i t y constant 2500shape constant ${ rad ius4} ${ rad ius4}${ radius4a} b lock iness constant ${block iness1} ${block iness2}f i x pts5 pts15 p a r t i c l e t e m p l a t e / superquadr ic 12799 atom type 1 dens i t y constant 2500shape constant ${ rad ius5} ${ rad ius5}${ radius5a} b lock iness constant ${block iness1} ${block iness2}f i x pts6 pts16 p a r t i c l e t e m p l a t e / superquadr ic 12809 atom type 1 dens i t y constant 2500shape constant ${ rad ius6} ${ rad ius6} ${ radius6a} b lock iness constant ${block iness1} ${block iness2}f i x pts7 pts17 p a r t i c l e t e m p l a t e / superquadr ic 12821 atom type 1 dens i t y constant 2500shape constant ${ rad ius7} ${ rad ius7} ${ radius7a} b lock iness constant ${block iness1} ${block iness2}f i x pts8 pts18 p a r t i c l e t e m p l a t e / superquadr ic 12823 atom type 1 dens i t y constant 2500shape constant ${ rad ius8} ${ rad ius8} ${ radius8a} b lock iness constant ${block iness1} ${block iness2}f i x pts9 pts19 p a r t i c l e t e m p l a t e / superquadr ic 12829 atom type 1 dens i t y constant 2500shape constant ${ rad ius9} ${ rad ius9} ${ radius9a} b lock iness constant ${block iness1} ${block iness2}f i x pts10 pts20 p a r t i c l e t e m p l a t e / superquadr ic 12941 atom type 1 dens i t y constant 2500shape constant ${ radius10} ${ radius10}${ radius10a} b lock iness constant ${block iness1} ${block iness2}compute 1 a l l p roper ty / atom rad iusv a r i a b l e d ia equal c 1 [ 1 ]f i x pdd1 a l l p a r t i c l e d i s t r i b u t i o n / d i s c r e t e / massbased 86028157 ${ n r a d i i } pts1 ${ f r ac1} pts2${ f r ac2} pts3 ${ f rac3} pts4 ${ f r ac4} pts5 ${ f rac5} pts6 ${ f rac6} pts7 ${ f r ac7} pts8 ${ f rac8}pts9 ${ f rac9} pts10 ${ f rac10}####################################################################################### Proper ty o f the f i x e d lower boundary###########################################f i x zwal ls2 a l l wa l l / gran model he r t z t a n g e n t i a l h i s t o r y r o l l i n g f r i c t i o n epsd2 sur face superquadr icp r i m i t i v e type 1 zplane 0################################################### I n s e r t i o n o f p a r t i c l e s ##################################################### DEACTIVATE FOR RESTARTf i x ins a l l i n s e r t / pack seed 86028121 d i s t r i b u t i o n t e m p l a t e pdd1 maxattempt 100 o r i e n t a t i o n randomve l constant 0 0 0 i n s e r t e v e r y 1 over lapcheck yes a l l i n yes p a r t i c l e s i n r e g i o n 10000 reg ion reg1group nve group reg ion regf i x i n t e g r a l l nve / superquadr ic in tegra t ion scheme 1###########( Grav i t y i s i n a c t i v e . p a r t i c l e s drop wi th a uni form v e l o c i t y )f i x g rav i a l l g r a v i t y 0 vec to r 0.0 0.0 −1.0112v a r i a b l e t ime equal step *${dt}run 1group atom2 reg ion reg1compute p1 atom2 contac t / atom # CN step 1compute p2 atom2 reduce sum c p1 # CN step 2### number o f p a r t i c l e s = 10000 ######v a r i a b l e CN equal c p2 /10000# Stresses output as kPav a r i a b l e a1 equal ( c s t [ 1 ] / vo l /1000)v a r i a b l e a2 equal ( c s t [ 2 ] / vo l /1000)v a r i a b l e a3 equal ( c s t [ 3 ] / vo l /1000)v a r i a b l e a4 equal ( c s t [ 4 ] / vo l /1000)v a r i a b l e a5 equal ( c s t [ 5 ] / vo l /1000)v a r i a b l e a6 equal ( c s t [ 6 ] / vo l /1000)v a r i a b l e a7 equal c s t / vo l /1000v a r i a b l e a8 equal c p5s h e l l mkdir meshplate mod aa As h e l l mkdir outputf i x ex t ra a l l p r i n t 500 ” ${ t ime} ${CN} ${d ia }” f i l e output / i n s e r t i o n . t x t################################## i n s e r t the f i r s t p a r t i c l e s so t h a t dump i s not empty##############################the rmo s ty le custom step atoms cpu spcpu t ime vo lthermo ${ thermofreq}thermo modify l o s t ignore norm no####################v a r i a b l e dumptime equal 100s h e l l mkdir meshplate mod aa2b A scale#dump dmp1 atom2 custom / v tk ${dumptime} meshplate mod aa2b A scale / dump par t i c les *. superq . v t k type x y z i d vx vy vz f x f y f z rad ius shapex shapey shapez quat1 quat2 quat3 quat4b lock iness1 b lock iness2#dump dmp3 a l l mesh / v t k ${dumptime} meshplate mod aa2b A scale / dump * . mesh . v tk s t ress area cad cad1#dump dmp1 atom2 custom / v tk 50 meshplate mod aa A / dump par t i c les * . superq . v t k type x y z i dvx vy vz f x f y f z rad ius shapex shapey shapez quat1 quat2 quat3 quat4 b lock iness1 b lock iness2#dump dmp1a pts11 custom / v tk 50 meshplate mod aa A / dump par t i c les 11p ts * . superq . v t k type x y z i d113vx vy vz f x f y f z rad ius shapex shapey shapez quat1 quat2 quat3 quat4 b lock iness1 b lock iness2#dump dmp1b pts12 custom / v tk 50 meshplate mod aa A / dump par t i c les 12p ts * . superq . v t k type x y z i dvx vy vz f x f y f z rad ius shapex shapey shapez quat1 quat2 quat3 quat4 b lock iness1 b lock iness2#dump dmp1c pts13 custom / v tk 50 meshplate mod aa A / dump par t i c les 13p ts * . superq . v t k type x y z i dvx vy vz f x f y f z rad ius shapex shapey shapez quat1 quat2 quat3 quat4 b lock iness1 b lock iness2#dump dmp1d pts14 custom / v tk 50 meshplate mod aa A / dump par t i c les 14p ts * . superq . v t k type x y z i dvx vy vz f x f y f z rad ius shapex shapey shapez quat1 quat2 quat3 quat4 b lock iness1 b lock iness2#dump dmp1e pts15 custom / v tk 50 meshplate mod aa A / dump par t i c les 15p ts * . superq . v t k type x y z i dvx vy vz f x f y f z rad ius shapex shapey shapez quat1 quat2 quat3 quat4 b lock iness1 b lock iness2#dump dmp1f pts16 custom / v tk 50 meshplate mod aa A / dump par t i c les 16p ts * . superq . v t k type x y z i dvx vy vz f x f y f z rad ius shapex shapey shapez quat1 quat2 quat3 quat4 b lock iness1 b lock iness2#dump dmp1g pts17 custom / v tk 50 meshplate mod aa A / dump par t i c les 17p ts * . superq . v t k type x y z i dvx vy vz f x f y f z rad ius shapex shapey shapez quat1 quat2 quat3 quat4 b lock iness1 b lock iness2#dump dmp1h pts18 custom / v tk 50 meshplate mod aa A / dump par t i c les 18p ts * . superq . v t k type x y z i dvx vy vz f x f y f z rad ius shapex shapey shapez quat1 quat2 quat3 quat4 b lock iness1 b lock iness2#dump dmp1i pts19 custom / v tk 50 meshplate mod aa A / dump par t i c les 19p ts * . superq . v t k type x y z i dvx vy vz f x f y f z rad ius shapex shapey shapez quat1 quat2 quat3 quat4 b lock iness1 b lock iness2#dump dmp1j pts20 custom / v tk 50 meshplate mod aa A / dump par t i c les 20p ts * . superq . v t k type x y z i dvx vy vz f x f y f z rad ius shapex shapey shapez quat1 quat2 quat3 quat4 b lock iness1 b lock iness2s h e l l mkdir restar t meshplate mod aa Ar e s t a r t ${ r e s ta r t Eve r y} restar t meshplate mod aa A / p1 restar t meshplate mod aa A / p2run 1000*${dt} upto # DEACTIVATE FOR RESTARTu n f i x insrun 1• COMPRESSION#processors 2 3 1atom sty le superquadr icatom modify map ar rayh a r d p a r t i c l e s yes # NOTE{ Had to add l i n e due to Y mod > 1e9 }newton o f fcommunicate s i n g l e ve l yesu n i t s s i############################################################ Del ta t###########################################################v a r i a b l e d t equal .0000002t imestep ${dt}114v a r i a b l e nStepsRun equal 300/${dt}v a r i a b l e wr i teOutEvery equal 0 .005/${dt}v a r i a b l e thermofreq equal 0 .01 / ${dt}v a r i a b l e res ta r t Eve r y equal 100r e a d r e s t a r t compression2 / p1############################################################ def ine box reg ion###########################################################v a r i a b l e x l equal 0v a r i a b l e xh equal 0.18v a r i a b l e y l equal 0v a r i a b l e yh equal 0.18v a r i a b l e z l equal 0v a r i a b l e zh equal 0.18######################## I n s e r t i o n reg ion ########################v a r i a b l e x l1 equal 0.01v a r i a b l e xh1 equal 0.17v a r i a b l e y l1 equal 0.01v a r i a b l e yh1 equal 0.17v a r i a b l e z l1 equal 0.01v a r i a b l e zh1 equal 0.17################################################### Creat ion o f the S imu la t ion & I n s e r t i o n box ###########################################################reg ion reg block ${ x l} ${xh} ${ y l} ${yh} ${ z l} ${zh} u n i t s box########################neighbor 0.002 b inneigh modi fy delay 0################################################115############################ Ma te r i a l P rope r t i es######################### c o e f f i c i e n t R e s t i t u t i o nv a r i a b l e cf ww equal 1.0v a r i a b l e cf pw equal 0.2v a r i a b l e c f pp equal 0.2# c o e f f i c i e n t F r i c t i o nv a r i a b l e cof pp equal 0.35 # p a r t i c l e − p a r t i c l e f r i c t i o n# c o e f f i c i e n t R o l l i n g F r i c t i o nv a r i a b l e corf ww equal 0.0v a r i a b l e cor f pw equal 0.0v a r i a b l e co r f pp equal 0.1# Young ’ s Modulusv a r i a b l e Yw equal 70e+9v a r i a b l e Yp equal 70e+9#####################v a r i a b l e b lock iness1 equal 2.0v a r i a b l e b lock iness2 equal 2.0v a r i a b l e ar equal 1.0############################################################### Ma te r i a l p r o p e r t i e s assignment f o r new p a i r s t y l e s#########################################################f i x m1 a l l p roper ty / g loba l youngsModulus peratomtype ${Yw}f i x m2 a l l p roper ty / g loba l poissonsRat io peratomtype 0.25f i x m3 a l l p roper ty / g loba l c o e f f i c i e n t R e s t i t u t i o n peratomtypepai r 1 ${c f pp}f i x m4 a l l p roper ty / g loba l c o e f f i c i e n t F r i c t i o n peratomtypepai r 1 ${cof pp}f i x m5 a l l p roper ty / g loba l c o e f f i c i e n t R o l l i n g F r i c t i o n peratomtypepai r 1 ${co r f pp}echo both######################Contact Law######################## Her tz ian w i thou t cohesionp a i r s t y l e gran model he r t z t a n g e n t i a l h i s t o r y r o l l i n g f r i c t i o n epsd2 sur face superquadr ic l i m i t F o r c e onp a i r c o e f f * *################################################################# P a r t i c l e Morphology ( Size , Shape and D i s t r i b u t i o n )#########################################################v a r i a b l e n r a d i i equal 10116v a r i a b l e rad ius1 equal 0.005223v a r i a b l e radius1a equal ${ rad ius1 }*${ar}v a r i a b l e rad ius2 equal 0.004547v a r i a b l e radius2a equal ${ rad ius2 }*${ar}v a r i a b l e rad ius3 equal 0.003958v a r i a b l e radius3a equal ${ rad ius3 }*${ar}v a r i a b l e rad ius4 equal 0.003446v a r i a b l e radius4a equal ${ rad ius4 }*${ar}v a r i a b l e rad ius5 equal 0.003000v a r i a b l e radius5a equal ${ rad ius5 }*${ar}v a r i a b l e rad ius6 equal 0.002612v a r i a b l e radius6a equal ${ rad ius6 }*${ar}v a r i a b l e rad ius7 equal 0.002274v a r i a b l e radius7a equal ${ rad ius7 }*${ar}v a r i a b l e rad ius8 equal 0.001979v a r i a b l e radius8a equal ${ rad ius8 }*${ar}v a r i a b l e rad ius9 equal 0.001723v a r i a b l e radius9a equal ${ rad ius9 }*${ar}v a r i a b l e radius10 equal 0.001500v a r i a b l e radius10a equal ${ radius10 }*${ar}v a r i a b l e ravg equal ( ${ rad ius5}+${ rad ius6 } ) / 2# v a r i a b l e ravg equal ${ rad ius5}v a r i a b l e f rac1 equal 0.1v a r i a b l e f rac2 equal 0.1v a r i a b l e f rac3 equal 0.1v a r i a b l e f rac4 equal 0.1v a r i a b l e f rac5 equal 0.1v a r i a b l e f rac6 equal 0.1v a r i a b l e f rac7 equal 0.1v a r i a b l e f rac8 equal 0.1v a r i a b l e f rac9 equal 0.1v a r i a b l e f rac10 equal 0.1##################################f i x pts1 pts11 p a r t i c l e t e m p l a t e / superquadr ic 10487 atom type 1 dens i t y constant 2500 shape constant${ rad ius1} ${ rad ius1} ${ radius1a} b lock iness constant ${block iness1} ${block iness2}f i x pts2 pts12 p a r t i c l e t e m p l a t e / superquadr ic 11887 atom type 1 dens i t y constant 2500 shape constant${ rad ius2} ${ rad ius2} ${ radius2a} b lock iness constant ${block iness1} ${block iness2}f i x pts3 pts13 p a r t i c l e t e m p l a t e / superquadr ic 11923 atom type 1 dens i t y constant 2500 shape constant117${ rad ius3} ${ rad ius3} ${ radius3a} b lock iness constant ${block iness1} ${block iness2}f i x pts4 pts14 p a r t i c l e t e m p l a t e / superquadr ic 12589 atom type 1 dens i t y constant 2500 shape constant${ rad ius4} ${ rad ius4} ${ radius4a} b lock iness constant ${block iness1} ${block iness2}f i x pts5 pts15 p a r t i c l e t e m p l a t e / superquadr ic 12799 atom type 1 dens i t y constant 2500 shape constant${ rad ius5} ${ rad ius5} ${ radius5a} b lock iness constant ${block iness1} ${block iness2}f i x pts6 pts16 p a r t i c l e t e m p l a t e / superquadr ic 12809 atom type 1 dens i t y constant 2500 shape constant${ rad ius6} ${ rad ius6} ${ radius6a} b lock iness constant ${block iness1} ${block iness2}f i x pts7 pts17 p a r t i c l e t e m p l a t e / superquadr ic 12821 atom type 1 dens i t y constant 2500 shape constant${ rad ius7} ${ rad ius7} ${ radius7a} b lock iness constant ${block iness1} ${block iness2}f i x pts8 pts18 p a r t i c l e t e m p l a t e / superquadr ic 12823 atom type 1 dens i t y constant 2500 shape constant${ rad ius8} ${ rad ius8} ${ radius8a} b lock iness constant ${block iness1} ${block iness2}f i x pts9 pts19 p a r t i c l e t e m p l a t e / superquadr ic 12829 atom type 1 dens i t y constant 2500 shape constant${ rad ius9} ${ rad ius9} ${ radius9a} b lock iness constant ${block iness1} ${block iness2}f i x pts10 pts20 p a r t i c l e t e m p l a t e / superquadr ic 12941 atom type 1 dens i t y constant 2500 shape constant${ radius10} ${ radius10} ${ radius10a} b lock iness constant ${block iness1} ${block iness2}#####################################f i x pdd1 a l l p a r t i c l e d i s t r i b u t i o n / d i s c r e t e / massbased 86028157 ${ n r a d i i } pts1 ${ f r ac1} pts2 ${ f rac2}pts3 ${ f rac3} pts4 ${ f r ac4} pts5 ${ f rac5} pts6 ${ f r ac6} pts7 ${ f rac7} pts8 ${ f r ac8} pts9 ${ f rac9}pts10 ${ f rac10}#######v a r i a b l e s t ress equal 10000v a r i a b l e fo rce equal ${s t ress }* l x * l y####################################################################################### Proper ty o f the f i x e d lower boundary###########################################compute p6 atom2 p a i r / gran / l o c a l pos fo rcecompute p7 atom2 proper ty / atom update on run end yes x y z i d vx vy vz f x f y f zcompute p7a atom2 reduce update on run end yes min c p7 [ 1 ]compute p7b atom2 reduce update on run end yes min c p7 [ 2 ]compute p7c atom2 reduce update on run end yes min c p7 [ 3 ]compute p7d atom2 reduce update on run end yes max c p7 [ 1 ]compute p7e atom2 reduce update on run end yes max c p7 [ 2 ]compute p7f atom2 reduce update on run end yes max c p7 [ 3 ]118#minv a r i a b l e a14 equal c p7av a r i a b l e a15 equal c p7bv a r i a b l e a16 equal c p7c#maxv a r i a b l e a17 equal c p7dv a r i a b l e a18 equal c p7ev a r i a b l e a19 equal c p7 fv a r i a b l e x d i f f equal c p7d −c p7av a r i a b l e x r a t i o equal ( c p7d −c p7a+${ ravg } ) / 0 .15# v a r i a b l e x r a t i o equal ( c p7e −c p7b ) / 0 . 1 5v a r i a b l e y d i f f equal c p7e −c p7bv a r i a b l e y r a t i o equal ( c p7e −c p7b+${ ravg } ) / 0 .15v a r i a b l e z d i f f equal c p7f −c p7cv a r i a b l e z r a t i o equal ( c p7f −c p7c ) / 0 . 1 5# v a r i a b l e sysvo l equal ( c p7d −c p7a ) * ( c p7e −c p7b ) * ( c p7f −c p7c )v a r i a b l e sysvo l equal l x * l y * ( c p7f −c p7c+v ravg )v a r i a b l e x r equal ( l z / ( c p7f −c p7c ) ) * ( l z / ( c p7d −c p7a ) ) * ( l z / ( c p7e −c p7b ) )############################################## Proper ty o f the f i x e d lower boundary###########################################compute p1 atom2 contac t / atom # CN step 1 # THIS IS SHOULD BE ALWAYS BEFORE RUN 1f i x cad a l l mesh / sur face f i l e mesh / p l a t e z h i g h . s t l type 1 scale 0.99f i x cad1 a l l mesh / sur face f i l e mesh / p la te z low . s t l type 1 scale 0.99run 1v a r i a b l e s t ress equal 10000v a r i a b l e area equal l x * l y########################################compute p5 atom2 kev a r i a b l e a8 equal c p5 # K i n e t i c energyv a r i a b l e p a r t v o l equal mass ( atom2 )/2500v a r i a b l e vo id equal ( l x * l y * ( c p7f −c p7c+${ ravg})−${p a r t v o l } ) / ${p a r t v o l}# v a r i a b l e vo id equal ( l x * l y * ( c p7f −c p7c+${ ravg }) −0.000616)/0.000616v a r i a b l e t ime equal step *${dt}compute p2 atom2 reduce sum c p1 # CN step 2v a r i a b l e CN equal c p2 /10000119##### Ca l cu l a t i on o f the s t ress tensor o f the assembly #########Samplev a r i a b l e a14min equal ${a14}#+${ ravg}*5v a r i a b l e a14max equal ${a17}#−${ ravg}*5v a r i a b l e a15min equal ${a15}#+${ ravg}*5v a r i a b l e a15max equal ${a18}#−${ ravg}*5v a r i a b l e a16min equal ${a16}#+${ ravg}*5v a r i a b l e a16max equal ${a19}#−${ ravg}*5reg ion sample block ${a14min} ${a14max} ${a15min} ${a15max} ${a16min} ${a16max} u n i t s boxgroup sample1 reg ion sample#Lower p l a tev a r i a b l e a16a equal c p7c+${ ravg}reg ion lowerp la te b lock ${a14} ${a17} ${a15} ${a18} ${a16} ${a16a} u n i t s boxgroup atom reg ion lowerp la teset group atom block iness ${block iness1} ${block iness1}##### Ca l cu l a t i on o f the s t ress tensor using the lower p l a t ecompute p7lp atom2 proper ty / atom update on run end yes x y z i d vx vy vz f x f y f zcompute p7clp atom2 reduce update on run end yes min c p7 lp [ 3 ]compute p7 f l p atom2 reduce update on run end yes max c p7 lp [ 3 ]compute p8lp atom2 reduce update on run end yes sum c p7 lp [ 1 0 ]compute p9lp atom2 reduce update on run end yes sum c p7 lp [ 9 ]compute p10lp atom2 reduce update on run end yes sum c p7 lp [ 8 ]v a r i a b l e a20 equal −c p8 lp * v x r / v area /1000v a r i a b l e a21 equal c p9 lp * v x r / v area /1000v a r i a b l e a22 equal c p10lp * v x r / v area /1000compute peratom atom2 s t ress / atom v i r i a lcomputep atom2 reduce sum c peratom [ 1 ] c peratom [ 2 ] c peratom [ 3 ] c peratom [ 4 ] c peratom [ 5 ] c peratom [ 6 ]compute pvol a l l pressure thermo temp v i r i a lv a r i a b l e sig11 equal −c p [ 1 ] * 1 / v sysvo lv a r i a b l e sig22 equal −c p [ 2 ] * 1 / v sysvo lv a r i a b l e sig33 equal −c p [ 3 ] * 1 / v sysvo lv a r i a b l e sig12 equal −c p [ 4 ] * 1 / v sysvo lv a r i a b l e sig13 equal −c p [ 5 ] * 1 / v sysvo lv a r i a b l e sig23 equal −c p [ 6 ] * 1 / v sysvo l120v a r i a b l e pbvol equal −c pvo l *1* vo l / v sysvo l /1000v a r i a b l e sig11a equal −c p [ 1 ] * 1 / v sysvo l /1000v a r i a b l e sig22a equal −c p [ 2 ] * 1 / v sysvo l /1000v a r i a b l e sig33a equal −c p [ 3 ] * 1 / v sysvo l /1000v a r i a b l e sig12a equal −c p [ 4 ] * 1 / v sysvo l /1000v a r i a b l e sig13a equal −c p [ 5 ] * 1 / v sysvo l /1000v a r i a b l e sig23a equal −c p [ 6 ] * 1 / v sysvo l /1000#Upper and lower p a r t i c l e p la tes# Shearing p la te a lso Upper p l a t ev a r i a b l e a19a equal c p7f −${ ravg}reg ion upperp la te b lock ${a14} ${a17} ${a15} ${a18} ${a19a} ${a19} u n i t s boxgroup atom1 reg ion upperp la terun 1# Inner p a r t i c l e sv a r i a b l e a14a equal c p7a+${ ravg}v a r i a b l e a15a equal c p7b+${ ravg}v a r i a b l e a17a equal c p7d −${ ravg}v a r i a b l e a18a equal c p7e −${ ravg}reg ion p a r t i c l e b lock ${a14a} ${a17a} ${a15a} ${a18a} ${a16a} ${a19a} u n i t s boxgroup a l la tom reg ion p a r t i c l ev a r i a b l e p l t v e l equal 0# f i x 123456 atom1 move l i n e a r 0.000000000000001 NULL NULL u n i t s box# Always use the f i x a f t e r ” run 1”f i x wa l l1 a l l wa l l / gran model he r t z t a n g e n t i a l h i s t o r y r o l l i n g f r i c t i o n epsd2 sur face superquadr icmesh n meshes 2 meshes cad cad1 #v a r i a b l e s t ress1 equal ${s t ress } /1s h e l l mkdir meshplate mod aa2b A scale 3v a r i a b l e dumptime equal 1000#dump dmp1 atom2 custom / v tk ${dumptime} meshplate mod aa2b A scale 3 / dump par t i c les * . superq . v t k type x y zi d vx vy vz f x f y f z rad ius shapex shapey shapez quat1 quat2 quat3 quat4 b lock iness1 block iness2121#dump dmp3 a l lmesh / v t k ${dumptime} meshplate mod aa2b A scale 3 / dump * . mesh . v tk s t ress area cad cad1#dump dmp4 atom2 l o c a l ${dumptime} meshplate mod aa2b A scale 3 / dump * . fo rce . t x t c p6 [ 1 ] c p6 [ 2 ] c p6 [ 3 ]c p6 [ 4 ] c p6 [ 5 ] c p6 [ 6 ] c p6 [ 7 ] c p6 [ 8 ] c p6 [ 9 ]v a r i a b l e alpha equal 6e−7v a r i a b l e ve lx equal −${alpha } * ( ${s t ress1}−v s ig11 ) *2v a r i a b l e ve ly equal −${alpha } * ( ${s t ress1}−v s ig22 ) *2v a r i a b l e ve lz equal −${alpha } * ( ${s t ress1}−v s ig33 )v a r i a b l e ve lza equal ${alpha } * ( ${s t ress1}−v s ig33 )the rmo s ty le custom step v s ig11 v s ig22 v s ig33 v pbvo l v ve l x v ve l y v ve l z v ve lza v a8 v vo id v CNthermo modify l o s t warn norm norun 0v a r i a b l e t o t a l s t e p equal 200000v a r i a b l e loop s tep equal 1000v a r i a b l e loop count max equal ${ t o t a l s t e p } / ${ l oop s tep}v a r i a b l e loop count loop ${ loop count max}s h e l l mkdir compression3r e s t a r t 100 compression3 / p1 compression3 / p2 #DEACTIVATE FOR RESTARTf i x i n t e g r a l l nve / superquadr ic in tegra t ion scheme 1l a b e l loopf i x 1234 atom move l i n e a r NULL NULL NULL u n i t s box # f i x i n g the lower p l a t ef i x p la te1 a l l move / mesh mesh cad l i n e a r 0 . 0 . ${ve lz} #−1#f i x p la te2 a l l move / mesh mesh cad1 l i n e a r 0 . 0 . ${velza} ## 1f i x 2 a l l deform 1 x ve l ${ve lx} y ve l ${ve ly} u n i t s box remap none#remap v#dump dmp4 atom2 l o c a l ${dumptime} meshplate mod aa2b A / dump * . fo rce . t x t c p6 [ 1 ] c p6 [ 2 ] c p6 [ 3 ] c p6 [ 4 ] c p6 [ 5 ]c p6 [ 6 ] c p6 [ 7 ] c p6 [ 8 ] c p6 [ 9 ] c p6 [ 1 0 ] c p6 [ 1 1 ] c p6 [ 1 2 ] c p6 [ 1 3 ] c p6 [ 1 4 ] c p6 [ 1 5 ] c p6 [ 1 6 ] c p6 [ 1 7 ] c p6 [ 1 8 ]c p6 [ 1 9 ] c p6 [ 2 0 ] c p6 [ 2 1 ] c p6 [ 2 2 ] c p6 [ 2 3 ] c p6 [ 2 4 ] c p6 [ 2 5 ] c p6 [ 2 6 ] c p6 [ 2 7 ] c p6 [ 2 8 ]c p6 [ 2 9 ] c p6 [ 3 0 ] c p6 [ 3 1 ]run ${ l oop s tep}p r i n t ” ${ t ime} ${CN} ${vo id}${a8} ${sysvo l} ${ve lx} ${ve ly} ${ve lz} ${velza} ${sig11a} ${sig22a} ${sig33a}${pbvol} ${a20} ${a21} ${a22}” append output / compress3 . t x tnext loop countjump SELF looprun 2e6• SHEARING#processors 2 2 1122atom modify map ar rayh a r d p a r t i c l e s yes # NOTE{ Had to add l i n e due to Y mod > 1e9 }newton o f fcommunicate s i n g l e ve l yesu n i t s s i############################################################ Del ta t###########################################################v a r i a b l e d t equal .0000002t imestep ${dt}v a r i a b l e nStepsRun equal 300/${dt}v a r i a b l e wr i teOutEvery equal 0 .005/${dt}v a r i a b l e thermofreq equal 0 .01 / ${dt}v a r i a b l e res ta r t Eve r y equal 100r e a d r e s t a r t compression6 / p1############################################################ def ine box reg ion###########################################################v a r i a b l e x l equal 0v a r i a b l e xh equal 0.15v a r i a b l e y l equal 0v a r i a b l e yh equal 0.15v a r i a b l e z l equal 0v a r i a b l e zh equal 0.15######################## I n s e r t i o n reg ion ########################v a r i a b l e x l1 equal 0.01v a r i a b l e xh1 equal 0.14v a r i a b l e y l1 equal 0.01v a r i a b l e yh1 equal 0.14v a r i a b l e z l1 equal 0.01v a r i a b l e zh1 equal 0.14################################################### Creat ion o f the S imu la t ion & I n s e r t i o n box###################################################reg ion reg block ${ x l} ${xh} ${ y l} ${yh} ${ z l} ${zh} u n i t s box####################123########################neighbor 0.002 b inneigh modi fy delay 0############################################################################ Ma te r i a l P rope r t i es######################### c o e f f i c i e n t R e s t i t u t i o nv a r i a b l e cf ww equal 1.0v a r i a b l e cf pw equal 0.2v a r i a b l e c f pp equal 0.2# c o e f f i c i e n t F r i c t i o nv a r i a b l e cof pp equal 0.5 # p a r t i c l e − p a r t i c l e f r i c t i o n# c o e f f i c i e n t R o l l i n g F r i c t i o nv a r i a b l e corf ww equal 0.0v a r i a b l e cor f pw equal 0.0v a r i a b l e co r f pp equal 0.1# Young ’ s Modulusv a r i a b l e Yw equal 70e+9v a r i a b l e Yp equal 70e+9#####################v a r i a b l e b lock iness1 equal 2.0v a r i a b l e b lock iness2 equal 2.0v a r i a b l e ar equal 1.0############################################################### Ma te r i a l p r o p e r t i e s assignment f o r new p a i r s t y l e s#########################################################f i x m1 a l l p roper ty / g loba l youngsModulus peratomtype ${Yw}f i x m2 a l l p roper ty / g loba l poissonsRat io peratomtype 0.25f i x m3 a l l p roper ty / g loba l c o e f f i c i e n t R e s t i t u t i o n peratomtypepai r 1 ${c f pp}f i x m4 a l l p roper ty / g loba l c o e f f i c i e n t F r i c t i o n peratomtypepai r 1 ${cof pp}f i x m5 a l l p roper ty / g loba l c o e f f i c i e n t R o l l i n g F r i c t i o n peratomtypepai r 1 ${co r f pp}echo both######################Contact Law######################## Her tz ian w i thou t cohesionp a i r s t y l e gran model he r t z t a n g e n t i a l h i s t o r y r o l l i n g f r i c t i o n epsd2 sur face superquadr icp a i r c o e f f * *124################################################################# P a r t i c l e Morphology ( Size , Shape and D i s t r i b u t i o n )#########################################################v a r i a b l e n r a d i i equal 10v a r i a b l e rad ius1 equal 0.005223v a r i a b l e radius1a equal ${ rad ius1 }*${ar}v a r i a b l e rad ius2 equal 0.004547v a r i a b l e radius2a equal ${ rad ius2 }*${ar}v a r i a b l e rad ius3 equal 0.003958v a r i a b l e radius3a equal ${ rad ius3 }*${ar}v a r i a b l e rad ius4 equal 0.003446v a r i a b l e radius4a equal ${ rad ius4 }*${ar}v a r i a b l e rad ius5 equal 0.003000v a r i a b l e radius5a equal ${ rad ius5 }*${ar}v a r i a b l e rad ius6 equal 0.002612v a r i a b l e radius6a equal ${ rad ius6 }*${ar}v a r i a b l e rad ius7 equal 0.002274v a r i a b l e radius7a equal ${ rad ius7 }*${ar}v a r i a b l e rad ius8 equal 0.001979v a r i a b l e radius8a equal ${ rad ius8 }*${ar}v a r i a b l e rad ius9 equal 0.001723v a r i a b l e radius9a equal ${ rad ius9 }*${ar}v a r i a b l e radius10 equal 0.001500v a r i a b l e radius10a equal ${ radius10 }*${ar}v a r i a b l e ravg equal ( ${ rad ius5}+${ rad ius6 } ) / 2# v a r i a b l e ravg equal ${ rad ius5}v a r i a b l e f rac1 equal 0.1v a r i a b l e f rac2 equal 0.1v a r i a b l e f rac3 equal 0.1v a r i a b l e f rac4 equal 0.1v a r i a b l e f rac5 equal 0.1v a r i a b l e f rac6 equal 0.1v a r i a b l e f rac7 equal 0.1v a r i a b l e f rac8 equal 0.1v a r i a b l e f rac9 equal 0.1v a r i a b l e f rac10 equal 0.1##################################125f i x pts1 pts11 p a r t i c l e t e m p l a t e / superquadr ic 10487 atom type 1 dens i t y constant 2500shape constant ${ rad ius1} ${ rad ius1} ${ radius1a} b lock iness constant ${block iness1} ${block iness2}f i x pts2 pts12 p a r t i c l e t e m p l a t e / superquadr ic 11887 atom type 1 dens i t y constant 2500shape constant ${ rad ius2} ${ rad ius2} ${ radius2a} b lock iness constant ${block iness1} ${block iness2}f i x pts3 pts13 p a r t i c l e t e m p l a t e / superquadr ic 11923 atom type 1 dens i t y constant 2500shape constant ${ rad ius3} ${ rad ius3} ${ radius3a} b lock iness constant ${block iness1} ${block iness2}f i x pts4 pts14 p a r t i c l e t e m p l a t e / superquadr ic 12589 atom type 1 dens i t y constant 2500shape constant ${ rad ius4} ${ rad ius4} ${ radius4a} b lock iness constant ${block iness1} ${block iness2}f i x pts5 pts15 p a r t i c l e t e m p l a t e / superquadr ic 12799 atom type 1 dens i t y constant 2500shape constant ${ rad ius5} ${ rad ius5} ${ radius5a} b lock iness constant ${block iness1} ${block iness2}f i x pts6 pts16 p a r t i c l e t e m p l a t e / superquadr ic 12809 atom type 1 dens i t y constant 2500shape constant ${ rad ius6} ${ rad ius6} ${ radius6a} b lock iness constant ${block iness1} ${block iness2}f i x pts7 pts17 p a r t i c l e t e m p l a t e / superquadr ic 12821 atom type 1 dens i t y constant 2500shape constant ${ rad ius7} ${ rad ius7} ${ radius7a} b lock iness constant ${block iness1} ${block iness2}f i x pts8 pts18 p a r t i c l e t e m p l a t e / superquadr ic 12823 atom type 1 dens i t y constant 2500shape constant ${ rad ius8} ${ rad ius8} ${ radius8a} b lock iness constant ${block iness1} ${block iness2}f i x pts9 pts19 p a r t i c l e t e m p l a t e / superquadr ic 12829 atom type 1 dens i t y constant 2500shape constant ${ rad ius9} ${ rad ius9} ${ radius9a} b lock iness constant ${block iness1} ${block iness2}f i x pts10 pts20 p a r t i c l e t e m p l a t e / superquadr ic 12941 atom type 1 dens i t y constant 2500shape constant ${ radius10} ${ radius10} ${ radius10a} b lock iness constant ${block iness1} ${block iness2}#####################################f i x pdd1 a l l p a r t i c l e d i s t r i b u t i o n / d i s c r e t e / massbased 86028157 ${ n r a d i i } pts1 ${ f r ac1} pts2 ${ f rac2}pts3 ${ f rac3} pts4 ${ f r ac4} pts5 ${ f rac5} pts6 ${ f r ac6} pts7 ${ f rac7} pts8 ${ f r ac8} pts9 ${ f rac9}pts10 ${ f rac10}#######v a r i a b l e s t ress equal 50000v a r i a b l e fo rce equal ${s t ress }* l x * l y# f i x p la te1 a l l mesh / sur face / s t ress / servo f i l e mesh / p l a t e z h i g h . s t l type 1 scale .999 com 0.075 0.075 0.1475c t r lPV fo rce ax is 0 . 0 . −1. t a r g e t v a l ${ f o rce} vel max 1 kp 5 # normal s t ress = 100 kPa# f i x p la te2 a l l mesh / sur face / s t ress / servo f i l e mesh / p la te z low . s t l type 1 scale .999 com 0.075 0.075 0.000c t r lPV fo rce ax is 0 . 0 . 1 . t a r g e t v a l ${ f o rce} vel max 0.065 kp 5 # normal s t ress = 100 kPa####################################################################################### Proper ty o f the f i x e d lower boundary###########################################compute p6 atom2 p a i r / gran / l o c a l i d de l t a pos ve l con tac tPo in t fo rce force normal f o r c e t a n g e n t i a l torquecompute p7 atom2 proper ty / atom update on run end yes x y z i d vx vy vz f x f y f zcompute p7a atom2 reduce update on run end yes min c p7 [ 1 ]compute p7b atom2 reduce update on run end yes min c p7 [ 2 ]126compute p7c atom2 reduce update on run end yes min c p7 [ 3 ]compute p7d atom2 reduce update on run end yes max c p7 [ 1 ]compute p7e atom2 reduce update on run end yes max c p7 [ 2 ]compute p7f atom2 reduce update on run end yes max c p7 [ 3 ]#minv a r i a b l e a14 equal c p7av a r i a b l e a15 equal c p7bv a r i a b l e a16 equal c p7c#maxv a r i a b l e a17 equal c p7dv a r i a b l e a18 equal c p7ev a r i a b l e a19 equal c p7 fv a r i a b l e x d i f f equal c p7d −c p7av a r i a b l e x r a t i o equal ( c p7d −c p7a+${ ravg } ) / 0 .15# v a r i a b l e x r a t i o equal ( c p7e −c p7b ) / 0 . 1 5v a r i a b l e y d i f f equal c p7e −c p7bv a r i a b l e y r a t i o equal ( c p7e −c p7b+${ ravg } ) / 0 .15v a r i a b l e z d i f f equal c p7f −c p7cv a r i a b l e z r a t i o equal ( c p7f −c p7c ) / 0 . 1 5# v a r i a b l e sysvo l equal ( c p7d −c p7a ) * ( c p7e −c p7b ) * ( c p7f −c p7c )v a r i a b l e sysvo l equal l x * l y * ( c p7f −c p7c+v ravg )v a r i a b l e x r equal ( l z / ( c p7f −c p7c ) ) * ( l z / ( c p7d −c p7a ) ) * ( l z / ( c p7e −c p7b ) )############################################## Proper ty o f the f i x e d lower boundary###########################################v a r i a b l e fo rce equal 969.577903917179compute p1 atom2 contac t / atom # CN step 1 # THIS IS SHOULD BE ALWAYS BEFORE RUN 1f i x cad a l l mesh / sur face f i l e mesh / p l a t e z h i g h . s t l type 1 scale 0.99f i x cad1 a l l mesh / sur face f i l e mesh / p la te z low . s t l type 1 scale 0.99run 1v a r i a b l e s t ress equal 50000v a r i a b l e area equal l x * l y########################################compute p5 atom2 kev a r i a b l e a8 equal c p5 # K i n e t i c energyv a r i a b l e p a r t v o l equal mass ( atom2 )/2500v a r i a b l e vo id equal ( l x * l y * ( c p7f −c p7c+${ ravg})−${p a r t v o l } ) / ${p a r t v o l}v a r i a b l e t ime equal step *${dt}compute p2 atom2 reduce sum c p1 # CN step 2127v a r i a b l e CN equal c p2 /10000#############################################################################################################################compute 2 atom2 reduce sum c p1 #” c p1 < 1 .0 ”v a r i a b l e p n a l l equal count ( atom2 )v a r i a b l e coun equal c 2the rmo s ty le custom step v counthermo modify l o s t warn norm norun 0v a r i a b l e eatom0 atom ” c p1 < 1 .0 ”group par t0 v a r i a b l e eatom0v a r i a b l e pn0 equal count ( par t0 )v a r i a b l e f rac0 equal ${pn0} /1v a r i a b l e eatom1 atom ” c p1 < 2 .0 ”group par t1 v a r i a b l e eatom1v a r i a b l e pn1 equal count ( par t1 )v a r i a b l e f rac1 equal ${pn1} /1v a r i a b l e eatom2 atom ” c p1 < 3 .0 ”group par t2 v a r i a b l e eatom2v a r i a b l e pn2 equal count ( par t2 )v a r i a b l e f rac2 equal ${pn2} /1v a r i a b l e eatom3 atom ” c p1 > 3 .0 ”group par t3 v a r i a b l e eatom3v a r i a b l e pn3 equal count ( par t3 )v a r i a b l e f rac3 equal ${pn3} /1######## zg c a l c u l a t i o n ###########v a r i a b l e eatzg atom ” c p1 > 0 .0 ”group par tzg v a r i a b l e eatzgcompute cont1g par tzg reduce sum c p1v a r i a b l e cont2g equal c cont1gv a r i a b l e png equal count ( par tzg )v a r i a b l e zg equal c cont1g / ${png}######## zm c a l c u l a t i o n ###########v a r i a b l e eatzm atom ” c p1 > 1 .0 ”128group partzm v a r i a b l e eatzmcompute cont1m partzm reduce sum c p1v a r i a b l e cont2m equal c cont1mv a r i a b l e pnm equal count ( partzm )v a r i a b l e zm equal c cont1m / ${pnm}######## zm3 c a l c u l a t i o n ###########v a r i a b l e eatzm3 atom ” c p1 > 2 .0 ”group partzm3 v a r i a b l e eatzm3compute cont1m3 partzm3 reduce sum c p1v a r i a b l e cont2m3 equal c cont1m3v a r i a b l e pnm3 equal count ( partzm3 )v a r i a b l e zm3 equal c cont1m / ${pnm3}######## zm4 c a l c u l a t i o n ###########v a r i a b l e eatzm4 atom ” c p1 > 3 .0 ”group partzm4 v a r i a b l e eatzm4compute cont1m4 partzm4 reduce sum c p1v a r i a b l e cont2m4 equal c cont1m4v a r i a b l e pnm4 equal count ( partzm4 )v a r i a b l e zm4 equal c cont1m / ${pnm4}######################################run 1##### Ca l cu l a t i on o f the s t ress tensor o f the assembly #########Samplev a r i a b l e a14min equal ${a14}#+${ ravg}*5v a r i a b l e a14max equal ${a17}#−${ ravg}*5v a r i a b l e a15min equal ${a15}#+${ ravg}*5v a r i a b l e a15max equal ${a18}#−${ ravg}*5v a r i a b l e a16min equal ${a16}#+${ ravg}*5v a r i a b l e a16max equal ${a19}#−${ ravg}*5reg ion sample block ${a14min} ${a14max} ${a15min} ${a15max} ${a16min} ${a16max} u n i t s boxgroup sample1 reg ion sample#Upper and lower p a r t i c l e p la tes#Lower p l a tev a r i a b l e a16a equal c p7c+${ ravg}reg ion lowerp la te b lock ${a14} ${a17} ${a15} ${a18} ${a16} ${a16a} u n i t s box129group atom reg ion lowerp la te# Shearing p la te a lso Upper p l a t ev a r i a b l e a19a equal c p7f −${ ravg}#−${ ravg}reg ion upperp la te b lock ${a14} ${a17} ${a15} ${a18} ${a19a} ${a19} u n i t s boxgroup atom1 reg ion upperp la tegroup atom3 sub t rac t atom2 atom atom1v a r i a b l e c4 equal count ( atom3 )compute peratom atom2 s t ress / atom v i r i a lcompute p atom2 reduce sum c peratom [ 1 ] c peratom [ 2 ] c peratom [ 3 ] c peratom [ 4 ] c peratom [ 5 ] c peratom [ 6 ]compute pvol a l l pressure thermo temp v i r i a lv a r i a b l e sig11 equal −c p [ 1 ] * 1 / v sysvo lv a r i a b l e sig22 equal −c p [ 2 ] * 1 / v sysvo lv a r i a b l e sig33 equal −c p [ 3 ] * 1 / v sysvo lv a r i a b l e sig12 equal −c p [ 4 ] * 1 / v sysvo lv a r i a b l e sig13 equal −c p [ 5 ] * 1 / v sysvo lv a r i a b l e sig23 equal −c p [ 6 ] * 1 / v sysvo lv a r i a b l e pbvol equal −c pvo l *1* vo l / v sysvo l /1000v a r i a b l e sig11a equal −c p [ 1 ] * 1 / v sysvo l /1000v a r i a b l e sig22a equal −c p [ 2 ] * 1 / v sysvo l /1000v a r i a b l e sig33a equal −c p [ 3 ] * 1 / v sysvo l /1000v a r i a b l e sig12a equal c p [ 4 ] * 1 / v sysvo l /1000v a r i a b l e sig13a equal c p [ 5 ] * 1 / v sysvo l /1000v a r i a b l e sig23a equal c p [ 6 ] * 1 / v sysvo l /1000##### Ca l cu l a t i on o f the s t ress tensor using the lower p l a t ecompute p7lp atom1 proper ty / atom update on run end yes x y z i d vx vy vz f x f y f zcompute p8lp atom1 reduce update on run end yes sum c p7 lp [ 1 0 ]compute p9lp atom1 reduce update on run end yes sum c p7 lp [ 9 ]compute p10lp atom1 reduce update on run end yes sum c p7 lp [ 8 ]v a r i a b l e a20 equal c p8 lp / v area /1000v a r i a b l e a21 equal c p9 lp / ${area}/1000v a r i a b l e a22 equal c p10lp / ${area}/1000130run 1# Inner p a r t i c l e sv a r i a b l e a14a equal c p7a+${ ravg}v a r i a b l e a15a equal c p7b+${ ravg}v a r i a b l e a17a equal c p7d −${ ravg}v a r i a b l e a18a equal c p7e −${ ravg}reg ion p a r t i c l e b lock ${a14a} ${a17a} ${a15a} ${a18a} ${a16a} ${a19a} u n i t s boxgroup a l la tom reg ion p a r t i c l ev a r i a b l e p l t v e l equal 0f i x wa l l1 a l l wa l l / gran model he r t z t a n g e n t i a l h i s t o r y r o l l i n g f r i c t i o n epsd2 sur facesuperquadr ic mesh n meshes 2 meshes cad cad1 #v a r i a b l e s t ress1 equal ${s t ress }/1000s h e l l mkdir shearing1 modv a r i a b l e dumptime equal 20000v a r i a b l e alpha equal 9e−6v a r i a b l e ve lx equal −${alpha } * ( ${s t ress1}−v s ig11 ) *2v a r i a b l e ve ly equal −${alpha } * ( ${s t ress1}−v s ig22 ) *2v a r i a b l e ve lz equal −${alpha } * ( ${s t ress1}−v s ig33 )v a r i a b l e ve lza equal ${alpha } * ( ${s t ress1}−v s ig33 )the rmo s ty le custom step v s ig11a v s ig22a v s ig33a v s ig23a v pbvo lv ve l x v ve l y v ve l z v ve lza v a8 v vo id v CNthermo modify l o s t warn norm nossssrun 0v a r i a b l e t o t a l s t e p equal 1000v a r i a b l e loop s tep equal 1000v a r i a b l e loop count max equal ${ t o t a l s t e p } / ${ l oop s tep}v a r i a b l e loop count loop ${ loop count max}v a r i a b l e c1 equal count ( atom1 )v a r i a b l e c2 equal count ( atom )v a r i a b l e c3 equal count ( atom2 )s h e l l mkdir shearing2ar e s t a r t 200 shearing2a / p1 shearing2a / p2 #DEACTIVATE FOR RESTARTf i x i n t e g r a l l nve / superquadr ic in tegra t ion scheme 1v a r i a b l e ve l equal 0.010f i x 1234 atom move l i n e a r NULL 0.00000000001 0.0000000001 u n i t s box # f i x i n g the lower p l a t e131f i x 123456 atom1 move l i n e a r NULL ${ve l} 0.00000000000001 u n i t s box# Always use the f i x a f t e r ” run 1”f i x ex t ra a l l p r i n t 200 ” ${ t ime} ${CN} ${vo id} ${a8} ${sysvo l} ${ve lx} ${ve ly} ${ve lz}${velza} ${sig11a} ${sig22a} ${sig33a} ${pbvol} ${sig23a} ${ve l} ${a20} ${a21} ${a22} ${zg}${zm} ${zm3} ${zm4} ${ f rac0} ${ f rac1} ${ f r ac2} ${ f r ac3 }” f i l e output / shear ing . t x trun ${nStepsRun}*4 upto pre no post yes every 200 ’ i f ” ${sig23a}>35” then ” q u i t ” ’132Bibliography[1] William Araujo and Christian Ledezma. Factors that affect liquefaction-induced lateralspreading in large subduction earthquakes. Applied Sciences, 10(18):6503, 2020.[2] Edwin L Harp Randall W Jibson. Inventory of landslides triggered by the 1994 northridge,california earthquake. 1995.[3] JD Bray, JP Stewart, MB Baturay, T Durgunoglu, A Onalp, RB Sancio, JP Stewart, D Ural,A Ansal, JB Bardet, et al. Damage patterns and foundation performance in adapazari.Earthquake Spectra, 16(1):163–189, 2000.[4] M Lincoln. The story behind this car in silt photo. nzraw.co.nz, 2020.[5] Timothy P McCrink, Cynthia L Pridmore, John C Tinsley III, Robert R Sickler, Scott JBrandenberg, and Jonathan P Stewart. Liquefaction and other ground failures in imperialcounty, california, from the april 4, 2010, el mayor-cucapah earthquake. Technical report,US Geological Survey, 2011.[6] Ryunosuke Kido and Yosuke Higo. Microscopic characteristics of partially saturated densesand and their link to macroscopic responses under triaxial compression conditions. ActaGeotechnica, 15(11):3055–3073, 2020.[7] GC Cho, J Dodds, and JC Santamarina. Particle shape effects on packing density. Stiffnessand Strength of Natural and Crushed Sands-Internal Report, Georgia Institute ofTechnology, 33pp, 2004.[8] Takeji Kokusho, Tadashi Hara, and Ryousuke Hiraoka. Undrained shear strength ofgranular soils with different particle gradations. Journal of Geotechnical andGeoenvironmental Engineering, 130(6):621–629, 2004.[9] Yogi P. Vaid and Hadi Chern, Jing C.and Tumi. Confinement pressure, grain angularity andliquefaction. Journal of Geotechnical Engineering, 111(10):1229–1235, 1985.[10] Alaa K. Ashmawy, V. Vinh Hoang, and Beena Sukumaran. Evaluating the influence ofparticle shape on liquefaction behavior using discrete element modeling. In The ThirteenthInternational Offshore and Polar Engineering Conference, Honolulu, Hawaii, USA, pages542–549. International Society of Offshore and Polar Engineers, 2003.133[11] Ruben R. Vargas, Kyohei Ueda, and Kazuaki Uemura. Influence of the relative density andK0 effects in the cyclic response of Ottawa F-65 sand-cyclic Torsional Hollow-Cylindershear tests for LEAP-ASIA-2019. Soil Dynamics and Earthquake Engineering,133:106111, 2020.[12] Torsten Wichtmann and Theodoros Triantafyllidis. An experimental database for thedevelopment, calibration and verification of constitutive models for sand with focus tocyclic loading: part i?tests with monotonic loading and stress cycles. Acta Geotechnica,11(4):739–761, 2016.[13] Misko Cubrinovski and Kenji Ishihara. Maximum and minimum void ratio characteristicsof sands. Soils and foundations, 42(6):65–78, 2002.[14] M Jefferies and K Been. Soil liquefaction: A steady state approach, 2016.[15] Junxing Zheng and Roman D Hryciw. Roundness and sphericity of soil particles inassemblies by computational geometry. Journal of Computing in Civil Engineering,30(6):04016021, 2016.[16] Jian Jia, Junxing Zheng, and Xinjian Sun. The critical state parameters of sands from theirimage-based intrinsic properties. Acta Geotechnica, pages 1–12, 2021.[17] Andrew Drescher, Ioannis Vardoulakis, and Chunhua Han. A biaxial apparatus for testingsoils. Geotechnical Testing Journal, 13(3):226–234, 1990.[18] H Bolton Seed and Kenneth L Lee. Liquefaction of saturated sands during cyclic loading.Journal of the Soil Mechanics and Foundations Division, 92(6):105–134, 1966.[19] William H Peacock and H Bolton Seed. Sand liquefaction under cyclic loading simpleshear conditions. Journal of the Soil Mechanics and Foundations Division, 94(3):689–708,1968.[20] Kenji Ishihara and Hiroaki Takatsu. Effects of overconsolidation and k0, conditions on theliquefaction characteristics of sands. Soils and Foundations, 19(4):59–68, 1979.[21] Kinya Miura, Seiichi Miura, and Shosuke Toki. Deformation behavior of anisotropic densesand under principal stress axes rotation. Soils and Foundations, 26(1):36–52, 1986.[22] Marwan Khashila, Mahmoud N Hussien, Mohamed Chekired, and Mourad Karray. On thedynamic soil behavior under triaxial and simple shear modes. International Journal ofGeomechanics, 21(8):04021134, 2021.[23] C Cappellaro, M Cubrinovski, G Chiaro, ME Stringer, JD Bray, and MF Riemer.Undrained cyclic direct simple shear testing of christchurch sandy soils. In 20th NZGSGeotechnical Symposium, pages 1–8, 2017.[24] Committee on Soil Dynamics of the Geotechnical Engineering Division. Definition ofterms related to liquefaction. Journal of the Geotechnical engineering division,104(9):1197–1200, 1978.134[25] Kenji Ishihara and Fumio Yamazaki. Cyclic simple shear tests on saturated sand inmulti-directional loading. Soils and Foundations, 20(1):45–59, 1980.[26] WD Liam Finn, Dennison J Pickering, and Peter L Bransby. Sand liquefaction in triaxialand simple shear tests. Journal of the soil mechanics and foundations division,97(4):639–659, 1971.[27] Mary Ellen Hynes and Richard S Olsen. Influence of confining stress on liquefactionresistance. Routledge, 2018.[28] R Kuerbis and YP Vaid. Sand sample preparation-the slurry deposition method. Soils andFoundations, 28(4):107–118, 1988.[29] ASTM. Standard test methods for determining the amount of material finer than 75-µm(no. 200) sieve in soils by washing. ASTM West Conshohocken, PA, 2014.[30] YP Vaid, JM Fisher, RH , and D Negussey. Particle gradation and liquefaction. Journal ofGeotechnical Engineering, 116(4):698–703, 1990.[31] YU¨KSEL Yilmaz and MURAT Mollamahmutoglu. Characterization of liquefactionsusceptibility of sands by means of extreme void ratios and/or void ratio range. Journal ofgeotechnical and geoenvironmental engineering, 135(12):1986–1990, 2009.[32] HBK Nguyen, MM Rahman, and AB Fourie. Effect of particle shape on constitutiverelation: Dem study. Journal of Geotechnical and Geoenvironmental Engineering,146(7):04020058, 2020.[33] Yoginder P Vaid, Jing C Chern, and Hadi Tumi. Confining pressure, grain angularity, andliquefaction. Journal of Geotechnical Engineering, 111(10):1229–1235, 1985.[34] Fatima M Barreiros, Paulo J Ferreira, and M Margarida Figueiredo. Calculating shapefactors from particle sizing data. Particle & particle systems characterization,13(6):368–373, 1996.[35] Simon J Blott and Kenneth Pye. Particle shape: a review and new methods ofcharacterization and classification. Sedimentology, 55(1):31–63, 2008.[36] Melissa R Cox and Muniram Budhu. A practical approach to grain shape quantification.Engineering Geology, 96(1-2):1–16, 2008.[37] Hakon Wadell. Volume, shape, and roundness of rock particles. The Journal of Geology,40(5):443–451, 1932.[38] LL Sloss. Sequences in the cratonic interior of north america. Geological Society ofAmerica Bulletin, 74(2):93–114, 1963.[39] P. C. Rouse´, R. J. Fannin, and D. A. Shuttle. Influence of roundness on the void ratio andstrength of uniform sand. Ge´otechnique, 58(3):227–231, 2015.135[40] MC Thomas, RJ Wiltshire, and AT Williams. The use of fourier descriptors in theclassification of particle shape. Sedimentology, 42(4):635–645, 1995.[41] Peter A. Cundall and Otto D. L. Strack. A discrete numerical model for granularassemblies. Ge´otechnique, 29(1):47–65, 1979.[42] MD Jiang, ZX Yang, D Barreto, and YH Xie. The influence of particle-size distribution oncritical state behavior of spherical and non-spherical particle assemblies. Granular Matter,20(4):1–15, 2018.[43] C Nouguier-Lehon, B Cambou, and E Vincens. Influence of particle shape and angularityon the behaviour of granular materials: a numerical analysis. International journal fornumerical and analytical methods in geomechanics, 27(14):1207–1226, 2003.[44] Shiwei Zhao, Nan Zhang, Xiaowen Zhou, and Lei Zhang. Particle shape effects on fabricof granular random packing. Powder technology, 310:175–186, 2017.[45] Behzad Soltanbeigi, Alexander Podlozhnyuk, Christoph Kloss, Stefan Pirker, Jin Y Ooi,and Stefanos-Aldo Papanicolopulos. Influence of various dem shape representationmethods on packing and shearing of granular assemblies. Granular Matter, 23(2):1–16,2021.[46] Jean-Francois Ferellec and Glenn R. McDowell. A method to model realistic particleshape and inertia in DEM. Granular Matter, 12(5):459–467, 2010.[47] Shiwei Zhao, T. Matthew Evans, and Xiaowen Zhou. Effects of curvature-related DEMcontact model on the macro-and micro-mechanical behaviours of granular soils.Ge´otechnique, 68(12):1085–1098, 2018.[48] Hoang Bao Khoi Nguyen, Md Mizanur Rahman, and ANDY B Fourie. How particle shapeaffects the critical state, triggering of instability and dilatancy of granular materials–resultsfrom a dem study. Ge´otechnique, pages 1–16, 2020.[49] R. Dyvik, S. Lacasse, T. Berre, and B. Raadim. Comparison of truly undrained andconstant volume direct simple shear tests. Ge´otechnique, 37(1):3–10, 1987.[50] Tang-Tat Ng and Ricardo Dobry. Numerical simulations of monotonic and cyclic loadingof granular soil. Journal of Geotechnical Engineering, 120(2):388–403, 1994.[51] Robert D Holtz, William D Kovacs, and Thomas C Sheahan. An introduction togeotechnical engineering, volume 733. Prentice-Hall Englewood Cliffs, 1981.[52] AW Skempton. Soil mechanics in relation to geology. Proceedings of the YorkshireGeological Society, 29(1):33–62, 1953.[53] JE Goodwin, Wendy Sage, and GP Tilly. Study of erosion by solid particles. Proceedingsof the Institution of Mechanical Engineers, 184(1):279–292, 1969.136[54] T Lambe. The permeability of fine-grained soils. In Symposium on Permeability of Soils.ASTM International, 1955.[55] Ralph H Kuerbis. The effect of gradation and fines content on the undrained loadingresponse of sand. PhD thesis, University of British Columbia, 1989.[56] Xiao Wei, Jun Yang, Yan-Guo Zhou, and Yunmin Chen. Influence of particle-size disparityon cyclic liquefaction resistance of silty sands. Ge´otechnique Letters, 10(2):155–161,2020.[57] I Alpan. The geotechnical properties of soils. Earth-Science Reviews, 6(1):5–49, 1970.[58] Mehmet Murat Monkul, Ehsan Etminan, and Aykut S¸enol. Influence of coefficient ofuniformity and base sand gradation on static liquefaction of loose sands with silt. SoilDynamics and Earthquake Engineering, 89:185–197, 2016.[59] Joseph Simon Zelasko. An investigation of the influences of particle size, size gradationand particle shape on the shear strength and packing behavior of quartziferous sands.Northwestern University, 1966.[60] Alessandro Simoni and Guy T Houlsby. The direct shear strength and dilatancy ofsand–gravel mixtures. Geotechnical & Geological Engineering, 24(3):523–549, 2006.[61] Alessandro Carrera, MR Coop, and Renato Lancellotta. Influence of grading on themechanical behaviour of stava tailings. Ge´otechnique, 61(11):935–946, 2011.[62] J Yang and XD Luo. The critical state friction angle of granular materials: does it dependon grading? Acta Geotechnica, 13(3):535–547, 2018.[63] Emilien Aze´ma, Sandra Linero, Nicolas Estrada, and Arcesio Lizcano. Shear strength andmicrostructure of polydisperse packings: The effect of size span and shape of particle sizedistribution. Physical Review E, 96(2):022902, 2017.[64] David Cantor, Emilien Aze´ma, Philippe Sornay, and Farhang Radjai. Rheology andstructure of polydisperse three-dimensional packings of spheres. Physical Review E,98(5):052910, 2018.[65] David Muir Wood and Kenichi Maeda. Changing grading of soil: effect on critical states.Acta Geotechnica, 3(1):3, 2008.[66] Fletcher M Wood, Jerry A Yamamuro, and Poul V Lade. Effect of depositional method onthe undrained response of silty sand. Canadian Geotechnical Journal, 45(11):1525–1537,2008.[67] Mitsutoshi Yoshimine and Reiko Koike. Liquefaction of clean sand with stratified structuredue to segregation of particle size. soils and Foundations, 45(4):89–98, 2005.[68] JD Frost and JY Park. A critical assessment of the moist tamping technique. GeotechnicalTesting Journal, 26(1):57–70, 2003.137[69] D ASTM. Standard test method for load controlled cyclic triaxial strength of soil. Annualbook of ASTM, 11, 2013.[70] Claudio Cappellaro, Misko Cubrinovski, Jonathan D Bray, Gabriele Chiaro, Michael FRiemer, and Mark E Stringer. Liquefaction resistance of christchurch sandy soils fromdirect simple shear tests. Soil Dynamics and Earthquake Engineering, 141:106489, 2021.[71] Md Mizanur Rahman and Sik-Cheung Robert Lo. Effect of sand gradation and fines typeon liquefaction behaviour of sand-fines mixture. In Geotechnical earthquake engineeringand soil dynamics IV, pages 1–11. 2008.[72] M Belkhatir, A Arab, T Schanz, H Missoum, and N Della. Laboratory study on theliquefaction resistance of sand-silt mixtures: effect of grading characteristics. GranularMatter, 13(5):599–609, 2011.[73] Geoffrey R Martin, H Bolton Seed, and WD Liam Finn. Effects of system compliance onliquefaction tests. Journal of the geotechnical engineering division, 104(4):463–479, 1978.[74] A. Tsomokos and V. N. Georgiannou. Effect of grain shape and angularity on theundrained response of fine sands. Canadian Geotechnical Journal, 47(5):539–551, 2010.[75] Zhen-Zhen Nong, Sung-Sik Park, and Dong-Eun Lee. Comparison of sand liquefaction incyclic triaxial and simple shear tests. Soils and Foundations, 2021.[76] Gye-Chun Cho, Jake Dodds, and J. Carlos Santamarina. Particle shape effects on packingdensity, stiffness, and strength: natural and crushed sands. Journal of Geotechnical andGeoenvironmental Engineering, 132(5):591–602, 2006.[77] Anna-Lena Persson. Image analysis of shape and size of fine aggregates. EngineeringGeology, 50(1-2):177–186, 1998.[78] Zubair Ahmed and Maxim Lebedev. Elastic properties of sands, part 1: Micro computedtomography image analysis of grain shapes and their relationship with microstructure.Geophysical Prospecting, 67(4-Rock Physics: from microstructure to seismicsignatures):723–744, 2019.[79] Robert M Koerner. Effect of particle characteristics on soil strength. Journal of the SoilMechanics and Foundations Division, 96(4):1221–1234, 1970.[80] Kunio Shinohara, Mikihiro Oida, and Boris Golman. Effect of particle shape on angle ofinternal friction by triaxial compression test. Powder technology, 107(1-2):131–136, 2000.[81] Kenneth L Lee and Iraj Farhoomand. Compressibility and crushing of granular soil inanisotropic triaxial compression. Canadian geotechnical journal, 4(1):68–86, 1967.[82] VN Georgiannou. The undrained response of sands with additions of particles of variousshapes and sizes. Ge´otechnique, 56(9):639–649, 2006.138[83] Soheil Ghadr and Arya Assadi-Langroudi. Effect of grain size and shape on undrainedbehaviour of sands. International Journal of Geosynthetics and Ground Engineering,5(3):1–9, 2019.[84] S Joseph Antony and Matthew R Kuhn. Influence of particle shape on granular contactsignatures and shear strength: new insights from simulations. International Journal ofSolids and Structures, 41(21):5863–5870, 2004.[85] Jonathan F. Hubler, Adda Athanasopoulos-Zekkos, and Dimitrios Zekkos. Monotonic,cyclic, and postcyclic simple shear response of three uniform Gravels in constant volumeconditions. Journal of Geotechnical and Geoenvironmental Engineering,143(9):04017043, 2017.[86] Shengjie Rui, Zhen Guo, Tongling Si, and Yujie Li. Effect of particle shape on theliquefaction resistance of calcareous sands. Soil Dynamics and Earthquake Engineering,137:106302, 2020.[87] Jun Ai, Jian-Fei Chen, J Michael Rotter, and Jin Y Ooi. Assessment of rolling resistancemodels in discrete element simulations. Powder Technology, 206(3):269–282, 2011.[88] John R Williams and Alex P Pentland. Superquadrics and modal dynamics for discreteelements in interactive design. Engineering Computations, 9(2):115–127, 1992.[89] R. Hart, Peter A. Cundall, and J. Lemos. Formulation of a three-dimensional distinctelement model—Part II. Mechanical calculations for motion and interaction of a systemcomposed of many polyhedral blocks. International journal of rock mechanics and miningsciences & Geomechanics abstracts, 25(3):117–125, 1988.[90] Keng-Wit Lim and Jose´ E Andrade. Granular element method for three-dimensionaldiscrete element calculations. International Journal for Numerical and Analytical Methodsin Geomechanics, 38(2):167–188, 2014.[91] Reid Kawamoto, Edward Ando`, Gioacchino Viggiani, and Jose´ E. Andrade. Level setdiscrete element method for three-dimensional computations with triaxial case study.Journal of the Mechanics and Physics of Solids, 91:1–13, 2016.[92] Matthew R. Kuhn, Hannah E. Renken, Austin D. Mixsell, and Steven L. Kramer.Investigation of cyclic liquefaction with discrete element simulations. Journal ofGeotechnical and Geoenvironmental Engineering, 140(12):04014075, 2014.[93] Tokio Morimoto, Masahide Otsubo, and Junichi Koseki. Microscopic investigation intoliquefaction resistance of pre-sheared sand: Effects of particle shape and initial anisotropy.Soils and Foundations, 61(2):335–351, 2021.[94] Alan H Barr. Superquadrics and angle-preserving transformations. IEEE Computergraphics and Applications, 1(1):11–23, 1981.139[95] Behzad Soltanbeigi, Alexander Podlozhnyuk, Stefanos-Aldo Papanicolopulos, ChristophKloss, Stefan Pirker, and Jin Y Ooi. Dem study of mechanical characteristics ofmulti-spherical and superquadric particles at micro and macro scales. Powder Technology,329:288–303, 2018.[96] Christoph Kloss, Christoph Goniva, Alice Hager, Stefan Amberger, and Stefan Pirker.Models, algorithms and validation for opensource DEM and CFD–DEM. Progress inComputational Fluid Dynamics, an International Journal, 12(2-3):140–152, 2012.[97] Alexander Podlozhnyuk, Stefan Pirker, and Christoph Kloss. Efficient implementation ofsuperquadric particles in Discrete Element Method within an open-source framework.Computational Particle Mechanics, 4(1):101–118, 2017.[98] Nikolai V Brilliantov, Frank Spahn, Jan-Martin Hertzsch, and Thorsten Po¨schel. Model forcollisions in granular gases. Physical Review E, 53(5):5382, 1996.[99] Raymond D Mindlin and HERBERT Deresiewicz. Elastic spheres in contact under varyingoblique forces. 1953.[100] Alberto Di Renzo and Francesco Paolo Di Maio. Comparison of contact-force models forthe simulation of collisions in DEM-based granular flow codes. Chemical engineeringscience, 59(3):525–541, 2004.[101] CM Wensrich and A Katterfeld. Rolling friction as a technique for modelling particleshape in dem. Powder Technology, 217:409–417, 2012.[102] Ellen M Rathje, Clint Dawson, Jamie E Padgett, Jean-Paul Pinelli, Dan Stanzione, AshleyAdair, Pedro Arduino, Scott J Brandenberg, Tim Cockerill, Charlie Dey, et al. Designsafe:new cyberinfrastructure for natural hazards engineering. Natural Hazards Review,18(3):06017001, 2017.[103] S. Yimsiri and Kenichi Soga. DEM analysis of soil fabric effects on behaviour of sand.Ge´otechnique, 60(6):483–495, 2015.[104] Jiangtao Wei and Gang Wang. Discrete-element method analysis of initial fabric effects onpre-and post-liquefaction behavior of sands. Ge´otechnique Letters, 7(2):161–166, 2017.[105] Ming Yang, Mahdi Taiebat, Patrick Mutabaruka, and Farhang Radjaı¨. Evolution of granularmaterials under isochoric cyclic simple shearing. Physical Review E, 103(3):032904, 2021.[106] Adnan Sufian, A. R. Russell, and Andrew J. Whittle. Anisotropy of contact networks ingranular media and its influence on mobilised internal friction. Ge´otechnique,67(12):1067–1080, 2017.[107] Ning Guo and Jidong Zhao. The signature of shear-induced anisotropy in granular media.Computers and Geotechnics, 47:1–15, 2014.140[108] Stephen K. Park and Keith W. Miller. Random number generators: good ones are hard tofind. Communications of the ACM, 31(10):1192–1201, 1988.[109] Jiangtao Wei and Gang Wang. Evolution of fabric anisotropy in cyclic liquefaction ofsands. Journal of Micromechanics and Molecular Physics, 1(03n04):1640005, 2016.[110] Gang Wang and Jiangtao Wei. Microstructure evolution of granular soils in cyclic mobilityand post-liquefaction process. Granular Matter, 18(3):1–13, 2016.[111] Lei Zhang and T. Matthew Evans. Boundary effects in discrete element method modelingof undrained cyclic triaxial and simple shear element tests. Granular Matter, 20(4):60,2018.[112] Kenji Ishihara, Fumio Tatsuoka, and Susumu Yasuda. Undrained deformation andliquefaction of sand under cyclic stresses. Soils and foundations, 15(1):29–44, 1975.[113] GDR MiDi gdrmidi@ polytech. univ-mrs. fr http://www. lmgc. univ-montp2. fr/MIDI/. Ondense granular flows. The European Physical Journal E, 14:341–365, 2004.[114] Eduardo L. Martin, Colin Thornton, and Stefano Utili. Micromechanical investigation ofliquefaction of granular media by cyclic 3D DEM tests. Ge´otechnique, 70(10):906–915,2020.[115] Roberto R Olivera Bonilla et al. Numerical simulations of undrained granular media.University of Waterloo Waterloo, ON, Canada, 2004.[116] Yasuhiko Okada and Hirotaka Ochiai. Coupling pore-water pressure with distinct elementmethod and steady state strengths in numerical triaxial compression tests under undrainedconditions. Landslides, 4(4):357–369, 2007.[117] Guang Liu, Guan Rong, Jun Peng, and Chuangbing Zhou. Numerical simulation onundrained triaxial behavior of saturated soil by a fluid coupled-dem model. Engineeringgeology, 193:256–266, 2015.[118] Matthew R Kuhn and Ali Daouadji. Simulation of undrained quasi-saturated soil with porepressure measurements using a discrete element (dem) algorithm. Soils and Foundations,60(5):1097–1111, 2020.[119] Catherine O’Sullivan. Particulate discrete element modelling: a geomechanics perspective.CRC Press, 2011.[120] Colin Thornton. Numerical simulations of deviatoric shear deformation of granular media.Ge´otechnique, 50(1):43–53, 2000.[121] David O Potyondy and PA Cundall. A bonded-particle model for rock. Internationaljournal of rock mechanics and mining sciences, 41(8):1329–1364, 2004.141[122] Xin Huang, Kevin J Hanley, Catherine O’Sullivan, and Chung Yee Kwok. Exploring theinfluence of interparticle friction on critical state behaviour using dem. InternationalJournal for Numerical and Analytical Methods in Geomechanics, 38(12):1276–1297, 2014.[123] Weimin Song, Fei Xu, Hao Wu, and Zihao Xu. Investigating the skeleton behaviors ofopen-graded friction course using discrete element method. Powder Technology,385:528–536, 2021.[124] Masanobu Oda. Fabric tensor for discontinuous geological materials. Soils andFoundations, 22(4):96–108, 1982.[125] HF Taylor, C O’Sullivan, T Shire, and WW Moinet. Influence of the coefficient ofuniformity on the size and frequency of constrictions in sand filters. Ge´otechnique,69(3):274–282, 2019.[126] Kate Calder. Statistical inference. New York: Holt, 1953.[127] Charles Voivret, Farhang Radjai, J-Y Delenne, and Moulay Saı¨d El Youssoufi. Multiscaleforce networks in highly polydisperse granular media. Physical review letters,102(17):178001, 2009.[128] Patrick Mutabaruka, Mahdi Taiebat, Roland J.-M. Pellenq, and Farhang Radjai. Effects ofsize polydispersity on random close-packed configurations of spherical particles. PhysicalReview E, 100(4):042906, 2019.[129] Jerry A Yamamuro and Poul V Lade. Steady-state concepts and static liquefaction of siltysands. Journal of geotechnical and geoenvironmental engineering, 124(9):868–877, 1998.[130] Shiwei Zhao and Jidong Zhao. A poly-superellipsoid-based approach on particlemorphology for dem modeling of granular media. International Journal for Numerical andAnalytical Methods in Geomechanics, 43(13):2147–2169, 2019.[131] Junxing Zheng, Hantao He, and Hossein Alimohammadi. Three-dimensional wadellroundness for particle angularity characterization of granular soils. Acta Geotechnica,16(1):133–149, 2021.[132] Ron Goldman. Curvature formulas for implicit curves and surfaces. Computer AidedGeometric Design, 22(7):632–658, 2005.[133] T. G. Sitharam, J. S. Vinod, and B. V. Ravishankar. Post-liquefaction undrained monotonicbehaviour of sands: experiments and DEM simulations. Ge´otechnique, 59(9):739–749,2015.[134] Daniela Dominica Porcino, Theodoros Triantafyllidis, Torsten Wichtmann, and GiuseppeTomasello. Application of critical state approach to liquefaction resistance of sand–siltmixtures under cyclic simple shear loading. Journal of Geotechnical andGeoenvironmental Engineering, 147(3):04020177, 2021.142[135] COR Abbireddy and CRI Clayton. Varying initial void ratios for dem simulations.Geotechnique, 60(6):497–502, 2010.[136] K Michael Salerno, Dan S Bolintineanu, Gary S Grest, Jeremy B Lechman, Steven JPlimpton, Ishan Srivastava, and Leonardo E Silbert. Effect of shape and friction on thepacking and flow of granular materials. Physical Review E, 98(5):050901, 2018.[137] Francesco Zamponi. Packings close and loose. Nature, 453(7195):606–607, 2008.[138] Giulia Guida, Itai Einav, Benjy Marks, and Francesca Casini. Linking micro grainsizepolydispersity to macro porosity. International Journal of Solids and Structures,187:75–84, 2020.[139] Andres R Barrero, William Oquendo, Mahdi Taiebat, and Arcesio Lizcano. Cyclicshearing response of granular material in the semi-fluidized regime. In GeotechnicalEarthquake Engineering and Soil Dynamics V: Numerical Modeling and Soil StructureInteraction, pages 100–107. American Society of Civil Engineers Reston, VA, 2018.[140] Dharma Wijewickreme, Somasundaram Sriskandakumar, and Peter Byrne. Cyclic loadingresponse of loose air-pluviated Fraser River sand for validation of numerical modelssimulating centrifuge tests. Canadian Geotechnical Journal, 42(2):550–561, 2005.[141] Reza Rezvani, Mohammad Amin Tutunchian, and Habib Shahnazari. Methods forestimating the post-cyclic settlement of hormuz island calcareous sand: an experimentalstudy. The European Physical Journal Plus, 135(3):1–20, 2020.[142] H Bolton Seed, IM Idriss, and Ignacio Arango. Evaluation of liquefaction potential usingfield performance data. Journal of geotechnical engineering, 109(3):458–482, 1983.[143] Peter K Robertson and CE Wride. Evaluating cyclic liquefaction potential using the conepenetration test. Canadian geotechnical journal, 35(3):442–459, 1998.[144] Kenji Ishihara. Liquefaction and flow failure during earthquakes. Ge´otechnique,43(3):351–451, 1993.[145] Colin Thornton. Granular dynamics, contact mechanics and particle system simulations. ADEM study. Particle Technology Series, 24, 2015.[146] Trushant S Majmudar and Robert P Behringer. Contact force measurements andstress-induced anisotropy in granular materials. Nature, 435(7045):1079–1082, 2005.[147] Ming Yang. Micro- and macromechanical modeling of granular materials under constantvolume cyclic shearing. PhD thesis, The University of British Columbia, Vancouver,Canada, 2020.[148] Farhang Radjai, Michel Jean, Jean-Jacques Moreau, and Ste´phane Roux. Forcedistributions in dense two-dimensional granular systems. Physical Review Letters,77(2):274–277, 1996.143[149] Nan Zhang and T Matthew Evans. Three dimensional discrete element method simulationsof interface shear. Soils and foundations, 58(4):941–956, 2018.[150] Robert Nau. Statistical forecasting: notes on regression and time series analysis. Stepwiseand All Possible Regressions. Available online: https://people. duke. edu/˜ rnau/regstep.htm (accessed on 2 May 2019), 2016.[151] Eloı¨se Marteau and Jose´ E Andrade. An experimental study of the effect of particle shapeon force transmission and mobilized strength of granular materials. Journal of AppliedMechanics, 88(11):111009, 2021.[152] Etienne Guyon, Jean-Yves Delenne, and Farhang Radjaı¨. Built on Sand: The Science ofGranular Materials. Mit Press, 2020.[153] Md Mizanur Rahman and TG Sitharam. Cyclic liquefaction screening of sand withnon-plastic fines: Critical state approach. Geoscience Frontiers, 11(2):429–438, 2020.[154] Xiaoqiang Gu, Jiachen Zhang, and Xin Huang. Dem analysis of monotonic and cyclicbehaviors of sand based on critical state soil mechanics framework. Computers andGeotechnics, 128:103787, 2020.[155] Jun Yang and H. Y. Sze. Cyclic behaviour and resistance of saturated sand undernon-symmetrical loading conditions. Ge´otechnique, 61(1):59–73, 2011.[156] Xin Huang, Kevin J. Hanley, Zixin Zhang, and Chung-yee Kwok. Structural degradation ofsands during cyclic liquefaction: Insight from DEM simulations. Computers andGeotechnics, 114:103139, 2019.[157] Kostya Shundyak, Martin van Hecke, and Wim van Saarloos. Force mobilization andgeneralized isostaticity in jammed packings of frictional grains. Physical Review E,75(1):010301, 2007.[158] Gary W Delaney and Paul W Cleary. The packing properties of superellipsoids. EPL(Europhysics Letters), 89(3):34002, 2010.[159] Lufeng Liu, Zhiyuan Yu, Weiwei Jin, Ye Yuan, and Shuixiang Li. Uniform and decoupledshape effects on the maximally dense random packings of hard superellipsoids. Powdertechnology, 338:67–78, 2018.[160] HY Sze and J Yang. Failure modes of sand in undrained cyclic loading: impact of samplepreparation. Journal of geotechnical and geoenvironmental engineering, 140(1):152–169,2014.[161] K Pan and ZX Yang. Effects of initial static shear on cyclic resistance and pore pressuregeneration of saturated sand. Acta Geotechnica, 13(2):473–487, 2018.[162] Dorthe Wildenschild and Adrian P Sheppard. X-ray imaging and analysis techniques forquantifying pore-scale structure and processes in subsurface porous medium systems.Advances in Water resources, 51:217–246, 2013.144[163] A Drescher and G De Josselin De Jong. Photoelastic verification of a mechanical modelfor the flow of a granular material. Journal of the Mechanics and Physics of Solids,20(5):337–340, 1972.[164] Gonzalo Castro. Liquefaction and cyclic mobility of saturated sands. Journal ofGeotechnical Engineering Division, 101(GT6):551–569, 1975.[165] Xin Huang, Chung-yee Kwok, Kevin J Hanley, and Zixin Zhang. Dem analysis of theonset of flow deformation of sands: linking monotonic and cyclic undrained behaviours.Acta Geotechnica, 13(5):1061–1074, 2018.[166] Mingjing Jiang, An Zhang, and Tao Li. Distinct element analysis of the microstructureevolution in granular soils under cyclic loading. Granular Matter, 21(2):1–16, 2019.145
Thesis/Dissertation
eng
Civil Engineering
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Taiebat, Mahdi
Effects of particle size distribution and particle shape on cyclic liquefaction response of granular materials
Text