## 2 Method
### Virtual Flux Method  
In this study, the virtual flux method [20; 21] was used to represent circular rigid particles on a Cartesian grid. This method has several advantages. It effectively captures physical quantities around objects, is easy to implement, operates faster than the commonly used immersed boundary method (IBM), and allows for arbitrary selection of interpolation accuracy based on the distance between particles [22; 23]. Figure 3 shows a schematic diagram of the virtual flux method.  
In the virtual flux method, a virtual distribution function is used at points near the virtual boundary points. This virtual distribution function is required to satisfy the boundary conditions at the virtual boundary. The no-slip boundary condition is applied to the velocity, while a zero pressure gradient in the normal direction is used as an approximate boundary condition for pressure, as indicated in Equations (23) and (24):  
\[\mathbf{u}_{\mathrm{vb}}=\mathbf{u}_{\mathrm{wall}}, \tag{23}\]  
\[\frac{\partial p}{\partial\mathbf{n}}=0. \tag{24}\]  
The distribution function at the virtual boundary is expressed using the following equation:  
\[f_{\mathrm{a}}(t,x_{\mathrm{vb}})=f_{\mathrm{a}}^{\mathrm{eq}}(t,x_{\mathrm{vb }})+\Big{\{}f_{\mathrm{a}}(t,x_{\mathrm{D}})-f_{\mathrm{a}}^{\mathrm{eq}}(t,x _{\mathrm{D}})\Big{\}}. \tag{25}\]  
Referring to Figure 3, the equilibrium distribution function at the virtual boundary point \(x_{\mathrm{vb}}\) can be obtained using the following equation:  
\[f_{\mathrm{a}}^{\mathrm{eq}}(t,x_{\mathrm{vb}})=\omega_{\mathrm{a}}\rho\Bigg{[} 1+\frac{3(\mathbf{e}_{\mathrm{a}}\cdot\mathbf{u}_{\mathrm{vb}})}{c^{2}}+\frac{9(\mathbf{e} _{\mathrm{a}}\cdot\mathbf{u}_{\mathrm{vb}})^{2}}{2c^{4}}-\frac{3\mathbf{u}_{\mathrm{vb }}^{2}}{2c^{2}}\Bigg{]}. \tag{26}\]  
where \(\mathbf{u}_{\mathrm{vb}}\) denotes the velocity vector at the virtual boundary point, while \(\mathbf{n}\) represents the normal vector to the virtual boundary. Additionally, the virtual distribution function \(f_{\mathrm{a}}^{*}(t,x_{\mathrm{E}})\) and the virtual equilibrium distribution function \(f_{\mathrm{a}}^{\mathrm{eq^{\prime}}}(t,x_{\mathrm{E}})\) at point E are expressed via linear extrapolation using the internal division ratios \(a\) and \(b\), as given in the following equations:  
\[f_{\mathrm{a}}^{*}(t,x_{\mathrm{E}})=\frac{a+b}{a}f_{\mathrm{a}}(t,x_{\mathrm{vb }})-\frac{b}{a}f_{\mathrm{a}}(t,x_{\mathrm{D}})\, \tag{27}\]  
Figure 3: Schematic of virtual flux method.  
\[f_{\alpha}^{\rm eq^{*}}(t,x_{\rm E})=\frac{a+b}{a}f_{\alpha}^{\rm eq}(t,x_{\rm vb })-\frac{b}{a}f_{\alpha}^{\rm eq}(t,x_{\rm D}) \tag{28}\]  
Because we used an incompressible formulation, the equilibrium function \(p_{\alpha}^{\rm eq}(t,x_{\rm vb})\) at the object surface is as follows:  
\[p_{\alpha}^{\rm eq}(t,x_{\rm vb})=\omega_{\alpha}\Bigg{[}p_{\rm vb}+\rho_{0} \Bigg{\{}(\vec{e}_{\alpha}\cdot\vec{u}_{\rm vb})+\frac{3(\vec{e}_{\alpha}\cdot \vec{u}_{\rm vb})^{2}}{2\epsilon^{2}}-\frac{\vec{u}_{\rm vb}^{2}}{2}\Bigg{\}} \Bigg{]}. \tag{29}\]  
where \(p_{\rm vb}\) represents the pressure at the virtual boundary. When the pressures \(p_{1}\) and \(p_{2}\) at points located at a distance \(h_{1}\) and \(h_{2}\) away from the boundary point in the normal direction to the wall are used, \(p_{\rm vb}\) can be expressed using the following equation:  
\[p_{\rm vb}=\frac{{h_{2}}^{2}p_{1}-{h_{1}}^{2}p_{2}}{{h_{2}}^{2}-{h_{1}}^{2}}. \tag{30}\]  
where \(p_{1}\) and \(p_{2}\) are calculated by interpolating the pressures at each of the four surrounding grid points. In this case, \(h_{1}\) and \(h_{2}\) were set to \(\sqrt{2}\) and \(2\sqrt{2}\) times the lattice size, respectively.