## 2 Numerical Method for Dataset Generation  
### Numerical Method  
LBM, as a mesoscopic method, can be considered as a numerical method independent of continuous kinetic theory. Compared with traditional macroscopic numerical methods, LBM has the advantages of easy coding implementation [46], good parallelism, easily handled geometrically complicated boundaries, and accurately describing the macroscopic motion phenomenon between gas-solid two phases at the microscopic level [46; 47]. Therefore, LBM is suitable for the polydisperse system of granular flow by considering its computational accuracy and efficiency of large complex fluid systems. So far, LBM is widely adopted as a solver to describe the characteristics of granular flow in which rigid granules are immersed in the continuum Newtonian fluid. In this context, we start with the Navier-Stokes equation, which is described by the lattice Boltzmann equation [48] as:  
\[f_{\alpha}(\mathbf{x}+\mathbf{e}_{\alpha}\delta_{t},t+\delta_{t})-f_{\alpha}(\mathbf{x},t )=-\frac{1}{\tau}\Big{[}f_{\alpha}(\mathbf{x},t)-f_{\alpha}^{eq}(\mathbf{x},t)\Big{]},\ \alpha=0,1,\ldots,n-1 \tag{1}\]  
where \(f_{\alpha}(\mathbf{x},t)\) is the distribution function with discrete velocity \(\mathbf{e}_{\alpha}\) at a position \(\mathbf{x}\) and a time \(t\), \(\delta t\) is the time step size, \(n\) is the total number of discrete velocities, and \(\tau\) is the relaxation time. \(f_{\alpha}^{eq}\) is the equilibrium distribution function, which depends on the density \(\rho\), velocity \(\mathbf{u}\), and temperature \(T\) of the fluid. \(f_{\alpha}^{eq}\)[48] can be described as:  
\[f_{\alpha}^{eq}=\omega_{\alpha}\rho\left(1+\frac{\mathbf{e}_{\alpha}\cdot\mathbf{u}}{c _{\rm s}^{2}}+\frac{\left(\mathbf{e}_{\alpha}\cdot\mathbf{u}\right)^{2}}{2c_{\rm s}^{4 }}-\frac{\mathbf{u}^{2}}{2c_{\rm s}^{2}}\right),\ \alpha=0,1,\ldots,n-1 \tag{2}\]  
where \(\omega_{\alpha}\) is the model-dependent weight coefficient, \(c_{\rm s}=RT\) (\(R\) is the gas constant) is the lattice sound velocity. The fluid density \(\rho\) and velocity \(\mathbf{u}\) can be derived as the zero-th and first order moments of \(f_{\alpha}\), respectively, as:  
\[\rho=\sum_{i=0}^{n-1}f_{\alpha}, \tag{3}\]  
\[\rho\mathbf{u}=\sum_{i=0}^{n-1}\mathbf{e}_{i}f_{i}, \tag{4}\]  
The fluid pressure is defined directly as \(p=c_{\rm s}^{2}\rho\). Through the Chapman-Enskog expansion, the Navier-Stokes equations can be obtained with the viscosity [49] given by:  
\[v=\rho c_{\rm s}^{2}\left(\tau-\frac{1}{2}\right)\delta_{t}. \tag{5}\]  
According to the definition of the Reynolds number \(Re=\rho UL/v\), the relaxation time \(\tau\)[50] is calculated by:  
\[\tau=\frac{UL}{c_{\rm s}^{2}\delta_{t}Re}+\frac{1}{2}, \tag{6}\]  
where \(U\) and \(L\) are the characteristic velocity and length, respectively. To obtain the drag force of a polydisperse system of granular flow, the structure of a polydisperse granular cluster in the flow field was randomly generated. Regarding the boundary condition, our model was referred to the LBM model of [PERSON] et al. [51]. The inlet boundary condition of the fluid field was set to the constant velocity, which was \(U=0.1\) m/s. At the outlet boundary, a stress-free condition was set, and the periodic boundary conditions were applied to the other two directions. For these settings, as shown in Figure 1, the flows horizontally from left to right with a steady state were gradually achieved, which means the flow field reaches the equilibrium state at each moment. Considering our previous study [23], a \(300\times 100\times 100\) grid was employed, and a single-relaxation-time collision operator in conjunction with a D3Q19 lattice was adopted  
Figure 1: The sketch of polydisperse system of granular flow. for the accuracy and stability of the model. In order to improve the computational efficiency of a polydisperse system, a parallel algorithm (MPI+LBM) was used to implement the numerical simulation in the flow field.