## 2 Method
### Governing Equation for Fluid  
In this study, the regularized lattice Boltzmann method was used for the governing equation of the fluid [18]. This method retains the simplicity and efficiency of the original lattice Boltzmann approach while ensuring computational stability at high Reynolds numbers. Furthermore, it considerably decreases memory usage during the analysis. For this study, we adopted one of the most widely used models, the D2Q9 model. The division of this model is shown in Figure 2. In this model, the space is segmented into square lattices, where microscopic particles at a lattice point \((i,j)\) either remain at the same point \((i,j)\) or move in the direction indicated by the arrows shown in Figure 2 over one time step. The discrete velocities of the particles are defined by the following equation:  
\[\mathbf{e}_{\alpha}=\left\{\begin{array}{cl}c(0,0)&(\alpha=0)\\ c\left(\cos\frac{(\alpha-1)\Pi}{2},\sin\frac{(\alpha-1)\Pi}{2}\right)&(\alpha= 1\text{--}4).\\ \sqrt{2}c\left(\cos\frac{(2\alpha-9)\Pi}{2},\sin\frac{(2\alpha-9)\Pi}{2}\right) &(\alpha=5\text{--}8)\end{array}\right. \tag{2}\]  
Figure 1: A schematic of the computational model showing circular rigid particles suspended in a two-dimensional channel with random distribution.  
Here, \(\alpha\) represents the tensor index and corresponds to the movement direction of virtual particles in the D2Q9 model, and \(c\) represents the advection velocity, defined as \(c=\delta x/\delta t\), where \(\delta t\) is the time step and \(\delta x\) is the spatial step. When the Boltzmann equation is discretized using the discrete velocity \(\mathbf{e}_{\alpha}\) and the collision relaxation time \(\tau\), the following equation is derived:  
\[f_{\alpha}(t+\delta t,\mathbf{x}+\mathbf{e}_{\alpha}\delta t)=f_{\alpha}(t,\mathbf{x})+ \frac{1}{\tau}\Big{(}f_{\alpha}^{\text{eq}}(t,\mathbf{x})-f_{\alpha}(t,\mathbf{x}) \Big{)}. \tag{3}\]  
where \(f_{\alpha}\) denotes the distribution function, while \(f_{\alpha}^{\text{eq}}\) represents the equilibrium distribution function. By applying the Chapman-Enskog expansion to the abovementioned equation, the Navier-Stokes equation can be derived. At this point, the relaxation time \(\tau\) is defined by the following equation:  
\[\tau=\frac{3\
u}{c\delta\mathbf{x}}+\frac{1}{2}. \tag{4}\]  
where \(\
u\) represents the kinematic viscosity. The equilibrium distribution function is expressed as follows:  
\[f_{\alpha}^{\text{eq}}=\omega_{\alpha}\rho\Bigg{\{}1+\frac{3(\  
The nonequilibrium part of the pressure distribution \(p_{\alpha}^{\rm neq}\) and the nonequilibrium part of the stress tensor \(\Pi_{ij}^{\rm neq}\) can be expressed as follows:  
\[p_{\alpha}^{\rm neq}(t,\mathbf{x})=p_{\alpha}(t,\mathbf{x})-p_{\alpha}^{\rm eq}(t,\mathbf{x}), \tag{11}\]  
\[\Pi_{ij}^{\rm neq}(t,\mathbf{x})=\Pi_{ij}(t,\mathbf{x})-\Pi_{ij}^{\rm eq}(t,\mathbf{x}). \tag{12}\]  
By applying the Chapman-Enskog expansion, the nonequilibrium part of the pressure distribution function can be transformed as follows:  
\[p_{\alpha}^{\rm neq}(t,\mathbf{x})\cong-\delta t\tau\omega_{\alpha}Q_{\alpha ij} \partial_{i}\rho u_{j}. \tag{13}\]  
where \(Q_{\alpha ij}\) is the normalized tensor, defined as \(Q_{\alpha ij}=e_{\alpha i}e_{\alpha j}-c_{\alpha}^{2}\partial_{ij}\). Furthermore, using Equations (9) and (12), the nonequilibrium part of the stress tensor can be expressed as follows:  
\[\Pi_{ij}^{\rm neq}(t,\mathbf{x})\cong-\delta tc_{\rm s}^{2}\tau\big{(}\partial_{ ij}\rho u_{j}+\partial_{j}\rho u_{i}\big{)}. \tag{14}\]  
From Equations (11) and (12), \(p_{\alpha}^{\rm neq}\) can be expressed as follows:  
\[p_{\alpha}^{\rm neq}(t,\mathbf{x})\cong\frac{\omega_{\alpha}}{2c_{\rm s}^{2}}Q_{ \alpha ij}\Pi_{ij}^{\rm neq}(t,\mathbf{x}), \tag{15}\]  
\[p_{\alpha}(\,t+\delta t\,,\,\mathbf{x}+\mathbf{e}_{\alpha}\delta t\,)=p_{\alpha}^{\rm eq }(t,\mathbf{x})+\bigg{(}1-\frac{1}{\tau}\bigg{)}p_{\alpha}^{\rm neq}(t,\mathbf{x}). \tag{16}\]