This paper introduces the Enhanced Minkowski Portal Refinement (EMPR) algorithm, a novel contact detection method for convex polygons in discrete element modeling (DEM). The key contributions are: (1) Extending the Minkowski Portal Refinement (MPR) algorithm, originally developed for collision detection in computer graphics, to accurately compute contact features between convex shapes. (2) Identifying and addressing an interior point dependency issue in MPR, which can lead to inconsistent final portals and inaccuracies in contact features and energy conservation. (3) Proposing the EMPR algorithm, which includes a criterion to detect potentially inaccurate final portals and a corrective scheme to ensure precise results. (4) Demonstrating that EMPR outperforms the combined GJK and EPA methods for contact detection in DEM by providing accurate contact states and features in a unified, simple, efficient, and easy-to-implement manner. Numerical examples with various polygonal shapes and particle numbers validate the algorithm's correctness, robustness, and effectiveness.

Efficiently generating realistic 3D particle assemblies with a desired particle-size distribution is a significant challenge, primarily due to complex internal and external geometric constraints. This study presents a method for densely packing spheres into arbitrary 3D shapes. By explicitly accounting for the boundary domain, the method substantially reduces gaps and voids at the surface of the final assembly. The proposed approach can process arbitrarily complex geometries imported from an STL file via a three-step procedure: (1) generating an initial particle packing, (2) compressing the initial packing, and (3) refilling any remaining space. Several examples are presented to demonstrate the method's high efficiency, low hardware requirement, and strong adaptability for complex geometries. Notably, the algorithm generated a packing of 1,084,791 spheres within a cylindrical domain at a density of 58.6% in just 1218.2 s using a 2GHz laptop.

Efficiently assembling densely packed discs in complex domains remains an unresolved challenge for discrete element simulations. This study develops an efficient dense disc packing method based on the advancing front method (AFM) with closed form version. The new method enables a unified and efficient generation of intricate packings incorporating not just interior boundaries with multiple cavities, but also convex or non-convex external boundaries with complex geometric shapes. Furthermore, the new method can eliminate potential discontinuity encountered in certain DEM assemblies, such as the one-sided lifting problem that arises during discs packing using AFM with open form version. Several examples with complex geometries are presented to showcase the versatility, capability and performance of the proposed method. The proposed approach not only offers an efficient means of generating DEM packings, but also empowers DEM for enhanced simulations with complex geometric boundaries. Specifically, it takes only 70.92 s to generate 675.3 million discs on a 2 GHz laptop.

Contact modelling for ellipses/ellipsoids within the discrete element framework has long been a challenging area. Currently, four methods, namely the intersection method, geometric potential method, midway method, and common-normal method, have been proposed to define contact features between ellipses. This study addresses two key objectives: first, to evaluate whether any of these methods are energy-conserving, and second, to propose an accurate, stable and efficient contact detection framework based on the common-normal method for 2D problems. In this framework, the Gilbert-Johnson-Keerthi (GJK) algorithm is used for contact detection between ellipses, while a hybrid expanding polytope algorithm (EPA)-Levenberg-Marquardt (LM) approach is employed to compute contact features. The proposed framework is verified through two numerical examples focusing on energy conservation. Its robustness and applicability are demonstrated by modelling an illustration example involving a large number of ellipses. Detailed computational performance and advantages of the proposed framework are also provided.

The high-resolution DEM-IMB-LBM model can accurately describe pore-scale fluid-solid interactions, but its potential for use in geotechnical engineering analysis has not been fully unleashed due to its prohibitive computational costs. To overcome this limitation, a message passing interface (MPI) parallel DEM-IMB-LBM framework is proposed aimed at enhancing computation efficiency. This framework utilises a static domain decomposition scheme, with the entire computation domain being decomposed into multiple subdomains according to predefined processors. A detailed parallel strategy is employed for both contact detection and hydrodynamic force calculation. In particular, a particle ID re-numbering scheme is proposed to handle particle transitions across sub-domain interfaces. Two benchmarks are conducted to validate the accuracy and overall performance of the proposed framework. Subsequently, the framework is applied to simulate scenarios involving multi-particle sedimentation and submarine landslides. The numerical examples effectively demonstrate the robustness and applicability of the MPI parallel DEM-IMB-LBM framework.

