Compressive properties of Min-mod-type limiters in modelling shockwave-containing flows

The long-ignored compressive properties of Min-mod type limiter is investigated in this manuscript by demonstrating its potential in numerically modelling shockwave-containing flows, especially in SWBLI (Shock Wave/Boundary Layer Interaction) problems. Theoretical studies were firstly performed based on Sweby’s TVD ( Total Variation Diminishing ) limiter region and Spekreijse’s monotonicity-preserving limiter region to indicate Min-mod type limiters’ compressive properties. The influence of limiters on the solution accuracy was evaluated using a hybrid-order analysis method based on the grid-independent study in three typical shockwave-containing flows. The conclusions are that, Min-mod type limiter can be utilized as a dissipative and/or compressive limiter, but depending on the reasonable value of the compression parameter. The compressive Min-mod limiter tends to be more attractive in modelling shockwave-containing flows as compared to other commonly preferred limiters because of its stable computational process and its high-resolution predictions. However, the compressive Min-mod limiter may suffer from its slightly poor convergence, as that observed in other commonly-accepted smooth limiters in modelling SWBLI problems.


Introduction
The shockwave/boundary layer interactions (SWBLI) have been greatly attractive in engineering and scientific research for more than fifty years [1][2][3][4][5]. To investigate the SWBLI problems [6][7][8], the high-order (at least third-order) shock-capturing schemes [9][10][11][12] are the general choice. However, the high-order schemes are less robust, difficult to code and practice, especially in solving SWBLI problems which contain numerous steep gradients in the form of shockwaves and relatively thin thickness of boundary layers. For industrial applications, there are very few routinely utilised working CFD codes embedding higher than second-order accuracy. The classical second-order upwind schemes [13][14][15][16] are always preferable in practical aerodynamics because of their simplicity, efficiency, flexibility and robustness in practice.
The main approaches of conducting second-order schemes for strong shockwave-containing flows are MUSCL and non-MUSCL approach [17][18][19][20]. Although the non-MUSCL-based schemes have shown to provide excellent accuracy, they suffer from unsound convergence properties due to the inherent non-differentiable conducting process. As compared to the non-MUSCL-based schemes, the MUSCL-based approach is more prevalent in engineering due to its inborn simplicity and efficiency in modelling complex flows. To achieve the second order accuracy, the MUSCL method [21] proposed by Van Leer is commonly adopted. The linear interpolation method used in conducting a second-order scheme can cause spurious oscillations in steep gradient regions [22]. Therefore, the MUSCL-based linear interpolation needs to be improved with a special manner, i.e. the limiter functions, when the flow solution contains shockwave or other discontinuities. The application of limiter function in conducting upwind schemes can suppress the numerical oscillations, achieve the monotonically converged solutions with the aimed second order accuracy. In this paper, the authors mainly focus on the classical second-order schemes, with the emphasis being laid on the influence of limiter functions on the accuracy and the convergence behaviour in numerical simulations of SWBLI problems.
Existing literature shows that extensive research has already been conducted on developing accurate and robust limiters. The most significant breakthrough in this field is Sweby's second-order TVD limiter region [23] based on the one-dimensional scalar conservation law as shown in Fig.(1a). However, Goodman and Leveque [24] have stated that TVD schemes in two-dimensional cases can only achieve first-order accuracy at most. To achieve second-order accuracy in multi-dimensional problems, Spekreijse [25] extended Sweby's TVD limiter region to monotonic limiter region as shown in Fig.(1b). Barth and Jespersen [26] derived a multi-dimensional limiter suitable for unstructured grids based on Spekreijse's monotonic limiter region. In the application of the limiter in the upwind scheme, the frequently encountered problem is that it may severely hamper the numerical solution convergence.
This phenomenon is even more pronounced for a non-differentiable limiter. Venkatakrishnan [27] devised a limiter function, in which a threshold parameter based on local cell size was added. And then, the solutions can converge to the steady state as expected, and in regions where numerical oscillations were below the selected threshold, the limiter could be effectively switched off. This modification to limiter is similar to that of van Albada et al. [28] in a different context via the problem of capturing smooth extrema without clipping. Kim and Kim [29] proposed a MLP (Multi-dimensional Limiting Process) method which shows the outstanding feature of controlling numerical oscillations in multi-space dimensions with very desirable properties in terms of accuracy, efficiency and robustness.
Yoon and Kim [30] modified the aforementioned MLP and refined it for three-dimensional applications without assuming local gradients, which made an excellent improvement on the solution accuracy, convergence, as well as the robustness for the steady/unsteady flows.  Although much effort has been invested in developing and testing different kinds of limiters, there is yet no widely accepted unambiguous conclusion on their attributes. Scott [31] first performed a systematical investigation on the limiter functions in MUSCL-based upwind schemes. Among of the investigated limiters, van Albada's limiter was considered as the most attractive one because of its less compressive properties and less limitation on time marching step, whereas Min-mod limiter was considered as the most dissipative one with the least accurate prediction of discontinuities. Scott's conclusion on the investigated limiters has been widely accepted in the current practice of CFD. However, the conclusion on Min-mod limiter is still incomplete. The authors, in this paper, argue that Min-mod type limiter has various forms and its properties are actually determined by a special compression parameter (as described later). In [31], Scott's discussions only referred to one situation of Min-mod type limiter's properties, and the others with different compression parameter were not discussed at all. After that, there has been no relevant and compelling literature to study this issue in depth and in detail. Although Min-mod limiter is considered to be the most dissipative limiter in general, it is widely accepted in CFD practice and embedded in various in-house CFD codes and commercial CFD software packages, and surprisingly demonstrating its capability of achieving acceptable numerical results with robustness. Very few researchers [32,33] are aware of this inconsistent behaviour of Min-mod type limiter of being a dissipative limiter with more accurate predictions of flow discontinuities. A thorough understanding of its underlying reasons has never been explored.
As noted by Scott [31], the limiter function in MUSCL-based upwind schemes plays important roles on the numerical solutions of the shockwave-containing flows. In Section 2 the in-house-developed CFD code is presented firstly, and the modifications by the authors to the MUSCL interpolation are presented in detail in Section 3 with the form of Spekreijse's primary modification to limiter functions. The properties of several frequently utilized limiters [34] as listed in Section 3 are discussed systematically based on Sweby's TVD limiter region and Spekreijse's monotonic limiter region. In Section 4 a hybrid-order estimator based on grid-independent theory is introduced, and it is used to perform the current investigation of limiter functions' properties. Following in Section 5, the numerical solutions from various limiters are discussed in detail. Finally, a general conclusion is drawn on the properties of commonly used limiter functions in the numerical simulation of the flows containing shockwaves, in particular on Min-mod type limiter's compressive properties which have never been referred to or adequately investigated in previous literature because of the historic overemphasis of it being the most dissipative limiter. One thing needs to be noticed is that, all the numerical test cases in this work are performed based on the calorically perfect gas model with a constant specific heat of 1.4.

