## 2 DEM-FDM Coupled Mechanism  
The coupling of the DEM and FDM was realized by adopting interface elements between the discrete and continuous domains. Interface elements transfer the load from the discrete domain to the continuous domain and transfer the nodal velocity from the continuous domain at the interface to the discrete domain. The interface between the discrete and continuous domains is divided into triangles. A contact diagram for the discrete and continuous domains is presented in Figure 2[26].  
Figure 2 represents a ball particle in contact with a triangular interface element. \(C\) refers to the contact point and \(CP\) refers to the position of the point closest to \(C\) on the triangular interface element. The positions of the three vertices of the triangular interface element are denoted as \(x_{i}\) (\(i=1,2,3\)). The areas of the three triangles created by connecting the three vertices of the triangle at \(CP\) are denoted as \(A_{i}\) (\(i=1,2,3\)). The total area of the three triangles is denoted as \(A=\sum\limits_{i=1}^{3}A_{i}\).  
The weighting factors \(\omega_{i}\) for each vertex are determined by dividing the triangular area opposite from the vertex by the total area of the larger triangle as \(\omega_{i}=A_{i}/A\). Therefore, the expression \(\sum\omega_{i}=1\) is satisfied, which ensures that the sum of the values extrapolated from \(CP\) to the vertices is equal to the value at \(CP\).  
Here \(r_{i}\) (\(i=1,2,3\)) refers to the vectors pointing from \(CP\) to the triangle vertices (i.e., \(r_{i}=x_{i}-CP\)). Let the forces applied at the contact point \(C\) be denoted as \(\overset{\rightarrow}{F}\). Because the contact point \(C\) on the ball particle and the point \(CP\) on the interface element may not be at the same location, the total moment of the ball particle acting on the interface element is calculated as  
\[\overset{\rightarrow}{M}=(C-CP)\times\overset{\rightarrow}{F}. \tag{1}\]  
Take the forces applied on the triangle vertices to be \(F_{i}\) (\(i=1,2,3\)), then the following equations can be obtained:  
\[\sum F_{i}=\overset{\rightarrow}{F}, \tag{2}\]  
\[\sum\overset{\rightarrow}{r}_{i}\times F_{i}=\overset{\rightarrow}{M}. \tag{3}\]  
In the local coordinate system, the x component coincides with the normal direction of the triangle and the y component coincides with the tangential direction of the triangle. Therefore, in the local coordinate system, the equivalent force system on the vertices of the triangle can be expressed as  
\[\sum\overset{\rightarrow}{F}_{i,x}=\overset{\rightarrow}{F}_{x}, \tag{4}\]  
\[\sum\overset{\rightarrow}{F}_{i,y}=\overset{\rightarrow}{F}_{y}, \tag{5}\]  
\[\sum\overset{\rightarrow}{F}_{i,z}=\overset{\rightarrow}{F}_{z}=0, \tag{6}\]  
\[\sum\left(\overset{\rightarrow}{r}_{i,y}\times\overset{\rightarrow}{F}_{i,z }-\overset{\rightarrow}{r}_{i,z}\times\overset{\rightarrow}{F}_{i,y}\right) =\overset{\rightarrow}{M}_{x}, \tag{7}\]  
\[\sum\left(\overset{\rightarrow}{r}_{i,z}\times\overset{\rightarrow}{F}_{i, x}-\overset{\rightarrow}{r}_{i,x}\times\overset{\rightarrow}{F}_{i,z}\right) =\overset{\rightarrow}{M}_{y}, \tag{8}\]  
\[\sum\left(\overset{\rightarrow}{r}_{i,x}\times\overset{\rightarrow}{F}_{i, y}-\overset{\rightarrow}{r}_{i,y}\times\overset{\rightarrow}{F}_{i,x}\right) =\overset{\rightarrow}{M}_{z}. \tag{9}\]  
Additionally, the following constraints are introduced to derive a specific solution to the equations above:  
\[\sum(\overset{\rightarrow}{r}_{i,z}\times\overset{\rightarrow}{F}_{i,z})=0. \tag{10}\]  
The equations above are the basic governing equations for the adopted DEM-FDM coupled method. During the simulation process, the node velocity, displacement, and load information are continuously transmitted between the continuous and discrete domain to realize coupling calculations.  
Figure 2: Contact diagram for the coupled model.