Multifield coupling is frequently encountered and also an active area of research in geotechnical engineering. In this work, a particle-resolved direct numerical simulation (PR-DNS) technique is extended to simulate particle-fluid interaction problems involving heat transfer at the grain level. In this extended technique, an immersed moving boundary (IMB) scheme is used to couple the discrete element method (DEM) and lattice Boltzmann method (LBM), while a recently proposed Dirichlet-type thermal boundary condition is also adapted to account for heat transfer between fluid phase and solid particles. The resulting DEM-IBM-LBM model is robust to simulate moving curved boundaries with constant temperature in thermal flows. To facilitate the understanding and implementation of this coupled model for non-isothermal problems, a complete list is given for the conversion of relevant physical variables to lattice units. Then, benchmark tests, including a single-particle sedimentation and a two-particle drafting-kissing-tumbling (DKT) simulation with heat transfer, are carried out to validate the accuracy of our coupled technique. To further investigate the role of heat transfer in particle-laden flows, two multiple-particle problems with heat transfer are performed. Numerical examples demonstrate that the proposed coupling model is a promising high-resolution approach for simulating the heat-particle-fluid coupling at the grain level.

In this work, a Minkowski difference‐based advancing front approach is proposed to generate convex and non‐circular particles in a predefined computational domain. Two specific algorithms are developed to handle the contact conformity of generated particles with the boundaries of the computational domain. The first, called the open form, is used to handle the smooth contact of generated particles with (external) boundaries, while the other, called the closed form, is proposed to handle the internal boundaries of a computational domain with a complex cavity. The Gilbert‐Johnson‐Keerthi (GJK) method is used to efficiently solve the contact detection between the newly generated particle at the front and existing particles. Furthermore, the problem of one‐sided particle lifting, which can cause some defects in the packing structure in existing advancing front methods during packing generation, is highlighted and an effective solution is developed. Several examples of increasing complexity are used to demonstrate the efficiency and applicability of the proposed packing generation approach. The numerical results show that the generated packing is not only more uniform, but also achieves a higher packing density than existing advancing front methods.

The coupled discrete element and lattice Boltzmann method using an immersed moving boundary scheme was extended to simulate methane hydrate exploitation involving mass transport and particle dissolution. In this coupled DEM‐IMB‐LBM model, a new Dirichlet‐type thermal boundary condition is extended to simulate moving curved boundaries with constant concentration. A novel periodic boundary including an efficient searching algorithm for particle contact is proposed to reduce the computational cost and boundary effect. Then this model is validated by two numerical examples: a circular particle with concentration convection‐diffusion moving in a horizontal channel and mass transport from a cylinder particle in a simple shear flow. The numerical results obtained from the proposed model agree well with previous studies. To further demonstrate the capacity of the proposed model, simulations of methane hydrate exploitation including two formations in marine sediments are carried out. The numerical results indicate that the coupled DEM‐IMB‐LBM is not only capable of simulating the dissolution of hydrate particles at the grain level, but also recover the sand erosion and migration process in a fundamental perspective during the methane hydrate exploitation process.

在计算域内生成密实圆形颗粒集合体是离散单元法模拟中的一个重要问题.颗粒生成算法中,波前法(包括closed型和open型)作为一种纯几何算法,由于其高效性得到了广泛应用.然而,波前法生成的颗粒集合体在域边界上存在较大的间隙,从而导致边界的非光滑性.为克服这一问题,将closed型波前法中处理外域边界方法拓展到open型波前法.该算法通过Netwon-Raphson迭代方法得到颗粒的中心坐标和半径,从而保证新生成的颗粒与域边界相切,最终消除了边界间隙和凹凸性.针对open型波前法在堆积过程中出现的右端抬升问题,给出了解决方法,消除了潜在的优势结构面,进一步提高了计算效率.结果表明:提出的closed型和open型算法对不同的内外域边界都有很好的适用性.相比其他颗粒生成算法,open型算法效率非常高,其效率分别是离散元模拟软件和closed型的20500倍和4.4倍.在考虑域边界情况下,open型波前法在2.3 GHz的笔记本电脑上生成410万密实圆形颗粒集合体只需0.9 s.

