Three-dimensional Discrete Element Method Simulations Research Articles

Three-dimensional discrete element simulations on pressure ridge formation

Abstract. This study presents the first three-dimensional discrete element method simulations of pressure ridge formation. Pressure ridges are an important feature of the sea-ice cover, as they contribute to the mechanical thickening of ice and likely limit the strength of sea ice at large scales. We validate the simulations against laboratory-scale experiments, confirming their accuracy in predicting ridging forces and ridge geometries. Then we demonstrate that Cauchy–Froude scaling applies for translating laboratory-scale results on ridging to full-scale scenarios. We show that non-simultaneous failure, where an ice floe fails at distinct locations across the ridge length, is required for an accurate representation of the ridging process. This process cannot be described by two-dimensional simulations. We also find a linear relationship between the ridging forces and the ice thickness, contrasting with earlier results in the literature obtained by two-dimensional simulations.

Fabric-based jamming phase diagram for frictional granular materials.

A jamming phase diagram maps the phase states of granular materials to their intensive properties such as shear stress and density (or packing fraction). We investigate how different phases in a jamming phase diagram of granular materials are related to their fabric structure via three-dimensional discrete element method simulations. Constant-volume quasi-static simple shear tests ensuring uniform shear strain field are conducted on bi-disperse spherical frictional particles. Specimens with different initial solid fractions are sheared until reaching steady state at a large shear strain (200%). The jamming threshold in terms of stress, non-rattler fraction, and coordination numbers (Z's) of different contact networks is discussed. The evolution of fabric anisotropy (F) of each contact network during shearing is also examined. By plotting the fabric data in the F-Z space, a unique critical fabric surface (CFS) becomes apparent across all specimens, irrespective of their initial phase states. Through the correlation of this CFS with fabric signals corresponding to jamming transitions, we introduce a novel jamming phase diagram in the fabric F-Z space, offering a convenient approach to distinguish the various phases of granular materials solely through the direct observation of geometrical arrangements of particles. This jamming phase diagram underscores the importance of the microstructure underlying the conventional jamming phenomenon and introduces a novel standpoint for interpreting the phase transitions of granular materials that have been exposed to processes such as compaction, shearing, and other complex loading histories.

The force and dynamic response of low-velocity projectile impact into 3D dense wet granular media

The presence of liquid bridges between particles in a wet granular media gives rise to capillary forces at the microscopic scale that dramatically alters the material's microstructure and macroscopic responses. Previous experimental studies show significant differences in the penetration depth of the projectile into assemblies of dry or wet particles. Still, there is a lack of fundamental understanding of the difference induced by the presence of liquid bridges in dynamic regimes. In this study, the contact model parameters of wet glass beads are calibrated employing the angle of repose tests and cylinder lifting tests, and a series of three-dimensional discrete element method simulations of spheres impact into wet granular packings are conducted. The dynamic characteristics of the projectile are obtained numerically and found to be in good agreement with the experimental data. The analyses indicate that the final penetration depth in the wet granular media is linearly related to the initial velocity and has a power function relationship with the impact energy. The generalized Poncelet law is suitable for wet granular materials, but the presence of liquid significantly affects the values of the parameters. There is a velocity “stagnation zone” at the initial stage of the wet granular impact, and the final penetration depth and the duration interaction time are smaller than that of the dry case. The difference can be attributed to the resistance evolution of the projectile exerted by particles and gives the parameter scaling law of the peak impact force. Further, a macro-micro transition technique is employed to characterize the velocity field, pressure field, and stress field inside the wet granular media and the spatial distribution is quantified based on stress theory.

Macro-micro mechanical study of principal stress rotation in granular materials using DEM simulations of hollow cylinder test

The dependency of sand behavior anisotropy on microstructure reveals the importance of macro-micro mechanical analysis of the principal stresses rotation effect on the strength-deformation response of granular material. In this study multi-scale analysis of the granular material's behavior using three-dimensional discrete element method simulations of the hollow cylinder test was presented. The effect of principal stresses rotation (α) on the macro-mechanical behavior of the soil was investigated from the micro-mechanical anisotropy aspect. It was observed that principal stress rotation affects the anisotropy coefficient and directional distribution of all micro-parameters which can change the soil's behavior according to loading conditions. The stress-force-fabric relationship was used to analyze the microstructure's stress state based on micro-parameters anisotropy coefficients. In addition to the good agreement between stress-force-fabric and general stress tensors, the decreasing trend of the maximum deviator stress and internal friction coefficient angle with increasing α can be observed for both stress tensors.

Liquefaction behavior of inherently anisotropic sand under cyclic simple shear: Insights from three-dimensional DEM simulations

Inherent fabric anisotropy has been widely recognized to have significant effects on the liquefaction behavior of sand. A discrete element method (DEM) study was conducted to investigate the liquefaction behavior of inherently anisotropic sand under cyclic simple shear conditions. In the DEM model, the under-compaction method was modified to generate uniform and inherently anisotropic specimens. A series of DEM simulations of undrained cyclic simple shear tests were performed to investigate the effects of bedding plane angles with a full spectrum ranging from −90° to +90°. The relative magnitudes of the three effective normal stress components were compared and analyzed. The internal structure evolution was quantified through a contact-normal-based fabric tensor that could describe the load-bearing structure of granular specimens. The simulation results show that the effect of bedding plane angle on the liquefaction resistance of anisotropic sand is not monotonic. For example, as the bedding plane angle of specimen increases from 0° to 90°, the liquefaction resistance first decreases and then increases, with the most detrimental scenario occurs at 45°. This could be explained by correlating the magnitudes of shear modulus with the discrepancy between the loading direction and the direction of fabric during each cycle.

Construction of motional phase maps for granular dampers

Harmonically vibrated granular media exhibit a variety of motional behaviours depending on amplitude, frequency, and vibration-to-gravity directional orientation. Motional behaviour defines the physical interactions of particles in the granular media and therefore the energy dissipation performance. A "phase map" that describes motional behaviour over broad ranges of frequency and amplitude is therefore a very useful tool in damper design. However, at present, identification of the operating motional conditions within the granular media has only been conducted by visual observation of the particles following a particle-level simulation or a specifically designed experiment. Because of this, design optimisation over a broad range of amplitude and frequency becomes costly. This paper aims to help reduce this cost through the development of approximate phase maps based on expected dissipative interactions of particles. Three-dimensional discrete element method simulations are conducted over a wide range of excitation intensities under two different vibration-to-gravity directional orientations (i.e., perpendicular, and parallel to the standard gravity direction) to allow the observation of as many motional phases as possible. The effect of particle size, volume filling ratio and particle shape on granular energy dissipation sources are also investigated.

Use of Combined Static and Dynamic Testing to Quantify the Participation of Particles in Stress Transmission

A number of research studies have recognized that not all particles in a sand specimen are active in stress transmission, particularly in the case of gap-graded soils. This has implications for the use of the global void ratio (e), and variables/parameters that depend upon e, to predict the influence of the state on the mechanical behavior of gap-graded soils. This study explores the possibility of comparing the shear stiffness determined in dynamic wave propagation tests (Gdyn) with stiffness values determined in small-strain probes (Gsta) to assess the extent to which all the particles in a specimen are actively engaged in stress transmission. The idea is initially developed using three-dimensional discrete element method simulations and considering ideal specimens. The practical application of this approach is then tested in a series of drained triaxial compression tests. The numerical studies considered both specimens with continuous, linear particle size distributions as well as gap-graded soils. The results show that the ratio Gsta/Gdyn may be associated with the ratio between the bulk density calculated considering only the stress transmitting particles, ρm, and the bulk density calculated considering all particles, ρ. The ratio ρm/ρ is determined by the proportion of inactive particles, which varies with the proportion by mass of finer particles in the soil (Ffiner) for gap-graded soils. The relationship between Gsta/Gdyn(e) and Ffiner enables a qualitative assessment of the proportion of inactive particles. The corresponding experimental test results show a similar but weaker correlation between the ratio of Gsta/Gdyn and Ffiner, this weaker correlation may be attributed to the differences between the simulations and the experimental conditions. However, despite the challenges with experimental implementation, the data presented here support the idea that it may be possible to qualitatively estimate the proportion by volume/mass of inactive particles in physical samples by comparing the stiffness results of static and dynamic testing.

Predicting the Timing of Catastrophic Failure During Triaxial Compression: Insights from Discrete Element Method Simulations

Examining the deformation of rocks during triaxial compression may provide insights into the precursory phase that leads to large earthquakes by revealing the components of the deformation field that evolve with predictable, systematic behavior preceding catastrophic failure. Here, we build three-dimensional discrete element method simulations of the triaxial compression of rock cores that include a variety of fault geometries in order to identify the components of the deformation field, including the velocity and strain components, that enable machine learning models to predict the timing of macroscopic failure. The results suggest that the velocity field provides more valuable information about the timing of macroscopic failure than the strain field, and in particular, the velocity component parallel to the maximum compression direction. The models also strongly depend on the second invariant of the strain deviator tensor, {J}_{2}, indicative of the shear strain. The importance of {J}_{2} on the model predictions increases with confining stress, consistent with laboratory observations that show a transition from tensile- to shear-dominated deformation with increasing confining stress. In contrast to expectations from previous machine learning analyses, none of the models strongly depend on components of the strain tensor indicative of dilation, such as the first invariant of the strain tensor. This difference may arise because the simulations host less dilative strain than the experimental rocks.

Comparison of clumps and rigid blocks in three-dimensional DEM simulations: curvature-based shape characterization

Clumps and rigid blocks are two established models in the discrete element method (DEM) for representing digitized shapes of granules. However, few studies compared their shape characterization effectiveness and computational efficiency in describing actual shapes with different accuracies. This research proposes a two-dimensional roundness index (RI2D) and concave index (CI) to quantify surface morphology using a curvature-based method. After collecting material shape information via laser scanning, the models were assessed through three-dimensional DEM simulations at the grain and specimen levels, with balls used as a reference. The findings confirm RI2D and CI as two reliable indicators of the degree to which clumps and rigid blocks vary at the angular scale. Based on RI2D, both clumps and rigid blocks can robustly and reasonably characterize irregular surfaces. However, further improving the accuracy of clumps while satisfying CI requires a significantly greater number of pebbles (spherical units); rigid blocks are incapable of capturing the concavity feature of real particles. RI2D and CI correlate well with the final state in column collapse tests and the shear strength in coupled DEM − FDM (finite difference method) triaxial tests. The difference in simulating granule mechanical behaviors between the two particle models is insignificant with varying levels of accuracy.

Unified penetration depth of low-velocity intruders into granular packings.

Penetration of intruders into granular packings is well described by separately considering the dry or wet case of granular environments in previous experiments and simulations; however, the unified description of such penetration depth in these two granular media remains elusive due to lacking clear explanations about its origins. Based on three-dimensional discrete element method simulations, we introduce a power-law fitting form of the final penetration depth of a spherical intruder with low velocity vertically penetrating into dry and wet granular packings, excellently expressed on a master curve as a power-law function of a dimensionless impact number that is defined as the square root of the ratio between the inertial stress of the intruder and the linear combination of the mean gravitational stress and the cohesive stress exerted on each grain in the packings, as a remarkable extension of the inertial number in dry granular flows. This scaling robustly provides physical insights inherent in the unified description of the material properties of granular packings and the impactor penetration conditions on the final penetration depth in the impact tests, providing evidence of impact properties in different disciplines and applications in science and engineering.

Lethe-DEM: an open-source parallel discrete element solver with load balancing

Approximately $${75}\%$$ of the raw material and $${50}\%$$ of the products in the chemical industry are granular materials. The discrete element method (DEM) provides detailed insights of phenomena at particle scale, and it is therefore often used for modeling granular materials. However, because DEM tracks the motion and contact of individual particles separately, its computational cost increases nonlinearly $$O(n_\mathrm{p}\log (n_\mathrm{p}))$$ – $$O(n_\mathrm{p}^2)$$ (depending on the algorithm) with the number of particles ( $$n_\mathrm{p}$$ ). In this article, we introduce a new open-source parallel DEM software with load balancing: Lethe-DEM. Lethe-DEM, a module of Lethe, consists of solvers for two-dimensional and three-dimensional DEM simulations. Load balancing allows Lethe-DEM to significantly increase the parallel efficiency by $$\approx {25}$$ – $${70}\%$$ depending on the granular simulation. We explain the fundamental modules of Lethe-DEM, its software architecture, and the governing equations. Furthermore, we verify Lethe-DEM with several tests including analytical solutions and comparison with other software. Comparisons with experiments in a flat-bottomed silo, wedge-shaped silo, and rotating drum validate Lethe-DEM. We investigate the strong and weak scaling of Lethe-DEM with $${1}\le n_\mathrm{c} \le {192}$$ and $${32}\le n_\mathrm{c} \le {320}$$ processes, respectively, with and without load balancing. The strong-scaling analysis is performed on the wedge-shaped silo and rotating drum simulations, while for the weak-scaling analysis, we use a dam-break simulation. The best scalability of Lethe-DEM is obtained in the range of $${5000}\le n_\mathrm{p}/n_\mathrm{c} \le {15{,}000}$$ . Finally, we demonstrate that large-scale simulations can be carried out with Lethe-DEM using the simulation of a three-dimensional cylindrical silo with $$n_\mathrm{p}={4.3}\times 10^6$$ on 320 cores.

A semi-empirical re-evaluation of the influence of state on elastic stiffness in granular materials

This study uses data acquired from three-dimensional discrete element method simulations to reconsider what measure of state can be used to predict stiffness in granular materials. A range of specimens with linear and gap-graded particle size distributions are considered and stiffness is measured using small amplitude strain probes. Analysis of the data firstly confirms that the void ratio, which is typically used as a measure of state in experimental soil mechanics, does not correlate well with shear stiffness. However, the empirical expressions developed by Hardin and his colleagues can capture variations in stiffness, provided an appropriate state variable is used. The study then highlights that the contribution of individual contacts to the overall stiffness is highly variable, depending on both the contact force transmitted and the particle size. Analyses explore how the stress transmission both within and between the different size fractions affects the overall stiffness. This heterogeneity in stiffness relates to the heterogeneity in the stress transmission amongst the different fractions. By considering the heterogeneity of stress distribution amongst different particle size fractions, a new semi-empirical stress-based state variable is proposed that provides insight into the factors that influence stiffness.

Three-dimensional DEM simulation of polydisperse particle flow in rolling mode rotating drum

The flow behavior of polydisperse particles in rolling mode rotating drum is important in industry, but it is still not clear. In this work, the discrete element method (DEM) is used to explore the flow of particles in the three-dimensional drum. After the accuracy of simulation verified by the experimental data from literature, the particle bed is divided into active layer and passive layer. The flow behavior of polydisperse particles was revealed by studying the particle velocity, residence time, mixing and axial dispersion in different regions in the drum. By exploring the parameters such as rotation speed, filling level and the shape and position of the lifter, some ideas are provided for improving the performance of the rolling mode rotating drum. In addition, the difference of flow behavior between polydisperse particles and monodisperse particles is also given, which can provide some reference for the industrial application of complex particle system.

The influence of particle size distribution on the stress distribution in granular materials

This study systematically explores the effect of the shape of the particle size distribution (PSD) on stress transmission in granular materials using three-dimensional discrete-element method simulations. Extending prior studies that have focused on bimodal mixtures of coarser and finer grains, a broad range of isotropically compressed specimens with spherical particles and linear, fractal and gap-graded particle size distributions is considered. Considering isotropic stress conditions, the nature of stress distribution was analysed by determining the mean effective particle stresses and considering the proportion of this stress transmitted by particles with different sizes. For gap-graded materials a contact-based perspective was adopted to consider the stress transmission both within and between the different size fractions. A clear correlation emerged between the cumulative distribution of particle sizes by volume and the cumulative distribution of particle sizes by mean effective stress for specimens with continuous PSDs. This correlation does not hold universally for gap-graded materials. In gap-graded materials the distribution of effective stress between the different size fractions depends upon the size ratio and the percentage of finer grains in the specimen. In contrast to specimens with continuous gradings, in gap-graded specimens the distribution of effective stress among the different size fractions exhibits a marked sensitivity to density. Basic network analysis is shown to provide a useful insight into effective stress transmission in bimodal gap-graded materials.

Three-dimensional Discrete Element Method simulation system of the interaction between irregular structure wheel and lunar soil simulant

The experimental execution, theoretical analysis, structural optimization and performance prediction of lunar rover wheels are subject to many limitations due to the particularity and complexity of the lunar environment. The numerical simulation has notable advantages in analyzing the interaction between wheels and lunar soil simulant and designing the lunar rover wheels applied to complex environments. In this paper, based on the theory of terramechanics, a three-dimensional (3D) dynamic system that simulated the interaction between irregular structure wheels and lunar soil simulant was established using the interface code and embedded FISH language of PFC3D. The rationality of the simulation system was verified by the wheel soil-bin test. Then, in order to improve and optimize the 3D dynamic simulation system, the interaction between wheel and lunar soil simulant was simulated in two different kinds of simulated lunar environments besides the flat ground, and the climbing property of the wheel was analyzed qualitatively. Finally, the mesoscopic changes of lunar soil stimulant particles were characterized by the color gradient variations, and the simulation results were processed in the IV program so that the distribution of lunar soil stimulant particles under the wheel could be observed more clearly. This research provided a reliable method for the study of the interaction between irregular wheel and loose lunar soil in complex environments and of designing the wheel structure.

Linking macroscopic failure with micromechanical processes in layered rocks: How layer orientation and roughness control macroscopic behavior

To constrain the impact of preexisting mechanical weaknesses on strain localization culminating in macroscopic shear failure, we simulate triaxial compression of layered sedimentary rock using three-dimensional discrete element method simulations. We develop a novel particle packing technique that builds layered rocks with preexisting weaknesses of varying orientations, roughness, and surface area available for slip. We quantify how the geomechanical behavior, characterized by internal friction coefficient, μ0, and failure strength, σF, vary as a function of layer orientation, θ, interface roughness, and total interface area. Failure of the simulated sedimentary rocks mirrors key observations from laboratory experiments on layered sedimentary rock, including minima σF and μ0 for layers oriented at 30° with respect to the maximum compressive stress, σ1, and maxima σF and μ0 for layers oriented near 0° and 90° to σ1. The largest changes in σF (66%) and μ0 (20%) occur in models with the smoothest interfaces and largest interface area. Within the parameter space tested, layer orientation exerts the most significant impact on σF and μ0. These simulations allow directly linking micromechanical processes observed within the models to macroscopic failure behavior. The spatial distributions of nucleating microfractures, and the rate and degree of strain localization onto preexisting weaknesses, rather than the host rock, are systematically linked to the distribution of failure strengths. Preexisting weakness orientation more strongly controls the degree and rate of strain localization than the imposed confining stress within the explored parameter space. Using the upper and lower limits of μ0 and σF obtained from the models, estimates of the Coulomb shear stress required for failure of intact rock within the upper seismogenic zone (7 km) indicates that a rotation of 30° of σ1 relative to the weakness orientation may reduce the shear stress required for failure by up to 100 MPa.

Three-dimensional discrete element method simulation of core disking

The phenomenon of core disking is commonly seen in deep drilling of highly stressed regions in the Earth’s crust. Given its close relationship with the in situ stress state, the presence and features of core disking can be used to interpret the stresses when traditional in situ stress measuring techniques are not available. The core disking process was simulated in this paper using the three-dimensional discrete element method software PFC3D (particle flow code). In particular, PFC3D is used to examine the evolution of fracture initiation, propagation and coalescence associated with core disking under various stress states. In this paper, four unresolved problems concerning core disking are investigated with a series of numerical simulations. These simulations also provide some verification of existing results by other researchers: (1) Core disking occurs when the maximum principal stress is about 6.5 times the tensile strength. (2) For most stress situations, core disking occurs from the outer surface, except for the thrust faulting stress regime, where the fractures were found to initiate from the inner part. (3) The anisotropy of the two horizontal principal stresses has an effect on the core disking morphology. (4) The thickness of core disk has a positive relationship with radial stress and a negative relationship with axial stresses.

Experimental validation study of 3D direct simple shear DEM simulations

Simple shear element tests can be used to examine numerous geotechnical problems; however, the cylindrical sample (NGI-type) direct simple shear (DSS) device has been criticized for its inability to apply uniform stresses and strains, as well as for its inability to fully define the stress state of the soil during shearing. Discrete element method (DEM) simulations offer researchers a means to explore the fundamental mechanisms driving the overall behavior of granular soil in simple shear and to improve the understanding of the DSS device itself. Here, three-dimensional DEM simulations of laminar NGI-type direct simple shear element tests and equivalent physical tests are compared to validate the numerical model. This study examines the sensitivity of the DEM simulation results to sample size, contact model and stiffness inputs, and ring wall boundary effects. Sample inhomogeneities are also considered by examining the radial and vertical void ratio distributions throughout the samples. Both the physical experiments and the DEM simulations presented herein indicate that the observed material response is highly sensitive to the particle size relative to the sample dimensions. The results show that samples with a small number of relatively large particles are very sensitive to small changes in packing; and thus, an exact match with the DEM simulation data cannot be expected. While increasing the number of particles greatly improved the agreement of the volumetric and stress–strain responses, the dense DEM samples were still initially much stiffer than the experimental results. This is most likely due to the fact that the inter-particle friction was artificially lowered, during the sample preparation for the DEM simulations, in order to increase the sample density.

Three-dimensional discrete element method simulation of the effect of bottom structure on solid flow in COREX shaft furnace

COREX is an industrially and commercially proven smelting reduction process. The shaft furnace (SF) for the pre-reduction of iron ore is one of the two major reactors of COREX. As one of the industrial systems using screw feeders, the burden descending behaviour in SF can be directly affected by the bottom structure such as the design of screw and associated charging container. In this work, a three-dimensional actual size model of COREX-3000 SF is established by means of the discrete element method. The effect of bottom structure, for example, the guiding cone and the design of screw, on solid flow is investigated. The results show that burdens are mainly drawn down from the first flight of the screw in the reference case, while the slow-moving particles located above the end of the screw. Placing a guiding cone and reducing the flight diameter of screw have a benefit for restraining the descending velocities in the central area of the SF and help to obtain a uniform flow pattern. An optimised case is also proposed. In the optimised case, a relatively uniform solid flow profile can be obtained and the uniformity of descending velocity along the radius is improved greatly. The findings of this work should be useful for the design, control and optimisation of the SF operation.

Local fluctuations and spatial correlations in granular flows under constant-volume quasistatic shear

We investigate the local fluctuations in dense granular media subjected to athermal, quasistatic shearing, based on three-dimensional discrete element method simulations. By shearing granular assemblies of different polydispersities under constant-volume constraint, we quantify the characteristics of local structures (in terms of local volume and local anisotropy) and local deformation (using local shear strain and nonaffine displacement). The distribution of the local volume in a granular medium is found unchanged during the entire shearing process, which indicates a constant temperaturelike compactivity for the material. The compactivity is not, however, equilibrated among different particle groups in a polydisperse assembly. The local structures of a disordered granular assembly are inherently anisotropic. The fluctuations in local anisotropy can be well captured by a gamma or mixed-gamma distribution function, which is also unchanged during the shear. The local anisotropic orientation evolves towards the coaxial direction of the stress anisotropy with shear. The deformation characteristics of a jammed granular medium have their origins in the structural amorphousness. The local shear strain field depicts clear shear transformation zones which act as plasticity carriers. The spatial correlation of the local shear strains exhibits a fourfold pattern which is stronger in the stress deviatoric planes than in the stress isotropic plane. The fluctuations of nonaffine displacement suggest an isotropic granular temperature and an isotropic spatial correlation independent of the stress state. Both the local strain and the nonaffine displacement exhibit a power-law decayed distribution with a long-range correlation. We further modify the shear-transformation-zone theory to predict the pressure-dependent constitutive behavior of a sheared granular material and compare its prediction with our simulation data.