## 1 Introduction  
With the development of science and technology, scholars are paying more and more attention to the flow mechanism. Numerical simulation can shorten the research cycle and save costs, which is one of the important ways to conduct flow mechanism research. For simulation methods such as large eddy simulation (LES) and direct numerical simulation (DNS), a large number of grids and longer computation time are often required to obtain very detailed calculation results, which goes against the advantages of numerical simulation. Therefore, high-order schemes have been studied and combined with a series of numerical simulation methods, enabling numerical simulations to capture more detailed physical details with fewer grids [1; 2]. In the process of solving the Navier-Stokes equations, scholars have observed that the solutions of the nonlinear equations exhibit discontinuities, regardless of the smoothness of the initial conditions. To address this issue, researchers have attempted to improve Weighted Essentially Non-Oscillatory (WENO) schemes, and it has been found that this approach can make numerical solving methods more efficient and accurate [3]. An essentially non-oscillatory (ENO) scheme was proposed in 1987, which achieved high-order accuracy and avoided oscillations in the solution. However, it was soon found to be unstable under certain unfavorable initial conditions [4]. Therefore, in the 1990s, the ENO scheme was improved and developed into the WENO scheme by [PERSON] et al. [5]. Subsequently, [PERSON] and [PERSON] [6] further improved the WENO scheme and proposed a general framework for the design of smoothness indicators and nonlinear weights by averaging \(k\) stencils to achieve a \(2k-1\) order of accuracy in smooth regions while maintaining a \(k\) order of accuracy at discontinuities. This led to the classical WENO schemes [7; 8; 9]. Although the classical WENO schemes have many advantages, they exhibit several drawbacks in engineering applications [10]:  
1. The computational cost is very high, and the calculation process is complex. 2. The optimal (linear) weights depend on the geometry of the mesh and may become negative in some cases, meaning they lack robustness. 3. The drawbacks become more pronounced with an increase in spatial dimension. To enhance the engineering application value of the WENO scheme, [PERSON] and [PERSON] [11] developed a new fifth-order finite-difference WENO-ZQ scheme in 2016, which was later extended to a finite-volume version in multiple dimensions [12]. This new class of WENO schemes uses a convex combination of a quartic polynomial with two linear polynomials on unequal-sized spatial stencils in one dimension and is extended to two dimensions in a dimension-by-dimension fashion. It is more efficient, simpler, and easily extendable to multiple dimensions. The associated linear weights are artificially set to be any random positive numbers, with the only requirement being that their sum equals one. This property ensures that the optimal (linear) weights do not depend on the geometry of the mesh, and the problem of negative linear weights observed in classical WENO-JS schemes [8] is mitigated. Additionally, in the context of solving steady-state problems using classical WENO schemes [6; 7], it is common for the residual to plateau at a level above machine zero, equivalent to the truncation error. This occurs even when there is minimal observable change in the physical variables over successive time iterations [13]. Numerous scholars have explored various methods to mitigate this issue. However, research indicates [13] that WENO-ZQ can converge to machine zero without requiring any modifications across a range of standard test cases. These encompass scenarios with strong shocks, contact discontinuities, rarefaction waves, and their interactions, and intricate wave dynamics interacting with computational boundaries. In recent years, numerous scholars have improved WENO-ZQ and applied it in various engineering applications.