Generation of a realistic 3D sand assembly using X‐ray micro‐computed tomography and spherical harmonic‐based principal component analysis

Three‐dimensional particle morphology is a significant problem in the discrete element modeling of granular sand. The major technical challenge is generating a realistic 3D sand assembly that is composed of a large number of random‐shaped particles containing essential morphological features of natural sands. Based on X‐ray micro‐computed tomography data collected from a series of image processing techniques, we used the spherical harmonics (SH) analysis to represent and reconstruct the multi‐scale features of real 3D particle morphologies. The SH analysis was extended to some highly complex particles with sharp corners and surface cavities. We then proposed a statistical approach for the generation of realistic particle assembly of a given type of sand based on the principle component analysis (PCA). The PCA aims to identify the major pattern of the coefficient matrix, which is made up of the SH coefficients of all the particles involved in the analysis. This approach takes into account the particle size effect on the variation of particle morphology, which is observed from the available results of micro‐computed tomography and QICPIC analyses of sand particle morphology. Using the aforementioned approach, two virtual sand samples were generated, whose statistics of morphological parameters were compared with those measured from real sand particles. The comparison shows that the proposed approach is capable of generating a realistic sand assembly that retains the major morphological features of the mother sand. Copyright © 2016 John Wiley & Sons, Ltd.     

A coupled 3‐dimensional bonded discrete element and lattice Boltzmann method for fluid‐solid coupling in cohesive geomaterials

This paper presents a 3D bonded discrete element and lattice Boltzmann method for resolving the fluid‐solid interaction involving complicated fluid‐particle coupling in geomaterials. In the coupled technique, the solid material is treated as an assembly of bonded and/or granular particles. A bond model accounting for strain softening in normal contact is incorporated into the discrete element method to simulate the mechanical behaviour of geomaterials, whilst the fluid flow is solved by the lattice Boltzmann method based on kinetic theory and statistical mechanics. To provide a bridge between theory and application, a 3D algorithm of immersed moving boundary scheme was proposed for resolving fluid‐particle interaction. To demonstrate the applicability and accuracy of this coupled method, a benchmark called quicksand, in which particles become fluidised under the driving of upward fluid flow, is first carried out. The critical hydraulic gradient obtained from the numerical results matches the theoretical value. Then, numerical investigation of the performance of granular filters generated according to the well‐acknowledged design criteria is given. It is found that the proposed 3D technique is promising, and the instantaneous migration of the protected soils can be readily observed. Numerical results prove that the filters which comply with the design criteria can effectively alleviate or eliminate the appearance of particle erosion in dams.     

Pore‐scale modeling of surface erosion in a particle bed

A fully coupled transient two‐dimensional model was employed to study fundamentals of flood‐induced surface erosion in a particle bed. The interaction of the liquid and solid phases is the key mechanism related to surface erosion. The solid phase was idealized at a particle scale by using the discrete element method. The fluid phase was modeled at a mesoscale level and solved using the lattice Boltzmann method. The fluid forces applied on the particles were calculated on the basis of the momentum the fluid exchanges with the particle. The proposed approach was used to model both single particles and particle beds subjected to Couette flow conditions. The behavior of both the single particle and the particle bed depended on particle diameter and surface shear fluid velocity. The conducted simulations show that the fluid flow profile penetrates the bed for a small distance. This penetration initiates sheet‐flow and surface erosion as the fluid interacts with particles. The effect of suppressing particle rotation on the fluid‐induced forces on the particle was also examined. Suppressing particle spinning may lead to underestimated erosion rate. Results of fluid and particle velocities were compared against experimental results and appeared to agree with the observed trends.Copyright © 2013 John Wiley & Sons, Ltd.     

Sensitivity of incipient particle motion to fluid flow penetration depth within a packed bed

A discrete element method is applied to a three‐dimensional analysis related to sediment entrainment on a micro‐scale. Sediment entrainment is the process by which a fluid medium accelerates particles from rest and advects them upward until they are either transported as bedload or suspended by the flow. Modelling of the entrainment process is a critically important aspect for studies of erosion, pollutant resuspension and transport, and formation of bedforms in environmental flows. Previous discrete element method studies of sediment entrainment have assumed the flow within the particle bed to be negligible and have only allowed for the motion of the topmost particles. At the same time, micro‐scale experimental studies indicate that there is a small slip of the fluid flow at the top of the bed, indicating the presence of non‐vanishing fluid velocity within the topmost bed layers. The current study demonstrates that the onset of particle incipient motion, which immediately precedes particle entrainment, is highly sensitive to this small fluid flow within the topmost bed layers. Using an exponential decay profile for the inner‐bed fluid flow, the discrete element method calculations are repeated with different fluid penetration depths within the bed for several small particle Reynolds numbers. For cases with slip velocity corresponding to that observed in previous experiments with natural sediment, the predicted particle velocity is found to be a few percent of the fluid velocity at the top of the viscous wall layer, which is a reasonable range of velocities for observation of incipient particle motion. This method for prescribing the fluid flow within the particle bed allows for the current discrete element method to be extended in future studies to the analysis of sediment entrainment under the influence of events such as turbulent bursting. Additionally, predictions for the slip velocities and fluid flow profile within the bed suggest the need for further experimental studies to provide the data necessary for additional improvement of the discrete element method models.     