Numerical Method
The two-dimensional governing equations of Navier-Stokes flows are presented in the numerical computation domain ( , )  as following: ( Previous literature showed that the computational results would be significantly influenced by the different intrinsic dissipation and dispersion properties of the spatial schemes, even with the same limiter function used [31]. A common belief is that the FVS-type scheme is robust in application, but with the loss of the accuracy in solutions due to the hefty numerical dissipation. To achieve the sufficient accuracy in solutions with FVS scheme, one has to use a more refined computational grid to decrease the computational dissipation caused by the insufficient accuracy of spatial discretisation. Roe's FDS scheme is another popular choice for solving flow governing equations with high resolution and high fidelity, but there is a possibility of violating the entropy condition, which may lead to significant difficulties in simulating high-speed complex flows. In Ref. [36] AUSMpw+ was proposed as an improved AUSM-type of scheme [37], which has the advantages of having not only the robustness of FVS schemes but also the high-resolution and high-accuracy of FDS schemes for CFD practical problems. Therefore, AUSMpw+ scheme is chosen to perform the aimed limiter function investigation in the current research.

Review of Limiter Functions
The original van Leer's MUSCL interpolation [21] is written as following: (2) Where L q and R q are primitive flow variables at left and right sides of the grid cell faces, and  is a parameter that determines the spatial accuracy. As  increases from -1 to 1/3, the solution accuracy would increase and 0  = and 13  = are the commonly adopted options in practice.
It has been proven that the linear interpolation formulas used to obtain second-order accurate algorithm, like Eq.(2), are not adequate since they can cause unphysical oscillations in solutions where discontinuities exist.
According to Spekreijse [25], Eq. (2) can be reformulated as: Where () r  is expressed as following： (4) and In computing numerical fluxes, limiter function will result in a reduction of accuracy in discontinuous regions. The first form of Min-mod limiter discussed by Sweby [19] is expressed as: Scott [31] and Anderson [38] employed different forms of Min-mod limiter in their work respectively, of which the characteristics are determined by a compression parameter. To make a completely transparent comparison between selected limiter functions, Min-mod limiter employed by Scott [31] and Anderson [38] is reformulated to have a similar form as that of Eq. (6).
The MUSCL interpolation formula used in [31,37] is presented below: where and  is the compression parameter given by In general, the maximum allowable value of  is adopted.  Eq.(8) will produce several forms of Min-mod-type limiter with different  and  : It is obvious that Min-mod type limiter shown in Eq.(10d) really satisfies monotonicity-preserving condition by Therefore, it would be an incomplete conclusion to take Min-mod limiter as the most dissipative limiter without investigating its compressive features. The compressive behaviours of Min-mod type limiter in practical applications will be further explored through three shockwave-containing flows in the following sections. One crucial aspect to be noted for Superbee limiter is that, it tends to turn smooth waves into square waves, and makes the gradient sharper in solving practical flow. The overly compressive nature of Superbee limiter in multiple dimensions may lead to stair-casing effects at flow discontinuities. Due to its impracticality in engineering, Superbee limiter will not be discussed in the following sections. Min-mod limiter of 2  = will also not be considered because of its Superbee-like features in the region of 1 r  as shown in Fig.(2b).

Grid Convergence Error Analysis Method
There has been consensus on the effect of limiters on the convergence properties and accuracy of numerical solutions. However, most of the existing literature is based on the conclusion of qualitative analysis, and it is difficult to provide quantitative analysis results for the characteristics of different limiters and their comparison due to the lack of effective error analysis tools. For SWBLI problem, the general numerical solutions based on upwind schemes can only achieve up to third-order accuracy due to the existence of shockwave discontinuity in the flow field, and even only the first order can be reached at the discontinuity. Roach [39,40] proposed a global numerical error analysis method based on Richardson extrapolation theory, but the effectivity on the discontinuity-containing flow field error analysis was less than expected. Roy et al. [41][42][43] demonstrated that, the first-order and second-order errors coexist in the numerical solutions of discontinuity-containing flows, and some flow quantities are non-monotonically converged with the grid refinement. The occurrence of non-monotonic grid convergence solution is due to the opposite sign of the first-and second-order errors. Therefore, Roy [41] proposed a hybrid-order error analysis method for this non-monotonic grid convergence phenomenon in discontinuity-containing flows. In this research, the same hybrid-order error analysis method by Roy is adopted to evaluate the performance of selected limiter functions.
The hybrid-order analysis method by Roy [41] is expressed as the following Eq.(12). (12) Where k f is the flow quantity on the k th level grid, and the densest grid is indicated as as the exact solution of flow quantities. The grid size scaling factor of 1 k + th grid to k th grid is defined as The detailed information of other symbols, such as By setting the grid scaling factor shown in Eq.(13) as a constant, the approximate solution of 1 g , 2 g and exact f can be obtained from Eq. (12), which is expressed as following: Where , The approximate solution of exact f shown in Eq. (16) is generally third-order accurate, and the spatial discretization error on the k th grid with respect to exact f can be expressed as following: Based on Roy's hybrid-order error analysis method, the first-and second-order errors can be given as follows, and the sum of them is presented at the same time in the following formulation:

NACA 0012 Airfoil
Anderson [38] studied the properties of Min-mod limiter of 4  = and the differentiable Albada-like limiter in two-dimensional transonic Euler flow around NACA 0012 airfoil. With the same 2D case in present work, a series of computational grids are generated to perform the grid-independence investigation. The free stream flow parameters are given as Although there is special interest in SWBLI problems in present work, it is a common way to verify the characteristics of a spatial scheme by performing Euler-based numerical simulations, and the same idea is adopted in verifying the limiter functions' properties in this test. Table 1 provides the detailed grid descriptions of the grid dimensions and the grid spacing. Figure 3 shows the drag coefficients calculated by different limiters on all grids, and the estimated grid-converged drag coefficients on 0 h = grid. It can also be concluded that, with grid refinement, the computed drag coefficients are grid-converged, and the grid-converged drag coefficients obtained by Eq.(16) for different limiter function are slightly different. The varying of grid-converged flow quantities indicates that, the numerical prediction accuracy is related to the intrinsic properties of the limiter function embedded in the numerical scheme.    Table 2, which indicates that the spatial discretization error from Min-mod limiter of 4  = is the least, while 1  = is the largest.    Another important aspect that may affect the computational solution accuracy is the iterative convergence error in numerical simulations. Venkatakrishnan [27] discussed the influences of limiters on converging to steady solutions of numerical simulation in great detail, pointed out that a non-differentiable limiter might more severely hamper the convergence process of a numerical simulation and result in less accurate solutions than that of a differentiable limiter. Figure 7 shows the residual convergence process with respect to iteration steps for all abovementioned limiters on the densest grid, and clearly demonstrates that the three differentiable limiters, i.e. van Leer limiter, Hemker-Koren limiter and van Albada limiter, exhibit excellent convergence performance by approaching to the machine zero level. However, as depicted in [38], the residuals of non-differentiable Min-mod type limiters fail to approach machine zero and exhibit numerical oscillations after a few orders of magnitude drops. Anderson et al. [38] also pointed out that the limit cycle oscillation in the iterative process for Min-mod limiter of 4  = was mainly due to its non-differentiable feature. Also, Ventakarishnan [27] concluded that the convergence behaviour is even worse in the case of non-differentiable limiter functions compared with that of smoothing and differentiable limiter functions.
The discussion mentioned above clearly demonstrates that conclusions from Fig.7 are consistent with those in [27] and [38], i.e. the non-differentiable limiter might severely deteriorate the convergence process of numerical simulations. The contrast, that the compressive Min-mod limiters with 3  = and 4  = produce larger spatial errors than that of the less compressive Hemker-Koren limiter, as shown in Fig.4, mainly comes from the poorly converged iterative residual due to the non-differentiable feature of Min-mod type limiter. The current study also demonstrates that, with a reasonable compression parameter  , Min-mod limiter could behave as a compressive limiter predicting the solution accurately with small spatial errors. Meanwhile, it might also become a dissipative limiter resulting in significant spatial errors, which correlate with the theoretical analysis shown in Fig.2.

Supersonic SWBLI on Flat Plate
Hakkinen et al. [44] conducted an experiment on the laminar boundary layer flow interacting with an incident shock on a flat plate, as illustrated in Fig.8. This experiment has been frequently used as a benchmark to verify various numerical algorithms. Here this same case is selected to perform the numerical evaluation of limiter functions.
The   Table 3 describes the detailed information of the computtional grids for the current study. Figure 9 presents the calculated drag coefficients by different limiter function on all grids, and the grid-converged drag values by Eq. (16) are also plotted. For different limiter function, the extrapolated grid-converged drag coefficients are slightly different, and the maximum difference between different limiters is about 1 count. The highest grid-converged drag is from Min-mod limiter of 4  = , while the lowest is from Hemker-Koren limiter.     Table 4.
(a) Skin friction coefficient

Fig. 14 Comparison of computed results
Detailed information about the separation regions for different limiters is shown in Table 4 and illustrated in Fig.(14a). For the distribution of surface pressure shown in Fig.(14b), the primary discrepancy of different limiters is at the impacting point, where there is an obvious pressure plateau, which is determined by the size of the separation bubble. The larger the separate region is, the flatter the plateau of the pressure distribution is. Figure 15 shows the pictures of SWBLI separation bubble predicted by different limiters and the streamline spectrum in the vicinity. It can be seen that, the stronger the dissipation of the limiter, the smaller the predicted separation bubble size. The least dissipative Min-mod limiter of 3, 4  = predicts the largest separation bubble size, which is consistent with the separation bubble information given in Table 4. It is also stated that Min-mod limiter of 3, 4  = is a compressive type limiter with less dissipation.

Hypersonic Flow about a 24-deg Compression Ramp
As presented schematically in Fig.16, the current test is a complex hypersonic case with the strong interactions of shockwaves, expansion waves, and boundary layer as well. As shown in Fig.16 Table 5 presents the detailed information of the computational grids. Figure 17 shows the comparison of the predicted drag coefficients on all grids with different limiters, which indicates the computations are grid-independent. The estimated grid-independent drag coefficients with hybrid-order extrapolation method for all limiters are also plotted in Fig.17. It can be concluded that, as the grid spacing approaches to zero, the predicted drag coefficients tend to converge to the grid-independent solutions for each limiter respectively, but the grid-independent values are different for different limiter. The variation of the extrapolated grid-converged solutions shown in the current test case, as well as in the aforementioned two test cases, indicates that, the grid-converging process of different limiter functions is determined by their inherent dissipation properties.  As noted by Dolling [1], Knight [2,4], Zheltovodovo [3,5] and Rudy [46], the flow field of hypersonic SWBLI over a ramp with large ramp angle is dominated greatly by the The distribution of pressure coefficients on the ramp plate is shown in Fig.(20a)  = obtains the smallest pressure peak value at the most upstream peak location. Therefore, the peak value decreases and peak location moves upstream as the dissipation of limiter functions, as shown in Fig.(20a).
For the skin friction coefficients shown in Fig.(20b),   Table 6 presents the detail information of the separation region predicted with different limiters, and Min-mod

Conclusions
The compressive properties of Min-mod type limiter are investigated theoretically and numerically in the current research. As a comparison, a series of commonly-used MUSCL-based limiter functions are assessed together. Three typical shockwave-containing flows, including the transonic flow about NACA 0012 airfoil and two high-speed laminar SWBLI problems, are performed. A hybrid-order spatial error estimator is introduced for the first time in public literature to perform the limiter function investigations. Some key conclusions can be drawn as follows: 1) The MUSCL-based interpolation commonly used to obtain formally second-order accurate scheme is reformulated according to Sweby's second-order TVD limiter region and Spekreijse's second-order monotonic limiter region. The properties of Min-mod type limiter are found to be determined by the compression parameter  . This notion complements the previous incomplete comments on Min-mod type limiter by demonstrating its compressive properties, not only the dissipative properties.
2) The hybrid-order spatial error estimator devised by Roy is introduced for the first time to perform the analysis and comparison of limiter functions' dissipative characteristics. Studies have shown that, the inherent dissipation of the limiter function has a significant impact on the numerical solution accuracy in shockwave-containing flows. The compressive limiter, like Min-mod limiter of 4  = , can predict the flow discontinuities more accurately than a dissipative limiters does.
3) It is very hard to obtain the well-converged solution in strong shockwave-containing flows, such as the selected supersonic/hypersonic SWBLI problems, no matter the limiter is differentiable or non-differentiable, and the differentiable/non-differentiable properties of limiters are no longer the indication of the iterative convergence to machine zero for strong shockwave-containing flows.
4) This work demonstrates that not only theoretically but also numerically Min-mod type limiter changes from the dissipative limiter to the compressive limiter by selecting the reasonable compression parameter. The compressive properties of Min-mod type limiter have been substantiated in current research through simulating strong shockwave-containing flows, resulting in high resolution of flow structures. This is significant for the further development of new numerical schemes via Min-mod-like procedure to simulate strong discontinuity-containing problems more accurately.