## 2 Methodology  
### Governing Equations for Inviscid Compressible Flows  
This study focuses on the simulation of inviscid compressible flows governed by the Euler equations. It can be expressed as follows based on its conservation law:  
\[\frac{\partial\mathbf{W}}{\partial t}+\frac{\partial\mathbf{F}}{\partial x}+\frac{ \partial\mathbf{G}}{\partial y}=0, \tag{1}\]\[\mathbf{W}=\begin{bmatrix}\rho\\ \rho u\\ \rho v\\ \rho E\end{bmatrix},\,\mathbf{F}=\begin{bmatrix}\rho u\\ \rho u^{2}+p\\ \rho uv\\ \rho uH\end{bmatrix},\,\mathbf{G}=\begin{bmatrix}\rho v\\ \rho vu\\ \rho v^{2}+p\\ \rho vH\end{bmatrix}. \tag{2}\]  
where \(\rho\) and \(p\) are the density and pressure of the flow field, \(u\) and \(v\) are the velocity components. \(E\) and \(H\) represent the total energy and enthalpy, expressed as follows:  
\[E=\tfrac{1}{\gamma-1}\tfrac{p}{\rho}+\tfrac{1}{2}\big{(}u^{2}+v^{2}\big{)}, \tag{3}\] \[H=\tfrac{p}{\rho}+E.\]  
where the \(\gamma\) is the specific heat ratio, taken as \(\gamma=1.4\) in air.  
The finite difference method is used to discretize Equation (1). The computational area is uniformly discretized to create a regular grid. \(\Delta x\) and \(\Delta y\) denote the size of the space. The finite difference scheme directly approximates spatial derivatives, constructing the conservative semi-discretization scheme.  
\[\frac{d\mathbf{W}_{ij}(t)}{dt}=-\frac{1}{\Delta x}\Big{(}\hat{\mathbf{F}}_{i+1/2,j}- \hat{\mathbf{F}}_{i-1/2,j}\Big{)}-\frac{1}{\Delta y}\Big{(}\hat{\mathbf{G}}_{i,j+1/2} -\hat{\mathbf{G}}_{i,j-1/2}\Big{)}, \tag{4}\]  
where \(\hat{\mathbf{F}}_{i+1/2,j}\) and \(\hat{\mathbf{G}}_{i,j+1/2}\) are the numerical fluxes at the \(x_{i+1/2,j}\) and \(y_{i,j+1/2}\) zone boundaries, respectively, in the following form:  
\[\begin{array}{l}\hat{\mathbf{F}}_{i+1/2,j}=\mathbf{F}\big{(}\mathbf{W}_{i-r,j},\cdots, \mathbf{W}_{i+s,j}\big{)},\\ \hat{\mathbf{G}}_{i,j+1/2}=\hat{\mathbf{G}}\big{(}\mathbf{W}_{i,j-r},\cdots,\mathbf{W}_{i,j+s} \big{)},\end{array} \tag{5}\]  
The equation satisfies two conditions: the numerical flux function matches the physical flux function, and each component is Lipschitz continuous.  
\[\begin{array}{l}\hat{\mathbf{F}}(W,\cdots,W)=\mathbf{F}(W),\\ \hat{\mathbf{G}}(W,\cdots,W)=\mathbf{G}(W).\end{array} \tag{6}\]  
The temporal discretization is achieved by transforming the system in Equation (4) into the following form  
\[\frac{d\mathbf{W}(t)}{dt}=\chi(\mathbf{W}), \tag{7}\]  
Then, we employ the TVD Range-Kutta methods developed by [PERSON] and [PERSON] (2018) to integrate Equation (7) forward in time. As shown in the following equation:  
\[\begin{array}{l}W^{(1)}=W^{n}+\Delta h\chi(W^{n}),\\ W^{(2)}=\tfrac{3}{4}W^{n}+\tfrac{1}{4}W^{(1)}+\tfrac{1}{4}\Delta h\chi(W^{(1) }),\\ W^{n+1}=\tfrac{1}{3}W^{n}+\tfrac{2}{3}W^{(2)}+\tfrac{2}{3}\Delta h\chi(W^{(2) }).\end{array} \tag{8}\]  
where \(\chi\) is the spatial discretization operator, \(\Delta h\) is the time step.  
The numerical flux in Equation (4) can be obtained by using Riemann solvers such as LF, ROE, and HLL (2018, 2019, 2020), as mentioned in the introduction. Alternatively, it may also be assessed using mesoscopic approaches like LBFS and RLBFS (2018, 2019). To achieve a high-resolution flow field and better capture shocks and discontinuities, a numerical algorithm with high accuracy, low dissipation, and greater stability should be constructed. These objectives can be achieved by combining the WENO reconstruction scheme, FDM, and URLBFS. Next, we review the RLBFS scheme based on the D1Q4 model, make improvements, propose the URLBFS numerical scheme, and then introduce the higher-order WENO5 and WENO7 schemes.