Refined modeling and movement characteristics analyses of irregularly shaped particles

Irregularly shaped (IRS) particles widely exist in many engineering and industrial fields. The macro physical and mechanical properties of the particle system are governed by the interaction between the particles in the system. The interaction between IRS particles is more complicated because of their complex geometric shape with extremely irregular and co‐existed concave and convex surfaces. These particles may interlock each other, making the sliding and friction of IRS particles more complex than that of particles with regular shape. In order to study the interaction of IRS particles more efficiently, a refined method of constructing discrete element model based on computed tomography scanning of IRS particles is proposed. Three parameters were introduced to control the accuracy and the number of packing spheres. Subsequently, the inertia tensor of the IRS particle model was optimized. Finally, laboratory and numerical open bottom cylinder tests were carried out to verify the refined modeling method. The influence of particle shape, particle position, and mesoscopic friction coefficient on the interaction of particles was also simulated. It is noteworthy that with the increase of mesoscopic friction coefficient, the fluidity of IRS particle assembly decreases, and intermittent limit equilibrium state may appear. Copyright © 2014 John Wiley & Sons, Ltd.     

Discrete element method simulations of the seismic response of shallow foundations including soil‐foundation‐structure interaction

A novel three‐dimensional particle‐based technique utilizing the discrete element method is proposed to analyze the seismic response of soil‐foundation‐structure systems. The proposed approach is employed to investigate the response of a single‐degree‐of‐freedom structure on a square spread footing founded on a dry granular deposit. The soil is idealized as a collection of spherical particles using discrete element method. The spread footing is modeled as a rigid block composed of clumped particles, and its motion is described by the resultant forces and moments acting upon it. The structure is modeled as a column made of particles that are either clumped to idealize a rigid structure or bonded to simulate a flexible structure of prescribed stiffness. Analysis is done in a fully coupled scheme in time domain while taking into account the effects of soil nonlinear behavior, the possible separation between foundation base and soil caused by rocking, the possible sliding of the footing, and the dynamic soil‐foundation interaction as well as the dynamic characteristics of the superstructure. High fidelity computational simulations comprising about half a million particles were conducted to examine the ability of the proposed technique to model the response of soil‐foundation‐structure systems. The computational approach is able to capture essential dynamic response patterns. The cyclic moment–rotation relationships at the base center point of the footing showed degradation of rotational stiffness by increasing the level of strain. Permanent deformations under the foundation continued to accumulate with the increase in number of loading cycles. Copyright © 2011 John Wiley & Sons, Ltd.     

3D finite element modeling of sand particle fracture based on in situ X‐Ray synchrotron imaging

Compressive loading of granular materials causes inter‐particle forces to develop and evolve into force chains that propagate through the granular body. At high‐applied compressive stresses, inter‐particle forces will be large enough to cause particle fracture, affecting the constitutive behavior of granular materials. The first step to modeling particle fracture within force chains in granular mass is to understand and model the fracture of a single particle using actual three‐dimensional (3D) particle shape. In this paper, the fracture mode of individual silica sand particles was captured using 3D x‐ray radiography and Synchrotron Micro‐computed Tomography (SMT) during in situ compression experiments. The SMT images were used to reconstruct particle surfaces through image processing techniques. Particle surface was then imported into Abaqus finite element (FE) software where the experimental loading setup was modeled using the extended finite element method (XFEM) where particle fracture was compared to experimental fracture mode viewed in radiograph images that were acquired during experimental loading. Load‐displacement relationships of the FE analysis were also compared with experimental measurements. 3D FE modeling of particle fracture offers an excellent tool to map stress distribution and monitors crack initiation and propagation within individual sand particles. Copyright © 2015 John Wiley & Sons, Ltd.     

Numerical study of inter‐particle bond failure by 3D discrete element method

