## 3 Discrete Element Method  
In this paper, particle motion and collision are simulated using the discrete element method. Lagrangian particle's motion follows [PERSON]'s second law in Equations (14) and (15). Particle translation is accelerated by hydrodynamic interaction force \(\text{F}_{f}\), particles collision force \(\text{F}_{c}\), and gravity force \(\text{G}\) every timestep. Particle rotational acceleration depends on torque \(\text{M}_{f}\) and \(\text{M}_{c}\), derived from \(\text{F}_{f}\) and \(\text{F}_{c}\).  
\[m\frac{d\text{U}_{i}}{dt}=\text{F}_{f}+\text{F}_{c}+\text{G} \tag{14}\]  
\[I\frac{d\text{O}_{i}}{dt}=\text{M}_{f}+\text{M}_{c} \tag{15}\]  
In Equations (14) and (15), \(\text{U}_{i}\) and \(\text{O}_{i}\) are translational and rotational velocity, respectively; \(m\) and \(l\) are mass and moment of inertia, respectively; \(\text{F}_{f}\) is hydrodynamic force on particle calculated by immersed moving boundary, which will be discussed in the next section. \(\text{x}_{i}\) and \(\text{x}_{j}\) are two particles' center vector, \(r_{i}\) and \(r_{j}\) are two particles' radius. \(\text{v}_{i}\) and \(\text{v}_{j}\) are two particles' velocity. The collision surface normal and tangential vector of particle \(i\) are \(\text{n}=\big{(}\text{x}_{i}-\text{x}_{j}\big{)}/\big{|}\text{x}_{j}-\text{x} _{i}\big{|}\) and \(t\), respectively. Collision force in the normal and tangential direction is \(\text{F}_{n}\) and \(\text{F}_{f}\).  
\(\text{F}_{n}\) and \(\text{F}_{t}\) are divided into linear elastic force and damping force. For \(\text{F}_{n}\), linear elastic force is related to overlap of collision particle \(\delta_{ij}=\big{(}\text{r}_{i}+r_{j}\big{)}-\big{|}\text{x}_{j}-\text{x}_{i} \big{|}\). Damping force is related to the particle's relative velocity \(\dot{\delta}_{ij}=\big{(}\text{v}_{i}-\text{v}_{j}\big{)}\cdot\text{n}\). The particle i normal force is given by  
\[\text{F}_{n}=\big{(}k_{n}\delta_{ij}-\gamma_{n}\dot{\delta}_{ij}\big{)}\cdot \text{n} \tag{16}\]  
where constant \(k_{n}\) is normal stiffness coefficient. constant \(\gamma_{n}\) is normal damping coefficient. Tangential force \(\text{F}_{t}\) is caused by the rough surface sliding between two particles. Previous time step tangential overlap is forgotten in the next time step. The relative tangential velocity is \(\dot{\delta}_{ij}^{t}=\big{(}\text{v}_{i}-\text{v}_{j}\big{)}\cdot\text{t}- \big{(}r_{i}\Omega_{i}+r_{j}\Omega_{j}\big{)}\times\text{n}\), and the tangential overlap is \(\delta_{ij}^{t}=\dot{\delta}_{ij}^{t}\cdot dt\). The tangential stiffness coefficient is \(\text{k}_{t}=\frac{2}{7}\text{k}_{n}\) and the tangential damping coefficient is \(r_{t}=\frac{2}{7}\gamma_{n}\)[26]. The particle i tangential force is given by Equation (17).  
\[\text{F}_{t}=\big{(}k_{t}\delta_{ij}^{t}-\gamma_{t}\dot{\delta}_{ij}^{t} \big{)}\cdot t \tag{17}\]  
The tangential force is limited by the [PERSON]'s friction law: \(|\text{F}_{t}|<\mu\text{F}_{n}\), thereby the final tangential force expression is given by  
\[\text{F}_{t}=\begin{cases}\text{sgn}(\text{F}_{t})\cdot\mu\text{F}_{n},&| \text{F}_{t}|\geq\mu\text{F}_{n}\\ \text{F}_{t},&|\text{F}_{t}|<\mu\text{F}_{n}\end{cases} \tag{18}\]  
where \(\mu\) is friction coefficient.