Grid Convergence Analysis for MUSCL-based Numerical Scheme in Shockwave-containing Flows

. This paper investigated the influence of limiter functions widely utilized in MUSCL-type (Monotone Upstream-centred Schemes for Conservation Laws) upwind numerical schemes on the solution accuracy of shockwave-containing flows. An incident shock interacting with laminar boundary layer developed on a flat plate was numerically simulated with the in-house developed code. A mixed-order grid convergence study was performed to assess the spatial errors of different limiters in simulating the selected shockwave-containing flow on flat plate. The conclusions are that, limiter functions implemented in the current in-house code play the critical roles in accurately predicting shockwave-containing flows. The mixed-order error estimator based on grid convergence study was proved to be applicable to evaluate the spatial errors of shockwave-containing flows, where the shock could reduce the nominal second- or third-order accuracy to first-order. The mixed-order estimator is conservative in the sense that the actual error is less than the error estimated, in the examined case.


Introduction
As the advances in computing power, there are great demands for the high confidence in modeling and simulating engineering problems. It has become a trend to employ high-order shock capturing schemes with at least a third-order spatial accuracy to perform investigations for shockwave-containing flows [1][2][3]. However, the traditional second-order upwind schemes are more preferable in practice due to their simplicity, flexibility and robustness. For strong discontinuitycontaining flows there are two different strategies to conduct second-order accurate schemes: MUSCL type and non-MUSCL type. Van Leer's MUSCL approach [4] is the widely accepted method to achieve a second-order upwind scheme, and the MUSCL interpolation process can be enhanced with the so-called limiter functions or limiters to prevent the overshoot or undershoot in numerical solutions if discontinuous solutions existing in the flow field. The purpose of utilising a limiter function is to suppress the non-physical oscillations and maintain monotonicity while retaining the expected second-order accuracy in numerical solutions.
For shockwave containing flows, the most common numerical methods are characteristic-based upwind methods that can achieve nominally third-order accuracy at most. However, the solution accuracy can reduce to first order at discontinuities due to the presence of shock waves. Previous work by Roy et al. [5][6][7] verified that, the use of the nominal second-order scheme in discontinuity-containing flows can result in the presence of both firstand second-order error in solutions. The non-monotonic grid convergence behaviour was found to occur when the first-and second-order error terms were of opposite signs, thus leading to error cancellation. When the solution variables behave in a non-monotonic fashion, a mixedorder error estimator by Roy [5] could be applied to assess the spatial errors of mixed-order solutions. This mixedorder error estimator was proven to provide conservative error estimate in the sense that the actual error is less than the estimated mixed-order error. Therefore, Roy's mixedorder error estimator is employed here to verify the solution accuracy of the selected limiter functions.
The limiter functions used in conducting second-order upwind schemes have great influence on the numerical resolutions of expansion fans, contact discontinuities and shock waves in flow fields. A brief description of the numerical method and the in-house code used in current study are provided In Section 2, and the formulation of MUSCL interpolation with limiter functions are also presented in details. A mixed-order error analysis method is introduced in Section 3 to prepare for grid-independent study and to investigate the influence of limiters on the accuracy of the solution of shockwave-containing flows. The influences of various limiters on accuracy, and grid convergence behaviours of numerical are demonstrated in Section 4. Finally, meaningful conclusions on the properties of the widely accepted limiters in shockwavecontaining flows are drawn.

Numerical Method
The two-dimensional compressible Navier-Stokes equations with thin shear layer approximation in the generalized coordinate ( , )   plane are given in conservation form as follows: Where is conservative vector, and are convective flux vectors, and are viscous flux vectors. An in-house CFD code named ATACF (Analysis Toolkits for All speed Compressible Flows), which is the updated version of ATTF (Analysis Toolkits for Transonic Flows) in [8], developed by Li, is applied to perform the simulation in the current research. In ATACF, the convective terms have several options to be discretized: central or upwind-biased methods, including central-difference scheme, flux-vector splitting scheme (FVS), flux-difference splitting schemes (FDS) and the series of AUSM (Advection Upstream Splitting Method) schemes. As compared with FVS and FDS schemes, the AUSMpw+ scheme [9] tends to be more attractive as an improvement of the AUSM scheme [10]. Therefore, the AUSMpw+ scheme is employed to perform the simulations to verify the properties of the relevant limiter functions in the present paper.