The cohesive‐frictional nature of cementitious geomaterials raises great interest in the discrete element method (DEM) simulation of their mechanical behavior, where a proper bond failure criterion is usually required. In this paper, the failure of bond material between two spheres was investigated numerically using DEM that can easily reproduce the failure process of brittle material. In the DEM simulations, a bonded‐grain system (composed of two particles and bond material in between) was discretized as a cylindrical assembly of very fine particles connecting two large end spheres. Then, the bonded‐grain system was subjected to compression/tension, shear, rolling and torsion loadings and their combinations until overall failure (peak state) was reached. Bonded‐grain systems with various sizes were employed to investigate bond geometry effects. The numerical results show that the compression strength is highly affected by bond geometry, with the tensile strength being dependent to a lesser degree. The shear, rolling and torsion strengths are all normal force dependent; i.e., with an increase in the normal force, these strengths first increase at a declining rate and then start to decrease upon the normal force exceeding a critical value. The combined actions of shear force, rolling moment and torque lead to a spherical failure envelope in a normalized loading space. The fitted bond geometry factors and bond failure envelopes obtained numerically in this three‐dimensional study are qualitatively consistent with those in previous two‐dimensional experiments. The obtained bond failure criterion can be incorporated into a future bond contact model. Copyright © 2015 John Wiley & Sons, Ltd.     

An efficient flexible membrane boundary condition for DEM simulation of axisymmetric element tests

The mechanical response of an assembly of particles depends on the applied boundary conditions. Robust calibration of numerical discrete systems to laboratory results is also a primary step in many studies of granular materials. In this study, a new membrane model was developed for simulating axisymmetric element tests. This membrane model uses a simple algorithm of an array of independently controlled walls and is computationally efficient. The effect of boundary flexibility on the system response was investigated by simulating a series of triaxial tests on dense and loose specimens. At the specimen scale, differences in shear strength and volume change of specimens were observed. It was shown that localization pattern depends on the applied boundary conditions. At the particle scale, particle‐membrane contact forces, coordination number, local void ratio, and anisotropy of fabric were all affected by the boundary flexibility.     

A numerical study of particle motion and two‐phase interaction in aeolian sand transport using a coupled large eddy simulation – discrete element method

Wind‐blown sand movement, considered as a particle‐laden two‐phase flow, was simulated by a new numerical code developed in the present study. The discrete element method was employed to model the contact force between sand particles. Large eddy simulation was used to solve the turbulent atmospheric boundary layer. Motions of sand particles were traced in the Lagrangian frame. Within the near‐surface region of the atmospheric boundary layer, interparticle collisions will significantly alter the velocity of sand. The sand phase is quite dense in this region, and its feedback force on fluid motion cannot be ignored. By considering the interparticle collision and two‐phase interaction, four‐way coupling was achieved in the numerical code. Profiles of sand velocity from the simulations were in good agreement with experimental measurements. The mass flux shows an exponential decay and is comparable to reported experimental and field measurements. The turbulence intensities and shear stress of sand particles were estimated from particle root‐mean‐square velocities. Distributions of slip velocity and feedback force were analysed to reveal the interactions between sand particles and the continuous fluid phase.     

Numerical studies of shear banding in interface shear tests using a new strain calculation method

Strain localization is closely associated with the stress–strain behaviour of an interphase system subject to quasi‐static direct interface shear, especially after peak stress state is reached. This behaviour is important because it is closely related to deformations experienced by geotechnical composite structures. This paper presents a study using two‐dimensional discrete element method (DEM) simulations on the strain localization of an idealized interphase system composed of densely packed spherical particles in contact with rough manufactured surfaces. The manufactured surface is made up of regular or irregular triangular asperities with varying slopes. A new simple method of strain calculation is used in this study to generate strain field inside a simulated direct interface shear box. This method accounts for particle rotation and captures strain localization features at high resolution. Results show that strain localization begins with the onset of non‐linear stress–strain behaviour. A distinct but discontinuous shear band emerges above the rough surface just before the peak stress state, which becomes more expansive and coherent with post‐peak strain softening. It is found that the shear bands developed by surfaces with smaller roughness are much thinner than those developed by surfaces with greater roughness. The maximum thickness of the intense shear zone is observed to be about 8–10 median particle diameters. The shear band orientations, which are mainly dominated by the rough boundary surface, are parallel with the zero extension direction, which are horizontally oriented. Published in 2007 by John Wiley & Sons, Ltd.     

冻结黏土三轴试验微观变形机理的离散元分析

