## 2 Methodology
### Evaluating Numerical Flux for Compressible Flow
#### 2.2.1 The RLBFS Scheme  
The key to the RLBFS scheme is to solve the numerical flux in different directions through LBFS based on the D1Q4 model. The model proposed by [PERSON] et al. [39] overcomes the problem of LBFS scheme performance degradation caused by too many user-specified parameters in other 1D models [35; 40]. This model mainly consists of equilibrium distribution functions \(g_{m}(m=1,2,3,4)\) and respective velocities \(d_{m}(m=1,2)\), as shown in Figure 1. Using the conservation form of moments, they can be given by:  
\[ \begin{split}& g_{1}=\frac{\rho(-d_{1}d_{2}^{2}-d_{2}^{2}u+d_{1}u^{2}+ d_{1}\varrho^{2}+u^{3}+3u\theta^{2})}{2d_{1}(d_{1}^{2}-d_{2}^{2})},\\ & g_{2}=\frac{\rho(-d_{1}d_{2}^{2}+d_{2}^{2}u+d_{1}u^{2}+d_{1} \varrho^{2}-u^{3}-3u\theta^{2})}{2d_{1}(d_{1}^{2}-d_{2}^{2})},\\ & g_{3}=\frac{\rho(d_{1}^{2}d_{2}+d_{1}^{2}u-d_{2}u^{2}-d_{2} \varrho^{2}-u^{3}-3u\theta^{2})}{2d_{2}(d_{1}^{2}-d_{2}^{2})},\\ & g_{4}=\frac{\rho(d_{1}^{2}d_{2}^{2}-d_{1}^{2}u-d_{2}u^{2}-d_{2} \varrho^{2}+u^{3}+3u\theta^{2})}{2d_{2}(d_{1}^{2}-d_{2}^{2})},\\ & d_{1}=\sqrt{u^{2}+3\theta^{2}-\sqrt{4u^{2}\theta^{2}+6\theta^ {4}}},\\ & d_{2}=\sqrt{u^{2}+3\theta^{2}+\sqrt{4u^{2}\theta^{2}+6\theta^ {4}}}.\end{split}  \tag{9}\]  
where \(\theta=\sqrt{Dp/\rho}\) is the peculiar velocity of particles, if it is a one-dimensional model, then \(D\) takes the value of \(1\). By solving the \(g_{m}(m=1,2,3,4)\) and \(d_{m}(m=1,2)\), the Euler equation can be restored using the physical conservation laws [35; 37]:  
\[ \begin{split}&\rho=\sum\limits_{i=1}^{4}g_{i},\\ &\rho u=\sum\limits_{i=1}^{4}g_{i}\xi_{i},\\ &\rho uu+p=\sum\limits_{i=1}^{4}g_{i}\xi_{i}\xi_{i},\\ &\rho E=\sum\limits_{i=1}^{4}g_{i}(\frac{1}{2}\xi_{i}\xi_{i}+e_{ p}),\\ &(\rho E+p)=\sum\limits_{i=1}^{4}g_{i}\bigg{(}\frac{1}{2}\xi_{i} \xi_{i}+e_{p}\bigg{)}\xi_{i}.\end{split}  \tag{10}\]  
where \(\xi_{i}\) represents the particle velocity in i-direction, \(\xi_{1}=d_{1},\xi_{2}=-d_{1},\xi_{3}=d_{2},\xi_{4}=-d_{2}\), \(e_{p}=[1-\frac{D}{2}(\gamma-1)]e\) and \(e=p/[(\gamma-1)\rho]\). It is important to note that when addressing multidimensional problems, the previously mentioned model is applied solely in the normal direction of the interface [31; 32], for instance, in the two-dimensional scenario depicted in Figure 2.  
Figure 1: Distribution of discrete lattice velocities for 1D model. Consider \((n_{x},n_{y})\) to represent the normal vector at the interface, we can derive a relationship as follows:  
\[ \begin{split}& U_{n}=n_{x}u+n_{y}v,\ U_{\rm{T}}=n_{x}v-n_{y}u\,\\ & u=n_{x}U_{n}-n_{y}U_{\rm{T}},\ v=n_{x}U_{\rm{T}}+n_{y}U_{\rm{H}}, \end{split}  \tag{11}\]  
Then, Equation (2) can be represented as:  
\[ \begin{split}\mathbf{W}&=\left[\begin{array}{c}\rho \\ \rho\big{(}n_{x}U_{\rm{H}}-n_{y}U_{\rm{T}}\big{)}\\ \rho\big{(}n_{x}U_{\rm{T}}+n_{y}U_{\rm{H}}\big{)}\\ \rho\big{(}U_{\rm{H}}^{2}/2+e\big{)}+\rho U_{\rm{T}}^{2}/2\end{array}\right],\\ &\mathbf{F}_{\rm{c}}=\left[\begin{array}{c}\rho U_{\rm{H}}\\ \big{(}\rho U_{\rm{H}}^{2}+p\big{)}n_{x}-\rho U_{\rm{H}}U_{\rm{T}}n_{y}\\ \rho U_{\rm{H}}^{2}+p\big{)}n_{y}+\rho U_{\rm{H}}U_{\rm{T}}n_{x}\\ \rho\big{(}U_{\rm{H}}^{2}/2+e\big{)}+p\big{)}U_{\rm{H}}+\rho U_{\rm{H}}U_{\rm{ T}}^{2}/2\end{array}\right].\end{split}  \tag{12}\]  
Assume the position at the interface is x = 0. The inviscid flux at the interface is given by:  
\[\mathbf{F}_{c,i+1/2}^{*}=\Big{[}\rho U_{\rm{H}}\quad\rho U_{\rm{H}}U_{\rm{H}}+ p\quad\Big{(}\rho\Big{(}\tfrac{1}{2}U_{\rm{H}}U_{\rm{H}}+e\Big{)}+p\Big{)}U_{ \rm{H}}\Big{]}^{T}=\sum_{i=1}^{4}\xi_{i}\mathbf{\varphi}_{\alpha}f_{i}(0,t), \tag{13}\]  
where \(\mathbf{\varphi}_{\alpha}=\Big{(}1,\xi_{i},\tfrac{1}{2}\xi_{i}^{2}+e_{p}\Big{)}^{ T}\), \(f_{i}(0,t)\) is the distribution function of equilibrium. And superscript * represents the numerical solution.