SummaryThe efficiency of solving equations plays an important role in implicit‐scheme discontinuous deformation analysis (DDA). A systematic investigation of six iterative methods, namely, symmetric successive over relaxation (SSOR), Jacobi (J), conjugate gradient (CG), and three preconditioned CG methods (ie, J‐PCG, block J‐PCG [BJ‐PCG], and SSOR‐PCG), for solving equations in three‐dimensional sphere DDA (SDDA) is conducted in this paper. Firstly, simultaneous equations of the SDDA and iterative formats of the six solvers are presented. Secondly, serial and OpenMP‐based parallel computing numerical tests are done on a 16‐core PC, the result of which shows that (a) for serial computing, the efficiency of the solvers is in this order: SSOR‐PCG > BJ‐PCG > J‐PCG > SSOR>J > CG, while for parallel computing, BJ‐PCG is the best solver; and (b) CG is not only the most sensitive to the ill‐condition of the equations but also the most time consuming under both serial and parallel computing. Thirdly, to estimate the effects of equation solvers acting on SDDA computations, an application example with 10 000 spheres and 200 000 calculation steps is simulated on this 16‐core PC using serial and parallel computing. The result shows that SSOR‐PCG is about six times faster than CG for serial computing, while BJ‐PCG is about four times faster than CG for parallel computing. On the other hand, the whole computation time using BJ‐PCG for parallel computing is 3.37 hours (ie, 0.061 s per step), which is about 36 times faster than CG for serial computing. Finally, some suggestions are given based on this investigation result.

A new C++ programming strategy with high modularization and good portability, and a novel data storage format for simultaneous equations with little computer memory consumption, no sorting operation, and simple addressing algorithm are proposed for the three-dimensional sphere discontinuous deformation analysis (3D SDDA) to overcome the shortcomings of existing computation programs. An object-oriented data structure for the 3D SDDA computing code that is highly modular and easily transplanted is designed. Then, to demonstrate the portability of the 3D SDDA computing code, two computation architectures are respectively constructed to form two independent computation programs for 3D SDDA. Finally, several benchmark tests are conducted to verify the correctness of the 3D SDDA model in the new computation program, and a 170,725-sphere landslide example is simulated on a desktop computer to demonstrate the capability of the new computation program in large-scale engineering applications. Comparison between the new and existing computation programs regarding computer memory and time consumed demonstrates the great advantages brought about by the new computation program.

Simulating large-scale problems are still challenging for discontinuous deformation analysis (DDA). To this end, this paper develops an efficient disk-based DDA (DDDA) model considering the efficiency in searching contact pairs and solving linear equations. First, an efficient contact search algorithm, namely the lattice search algorithm (LSA), is proposed, and efficiency tests of the LSA and direct search algorithm (DSA) demonstrate the high efficiency of the LSA. Second, three equations solvers, namely Jacobi iterative method (J), conjugate gradient method (CG), and preconditioned CG (PCG), are adopted to respectively solve the equations of the DDDA, and efficiency tests of these solvers show that the best solver is the PCG and the J is unsuited to solving the equations when using large penalty spring stiffness. Finally, a landslide simulation which includes 30,000 disks, 5,990 line segments, and 180,000 calculation steps is conducted, the result of which shows that: (1) up to 41.1 h, which is 99.2% of the time consumed in contact search using the DSA, is reduced by using the LSA; (2) up to 18.22 h, which is 69.8% of the time consumed in solving the equations using the CG, is reduced by using the PCG as the equation solver; (3) the time consumed in the simulation is 68.72 h when using the DSA to search contact pairs and using the CG to solve the equations; while the consumed time reduces to 9.4 h when using the LSA to search contact pairs and using the PCG to solve the equations, whose reduction proportion is 86.3%. The simulation indicates the large-scale computation capacity and further application in engineering with the DDDA.