基于三维颗粒离散单元法,赋予颗粒相应的细观参数,并采用黏结发生在接触颗粒间有限范围内的模型来考虑冻土颗粒中冰的胶结作用,建立了冻结黏土三维离散元数值模型.在相同围压、不同温度和相同温度、不同围压下对冻结黏土的室内三轴试验进行数值模拟,对比了数值试验与室内测试的应力-应变曲线,两者吻合较好.数值模拟结果表明:围压增大会使得接触黏结逐渐失效,在剪切带中胶结冰的破坏区域将增大,而温度的降低则会产生相反结果,这些微观变化都将对冻结黏土的宏观力学变形产生较大影响,同时,细观参数对温度的依赖性也很明显.冻结黏土三轴试验微观变形离散元模拟思路及方法可为今后运用离散单元法研究冻土力学行为提供一定的参考.     

Extended CFD–DEM for free‐surface flow with multi‐size granules

Computational fluid dynamics and discrete element method (CFD–DEM) is extended with the volume of fluid (VOF) method to model free‐surface flows. The fluid is described on coarse CFD grids by solving locally averaged Navier–Stokes equations, and particles are modelled individually in DEM. Fluid–particle interactions are achieved by exchanging information between DEM and CFD. An advection equation is applied to solve the phase fraction of liquid, in the spirit of VOF, to capture the dynamics of free fluid surface. It also allows inter‐phase volume replacements between the fluid and solid particles. Further, as the size ratio (SR) of fluid cell to particle diameter is limited (i.e. no less than 4) in coarse‐grid CFD–DEM, a porous sphere method is adopted to permit a wider range of particle size without sacrificing the resolution of fluid grids. It makes use of more fluid cells to calculate local porosities. The developed solver (cfdemSolverVOF) is validated in different cases. A dam break case validates the CFD‐component and VOF‐component. Particle sedimentation tests validate the CFD–DEM interaction at various Reynolds numbers. Water‐level rising tests validate the volume exchange among phases. The porous sphere model is validated in both static and dynamic situations. Sensitivity analyses show that the SR can be reduced to 1 using the porous sphere approach, with the accuracy of analyses maintained. This allows more details of the fluid phase to be revealed in the analyses and enhances the applicability of the proposed model to geotechnical problems, where a highly dynamic fluid velocity and a wide range of particle sizes are encountered. Copyright © 2015 John Wiley & Sons, Ltd.     

Understanding heavy mineral enrichment using a three‐dimensional numerical model

Layered deposits of relatively light and heavy minerals can be found in many aquatic environments. Quantification of the physical processes which lead to the fine‐scale layering of these deposits is often limited with flumes or in situ field experiments. Therefore, the following research questions were addressed: (i) how can selective grain entrainment be numerically simulated and quantified; (ii) how does a mixed bed turn into a fully layered bed; and (iii) is there any relation between heavy mineral content and bed stability? Herein, a three‐dimensional numerical model was used as an alternative measure to study the fine‐scale process of density segregation during transport. The three‐dimensional model simulates particle transport in water by combining a turbulence‐resolving large eddy simulation with a discrete element model prescribing the motion of individual grains. The granular bed of 0·004 m in height consisted of 200 000 spherical particles (D50 = 500 μ m). Five suites of experiments were designed in which the concentration ratio of heavy (5000 kg m⁻³) to light particles (i.e. 2560 kg m⁻³) was increased from 6%, 15%, 35%, 60% to 80%. All beds were tested for 10 sec at a predefined flow speed of 0·3 m sec⁻¹. Analysis of the particle behaviour in the interior of the beds showed that the lighter particles segregated from the heavy particles with increasing time. The latter accumulated at the bottom of the domain, forming a layer, whereas the lighter particles were transported over the layer forming sweeps. Particles below the heavy particle layer indicated that the layer was able to armour the particles below. Consequentially, enrichment of heavy minerals in a layer is controlled by the segregation of a heavy mineral fraction from the light counterpart, which enhances current understanding of heavy mineral placer formation.     

Discrete simulations of a triaxial compression test for sand by DEM

A quasi‐static homogeneous drained triaxial compression test on cohesionless sand under constant lateral pressure was simulated using a three‐dimensional discrete element method. Grains were modelled by means of particle clusters composed of rigid spheres or spheres with contact moments imitating irregular particle shapes. Attention was paid to the effect of initial void ratio and grain shape mixture on the shear strength, volume changes, force chains, kinetic, elastic and dissipated energies. In addition, the effect of the mean grain size, grain size distribution, grain size range, specimen size and roughness and stiffness of boundaries was numerically analysed in initially dense sand. Some numerical results were compared with available experimental results. Copyright © 2014 John Wiley & Sons, Ltd.     

Analysis of soil‐structural interface behavior using three‐dimensional DEM simulations