Limiter Functions
A brief review of the original MUSCL interpolation equations proposed by van Leer [4] is presented as follows:     1 2 To suppress these unphysical oscillations in solutions, monotonicity-preserving principles are invoked in MUSCL interpolation stage by introducing limiter functions. Following the treatment of Spekreijse [12], the interpolation formula shown in Eq.(2) is rewritten as:   , Eq. (5d) becomes as following:

Grid Convergence Analysis Method
The mixed-order extrapolation method by Roy [5] is presented as in Eq. (7), where three discrete solutions on three different grids from fine to coarse are needed: Where 32 In general, the above estimation of the solution exact f is of third-order accuracy, and the estimated spatial error based on Eq.(11) on each grid level is calculated using the following relationship: As stated by Roy [5], if multiplied by a safety factor, for example F s = 3.0, Eq.(12) will produce the similar error estimate with that of the GCI (Grid Convergence Index) method proposed by Roache [13,14] when the generalized Richardson Extrapolation method is employed to determine exact f in Eq. (12) and when the estimated GCI error is less than 10-20%. Thus Eq.(12) is employed in current research to evaluate the spatial error of the grid convergence study.
The normalized magnitudes of the first-and secondorder error terms, as well as their sum, are formulated as:

Results and Discussions
An incident shock interacting with laminar boundary layer developed on a flat plate, as shown in Fig.1 is studied. This case study corresponds to the experiments by Hakkinen et. al. in [15], which is often used to validate the efficiency and accuracy of numerical methods. Three computational grids are generated in a rectangle domain with length 2 starting from the leading edge of the plate and height starting from the surface of the plate. The reference length measured from the leading edge of the plate to the impingement point on the plate is . The detailed grid information used for this computation is shown in Table.   . Table.1 shows the detailed grid information used in current simulations including grid dimensions used.   Fig.2 shows the predicted drag coefficients, as well as the extrapolated grid-independent drag coefficients with different limiters. It is evident that, as the grid is successively refined, the numerically predicted drag coefficients for each limiter converged to their gridindependent values respectively. Fig.2 also shows that the Min-mod limiter with 4   predicted the highest gridindependent drag coefficient, while the Hemker-Koren limiter predicts the lowest drag coefficient.   show the nonmonotonic convergence behaviour due to the cancellation of first-and second-order terms with opposite signs. However, the grid-convergence errors for all limiters approach to formally second-order range asymptotically on coarser grids, while to observed first-order range on finer grids. The spatial errors of the solution with the employment of limiter functions in the numerical scheme are of the mixed-order accuracy in current simulation. Fig.4 shows the spatial error levels with different limiters.
It is evident that the Min-mod limiter with 1   predicts the highest spatial error, while 4   possesses the lowest spatial error.  The separation bubbles due to the interactions of the shock wave with the boundary layer predicted by different limiters are shown in Fig.5, where the streamlines around the separation bubble are also plotted. The changes of the boundary layer along the entire flat plate can also be seen from Fig.5, especially in the region of separation bubble. As the dissipative characteristic of limiter gets stronger, the predicted size of separation bubble becomes smaller. The most dissipative Min-mod limiter with 1   predicts the smallest separate region.
The sizes of separation region predicted using the Minmod limiters of 3, 4   are larger than those predicted by implementing van Leer limiter, van Albada limiter and Hemker-Koren limiter.
From the above discussions, it can be concluded that, for discontinuity-containing flows, the spatial error, i.e. the solution accuracy, is mainly determined by the dissipation level of the limiters adopted: the less dissipation for the limiter, the more accurate the solution.

Conclusion
The properties of several frequently used MUSCL-type limiters, including van Leer limiter, van Albada limiter, and Hemker-Koren limiter functions and Min-mod limiter are assessed together. Grid independent studies based on the mixed-order spatial error estimator are conducted to investigate the influence of limiter's dissipation/compression level on the solution accuracy for the first time in available public publications. Based on the mixed-order error estimating method, the gridindependent solutions are of first-order accuracy in shockwave-containing flows and of second-order accuracy on the coarser grid. The spatial errors are of different levels with different limiters which possess different intrinsic dissipative/compressive properties. These intrinsic dissipative/compressive features play an important role in the solution accuracy in predicting strong shockwave-containing flows.