Numerical investigation of force transmission in granular media using discrete element method

Nội dung text: Numerical investigation of force transmission in granular media using discrete element method

1. Vietnam Journal of Mechanics, VAST, Vol.42, No. 2 (2020), pp. 153 – 171 DOI: NUMERICAL INVESTIGATION OF FORCE TRANSMISSION IN GRANULAR MEDIA USING DISCRETE ELEMENT METHOD Thong Chung Nguyen1, Lu Minh Le1, Hai-Bang Ly2, Tien-Thinh Le3,∗ 1Vietnam National University of Agriculture, Hanoi, Vietnam 2University of Transport Technology, Hanoi, Vietnam 3Duy Tan University, Da Nang, Vietnam ∗E-mail: letienthinh@duytan.edu.vn Received: 19 January 2020 / Published online: 10 May 2020 Abstract. In this paper, a numerical Discrete Element Method (DEM) model was cali- brated to investigate the transmission of force in granular media. To this aim, DEM sim- ulation was performed for reproducing the behavior of a given granular material under uniform compression. The DEM model was validated by comparing the obtained shear stress/normal stress ratio with results published in the available literature. The network of contact forces was then computed, showing the arrangement of the material microstruc- ture under applied loading. The number and distribution of the contacts force were also examined statistically, showing that the macroscopic behavior of the granular medium highly depended on the force chain network. The DEM model could be useful in explor- ing the mechanical response of granular materials under different loadings and boundary conditions. Keywords: granular mechanics, discrete element method, force chain, compression test. 1. INTRODUCTION A granular medium is composed of separate particles that move without depen- dence and interact with other particles via contact points [1]. Typical granular materials could be found in civil engineering, such as geotechnical engineering, mining or energy production, chemical, pharmaceutical, and agricultural industries [2–4]. Research and development of machinery/device for processing granular materials have been consid- erably increased over the past ten years, requiring above all a good knowledge of inter- actions between particulate systems itself and with machine parts [5]. For instance, the coefficient of friction has been introduced, measured to characterize the dissipation of en- ergy when the particles collide [6]. These particulate interactions have been investigated for many years using analytical, semi-analytical, or experimental approaches [3,7,8]. De- spite all the efforts, it is not always possible to carry out a large number of configurations c 2020 Vietnam Academy of Science and Technology
2. 154 Thong Chung Nguyen, Lu Minh Le, Hai-Bang Ly, Tien-Thinh Le taking into account all the possible parameters [6]. Moreover, experimental works might not have the required ability to investigate the local interactions, particularly in terms of transmission of stress, collapse of force chain under deformation and so on [9]. It clearly showed that a more robust manner is thus required for better understanding and charac- terizing the mechanical properties of granular materials [10]. From a numerical simulation point of view, the mechanics of granular media can be modeled by either continuum [11–13] or discrete [14–16] approaches. More precisely, in a discrete approach, the Discrete Element Method (DEM) has been primarily employed to simulate granular materials [10, 17]. As an example, Than et al. [18] have developed a DEM model for investigating the plastic response of wet granular material under com- pression. Also, based on DEM technique, Xie et al. [19] have pointed out the influence of interlayer on the strength and deformation of layered rock specimens in uniaxial tests. In another study, Tran et al. [2] have employed DEM algorithm to simulate the behavior of concrete under triaxial loading. Xu et al. [20] have proposed a comparison between DEM simulation and experiments while investigating the mechanical behavior of sea ice. Lom- men et al. [17] have studied the relationship between particle stiffness and bulk material behavior in a numerical simulation context. Furthermore, the combination of DEM and other numerical techniques has been performed by Dratt and Katterfeld [21]. The authors have combined DEM with Finite Element Method (FEM) for investigating the dynamic deformation of machine parts in contact with particle flow. Besides, Zhou et al. [22] have combined DEM with Computational Fluid Dynamics (CFD) for modeling granular flow in hydraulic conveyor. So far, studies involving DEM technique could strongly improve the investigation of mechanical properties of particulate systems by enabling an access to the local behavior in a granular media. Such numerical simulation technique could also save time and cost compared with complex experiments in the design and development of machinery involving particulate systems. In this study, DEM model was developed for investigating the transmission of stress in granular media under the compression force. To this aim, the following steps were adopted as a methodology. First, a set of DEM parameters for the granular media was collected in the available literature, involving dimensional, gravimetric, mechanical, and interaction properties. Precisely, the DEM parameters were the size distribution, shape, mass density, Young’s modulus, Poisson’s ratio, shear modulus, coefficient of static fric- tion, coefficient of rolling friction and coefficient of restitution. In a second step, a com- pression test was designed and performed using DEM simulations. Simultaneously, lo- cal mechanical information of particles was recorded, including the stress, force chain transmission and so on. The obtained results allowed exploring the ability of DEM tech- nique in a mechanical context. Moreover, the features of DEM method were exposed to monitoring and analyzing the displacements and forces of all particles in the considered granular media.
3. Numerical investigation of force transmission in granular media using discrete element method 155 2. MATERIALS AND METHODS 2.1. Brief introduction to DEM DEM was developed based on the simulation of the motion of separate particles in a granular medium [23]. Such motion is determined by solving Newton’s translational and rotational equations of motion for individual particles. The translational equation of motion is given as below [24] dvi mi = ∑ Fij + mig , (1) dt j where mi is the mass of particle i, vi is the velocity, t is the time, Fij is the force of contact acting on the particle i from the particle j, and g is the gravity. The rotational equation of the motion is expressed as follow [23] dωi Ii = ∑ Tij , (2) dt j where Ii is the moment of inertia, ωi is the angular velocity, and Tij is the torque acting on the particle i from the particle j. In a DEM model, the contact force is commonly modeled by spring, dashpot, and frictional slider [25, 26]. One of the most used contact models is the Hertz–Mindlin model [27], involving various parameters such as Young’s modulus, Poisson’s ratio, shear modulus, coefficient of static friction, coefficient of rolling friction and coefficient of restitution [28]. These coefficients, relating the relationships between particle/particle and particle/wall, were introduced to characterize the loss of energy when the particles interact. Based on this principle, DEM simulation could reflect the interactions occurring inside the granular media [18]. Underlying assumptions of DEM model include isotropy and elasticity of the considered particles. On the other hand, the spherical element is the fundamental element in a DEM model. The description of DEM model is well documented in Lommen et al. [17] and Xie et al. [19]. One of the first applications of DEM was carried out by Cundall and Strack for investigating the mechanics of rock and soil [1]. Recently, the fast growth of computa- tional capacity makes it more and more practical to employ numerical methods for solv- ing engineering problems [16]. To date, many works using DEM technique for investi- gating the mechanical properties of granular materials have been published [2,20,29–31]. 2.2. Description of compression test The compression test used in this study is schematized in Fig.1. Granular material with characteristics introduced in Tab.1 was filled into a box container of 400 × 100 × 300 mm. The initial height of the granular medium was 280 mm, exhibiting more than 47.000 particles. At the top of the container, a compression plate is placed. The latter can move freely along the vertical direction (z-axis). A confinement force is exerted to the compression plate, which compresses the granular medium uniformly under a con- stant loading. Such compression force is a constant normal one applying to the particles,
4. Nguyen Chung Thong, Le Minh Lu, Ly Hai Bang and Le Tien Thinh 3 On the other hand, the spherical element is the fundamental element in a DEM model. The description of DEM model is well documented in Lommen et al. [17] and Xie et al. [19]. One of the first applications of DEM was carried out by Cundall and Strack for investigating the mechanics of rock and soil [1]. Recently, the fast growth of computational capacity makes it more and more practical to employ numerical methods for solving engineering problems [16]. To date, many works using DEM technique for investigating the mechanical properties of granular materials have been published [2,20,29–31]. 2.2. Description of compression test 156The compression Thong test Chung used Nguyen, in this Lu study Minh is Le, schematized Hai-Bang Ly, in Tien-Thinh Fig. 1. Granular Le material with characteristics introduced in Table 1 was filled into a box container of 400x100x300 mm. The initial height of the granular medium was 280 mm, exhibiting more than 47.000 particles. At the top of the whereas thecontainer force, a compression acting on plate the is upperplaced. The part latter in can the movex -directionfreely along the is vertical perpendicular direction (z-axis) to. the nor- mal force,A which confinement was force previously is exerted to mentioned. the compression Such plate, a which design compresses of the the test granular allows medium characteriz- uniformly under a constant loading. Such compression force is a constant normal one applying to the ing the transmissionparticles, whereas of the force force inacting the on granularthe upper part medium in the x-direction locally, is perpendicular under compression to the normal using a numericalforce DEM, which approach. was previously mentioned. Such a design of the test allows characterizing the transmission of force in the granular medium locally, under compression using a numerical DEM approach. Fig. 1. DesignFig. 1. Design of compressionof compression test test in this in study. this study 2.3. DEM input parameters 2.3. DEM inputIn this parameters study, the mechanical behavior of agricultural granular materials was investigated, such as dry soybean grains (Glycine max variety, moisture content lower than 10%) to develop and design the In thisseeding study, machine the. The mechanical microscopic parameters behavior of soybean of agricultural particles are commonly granular represented materials based on was inves- tigated, suchfour categories as dry, as soybean in the following. grains (Glycine max variety, moisture content lower than 10%) to developThe first and category design includes the gravimetric seeding properties machine. such as The the true microscopic density. The second parameters category of soy- includes dimensional properties, especially size (i.e., equivalent diameter) and shape. The third category bean particlesincludes are mechanical commonly propertie representeds, such as shear modulus, based Young’s on four modulus categories,, and Poisson’s as inratio. the The following. last The firstcategory category includes the includes interaction properties gravimetric, such as propertiesfriction (coefficient such of static as thefriction true particle/particle density. The sec- ond categoryand particle/wall, includes dimensional coefficient of rolling properties, friction particle/particle especially size and particle/wall), (i.e., equivalent restitution diameter) (coefficient of restitution particle/particle, and particle/wall). It should be noticed that the calibration of and shape. The third category includes mechanical properties, such as shear modulus, Young’s modulus, and Poisson’s ratio. The last category includes the interaction prop- erties, such as friction (coefficient of static friction particle/particle and particle/wall, coefficient of rolling friction particle/particle and particle/wall), restitution (coefficient of restitution particle/particle, and particle/wall). It should be noticed that the calibra- tion of all microscopic parameters for soybean grains is not an easy task 32]. Thus, in this study, the microscopic parameters (i.e., DEM input parameters) of particles were taken from the available literature of Ghodki et al. [32], as it was reported for the same variety of soybean. Moreover, Ghodki et al. [32] have admitted a single sphere modeling for the shape of particles, which allowed reducing the computational time considerably com- pared to multi-spheres or superquadric approaches [33]. It should be noticed that such single sphere modeling was selected based on the shape characterization of the consid- ered particles [32]. In this study, the LIGGGHTS R code (stand for Open Source Discrete Element Method Particle Simulation) was used for the DEM simulations [34]. A no-cohesion nonlinear Hertz–Mindlin model was used for simulating the contact between particle-particle and
5. Numerical investigation of force transmission in granular media using discrete element method 157 particle-wall, as recommended by various works, such as Raji et al. [25], or Horabik et al. [35]. Tab.1 indicates the details of DEM simulation performed in this study, including the DEM input parameters collected from the available literature [32]. The simulations Table 1. Parameters of DEM simulations in this study Parameter Description and value Unit Sliding friction: Hertz-Mindlin Contact model Rolling friction: constant directional torque Cohesion: none Gravity 9.81 m/s2 Particle shape model Spherical Time step 1e-5 s Particle size 6.24 mm True density of particles 1220 kg/m3 Young’s modulus of particles 50 MPa Poisson’s ratio of particles 0.26 Shear modulus of particles 19.84 MPa Young’s modulus of wall 3000 MPa Poisson’s ratio of wall 0.37 Shear modulus of wall 1095 MPa Coefficient of static friction particle/particle 0.26 Coefficient of static friction particle/wall 0.30 Coefficient of restitution particle/particle 0.17 Coefficient of restitution particle/wall 0.35 Coefficient of rolling friction particle/particle 0.08 Coefficient of rolling friction particle/wall 0.08 Length of container 400 mm Width of container 100 mm Number of particles 47362 Initial fill height 280 mm Final height 200 mm Mesh of wall Triangular type (STL) Number of elements (container and plate) 15604 Element area Average: 7.06e-5 m2 Minimum angle Average: 54.15 ˚ Aspect ratio Average: 1.05 Velocity of compression plate 10−1 m/s
6. 158Nguyen Thong Chung Chung Thong, Nguyen, Le LuMinh Minh Lu, Le, Ly Hai-Bang Hai Bang Ly, Tien-Thinh and Le Tien Le Thinh 5 were performedNumberNguyen of Nguyenelements using Chung Chung a (container LenovoThong, Thong, Le ThinkPadand MinhLe Minhplate) Lu, Lu,Ly L420 HaiLy Hai Bang15604 Intel Bang and Coreand Le TienLe i5-2520M Tien Thinh Thinh 2.50 GHz, 5 8 Gb5 of RAM, whereas the post-treatments were performed by using Matlab R2018a [36] and Element area Average: 7.06e-5 m2 Paraview 5.4.1 [37]. NumberMinimumNumber of elements of angle elements (container (container and andplate) plate) 1560415604 Average: 54.15 ° ElementInElement order area toarea ensure the relevance of theAverage: selectedAverage: 7.06e 7.06e set-5 of-5 DEM inputm2 parameters, m2 in- dicatedAspect in Tab. ratio1, size characterization and siloAverage: discharge 1.05 tests were performed. More MinimumMinimum angle angle Average:Average:- 154.15 54.15 ° ° precisely,Velocity the sizeof compression characterization plate allowed obtaining10 a particle size distributionm/s of par- AspectAspect ratio ratio Average:Average: 1.05 1.05 ticlesVelocity (for Velocity generating of compression of compression particle plate plate diameter in DEM10-110 - simulations),1 whereas them/sm/s silo discharge test allowed checking the efficiency of friction coefficients (i.e., static and rolling frictions In order to ensure the relevance of the selected set of DEM input parameters, indicated in Table particle/particle).1, size characterization Brief and details silo discharge of these test twos were investigations performed. More are precisely, following. the size characterization allowedTheIn orderIn size obtainingorder to characterizationensure to ensure a the particle relevancethe relevance size of distribution the particlesof theselected selected of set was particles ofset conductedDEM of DEM (forinput input generatingparameters using parametersa , particleindicated home-made, indicated diameter in Table in Table imaging in DEM 2 1,platform sizesimulations),1, size characterization characterization (4.42 whereas MP/cm and the andsilo silo silopixeldischarge discharge discharge density test stest were Fujifilms allowed were performed performed checking X-E2S. More. More camera theprecisely, efficiencyprecisely, with the sizethe of a Fujinonsizefriction characterization characterization coefficients XF18-55mm (i.e., allowedF2.8-4staticallowed Robtainingand LM obtainingrolling OIS a friction particle lens), a particle ssizeallowed particle/particle) size distribution distribution obtaining of. Brief particles of an particles details image (for of(for generating resolutionthese generating two particle investigations particleof 16 diameter pixels diameter are in per following. DEM in mm DEM (rec- simulations),simulations), whereas whereas the silothe silodischarge discharge test testallowed allowed checking checking the efficiencythe efficiency of friction of friction coefficients coefficients (i.e. (i.e., , staticommendedstatic and Theandrolling rollingsize for friction characterizationcharacterizing frictions particle/particle)s particle/particle) of particlesparticles. Brief. Brief wasd greateretails d conductedetails of these thanof these twousing 3 two mminvestigations a investigations home in size-made are [38 imagingarefollowing.]). following. Soybean platform grains (4.42 wereMP/cm² randomly pixel density selected Fujifilm for capturingX-E2S camera images with a (about Fujinon 900 XF18 grains-55mm were F2.8 tested).-4 R LM Fig.OIS lens)2(a), showsThe theThe size rawsize characterization characterization image, whereas of p ofarticles particles Fig. was 2(b) was conducted presentsconducted using using the a home processed a home-made-made imaging binary imaging platform image platform (4.42 indicating (4.42 MP/cm²allowedMP/cm² pixel obtaining pixel density density anFujifilm imageFujifilm X resolution-E2S X-E2S camera camera of 16with pixelswith a Fujinon a per Fujinon mm XF18 (rXF18ecommended-55mm-55mm F2.8 F2.8- 4for R-4 characterizingLM R LM OIS OIS lens) lens), particles, allowedthegreaterallowed equivalent obtaining than obtaining 3 mman diameter image an in image size resolution [39] ofresolution each). Soybean of particle.16 of pixels16 grainspixels per The permmwere equivalentmm (r randomlyecommended (recommended diameterselected for forcharacterizing for characterizing was capturing computed particles image particles s based (about greateron9greater00 the grainsthan obtained than 3 weremm 3 mm in areatested size in size of[39]). Fig the[39]). .Soybean particle.2).a Soybeanshows grains the Usinggrains rawwere were 900image randomly equivalentrandomly, whereas selected selected Fig diameters, for. 2b forcapturing presents capturing the image the particleimage processeds (abouts (about size binary dis- 9tribution00image 9grains00 grains indicating were is shownwere tested testedthe). in equivalentFig) Fig .Fig 2a . 2(c) shows2a diametershows, exhibiting the theraw of raw image each animage ,particle averagewhereas, whereas. The Fig of .Figequivalent 6.332b. presents2b mm presents anddiameter the theprocessed a standard processed was computed binary deviationbinary based imageofonimage 0.46 theindicating mm.obtainedindicating Itthe is areaequivalentthe seen equivalent of the that diameter particle. thediameter average of Using each of each particle 900 particle particle equivalent. The. diameterThe equivalent equivalent diameters, diameter obtained diameter the was particle bywas computed image computed size baseddistribution analysis based in is onthis showntheon study theobtained inobtained Fig was .area 2c very,area exhibitingof theof close theparticle. particle.an to average theUsing Using result 900of 6.33900 equivalent obtained equivalentmm and diameters, bya diameters,standard Ghodki the deviation theparticle et particle al. [size of32 size0.46] distribution for distributionmm. the It same isis seen is soy- that shownshown in Fig in .Fig 2c,. exhibiting2c, exhibiting an average an average of 6.33 of 6.33 mm mm and anda standard a standard deviation deviation of 0.46 of 0.46 mm. mm. It is It seen is seen that that beanthe averagevariety particle (i.e., 6.24 diameter mm). obtainedFinally, the by image particle analysis size distribution in this study was very used close for generating to the result theobtained the average average by particle Ghodki particle diameter diameteret al. obtained[33] obtained for the by image bysame image soybean analysis analysis invariety this in this study(i.e. study, was6.24 was verymm). very close Finally, close to the to the the result particle result size obtainedparticledistributionobtained by diameter Ghodkiby was Ghodki usedet inal. et DEMfor [33]al. generating[33] for simulations. thefor thesame particle same soybean soybean diameter variety variety in (i.e. DEM (i.e., 6.24 ,simulations. 6.24 mm). mm). Finally, Finally, the theparticle particle size size distributiondistribution was was used used for generatingfor generating particle particle diameter diameter in DEM in DEM simulations. simulations. (a) (b) (c) Fig.Fig. Fig.2. Size 2 2. . SizeSize characterization characterization of particles of of particles particles in this in in thisstudy: this study: study:(a) raw (a) (a) rawimage, raw image, image,(b) processed(b) processed(b) processed image image with image with an equivalent anwith equivalent an equivalent diameterdiameter of each ofof each eachparticle particle particle, and, and(c), and particle (c) (c) particle particle size sizedistribution size distribution distribution from from image from image analysis image analysis. analysis. . Fig. 2 . Size characterization of particles in this study: (a) raw image, (b) processed image with an equivalentRegardingRegardingRegarding diameterthe dischargethe dischargedischarge of eachtest, test, test,a particle, flat a a-flatbottomed flat-bottomed- andbottomed (c)rectangular particlerectangular rectangular silo size siloof distribution silo160 of 160 andof 160and 100 and100 frommm mm100 of image length ofmm length of analysisand length and and widthwidthwidth, respectively, ,respectively respectively, together,, togethertogether with with witha circular a acircular circular orifice orifice orifice of 50 of mm of50 50mm of mm diameter of diameterof diameter, was, was prepared, wasprepared prepared. A kg. A of kg. soybe Aof kgsoybe anof soybean an particles was randomly selected and filled into the silo, exhibiting a fill height of 100 mm. In the DEM particlesparticles was was randomly randomly selected selected and filledand filled into the into silo, the exhibiting silo, exhibiting a fill height a fill ofheight 100 mm. of 100 In themm. DEM In the DEM simulation,simulation,Regarding the thesame same the procedure dischargeprocedure was was test,applied. applied. a flat-bottomedAs Asfriction friction plays plays rectangularthe themost most critical critical silo role of role in 160 thein and therheology rheology 100 mm of simulation, the same procedure was applied. As friction plays the most critical role in the rheology behaviorlengthbehavior and of granular of width, granular materials respectively, materials [32] [32], the together, the efficiency efficiency with of the of a circular the selected selected coefficients orifice coefficients of 50 of mmstatic of static of and diameter, and rolling rolling was frictionprepared.behaviorfrictions particle/particles particle/particle of A granular kg of soybean( see materials ( Tablesee Table 1) particles[32] were 1), were thechecked checked wasefficiency based randomly based on of thison the this test. selected test. To thiTo scoefficients thi andaims aim, in filled ,DEM in DEM of intosimulation static simulation the and silo,, rolling, ex- frictions particle/particle (see Table 1) were checked based on this test. To this aim, in DEM simulation, thehibiting thecoefficient coefficient a fill of height static of static friction of 100friction wasmm. was varied In varied the in DEM ain 0. a18 0. simulation,-0.183-40. range34 range with the with a same step a step of procedure 0 of.0 4 0,.0 whereas4, whereas was the applied. the the coefficient of static friction was varied in a 0.18-0.34 range with a step of 0.04, whereas the coefficientAscoefficient friction of rolling playsof rolling friction the friction most was was varied critical varied between role betweenin 0.0 the5 0.0 - 0.15 rheology -4 0.1 with4 with a step behaviora step of 0.03 of 0.03. A of macroscopic. A granular macroscopic property materials property, ,[ 39], coefficientthe final mass of rolling retained friction in the silo was after varied discharged between, was 0.0 chosen5 - 0.14 to with make a stepcomparisons of 0.03. Abetween macroscopic experiment property , thethe final efficiency mass retained of thein the selected silo after discharged coefficients, was of chosen static to make and comparisons rolling frictions between particle/particle experiment andthe andDEM final DEM simulations. mass simulations. retained in the silo after discharged, was chosen to make comparisons between experiment and DEM simulations.
7. Numerical investigation of force transmission in granular media using discrete element method 159 (see Tab.1) were checked based on this test. To this aim, in DEM simulation, the coef- ficient of static friction was varied in a 0.18–0.34 range with a step of 0.04, whereas the coefficient of rolling friction was varied between 0.05–0.14 with a step of 0.03. A macro- scopic property, the final mass retained in the silo after discharged, was chosen to make comparisons between experiment and DEM simulations. 6 3.Nguyen RESULTS Chung Thong, AND Le Minh DISCUSSIONS Lu, Ly Hai Bang and Le Tien Thinh 3.1. Validation6 of numericalNguyen model Chung Thong, Le Minh Lu, Ly Hai Bang and Le Tien Thinh 3. RESULTS AND DISCUSSIONS In this3.1. section, Validation the of numerical numerical3. modelRESULTS DEM model AND DISCUSSIONS is compared with experimental work in the literature to evaluate the effectiveness of the model. Fig. 3(a) presents the initial assembly 3.1. ValidationIn this section, of numerical the numerical model DEM model is compared with experimental work in the literature of particlesto evaluate inIn thethis the section, box effectiveness container,the numerical of the DEM describedmodel. model Fig. is3a compared inpresents Section thewith initial experimental 2.2, assembly whereas work of particles in Fig. the literature 3(b)in the bshows ox the initial forcetocontainer, evaluate chain describedthe network effectiveness in Section of of the the2.2, model. medium,whereas Fig Fig. 3a. 3b aspresents shows well thethe as initialinitial a visualization force assembly chain of network particles of of in the the medium compressionbox , plate andcontainer,as its well triangular as described a visualization in mesh. Section of the 2.2, compression whereas Fig plate. 3b shows and its the triangular initial force mesh chain. network of the medium, as well as a visualization of the compression plate and its triangular mesh. (a) (b) Fig. 3. Visualization of: (a) particle assembly at initial configuration and (b) initial force chain network and Fig. 3. Visualization of: (a) particlecompression assembly at initialplate withconfiguration triangular and mesh (b). initial force chain network and Fig. 3. Visualization of: (a) particlecompression assembly plate atwith initial triangular configuration mesh. and (b) initial force chain As recommendednetwork by and various compression works in the plate literature with [10,40] triangular, the coefficient mesh of static friction particle/particleAs recommended is characterized by various by worksthe ratio in of the the literature shear stress [10,40] to the, the normal coefficient stress , ofwhile static the friction granular particle/particlematerial is subjected is characterized to loading. by In the this ratio case of of the compression, shear stress the to thestresses normal in thestress x-,axis while (tangential the granular to the As recommendedmaterialdirection is ofsubjected compression) byto loading. various and Inz-axis this works (normalcase of compression, into the the direction literature the stressesof compression) [ 10in the, 40 x-],axis were the (tangential calculated coefficient to by the the of static directioncorresponding of compression) wall reaction and forces z-axis. More (normal prec toisely, the sdirectionuch reactions of compression) were calculated were based calculated on the by reaction the friction particle/particlecorrespondingforces in each wall triangular reaction is mesh characterizedforces element. More in prec contactisely, by with such the particles reactions ratio [40] were of. Fig calculated thes. 4a shear and based 4b show stress on the the reaction evolution to the normal stress, whileforcesof normal thein each stress granular triangular and shear mesh material stress element over in iselapsed contact subjected time with. Theparticles tocomparison loading. [40]. Fig betweens. 4a In and the this 4b shear show case stressthe evolution of to compression,normal the stressesofstress normal in ratio the stress andx-axis andthe workshear (tangential ofstress Ghodki over etelapsed to al. the[33] time directionfor. Thethe considered comparison of compression) granular between material the shear is and stress shownz to-axis in normal Fig. (normal4c . to stressIn Ghodki ratio and et al.the [33] work, the of Ghodkiinter-particle et al. [33]friction for thecoefficient considered of granular0.26 was material calibrated is shown by combining in Fig. 4c .the the directionInexperimental Ghodki of compression) et al.angle [33] of, therepose inter weretest-particle and calculatedDEM friction simulation coefficient by (calibration the of 0.26 corresponding was result calibrated was indicated by wallcombining in Section reaction the 3.2 forces. More precisely,experimentalin Ghodki such et angle al. [33] reactions of )repose. As can test werebe and seen DEM calculatedin Fig. simulation 4c, the normal based(calibration stress on result on the the was reaction wall indicated starts forcestoin increaseSection in 3.2when each trian- in Ghodki et al. [33]). As can be seen in Fig. 4c, the normal stress on the wall starts to increase when gular meshthe elementcompression in plate contact contacts withwith the particles particle assembly [40] Figs.A small 4(a) overshoot and is4(b) also observedshow the, due evolutionto thethe compression first interactions plate between contacts particle with thes and particle compr assemblyession plate A small The compression overshoot is platealso observedis vertically, due moved to of normalthein stress orderfirst interactions to andcompress shear between the granular stress particle material overs and compr under elapsedession constant plate. time. velocity. The Thecompression As for comparison discrete plate element is vertically betweens, the moved particles the shear stress to normalinarranged order to incompress stress order to ratio the respond granular and to thematerial the loading. work under Finally, ofconstant Ghodki the velocity.granular et Asmedium al. for [32 discrete reach] for eselement the a convergence considereds, the particles in both granular arrangedthe normal in order and shear to respond stress. to Such the loading. convergence Finally, exhibits the granular the equilibrium medium reach of thees agranular convergence medium in both under material isthe shownnormal and in shear Fig. stress. 4(c) Such. In convergence Ghodki et exhibits al. [32 the], equilibrium the inter-particle of the granular friction medium under coefficient of constant loading. As shown in Fig. 4c, the ratio of shear stress to normal stress at equilibrium state under constant loading. As shown in Fig. 4c, the ratio of shear stress to normal stress at equilibrium state under 0.26 was calibratedconstant loading by is combininghighly correlated the compared experimental with the work angle of Ghodki of repose et al. [33] test for the and considered DEM simula- constant loading is highly correlated compared with the work of Ghodki et al. [33] for the considered granular material, showing a high effectiveness of the proposed numerical DEM model. tion (calibrationgranular material, result showing was indicateda high effectiveness in Section of the proposed 3.2 in numerical Ghodki DEM et model al. [.32 ]). As can be seen
8. 160 Thong Chung Nguyen, Lu Minh Le, Hai-Bang Ly, Tien-Thinh Le in Fig. 4(c), the normal stress on the wall starts to increase when the compression plate contacts with the particle assembly. A small overshootNguyen is Chung also Thong, observed, Le Minh Lu, due Ly Hai to theBang firstand Le Tien Thinh 7 interactions between particles and compression plate. The compression plate is vertically moved in order to compress the granular material under constant velocity. As for discrete elements, the particles arranged in order to respond to the loading. Finally, the granular medium reaches a convergence in both the normal and shear stress. Such convergence exhibits the equilibrium of the granular medium under constant loading. As shown in Fig. 4(c), the ratio of shear stress to normal stress at equilibrium state under constant loading is highly correlated compared with the work of Ghodki et al. [32] for the con- sidered granular material, showing a high effectiveness of the proposed numerical DEM model.NguyenNguyen Chung Chung Thong, Thong, Le MinhLe Minh Lu, Lu,Ly HaiLy Hai Bang Bang and and Le TienLe Tien Thinh Thinh 7 7 Fig. 4. Evaluation of: (a) normal stress, (b) shear stress, and (c) shear stress / normal stress ratio over time. (a) Normal (b) Tangential (c) Ratio In addition, the results of the silo discharge test are presented in Fig. 5. Visualization of discharge flow at different colored layers in a slice view mode is presented in Fig. 5a, showing the retention zone Fig. 4. Evaluation of: (a) normal stress,in the (b)flat shear-bottomed stress,silo. andFig. 5b (c) shows shearthestress/normal difference Δm between mass retained in the silo from DEM simulations and experiment, in function of the friction coefficients particle/particle. It is shown that the stressdifference ratio over Δm timecould vary between 2 and 30 g. The mass retained in the experiment was 176.7 g. It is seen that the couple of (0.26, 0.08) allowed obtaining the smallest value of Δm (2.1 g). Thus, the 8 8 NguyenNguyen Chung Chung Thong, Thong, Le Minh Le Lu,Minh Ly Lu,Hai Ly Bang Hai and Bang Le andTien LeThinh Tien Thinh efficiency of the selected friction coefficients was confirmed, allowed having more confident results. Fig.Fig. 4. E 4valuation. Evaluation of: (a)of: normal(a) normal stress, stress, (b) shear(b) shear stress stress, and, and(c) shear(c) shear stress stress / normal / normal stress stress ratio ratio over over time time. . In addition,In addition, the theresultsresults of the of thesilo silo discharge discharge test testare arepresented presented in Fig in .Fig5. .Visualization5. Visualization of discharge of discharge flowflowat differentat different colored colored layers layersin ain slice a slice view view mode mode is presented is presented in Fig in .Fig5a., 5ashowing, showing the the retention retention zone zone in thein theflat-flatbottomed-bottomedsilosilo. Fig. .Fig5b. shows5b showsthethe difference differenceΔmΔm between between mass mass retained retained in the in thesilo silo from from DEM DEM simulationssimulations and and experiment experiment, in ,function in function of the of thefrictionfriction coefficients coefficients particle/pa particle/particlerticle. It .isIt shown is shown that that the the differencedifference Δm Δmcould could vary vary between between 2 and 2 and 30 g.30 The g. The mass mass retained retained in thein theexperimentexperiment was was 176.7 176.7 g. Itg .is It is seenseen that thatthe the couple couple of (0.26, of (0.26, 0.08) 0.08) allowed allowed obtaining obtaining the thesmallest smallest value value of Δm of Δm (2.1 (2.1 g). g) Thus,. Thus, the the efficiencyefficiency of the of theselected selected friction friction coefficients coefficients was was confirmed confirmed, allowed, allowed hav havingingmoremore confident confident results results. . (a) (a) (b) (b) Fig. 5. Results of silo discharge(a) test: (a) visualization of particle fl ow at different colored layers(b) and retention Fig. 5. Results of silo discharge test: (a) visualization of particle flow at different colored layers and retention Fig. 5.zone, Results (b) evolution of silo of Δm discharge in function test: of the (a)coefficient visualization of static friction of particle particle/particle flow and at coefficient different of colored layers zone, (b) evolution of Δm in functionrolling of friction the coefficient particle/particle of static. friction particle/particle and coefficient of and retention zone, (b) evolutionrolling of ∆frictionm in particle/particle function of the. coefficient of static friction 3.2. Investigation of transmission of force 3.2. Investigationparticle/particle of transmission and coefficient of force of rolling friction particle/particle In this section, the numerical DEM model was used to investigate the transmission of force in the granular mediumIn this section, under compression. the numerical Fig DEM. 6 shows model the was evolution used to investigate of particle the velocity transmission (z-velocity, of force x- in the velocitygranular, and mediumvelocity magnitude under compression., respectively Fig) .at 6 differentshows the positions evolution of the of compression particle velocity plate. (zFig-velocity,. 7 x- presentsvelocity the, correspondingand velocity magnitudeconfigurations, respectively of the granular) at different medium positions, including of the the force compression chain network plate. Fig. 7 (z-direction,presents thex-direction corresponding, and force configurations magnitude, respectively of the granular). It is medium seen that, including the most themoving force particles chain network are (ztho-direction,se in contact x- directionwith the compression, and force magnitude plate. As the, respectively compression). plateIt is seenwas moved that the in mostthe z- direction,moving particles the arevelocity those in in z -contactdirection with was the dominant compression compared plate. to Asother the directions. compression plate was moved in the z-direction, the velocity in z-direction was dominant compared to other directions.
9. Numerical investigation of force transmission in granular media using discrete element method 161 In addition, the results of the silo discharge test are presented in Fig.5. Visualization of discharge flow at different colored layers in a slice view mode is presented in Fig. 5(a), showing the retention zone in the flat-bottomed silo. Fig. 5(b) shows the difference ∆m between mass retained in the silo from DEM simulations and experiment, in function of the friction coefficients particle/particle. It is shown that the difference ∆m could vary between 2 and 30 g. The mass retained in the experiment was 176.7 g. It is seen that the couple of (0.26, 0.08) allowed obtaining the smallest value of ∆m (2.1 g). Thus, the effi- ciency of the selected friction coefficients was confirmed, allowed having more confident results. 3.2. Investigation of transmission of force In this section, the numerical DEM model was used to investigate the transmission of force in the granular medium under compression. Fig.6 shows the evolution of particle velocity (z-velocity, x-velocity, and velocity magnitude, respectively) at different posi- tions of the compression plate. Fig.7 presents the corresponding configurations of the granular medium, including the force chain network (z-direction, x-direction, and force magnitude, respectively). It is seen that the most moving particles are those in contact with the compression plate. As the compression plate was moved in the z-direction, the velocity in z-directionNguyen was Chung dominant Thong, Le Minh compared Lu, Ly Hai to Bang other and directions.Le Tien Thinh 9 Fig. 6. Visualization of the velocity field of particles in the granular medium at different positions of the Fig. 6. Visualizationcompression plate of. The the colorbar velocity was adapted field offor particleseach velocity in field the in granular order to explore medium the most at appropriate different positions vision effect. of the compression plate. The colorbar was adapted for each velocity field in order to explore the Regarding the force chain most network appropriate (Fig. 7), at vision initial configuration effect (without loading from compression plate), the force chains with low amplitude were created at the bottom of the granular medium, showing the influence of the weight of particles at the top level. However, at initial configuration, the force chains generally had no specific orientations, i.e., the contact forces were uniformly distributed in the medium. When the compression plate contacts with the medium at heights of 270, 260, and 230 mm, the force chains were progressively created, also in increasing amplitude. The contact forces in the z-direction were significant compared to those in the x-direction. This is also proved when regarding the velocity field (Fig. 6). This exciting result showed how the compression forces were transmitted through the particulate system. The orientations of force chains are mainly parallel to the vertical axis, which is the direction of the compression loading. The transmission network also provides information on the structural arrangement, related to the change of the microstructure to respond to the loading.
10. 162 Thong Chung Nguyen, Lu Minh Le, Hai-Bang Ly, Tien-Thinh Le Regarding the force chain network (Fig.7), at initial configuration (without loading from compression plate), the force chains with low amplitude were created at the bottom of the granular medium, showing the influence of the weight of particles at the top level. However, at initial configuration, the force chains generally had no specific orientations, i.e., the contact forces were uniformly distributed in the medium. When the compression plate contacts with the medium at heights of 270, 260, and 230 mm, the force chains were progressively created, also in increasing amplitude. The contact forces in the z-direction were significant compared to those in the x-direction. This is also proved when regarding the velocity field (Fig.6). This exciting result showed how the compression forces were transmitted through the particulate system. The orientations of force chains are mainly parallel to the vertical axis, which is the direction of the compression loading. The trans- mission network also provides information on the structural arrangement, related to the change of10 the microstructureNguyen to Chung respond Thong, Le to Minh the Lu, loading. Ly Hai Bang and Le Tien Thinh Fig. 7. Visualization of force chain network in the granular medium at different positions of compression plate. Fig. 7. VisualizationThe colorbar of force was adapted chain f networkor each case in order the to granular explore the mediummost appropriate at different vision effect positions. of com- pression plate. The colorbar was adapted for each case in order to explore the most appropriate vision effect Fig. 8(a) presents the increase of number of contact forces in function of fill height, normalized to the number of contact force at initial configuration (i.e., 100%), whereas Fig. 8(b) shows the evolution of the number of contact forces in function of elapsed time. Fig. 8. Evaluation of the number of contact forces in function of (a) fill height and (b) elapsed time. Fig. 8a presents the increase of number of contact forces in function of fill height, normalized to the number of contact force at initial configuration (i.e., 100%), whereas Fig. 8b shows the evolution of the number of contact forces in function of elapsed time. It is seen that the number of contact forces
11. 10 Nguyen Chung Thong, Le Minh Lu, Ly Hai Bang and Le Tien Thinh 10 Nguyen Chung Thong, Le Minh Lu, Ly Hai Bang and Le Tien Thinh Numerical investigation of force transmission in granular media using discrete element method 163 It is seen that the number of contact forces linearly increased, as expected, because of the uniform movement of the compression plate. At the equilibrium state, the number of contact forces was increased by about 170% and remained a horizontal asymptotic, as seen in Fig. 8(b). The average and standard deviation values of the probability density distribution of the force chain network are also presented in Fig. 9(a), and 9(b), respec- tively, in function of fill height. As the compression is in the z-direction, the mean value of the contact force in the z-direction was the highest. However, contact forces exhibit ap- proximately the similar standard deviation values in all the directions. Finally, Fig. 9(c) presents the statistical distribution of the magnitude of contact forces, including their av- erage and standard deviation. Statistically, the contact forces increase in both amplitude Fig.and 7.standard Visualization deviation. of force chain network in the granular medium at different positions of compression plate. Fig.The 7. Visualization colorbar was of adapted force chain for each network case in in the order granular to explore medium the at most different appropriate positions vision of compression effect. plate. The colorbar was adapted for each case in order to explore the most appropriate vision effect. Fig. 8Fig Evaluation 8. Evaluation of the of number the number of contact of contact forces forces in functionin function of of(a) (a)fill fill height height and and (b) (b) elapsed elapsed time time. . (a) (b) Fig. 8Faigpresents. 8a presents the increasethe increase of number of number of contact of contact forces forces in infunction function of of fill fill height, height, normalized normalized toto theFig. numberthe 8. number Evaluation of contact of contact of force the force number at initial at initial of configuration contact configuration forces (i.e. in(i.e., function100%), 100%), whereas, of whereas (a) fill Fig. heightFig. 8b 8bshows andshows (b) the the elapsed evolution evolution time of of the numberthe number of contact of contact forces forces in function in function of elapsedof elapsed time time. It. isIt isseen seen that that the the number number of of contact contact forcesforces The force vectors presented in Fig.9 were employed to calculate the total force ex- erted on each particle in the system. At the maximum compression point of 200 mm of height, Fig. 10(a) presents the histogram of the number of particles in contact with a given particle, whereas the histogram of total force exerted on all the particles is shown in Fig. 10(b). The total force exerted on a given particle was calculated by the sum of all the contact forces of its surrounding particles. It can be noticed that at the maximum compression point of 200 mm, each particle was exposed to an average of 8 surrounding particles, whereas the average of total force exerted was 33 N, with a standard deviation of about 13 N. In Appendix, the measurement of the critical breakage force that causes the soybean grains to crack is presented. Results showed that the critical compressive force was in the range of 50-70 N (approximately 1.5 mm of particle deformation). Based on the re- sults obtained (Fig. 10(b)) and the measurement of breakage force, it can be seen that the 200 mm compression point was a critical limit for maintaining the bond between parti- cles. If the compression increased further, the total force exerted on particles would also increase, leading to the destruction of particle bonds.
12. Nguyen Chung Thong, Le Minh Lu, Ly Hai Bang and Le Tien Thinh 11 Nguyen Chung Thong, Le Minh Lu, Ly Hai Bang and Le Tien Thinh 11 linearly increased, as expected, because of the uniform movement of the compression plate. At the equilibrium state, the number of contact forces was increased by about 170% and remained a horizontal linearly increased, as expected, because of the uniform movement of the compression plate. At the asymptotic, as seen in Fig. 8b. Theequilibrium average andstate standard, the number deviation of contact values forces of wasthe increased probability by about density 170% and remained a horizontal distribution of the force chain networkasymptotic are also, as presented seen in Fig. in Fi 8bg 9aThe, and average9b, respectively, and standard in deviation function values of of the probability density fill height. As the compression is indistributionthe z-direction, of the theforce mean chain value network of the arecontact also presented force in inthe Fi zg-.direction9a, and 9b , respectively, in function of was the highest. However, contact forcesfill height. exhibit As the approximately compression isthe in similarthe z-direction, standard the deviation mean value values of the incontact force in the z-direction all the directions. Finally, Fig. 9c presentswas the highest.the statistical However, distribution contact forces of the exhibit magnitude approximately of contact the forces, similar standard deviation values in including their average and standardall the deviation.directions. Finall Statistically,y, Fig. 9c the presents contact the statistical forces increase distribution in bothof the magnitude of contact forces, amplitude and standard deviation. including their average and standard deviation. Statistically, the contact forces increase in both amplitude and standard deviation. Nguyen Chung Thong, Le Minh Lu, Ly Hai Bang and Le Tien Thinh 11 linearly increased, as expected, because of the uniform movement of the compression plate. At the equilibrium state, the number of contact forces was increased by about 170% and remained a horizontal asymptotic, as seen in Fig. 8b. The average and standard deviation values of the probability density distribution of the force chain network are also presented in Fig. 9a, and 9b, respectively, in function of fill height. As the compression is in the z-direction, the mean value of the contact force in the z-direction was the highest. However, contact forces exhibit approximately the similar standard deviation values in all the directions. Finall164y, Fig. 9c presents the statistical Thong distribution Chung Nguyen,of the magnitude Lu Minh of Le, contact Hai-Bang forces, Ly, Tien-Thinh Le including their average and standard deviation. Statistically, the contact forces increase in both amplitude and standard deviation. 12 12 NguyenNguyen Chung Chung Thong, Thong, Le Le Minh Minh Lu, Lu, Ly Ly Hai Hai Bang Bang and and Le Le Tien Tien Thinh Thinh Fig.Fig. 9. Evaluation9. Evaluation of ofcontact contact force force at differentat different fill fill heights heights: (a): (a) average average value, value, (b) (b) standard standard deviation deviation value value, and, and (a) (c)(c) contact contact force force magnitude magnitude. . (b) TheThe force force vector vectors presenteds presented in inFig. Fig. 9 were9 were employed employed to tocalculate calculate the the total total force force exerted exerted on on each each particleparticle in inthe the system. system.AtAt the the maximum maximum compression compression point point of of200 200 mm mm of ofheight height, Fig., Fig. 10a 10a presents presents the the histogramhistogram of ofthe thenumbernumber of ofparticles particles in incontact contact with with a giva givenen particle particle, whereas, whereas the thehistogramhistogram of of total total forceforce exerted exerted on onallall the the particle particles iss isshown shown in inFig Fig. 10. 10b. bTh. The totale total force force exerted exerted on on a givena given particle particle was was calculatedcalculated by by the thesumsum of ofall all the thecontactcontact forces forces of ofits its surrounding surrounding particles particles. It. Itcan can be be noticed noticed that that at atthe the maximummaximum compression compression point point of of200 200 mm, mm, each each particle particle wa was exposeds exposed to toan an average average of of8 8surrounding surrounding particles,particles, whereas whereas the theaverageaverage of oftotal total force force exerted exertedwaswas 33 33 N, N, with with a standarda standard deviation deviationofof about about 13 13 N. N. In Inthe theAppendix,Appendix, the the measurement measurement of ofthe the critical critical breakage breakage force force that that causes causes the the soybean soybean grains grains to to crackcrack is ispresent presented.ed. Results Results showed showed that that the the critical critical compressive compressive force force wa was ins in the the range range of of 50 50-70-70 N N (approximately(approximately 1.5 1.5 mm mm of of particle particle deformation). deformation). Based Based on on the the resul results ts obtained obtained(Fig.(Fig. 10b 10b) and) and the the measurementmeasurement of ofbreakage breakage force, force, it canit can be be seen seen that that (c)the the 200 200 mm mm compression compression point point wa was as criticala critical limit limit forfor maintaining maintaining the the bond bond between between particles. particles. If Ifthe thecompressioncompression increased increasedfurther,further, the the total total force force exerted exerted onon particlesFig. particles 9. Evaluation w ouldwould also also of increase contactincrease, forcelead, leading ating differentto tothe thedestructiondestruction fill heights: of of (a)particle particle average bonds bonds value,. . (b) standard deviation value, and (c) contact force magnitude Fig.Fig. 10 .E10.valuation Evaluation of ofcontact contact(a) force force at Heightat Height = 200= 200 mm mm: (a): (a) number number of ofparticles particles in(b) incontact contact, (b), (b) distribution distribution of of totaltotal force force exerted exerted on on particle particles. s. Fig. 10. Evaluation of contact force at Height = 200 mm: (a) number of particles in contact, 3.33 3Discussion. Discussions s (b) distribution of total force exerted on particles TheThe main main findings findings of ofthis this work work could could be be summarized summarized as asthe the following followings: s: • • A Acalibration calibration procedure procedure of ofa numericala numerical DEM DEM model model for for a granulara granular assembly assembly was was presentedpresented and and validated validated with with experimental experimental data; data; • • TheThe DEM DEM model model was was developed, developed, allowed allowed investigating investigating the the transm transmissionission of of force force in inthe the granular granular medium medium under under compression compressionloading;loading; • • TheThe force force chain chain network network was was statistically statistically quantified, quantified, showing showing the the response response of of the the granulargranular medium medium under under loading loading. . Indeed,Indeed, for for a givena given granular granular medium medium under under loading, loading, the the force force chain chain network network could could be be considered considered as asa loada load-bearing-bearing system system [41] [41]. The. The force force chain chain network network allows allows the the granular granular material material to toadapt adapt itself itself in in orderorder to tosupport support different different loadings loadings and and boundary boundary conditions conditions. From. From a materiala material point point of of view, view, the the force force chainchain network network characterizes characterizes the the microstructure microstructure of ofthe the granular granular material. material. It Ithas has been been pointed pointed out out that that the the granulargranular material material can can change change its its microstructure microstructure in inorder order to tobear bear the the given given loadings loadings [42] [42]. Moreover,. Moreover, the the weakweak network network of ofparticles particles surrounds surrounds the the force force chain chain network network has has been been reported reported as as the theprincipalprincipal energy energy dissipationdissipation source source by by various various works works in inthe the literature literature [43] [43]. Therefore,. Therefore, the the results results of of this this study study confirm confirm the the
13. Numerical investigation of force transmission in granular media using discrete element method 165 3.3. Discussions The main findings of this work could be summarized as the followings: - A calibration procedure of a numerical DEM model for a granular assembly was presented and validated with experimental data; - The DEM model was developed, allowed investigating the transmission of force in the granular medium under compression loading; - The force chain network was statistically quantified, showing the response of the granular medium under loading. Indeed, for a given granular medium under loading, the force chain network could be considered as a load-bearing system [41]. The force chain network allows the granular material to adapt itself in order to support different loadings and boundary conditions. From a material point of view, the force chain network characterizes the microstructure of the granular material. It has been pointed out that the granular material can change its microstructure in order to bear the given loadings [42]. Moreover, the weak network of particles surrounds the force chain network has been reported as the principal energy dissipation source by various works in the literature [43]. Therefore, the results of this study confirm the role of the force chain network in granular mechanics as such network governs the mechanical response of the materials. In order to explore the influence of microscopic particle parameters on the force chain network, especially from friction point of view, DEM simulation for compression test was repeated in changing the value of the coefficient of static friction particle/particle and co- efficient of rolling friction particle/particle, respectively. The results of such a sensitivity analysis are presented in Fig. 11. Figs. 11(a) and 11(b) show the deviation of the num- ber and magnitude of contact forces in function of the deviation of coefficient of static friction particle/particle, respectively (coefficient of static friction particle/particle var- ied in the range of [0.20, 0.26 and 0.32]). On the other hand, Figs. 11(c) and 11(d) show the deviation of the number and magnitude of contact forces in function of the deviation of the coefficient of rolling friction particle/particle, respectively (coefficient of rolling friction particle/particle varied in the range of [0.05, 0.08 and 0.11]). It could be ob- served that changing the value of microscopic parameters of particles modified the force chain network in terms of both number and magnitude of contact forces, with different amplitudes at different heights. A higher friction coefficient allowed obtaining a more significant number and magnitude of contact forces, which was also observed in differ- ent studies in the literature [44, 45]. Besides, it can be seen that the coefficient of static friction particle/particle exhibited a more critical role than the rolling friction coefficient in the compression problem, especially in terms of the magnitude of contact forces. The reason is that the rolling resistance is mainly considered in dynamic impact problems, as pointed out by Zhang et al. [46]. Certainly, further investigations should be carried out in order to explore the fail- ure of the granular materials, i.e., buckling of force chains. Such an investigation could indicate the stability of the force chain network and which parameters that the stabil- ity depends on. On the other hand, the relationship between the force chain network and the energy of the particulate system should also be examined to clarify the energy
14. NguyenNguyen Chung Chung Thong, Thong, Le Le Minh Minh Lu, Lu, Ly Ly Hai Hai Bang Bang and and Le Le Tien Tien Thinh Thinh 1313 NguyenNguyen Chung Chung Thong, Thong, Le Le Minh Minh Lu, Lu, Ly Ly Hai Hai Bang Bang and and Le Le Tien Tien Thinh Thinh 1313 rolerole of ofthe the force force chain chain network network in ingranular granular mechanics mechanicsasas such such network network governs governs the the mechanical mechanical response response ofroleroleof the ofthe ofmaterials the materialsthe force force. chain. chain network network in in granular granular mechanics mechanicsasas such such network network governs governs the the mechanical mechanical response response of the materials. of the materialsInIn order order to. toexplore explore the the influence influence of ofmicroscopic microscopic particle particle parameters parameters on on the the force force chain chain network network, , especiallyespeciallyInIn order fromorder from to friction to explorefriction explore point thepoint the influenceof influenceof view, view,DEM ofDEM ofmicroscopic microscopic simulation simulation particle for particle for compression compression parameters parameters test ontest onwas the was the forcerepeated repeatedforce chain chain in inchanging network changing network, , theespeciallyespeciallythe value value from of from ofthe frictionthe frictioncoefficientcoefficient point point of of of view, of view, static staticDEM DEM friction friction simulation simulation particle/particle particle/particle for for compression compression and and coefficient test coefficient test was was repeated repeated of of rolling rolling in inchanging frictionchanging friction particle/particle,thetheparticle/particle, value value of of the therespectively. respectively.coefficientcoefficient The of The of static results staticresults friction of friction of such such particle/particle a particle/particlesensitivitya sensitivity analysis analysisand and coefficientare coefficientare presented presented of of in rolling inFig rolling Fig. 1.1 friction1 . 1 frictionFig. Figs. s. 1particle/particle,1particle/particle,1a1 aand and 1 1 1b1 b show show respectively. respectively. the the deviation deviation The The ofresults ofresultsthe thenumber ofnumber of such such anda andsensitivitya sensitivity magnitude magnitude analysis analysis of of contact contactare are presented forces presented forces in in function inFig function Fig. 1.1 1.of 1Fig . of theFig s.the s. deviation111deviationa1 aand and 1 of1 1 b1 of b coefficient show coefficientshow the the deviation of deviation of static static of friction of frictionthe thenumber particle/particlenumber particle/particle and and magnitude magnitude, respectively, respectively of of contact contact (coefficient ( coefficient forces forces in of in function of staticfunction static friction of friction ofthe the particle/particledeviationdeviationparticle/particle of of coefficient coefficient varied varied in inofthe of the static rang static rang e frictionof e frictionof [0.20, [0.20, particle/particle 0.26particle/particle 0.26 and and 0.32 0.32]), .] respectively)On, . respectivelyOn the the other other (hand,coefficient (hand,coefficient Fig Figs. s.1 of1 1c of1 staticandc and static 1 1 1 frictiond 1 frictionshowd show theparticle/particleparticle/particlethe deviation deviation of varied of thevaried thenumber innumber in the the rang andrang ande ofe magnitude of magnitude[0.20, [0.20, 0.26 0.26 of and of contactand contact0.32 0.32] ) forces.] On) forces. On the the in other in functionother function hand, hand, ofFig ofFigthes. thes.11deviation1cdeviation1 andc and 11 1 d1 ofshow d of showthe the coefficientthethecoefficient deviation deviation of of of of rollingthe rollingthenumbernumber friction friction and and magnitude particle/particle, magnitude particle/particle, of of contact contact respectively respectively forces forces in in (coefficient function (coefficient function of ofthe of the ofdeviation rollingdeviationrolling friction of friction ofthe the particle/particlecoefficientcoefficientparticle/particle of of varied rolling varied rolling in in frictionthe frictionthe range range particle/particle,of particle/particle, of [0.05, [0.05, 0.08 0.08 and and respectively 0.11]). respectively 0.11]). It Itcould could (coefficient (coefficient be be observed observed ofof that rollingthatrolling changing changing friction friction the the valueparticle/particleparticle/particlevalue of of microscopic microscopic varied varied parameters inparameters in the the range range of of particlesof particlesof [0.05, [0.05, modified 0.08 modified 0.08 and and the 0.11]). the 0.11]).force force It chain It couldchain could network networkbe be observed observed in interms terms that that of ofchangingboth changingboth number number the the andvaluevalueand magnitude of magnitude of microscopic microscopic of of contact parameters contact parameters forces, forces, of of particleswith particleswith different different modified modified amplitudes amplitudes the the force force atchain atchain different different network network heights. heights.in interms terms A A of higher of higherboth both number friction number friction coefficientandandcoefficient magnitude magnitude allow allow ofed of ed contactobtaining contactobtaining forces, forces,a morea more with withsignificant significant different different number amplitudesnumber amplitudes and and magnitude at magnitude at different different of of heights.contact heights.contact forcesA forces A higher higher, which, which friction friction was was alsocoefficientcoefficientalso observed observed allow allow in indifferented eddifferentobtainingobtaining studies studies a amore inmore inthe significantthe literaturesignificant literature number[44,45] number[44,45]. andBesides,. andBesides, magnitude magnitude it icant can beof be ofseencontact seencontact that that forces theforces thecoefficient, coefficientwhich, which was of was of staticalsoalsostatic observed frictionobserved friction particle/particlein inparticle/particle different different studies studiesexhibitedexhibited in in the the literature a literature morea more critical [44,45] critical [44,45] role. Besides,role. Besides,thanthan the i tthe canitrolling canrolling be be seen friction seen friction that that thecoeffic coefficthecoefficientcoefficientientient in inthe of the of compressionstaticstaticcompression friction friction problem,particle/particle problem,particle/particle especially especiallyexhibitedexhibited in interms terms a morea of more ofthe thecritical magnitudecriticalmagnitude role rolethan of thanof contact thecontact therollingrolling forces. forces. friction friction The The reasoncoeffic reasoncoeffic ientis ientisthat inthat inthethe the 166 Thong Chung Nguyen, Lu Minh Le, Hai-Bang Ly, Tien-Thinh Le rcompressionollingcompressionrolling resistance resistance problem, problem, is ism ainm especiallyainly especiallylyconconsideredsidered in in terms interms indynamic dynamic of ofthe themagnitude impactmagnitude impact problems problems of of contact contact, as, as pointed forces.pointed forces. out The out Theby byreason Zhang reason Zhang iset isetthatal. al.that [46] the[46] the. . rollingrolling resistance resistance is ism mainainlylyconconsideredsidered in indynamic dynamic impact impact problems problems, as, as pointed pointed out out by by Zhang Zhang et etal. al. [46] [46]. . (a) (b) Fig.Fig. 11 1. 1Evaluation. Evaluation of ofthethenumbernumber of ofcontact contact force forces ins infunction function of: of (a): (a) coefficient coefficient of ofstatic static friction friction particle/particle,particle/particle,Fig.Fig. 11 1. 1Evaluation .(c) Evaluation (c) coefficient coefficient of of the ofthe number ofrollingnumber rolling frictof offrictcontaction contaction particle/particle; particle/particle;force forces ins infunction function evaluation evaluation of :of (a): (a) coefficientof coefficient ofthe themagnitudemagnitude of ofstatic static of frictionofcontact friction contact forces forces (c) (d) particle/particle,particle/particle, (c) (c) coefficient coefficient of of rolling rolling frict frictionion particle/particle; particle/particle; evaluation evaluation of ofthe themagnitudemagnitude of ofcontact contact forces forces Fig. 11. Evaluation of the number of contact forces in function of: (a) coefficient of static friction particle/particle, (c) coefficient of rolling friction particle/particle; evaluation of the magnitude of contact forces in function of: (b) the coefficient of static friction particle/particle, (d) coefficient of rolling friction particle/particle dissipation. Finally, more complex mechanical tests should be conducted, which could establish the relationship between the microscale parameters and macroscopic responses of the granular media. 4. CONCLUSIONS In this work, numerical simulation of the compression test for a granular medium has been presented to investigate the force transmission. The numerical DEM model was calibrated, taking into account the microscale parameters of the granular medium, especially the particle size distribution, mechanical properties, coefficient of static fric- tion, coefficient of rolling friction and coefficient of restitution. The DEM model was
15. Numerical investigation of force transmission in granular media using discrete element method 167 compared with experimental data, showing a good capability of simulation. The force chain network of the granular medium was calculated, showing the crucial role of such a network on the macroscopic mechanical response of the medium. Contact forces were analyzed statistically, presenting its evaluation in function of the applied load. Results showed that the more the compression applied, the higher number of contacts forces was created, along with an increase of amplitude and standard deviation. In addition, further works should be conducted and applied for different granular materials and also under different loading conditions. Moving walls with constant lateral confinement should be investigated in further studies as they can affect the settlement of the sample in the z- direction and the static/dynamic states of the sample at the time step. Finally, the failure of granular materials should be investigated regarding the collapse of the force chain network at the micro-level. ACKNOWLEDGMENT The authors gratefully acknowledge the Vietnam National University of Agriculture for supporting this research. REFERENCES [1] P. A. Cundall and O. D. L. Strack. A discrete numerical model for granular assemblies. Geotechnique, 29, (1), (1979), pp. 47–65. [2] V. T. Tran, F.-V. Donze,´ and P. Marin. A discrete element model of concrete un- der high triaxial loading. Cement and Concrete Composites, 33, (9), (2011), pp. 936–948. [3] Z. Kaliniewicz and Z. Zuk.˙ A relationship between friction plate roughness and the exter- nal friction angle of wheat kernels. International Journal of Food Properties, 20, (sup3), (2017), pp. S2409–S2417. [4] N. T. Cuong, B. H. Ha, and R. Fukagawa. Failure mechanism of two-dimensional granular columns: Numerical simulation and experiments. Vietnam Journal of Mechanics, 37, (4), (2015), pp. 239–250. [5] H. Fu, C. Jin, and J. Yu. The DEM-based digital design platform for agricultural machin- ery—AgriDEM. In International Conference on Discrete Element Methods, Springer, Springer, (2016), pp. 1253–1263. [6] J. Horabik and M. Molenda. Force and contact area of wheat grain in friction. Journal of Agricultural Engineering Research, 41, (1), (1988), pp. 33–42. 8634(88)90201-6. [7] E. P. Montella,` M. Toraldo, B. Chareyre, and L. Sibille. Localized fluidization in gran- ular materials: Theoretical and numerical study. Physical Review E, 94, (5), (2016). [8] D. T. Huan, T. H. Quoc, and T. M. Tu. Free vibration analysis of functionally graded shell panels with various geometric shapes in thermal environment. Vietnam Journal of Mechanics, 40, (3), (2018), pp. 199–215. [9] F. Nicot, H. Xiong, A. Wautier, J. Lerbet, and F. Darve. Force chain collapse as grain column buckling in granular materials. Granular Matter, 19, (2), (2017).
16. 168 Thong Chung Nguyen, Lu Minh Le, Hai-Bang Ly, Tien-Thinh Le [10] Y. Zhou, H. Wang, B. Zhou, and J. Li. DEM-aided direct shear testing of gran- ular sands incorporating realistic particle shape. Granular Matter, 20, (3), (2018). [11] Z. Deng, P. B. Umbanhowar, J. M. Ottino, and R. M. Lueptow. Continuum modelling of segregating tridisperse granular chute flow. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 474, (2211), (2018). [12] S. Dunatunga and K. Kamrin. Continuum modelling and simulation of granular flows through their many phases. Journal of Fluid Mechanics, 779, (2015), pp. 483–513. [13] S. Bidier and W. Ehlers. Particle simulation of granular media and homogeni- sation towards continuum quantities. PAMM, 13, (1), (2013), pp. 575–576. [14] M. Tolomeo, V. Richefeu, G. Combe, J. N. Roux, and G. Viggiani. An assessment of discrete element approaches to infer intergranular forces from experiments on 2D granular media. International Journal of Solids and Structures, 187, (2020), pp. 48–57. [15] J. Duriez and R. Wan. Stress in wet granular media with interfaces via homogeniza- tion and discrete element approaches. Journal of Engineering Mechanics, 142, (12), (2016). [16] M. Servin, D. Wang, C. Lacoursiere,` and K. Bodin. Examining the smooth and nonsmooth discrete element approaches to granular matter. International Journal for Numerical Methods in Engineering, 97, (12), (2014), pp. 878–902. [17] S. Lommen, D. Schott, and G. Lodewijks. DEM speedup: Stiffness effects on behavior of bulk material. Particuology, 12, (2014), pp. 107–112. [18] V. D. Than, S. Khamseh, A. M. Tang, J.-M. Pereira, F. Chevoir, and J.-N. Roux. Basic mechan- ical properties of wet granular materials: a DEM study. Journal of Engineering Mechanics, 143, (1), (2017). [19] L. Xie, P. Jin, T.-C. Su, X. Li, and Z. Liang. Numerical simulation of uniaxial compression tests on layered rock specimens using the discrete element method. Computational Particle Mechanics, (2019), pp. 1–10. [20] Z. Xu, A. M. Tartakovsky, and W. Pan. Discrete-element model for the in- teraction between ocean waves and sea ice. Physical Review E, 85, (1), (2012). [21] M. Dratt and A. Katterfeld. Coupling of FEM and DEM simulations to con- sider dynamic deformations under particle load. Granular Matter, 19, (3), (2017). [22] M. Zhou, S. Wang, S. Kuang, K. Luo, J. Fan, and A. Yu. CFD-DEM modelling of hydraulic conveying of solid particles in a vertical pipe. Powder Technology, 354, (2019), pp. 893–905. [23] M. Kobayakawa, S. Miyai, T. Tsuji, and T. Tanaka. Interaction between dry gran- ular materials and an inclined plate (comparison between large-scale DEM sim- ulation and three-dimensional wedge model). Journal of Terramechanics, (2019). [24] C. Gonzalez-Montellano,´ J. M. Fuentes, E. Ayuga-Tellez,´ and F. Ayuga. Determi- nation of the mechanical properties of maize grains and olives required for use in DEM simulations. Journal of Food Engineering, 111, (4), (2012), pp. 553–562.
17. Numerical investigation of force transmission in granular media using discrete element method 169 [25] A. O. Raji and J. F. Favier. Model for the deformation in agricultural and food particu- late materials under bulk compressive loading using discrete element method. I: Theory, model development and validation. Journal of Food Engineering, 64, (3), (2004), pp. 359–371. [26] V.-T. Nguyen and M.-Q. Le. Atomistic simulation of the uniaxial compression of black phosphorene nanotubes. Vietnam Journal of Mechanics, 40, (3), (2018), pp. 243–250. [27] K. L. Johnson, K. Kendall, and a. Roberts. Surface energy and the contact of elastic solids. In Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, The Royal Society London, Vol. 324, (1971), pp. 301–313, [28] A. Di Renzo and F. P. Di Maio. Comparison of contact-force models for the simulation of collisions in DEM-based granular flow codes. Chemical Engineering Science, 59, (3), (2004), pp. 525–541. [29] K. F. Malone and B. H. Xu. Determination of contact parameters for discrete ele- ment method simulations of granular systems. Particuology, 6, (6), (2008), pp. 521–528. [30] J. Ai, J.-F. Chen, J. M. Rotter, and J. Y. Ooi. Assessment of rolling resistance mod- els in discrete element simulations. Powder Technology, 206, (3), (2011), pp. 269–282. [31] J. M. Boac, R. P. K. Ambrose, M. E. Casada, R. G. Maghirang, and D. E. Maier. Applications of discrete element method in modeling of grain postharvest operations. Food Engineering Reviews, 6, (4), (2014), pp. 128–149. [32] B. M. Ghodki, M. Patel, R. Namdeo, and G. Carpenter. Calibration of discrete element model parameters: soybeans. Computational Particle Mechanics, 6, (1), (2019), pp. 3–10. [33] T. Xu, J. Yu, Y. Yu, and Y. Wang. A modelling and verification approach for soybean seed particles using the discrete element method. Advanced Powder Technology, 29, (12), (2018), pp. 3274–3290. [34] C. Kloss, C. Goniva, A. Hager, S. Amberger, and S. Pirker. Models, algorithms and validation for opensource DEM and CFD–DEM. Progress in Computational Fluid Dynamics, an Interna- tional Journal, 12, (2-3), (2012), pp. 140–152. [35] J. Horabik and M. Molenda. Parameters and contact models for DEM simulations of agri- cultural granular materials: A review. Biosystems Engineering, 147, (2016), pp. 206–225. [36] Matlab. Natick, MA, USA, (2018). The MathWorks. [37] U. Ayachit. The ParaView guide: A parallel visualization application, ParaView 4.3. Kitware, In- corporated, (2015). [38] D. O. C. Souza and F. C. Menegalli. Image analysis: Statistical study of particle size distribution and shape characterization. Powder Technology, 214, (1), (2011), pp. 57–63. [39] C. J. Coetzee and D. N. J. Els. Calibration of granular material parameters for DEM modelling and numerical verification by blade–granular material interaction. Journal of Terramechanics, 46, (1), (2009), pp. 15–26. [40] T. A. H. Simons, R. Weiler, S. Strege, S. Bensmann, M. Schilling, and A. Kwade. A ring shear tester as calibration experiment for DEM simulations in agitated mixers—a sensitivity study. Procedia Engineering, 102, (1), (2015), pp. 741–748.
18. 170 Thong Chung Nguyen, Lu Minh Le, Hai-Bang Ly, Tien-Thinh Le [41] W. Wu, G. Ma, W. Zhou, D. Wang, and X. Chang. Force transmission and anisotropic charac- teristics of sheared granular materials with rolling resistance. Granular Matter, 21, (4), (2019). [42] A. Salazar, E. Saez,´ and G. Pardo. Modeling the direct shear test of a coarse sand using the 3D Discrete Element Method with a rolling friction model. Computers and Geotechnics, 67, (2015), pp. 83–93. [43] A. Tordesillas, C. A. H. Steer, and D. M. Walker. Force chain and contact cycle evolution in a dense granular material under shallow penetration. Nonlinear Processes in Geophysics, 21, (2), (2014), pp. 505–519. [44] T.-X. Xiu, W. Wang, K. Liu, Z.-Y. Wang, and D.-Z. Wei. Characteristics of force chains in fric- tional interface during abrasive flow machining based on discrete element method. Advances in Manufacturing, 6, (4), (2018), pp. 355–375. [45] W. Wang, W. Gu, and K. Liu. Force chain evolution and force characteristics of shearing granular media in Taylor-Couette geometry by DEM. Tribology Transactions, 58, (2), (2015), pp. 197–206. [46] L. Zhang, N. G. H. Nguyen, S. Lambert, F. Nicot, F. Prunier, and I. Djeran- Maigre. The role of force chains in granular materials: from statics to dynamics. Eu- ropean Journal of Environmental and Civil Engineering, 21, (7-8), (2017), pp. 874–895.
19. Numerical investigation of force transmission in granular media using discrete element method 171 APPENDIX Critical breakage force In this work, the critical breakage force of soybean particles was measured using F.S. 20,000 kN testing machine, available at the Strength of materials laboratory, Faculty of En- gineering, Vietnam National University of Agriculture. In each test, three particles were positioned in the machine to form an equilateral triangle, as showing in Fig. A1(a). Three tests were finally conducted, as shown in Fig. A1(b) for displacement - compression force curves. The breakage of particles was observed at 50-70 N for most of the particles (i.e., a shortening higher than 1.5 mm compared to the average particle diameter of 6.24 mm). It should be noticed that the critical value should also be selected based on the germination 15 rate of particlesNguyen afterChung being Thong, deformed Le Minh (i.e., Lu, higher Ly Hai than Bang 85-90%). and Le Tien Consequently, Thinh 50 N was finally chosenNguyen as a critical Chung breakageThong, Le forceMinh Lu, of soybeanLy Hai Bang particles. and Le Tien Thinh 15 (a) (b) Fig. A1. Measurement of critical breakage force of particles: (a) compression test, (b) displacement - Fig.Fig. A1 A1. .Measurement Measurement of ofcritical critical breakage breakagecompression force force of ofparticles:force particles: curve. (a) (a)compression compression test, test, (b) (b)displacement displacement - compression- compression force force curve. curve REFERENCES [1] P.A. Cundall, O.D.L. Strack, A discrete numericalREFERENCES model for granular assemblies, Géotechnique. 29 (1979) 47–65. [1] P.A. Cundall, O.D.L. Strack, A discrete numerical model for granular assemblies, Géotechnique. 29 (1979) 47[2]– 65.V.T. Tran, F.-V. Donzé, P. Marin, A discrete element model of concrete under high triaxial loading, Cement and Concrete Composites. 33 (2011) 936–948. [2] V.T. Tran, F.-V. Donzé, P. Marin, A discrete element model of concrete under high triaxial loading, Cement [3] Z. Kaliniewicz, Z. Żuk, A relationship between friction plate roughness and the external friction angle of and Concretewheat Composites. kernels, International 33 (2011) 936 Journal–948. of Food Properties. 20 (2017)mp.2011.01.003. S2409–S2417. [3] Z. Kaliniewicz, Z. Żuk, A relationship between friction plate roughness and the external friction angle of wheat[4] N.T. kernels, Cuong, B.H. International Ha, R. Fukagawa, Journal Failure mechanism of Food of two Properties.-dimensional granular 20 (2017) columns: Numerical S2409–S2417. and experiments, Vietnam Journal of Mechanics. 37 (2015) 239–250. [4] N.T. Cuong, B.H. Ha, R. Fukagawa, Failure mechanism of two-dimensional granular columns: Numerical simulation[5] H. Fu, C. and Jin, J. Yu, experiments, The DEM- Based Vietnam Digital Design Journal Platform of for Agricultural Mechanics. Machinery 37 — (2015)AgriDEM, 239 in: –250. Li, Y. Feng, G. Mustoe-7136/37/4/5844. (Eds.), Proceedings of the 7th International Conference on Discrete Element Methods, Springer, Singapore, 2017: pp. 1253–1263. [5] H.[6] Fu, J.C. Horabik, Jin, J. Yu, M. The Molenda, DEM - ForceBased and Digital contact Design area Platform of wheat for grain Agricultural in friction, Machinery Journal of — AgriculturalAgriDEM, in: X. Li, EngineeringY. Feng, G. Research. Mustoe 41 (Eds.), (1988) Proceedings33–42. of the 7th International-8634(88)90201 Conference-6. on Discrete Element Methods, Springer, Singapore, 2017: pp. 1253–1263. [7] E.P. Montellà, M. Toraldo, B. Chareyre, L. Sibille, Localized fluidization in granular materials: Theoretical [6] J. Horabik,and numerical M. Molenda, study, Phys. Force Rev. and E. 94 contact (2016) 052905. area of wheat grain in friction, Journal of Agricultural Engineering[8] D.T. Huan, Research. T.H. Quoc, 41 (1988) T.M. Tu,33– Free42. analysis of functionally-8634(88)90201 graded shell panels-6. with various [7] E.P. Montellà,geometric M. shapes Toraldo, in B. thermal Chareyre, environment, L. Sibil le, Vietnam Localized Journal fluidization of Mechanics. in granular 40 materials: (2018) 199 Theoretical–215. and numerical study, Phys. Rev. E. 94 (2016) 052905. [9] F. Nicot, H. Xiong, A. Wautier, J. Lerbet, F. Darve, Force chain collapse as grain column buckling in [8] D.T. Huan,granula T.H.r materials, Quoc, GranularT.M. Tu, Matter. Free 19vibration (2017) 18. analysis of functionally graded-017 shell-0702 panels-0. with various geometric shapes in thermal environment, Vietnam Journal of Mechanics. 40 (2018) 199–215. [10] Y. Zhou, H. Wang, B. Zhou,-7136/10776. J. Li, DEM- aided direct shear testing of granular sands incorporating realistic particle shape, Granular Matter. 20 (2018) 55. [9] F.[11] Nicot, Z. Deng, H. Xiong, P.B. Umbanhowar, A. Wautier, J.M. J. Lerbet,Ottino, R.M. F. Darve, Lueptow, Force Continuum chain collapsemodelling as of grainsegregating column tridisperse buckling in granulagranularr materials, chute Granularflow, Proceedings Matter. of 19 the (2017) Royal Society18. A: Mathematical, Physical and Engineering-017-0702- 0.Sciences. 474 (2018) 20170384. [10] Y. Zhou, H. Wang, B. Zhou, J. Li, DEM-aided direct shear testing of granular sands incorporating realistic particle shape, Granular Matter. 20 (2018) 55. [11] Z. Deng, P.B. Umbanhowar, J.M. Ottino, R.M. Lueptow, Continuum modelling of segregating tridisperse granular chute flow, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences. 474 (2018) 20170384.