The shear behavior at the interface between the soil and a structure is investigated at the macroscale and particle‐scale levels using a 3‐dimensional discrete element method (DEM). The macroscopic mechanical properties and microscopic quantities affected by the normalized interface roughness and the loading parameters are analyzed. The macro‐response shows that the shear strength of the interface increases as the normalized roughness of the interface increases, and stress softening and dilatancy of the soil material are observed in the tests that feature rough interfaces. The particle‐scale analysis illustrates that a localized band characterized by intense shear deformation emerges from the contact plane and gradually expands as shearing progresses before stabilizing at the residual stress state. The thickness of the localized band is affected by the normalized roughness of the interface and the normal stress, which ranges between 4 and 5 times that of the median grain diameter. A thicker localized band is formed when the soil has a rough shearing interface. After the localized band appears, the granular material structuralizes into 2 regions: the interface zone and the upper zone. The mechanical behavior in the interface zone is representative of the interface according to the local average stress analysis. Certain microscopic quantities in the interface zone are analyzed, including the coordination number and the material fabric. Shear at the interface creates an anisotropic material fabric and leads to the rotation of the major principal stress.     

Investigation of slope instability induced by seepage and erosion by a particle method

A novel particle based Bluff Morphology Model (BMM) developed by the authors is extended in this paper to investigate the effect of two dimensional seepage on the stability and collapse of soil slopes and levees. To incorporate the seepage in the model, Darcy’s law is applied to the interactions among neighbouring soil particles and ghost particles are introduced along the enclosed soil boundary so that no fluid crosses the boundary. The contribution of partially saturated soils and matric suction, as well as the change in hydraulic conductivity due to seepage, are predicted well by the present model. The predicted time evolution of slope stability and seepage induced collapse are in reasonable agreement with the experimental results for homogeneous non-cohesive sand and multiple layered cohesive soils. Rapid drawdown over a sand soil is also investigated, and the location and time of the levee collapse occurrence are well captured. A toe erosion model is incorporated in the BMM model, and the location and quantity of erosion from lateral seepage flow is well predicted. The interplay of erosion, seepage and slope instability is examined.     

Discrete element modelling of a granular platform supported by piles in soft soil – Validation on a small scale model test and comparison to a numerical analysis in a continuum

This paper focuses on the mechanisms taking place in a granular platform supported by piles in soft soil. Several modelling approaches were explored. A two-dimensional small scale model test using the Taylor–Schneebeli soil analogue was first developed and the experimental results were compared to a discrete element model using the particle code PFC. The validation of this numerical approach allowed the parametric study to be extended numerically. Parametric studies were also performed on continuum model using the finite-difference code FLAC. Comparison of the parametric studies performed on each modelling approach underlined some differences and lead to a consideration on the macro- and micromechanical parameters.     

Bond rolling resistance and its effect on yielding of bonded granulates by DEM analyses

A discrete element modelling of bonded granulates and investigation on the bond effect on their behaviour are very important to geomechanics. This paper presents a two‐dimensional (2‐D) discrete element theory for bonded granulates with bond rolling resistance and provides a numerical investigation into the effect of bond rolling resistance on the yielding of bonded granulates. The model consists of mechanical contact models and equations governing the motion of bonded particles. The key point of the theory is that the assumption in the original bond contact model previously proposed by the authors (55th CSCE‐ASCE Conference, Hamilton, Ont., Canada, 2002; 313–320; J. Eng. Mech. (ASCE) 2005; 131 (11):1209–1213) that bonded particles are in contact at discrete points, is here replaced by a more reliable assumption that bonded particles are in contact over a width. By making the idealization that the bond contact width is continuously distributed with the normal/tangential basic elements (BE) (each BE is composed of spring, dashpot, bond, slider or divider), we establish a bond rolling contact model together with bond normal/tangential contact models, and also relate the governing equations to local equilibrium. Only one physical parameter β needs to be introduced in the theory in comparison to the original bond discrete element model. The model has been implemented into a 2‐D distinct element method code, NS2D. Using the NS2D, a total of 86 1‐D, constant stress ratio, and biaxial compressions tests have been carried out on the bonded granular samples of different densities, bonding strengths and rolling resistances. The numerical results show that: (i) the new theory predicts a larger internal friction angle, a larger yielding stress, more brittle behaviour and larger final broken contact ratio than the original bond model; (ii) the yielding stress increases nonlinearly with the increasing value of β, and (iii) the first‐yield curve (initiation of bond breakage), which define a zone of none bond breakage and which shape and size are affected by the material density, is amplified by the bond rolling resistance in analogous to that predicted by the original bond model. Copyright © 2006 John Wiley & Sons, Ltd.