## 2 Numerical Methods
### Virtual Flux Method  
The virtual flux method [22; 23] was applied to describe the curved boundary of particle shape on a Cartesian grid. Suspension flows including solid particle can be solved by reflecting the distribution functions calculated using the virtual flux method. One of the advantages of this method is simple algorism with no iterative calculations. The distribution function on the particle surface can be obtained by simply considering its hydrodynamic conditions. These procedures need no iterations. Compared to the immersed boundary method (IBM) [24], which is popularly coupled with LBM, the virtual flux method has superior spatial resolution without adding external force to the body surfaces to describe the arbitrary shape on the Cartesian grid.  
Figure 3 shows the virtual boundary on the Cartesian grid. For example, the distribution function streams to nine directions in the D2Q9 model; however, the distribution function at grid point 3 cannot propagate to grid point 1 (described as direction \(\alpha\)) because the virtual boundary point \(\rm w\) exists at the position where \(a\):\(b\) divides the space between grid points 1 and 3.  
A brief description of the calculation procedure to obtain the distribution function at grid point \(1/\pi_{\rm a,\,1}\) is as follows: First, the physical quantities on the virtual boundary (point w) were obtained by imposing both the no-slip boundary condition (\(\mathbf{u}_{\rm w}=\mathbf{u}_{\rm p}\)) and pressure boundary condition on the particle surface (\(\partial p/\partial\mathbf{n}=0\)). Pressure \(p_{\rm P}\) and \(p_{\rm Q}\) at point P and Q, which were located in the normal direction to the wall surface from the boundary point, were interpolated by using data neighboring four points. Pressure at grid point w \(p_{\rm w}\) was obtained by extrapolating pressure \(p_{\rm P}\) and \(p_{\rm Q}\) as  
\[p_{\rm w}=\frac{h_{\rm 2}{}^{2}p_{\rm P}-h_{\rm 1}{}^{2}p_{\rm Q}}{h_{\rm 2}{}^{2} -h_{\rm 1}^{2}}, \tag{26}\]where distances of \(h_{1}\) and \(h_{2}\) in this study were \(\sqrt{2}\), and \(2\sqrt{2}\), respectively. The equilibrium distribution function \(f^{\text{eq}}_{\alpha,\,\text{w}}\) on the virtual boundary was then obtained, and the distribution function \(f_{\alpha,\,\text{w}}\) was  
\[f_{\alpha,\,\text{w}}=f^{\text{eq}}_{\alpha,\,\text{w}}(u_{\text{w}},\,\,p_{ \text{w}})+(f_{\alpha,1}-f^{\text{eq}}_{\alpha,1}). \tag{27}\]  
Then, the virtual distribution function \(f^{*}_{\alpha,\,3}\) and the equilibrium distribution function \(f^{\text{eq}*}_{\alpha,\,3}\) at grid point 3 were obtained by extrapolation using these values at grid point 1 as follows:  
\[f^{*}_{\alpha,\,3}=\frac{a+b}{a}f_{\alpha,\,\text{w}}-\frac{b}{a}f_{\alpha,\,1}, \tag{28}\]  
\[f^{\text{eq}*}_{\alpha,\,3}=\frac{a+b}{a}f^{\text{eq}}_{\alpha,\,\text{w}}- \frac{b}{a}f^{\text{eq}}_{\alpha,\,1}. \tag{29}\]  
Finally, the distribution function at grid point 1 for the next time step considering the virtual boundary using \(f^{*}_{\alpha,\,3}\) and \(f^{\text{eq}*}_{\alpha,\,3}\) can be obtained as  
\[f_{\alpha,1}=f^{*}_{\alpha,\,3}+\frac{1}{\tau}\Big{(}f^{\text{eq}*}_{\alpha, \,3}-f^{*}_{\alpha,\,3}\Big{)}. \tag{30}\]