In order to investigate rock mechanical behavior under cyclic loading and unloading conditions, the clump parallel-bond model (CPBM) combined with the newly-proposed loading procedure was used in numerical simulations and the simulated results were compared with those of the experiment. Then the intact rock under different confining pressures and pre-cracked rock with different flaw inclination angles are also investigated in detail. The related numerical results have demonstrated that: (1) the CPBM, combined with the use of the newly-proposed loading procedure, can be used to modelling the intact rock and pre-cracked rock mechanical behavior under cycle loading and unloading conditions. The simulated results are in good agreement with the experimental results; (2) the number of microscopic cracks increases with the cycle increases during the loading process, while few microscopic cracks are generated during the unloading process. The higher the confining pressure, the less damage the rock has accumulated during the cyclic loading process; (3) the inclination of the pre-existing flaw has a significant effect on the macroscopic fracture pattern.

The particle simulation method is used to study the effects of loading waveforms (i.e. square, sinusoidal and triangle waveforms) on rock damage at mesoscopic scale. Then some influencing factors on rock damage at the mesoscopic scale, such as loading frequency, stress amplitude, mean stress, confining pressure and loading sequence, are also investigated with sinusoidal waveform in detail. The related numerical results have demonstrated that: 1) the loading waveform has a certain effect on rock failure processes. The square waveform has the most damage within these waveforms, while the triangle waveform has less damage than sinusoidal waveform. In each cycle, the number of microscopic cracks increases in the loading stage, while it keeps nearly constant in the unloading stage. 2) The loading frequency, stress amplitude, mean stress, confining pressure and loading sequence have considerable effects on rock damage subjected to cyclic loading. The higher the loading frequency, stress amplitude and mean stress, the greater the damage the rock accumulated; in contrast, the lower the confining pressure, the greater the damage the rock has accumulated. 3) There is a threshold value of mean stress and stress amplitude, below which no further damage accumulated after the first few cycle loadings. 4) The high-to-low loading sequence has more damage than the low-to-high loading sequence, suggesting that the rock damage is loading-path dependent.

Purpose -The main purpose of this paper is to present a comprehensive upscale theory of the thermo-mechanical coupling particle simulation for three-dimensional (3D) large-scale non-isothermal problems, so that a small 3D length-scale particle model can exactly reproduce the same mechanical and thermal results with that of a large 3D length-scale one.Design/methodology/approach -The objective is achieved by following the scaling methodology proposed by Feng and Owen (2014).Findings -After four basic physical quantities and their similarity-ratios are chosen, the derived quantities and its similarity-ratios can be derived from its dimensions. As the proposed comprehensive 3D upscale theory contains five similarity criteria, it reveals the intrinsic relationship between the particle-simulation solution obtained from a small 3D length-scale (e.g. a laboratory length-scale) model and that obtained from a large 3D length-scale (e.g. a geological length-scale) one. The scale invariance of the 3D interaction law in the thermo-mechanical coupled particle model is examined. The proposed 3D upscale theory is tested through two typical examples. Finally, a practical application example of 3D transient heat flow in a solid with constant heat flux is given to illustrate the performance of the proposed 3D upscale theory in the thermo-mechanical coupling particle simulation of 3D large-scale non-isothermal problems. Both the benchmark tests and application example are provided to demonstrate the correctness and usefulness of the proposed 3D upscale theory for simulating 3D non-isothermal problems using the particle simulation method.Originality/value -The paper provides some important theoretical guidance to modeling 3D large-scale non-isothermal problems at both the engineering length-scale (i.e. the meter-scale) and the geological lengthscale (i.e. the kilometer-scale) using the particle simulation method directly.

Coal and gas outburst is a very complex dynamic disaster in underground coal mining process. In this paper，a numerical model was established based on the theory of particle flow to simulate the development of micro-cracks，displacement，force and velocity fields，and to investigate the microscopic mechanism of gas pressure and layer-stiffness ratio. The simulation results show that the coal and gas outburst is a relatively quick process. The shear cracks were mainly concentrated in the front tip of outburst，the tensile cracks occurred deep inside the coal. The gas outburst has great influence on the coal and rock damage. When the gas pressure is relatively small，the shear cracks occur at the front and the tensile crack reached deeper. When the gas pressure is relatively large，the shear cracks and tensile cracks penetrate into the same depth. The damage of coal and rock and the shear crack ratio increases. When the layer-stiffness ratios is not the same，the speed and shape of crack propagation are not the same.

