# A Simulation Method for Layered Filling of Grain Piles Based on the Discrete Element Method
## 2 Materials and Methods
### Layered Filling Method  
In DEM simulations, a significant amount of computational time is dedicated to particle collision detection [26]. The objective of this algorithm is to reduce collision detection and optimize computation time. Drawing upon the Janssen effect [27; 28; 29], when a grain pile reaches a certain height, the influence of further filling at that height on the lower part of the pile can be disregarded. As a result, the need for contact detection during the calculation process can be reduced. The removed section of the pile is referred to as a sub-pile, while the remaining section is referred to as the main pile. The entire Implementation process is shown in Figure 1.  
The typical DEM simulation loop can be summarized as follows: collision detection, solving interactions, updating, and miscellaneous tasks [30]. The LFM, as depicted in Figure 2, operates between the updating and miscellaneous tasks. The specific steps of the LFM are outlined below:  
1. At intervals, the LFM traverses all particle properties and labels particles that meet kinetic energy conditions. It marks particles that meet specified criteria (such as Kinetic energy);
2. Once all particles within a certain range surrounding a labeled particle have been labeled, this particle is considered to be stable; \begin{table}
\begin{tabular}{l c} \hline \hline
**Parameters** & **Equations** \\ \hline Equivalent contact radius & \(\frac{1}{K}=\frac{1}{K}+\frac{1}{K^{\prime}}\) \\ Equivalent Young’s modulus & \(\frac{1}{K^{\prime}}=\frac{1-\frac{1}{K^{\prime}}}{E_{i}}+\frac{1-\frac{1}{K^ {\prime}}}{e^{\prime}_{i}}\) \\ Equivalent mass & \(\frac{1}{m^{\prime}}=\frac{1}{m_{i}}+\frac{1}{m_{i}}\) \\ Normal stiffness & \(S_{n}=2E^{*}\sqrt{R^{*}a}n\) \\ Damping coefficient & \(\beta=\frac{lue}{\sqrt{lu^{2}+\pi^{2}}}\) \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Fundamental equations. Figure 1: Implementation process of LFM. 3. The highest-positioned particle is located, and particles below this position are separated to form a sub-pile. The base of the pile consists of a static particle surface formed by the upper layer of particles from the former sub-pile;
4. After the piling process, the sub-piles are merged with the main pile and gravity settling is allowed until the imbalance forces meet the predetermined requirements. Figure 2: LFM simulation loop. Kinetic energy can reflect the motion state of particles. In this study, kinetic energy is selected as the criterion for determining particle stability. Kinetic energy can reflect the motion state of particles to some extent. Subsequently, the particles within a specific range surrounding the labeled particles need to be examined. To accomplish this, a simplified search algorithm proposed by [PERSON] is utilized [31], which leverages GPU acceleration for faster computations. The following illustrates how the algorithm divides the computational domain into multiple buckets with defined widths based on Equation (5):  
\[l_{b}=c_{b}\times r_{max} \tag{5}\]  
where \(c_{b}\) is a coefficient of bucket width determined by particle movement and \(r_{max}\) is the radius of the maximum particle. Figure 3 illustrates the process of determining the state of particles. The solid wire-frames represent particles allocated to different block areas. When searching for neighboring particles, the search group is limited to particles within nearby bins, thereby enhancing search efficiency. For example, to determine whether the red particle in Figure 3a is in a stable state, it is necessary to determine whether the particle itself and the particles around it that can affect its state are also in a stable or quasi-stable state. The involved block areas are the block areas where the particle is located, as shown in Figure 3a, and the other 8 block areas surrounding it.