## 2 Numerical Methods
### Governing Equation for Fluid Part  
The regularized lattice Boltzmann method (RLBM) [18] was used as a governing equation for the fluid part. This method is based on the lattice Boltzmann method (LBM). RLBM maintains simplicity, efficiency, and good scalability for parallel computing of LBM; in addition, it can reduce sufficient memory usage. The distribution function \(f_{\alpha}\) for the lattice Boltzmann equation can be written as  
\[f_{\alpha}=w_{\alpha}\Big{(}a_{0}+b_{i}e_{ai}+c_{ij}e_{ai}e_{aj}+d_{ijk}e_{ai}e _{aj}e_{ak}+\cdots\Big{)}. \tag{2}\]  
The expansion of the distribution function up to second-order moment is equivalent to incompressible Navier-Stokes equations; therefore, the distribution function is approximated as  
\[f_{\alpha}\ \approx\ w_{\alpha}\big{(}a_{0}+b_{i}e_{ai}+c_{ij}e_{ai}e_{aj} \big{)}, \tag{3}\]  
where \(e_{\alpha}\) is the discrete velocity vector and subscript \(\alpha\) represents velocity direction. The distribution functions propagate with velocities restricted to a finite set of discrete vectors, which is modeled in the D\(n\)Q\(m\) model (\(m\) denotes discrete velocities in \(n\) dimensions). The D2Q9 model [19] was applied in our simulations. This model is one of the most common lattice speed models for two-dimensional lattice Boltzmann simulations and it has the following nine discrete lattice velocity vectors:  
\[e_{\alpha}=\left\{ \begin{array}{ll}c(0,0)&(\alpha=0)\\ c\Big{(}\cos\frac{(\alpha-1)\pi}{2},\ \sin\frac{(\alpha-1)\pi}{2}\Big{)}&( \alpha=1\text{--}4)\\ \sqrt{2}c\ (\cos\frac{(2\alpha-9)\pi}{4},\ \sin\frac{(2\alpha-9)\pi}{4} \Big{)}&(\alpha=5\text{--}8)\end{array} \right. \tag{4}\]  
In this model, weight coefficients \(w_{\alpha}\) are defined as follows:  
\[w_{\alpha}=\left\{ \begin{array}{ll}4/9\ \ \ (\alpha=0)\\ 1/9\ \ \ (\alpha=1\text{--}4)\\ 1/36\ \ (\alpha=5\text{--}8)\end{array} \right., \tag{5}\]  
The parameters \(a_{0}\), \(b_{i}\), and \(c_{ij}\) are determined to satisfy the following relationships:  
\[\sum_{\alpha}f_{\alpha}=\rho, \tag{6}\]  
\[\sum_{\alpha}f_{\alpha}e_{\alpha i}=\rho u_{i}, \tag{7}\]  
\[\sum_{\alpha}e_{ai}e_{aj}\ f_{\alpha}\ =\ \frac{c^{2}}{3}\rho\delta_{ij}+\rho u_{i}u_{ j}+\Pi I_{ij}^{\text{neq}}, \tag{8}\]  
where \(\rho\) is the fluid density, \(u\) is the fluid velocity vector, \(c\) is the lattice speed defined as \(c=\delta x/\delta t\), \(\delta_{ij}\) is the Kronecker delta, and \(\Pi I_{ij}^{\text{neq}}\) is the non-equilibrium part of the stress tensor. The lattice units of \(\ \delta x=1\) and \(\delta t=1\) were used in this study. The parameters \(a_{0}\), \(b_{i}\), and \(c_{ij}\) are expressed by \(\rho\), \(u\), and \(\Pi I_{ij}^{\text{neq}}\). The distribution function \(f_{\alpha}\) can therefore be written as  
\[f_{\alpha}=w_{\alpha}\rho\left[1+\frac{3(e_{\alpha}u_{i})}{c^{2}}+ \frac{9(e_{\alpha}u_{i})^{2}}{2c^{4}}-\frac{3(u_{i}u_{j})}{2c^{2}}\right] \tag{9}\] \[\qquad\qquad+\frac{9w_{\alpha}}{2c^{2}}\left(\frac{e_{ai}e_{aj}}{ c^{2}}-\frac{1}{3}\delta_{ij}\right)\Pi I_{ij}^{\text{neq}}.\]The first term on the right-hand side is equal to Maxwell equilibrium function \(f_{\alpha}^{\text{eq}}\). The second term on the right-hand side is described as \(f_{\alpha}^{\text{1}}\), and the distribution function \(f_{\alpha}\) is expressed by  
\[f_{\alpha} =f_{\alpha}^{\text{eq}}+f_{\alpha}^{\text{1}}.\] (10) The time evolution equation for the regularized lattice Boltzmann equation is then \[f_{\alpha}\big{(}t+\delta t,\,\mathbf{x}+\mathbf{e}_{\alpha}\delta t\big{)} =f_{\alpha}^{\text{eq}}(t,\,\mathbf{x})+\left(1-\frac{1}{\tau}\right)f_{ \alpha}^{\text{1}}(t,\,\mathbf{x}),\] (11) where \(\tau\) is single relaxation time defined using kinetic viscosity of fluid \(\
u\) as follows: \[\tau =\frac{3\
u}{c\delta x}+\frac{1}{2}.\] (12) The single relaxation time \(\tau\) was set to 0.88 for all simulations. In addition, local pressure distribution function \(p_{\alpha}\) introduced by [PERSON] et al. [20] was applied to simulate incompressible flow.