A thermo-mechanical coupled particle model for simulation of thermally-induced rock damage based on the particle simulation method was proposed. The simulation results of three verification examples, for which the analytical solutions are available, demonstrate the correctness and usefulness of the thermo-mechanical coupled particle model. This model is applied to simulating an application example with two cases: one is temperature-independent elastic modulus and strength, while the other is temperature-dependent elastic modulus and strength. The related simulation results demonstrate that microscopic crack initiation and propagation process with consideration of temperature-independent and temperature-dependent elastic modulus and strength are different and therefore, the corresponding macroscopic failure patterns of rock are also different. On the contrary, considering the temperature-dependent elastic modulus and strength has no or little effect on the heating conduction behavior. Numerical results, which are obtained by using the proposed model with temperature-dependent elastic modulus and strength, agree well with the experimental results. This also reveals that the rock subjected to heating experiences much more cracking than the rock subjected to cooling.

Purpose – The purpose of this paper is to present an upscale theory of the thermal-mechanical coupling particle simulation for non-isothermal problems in two-dimensional quasi-static system, under which a small length-scale particle model can exactly reproduce the same mechanical and thermal results with that of a large length-scale one. Design/methodology/approach – The objective is achieved by extending the upscale theory of particle simulation for two-dimensional quasi-static problems from an isothermal system to a non-isothermal one. Findings – Five similarity criteria, namely geometric, material (mechanical and thermal) properties, gravity acceleration, (mechanical and thermal) time steps, thermal initial and boundary conditions (Dirichlet/Neumann boundary conditions), under which a small-length-scale particle model can exactly reproduce both the mechanical and thermal behavior with that of a large length-scale model for non-isothermal problems in a two-dimensional quasi-static system are proposed. Furthermore, to test the proposed upscale theory, two typical examples subjected to different thermal boundary conditions are simulated using two particle models of different length scale. Originality/value – The paper provides some important theoretical guidances to modeling thermal-mechanical coupled problems at both the engineering length scale (i.e. the meter scale) and the geological length scale (i.e. the kilometer scale) using the particle simulation method directly. The related simulation results from two typical examples of significantly different length scales (i.e. a meter scale and a kilometer scale) have demonstrated the usefulness and correctness of the proposed upscale theory for simulating non-isothermal problems in two-dimensional quasi-static system.

A numerical model based on the lattice Boltzmann method is presented to investigate the viscous fingering phenomena of miscible displacement processes in porous media, which involves the fluid flow, heat transfer and mass transport. Especially, temperature- and concentration-dependent pore-fluid viscosity is considered. A complete list is derived and given for the conversion of relevant physical variables to lattice units to facilitate the understanding and implementation of the coupled problems involving fluid flow, heat transfer and mass transport using the LBM. To demonstrate the proposed model capacity, two different complex geometry microstructures using high resolution micro-computed tomography (micro-CT) images of core sample have been obtained and incorporated as computation geometries for modeling miscible displacement processes in porous media. The viscous fingering phenomena of miscible displacement processes are simulated in two different cases, namely in a channel and a porous medium respectively. Some influencing factors on the miscible displacement process, such as the pore-scale microstructure, Le number and Re number, are studied in great detail. The related simulation results have demonstrated that: (1) the existence of the pore-scale microstructure can have a significant effect on the front morphologies and front propagation speed in the miscible displacement process; (2) as the Le number increases, the fluid front and thermal front evolve differently, with the thermal front being less unstable due to more diffusion; (3) a larger Re number can lead to an increase in the propagation speed of the front.

Based on the particle simulation method, a thermo-mechanical coupling particle model is proposed for simulating thermally-induced rock damage. In this model, rock material is simulated as an assembly of particles, which are connected to each other through their bonds, in the case of simulating mechanical deformation, but connected to each other through thermal pipes in the case of simulating heat conduction. The main advantages of using this model are that: (1) microscopic parameters of this model can be directly determined from the related macroscopic ones; (2) the temperature-dependent elastic modulus and strength are considered in an explicit manner, so that thermally-induced rock damage can be realistically simulated in a thermo-mechanical coupling problem. The related simulation results from an application example have demonstrated that: (1) the proposed model can produce similar behaviors to those observed in experiments; (2) the final failure is initiated from the outer surface of the testing sample and propagates toward the borehole; (3) microscopic crack initiation and propagation processes can be reasonably simulated at the cooling stage.