A generalized adaptive central-upwind scheme for compressible flow simulations and preventing spurious vortices
Amareshwara Sainadh Chamarthi





AA_address]Amareshwara Sainadh Chamarthi [cor1]Corresponding author. 
 
E-mail address: sainath@caltech.edu (Amareshwara Sainadh  Ch.).
[AA_address]Division of Engineering and Applied Science, California Institute of Technology, Pasadena, CA, USA

§ ABSTRACT 
This work introduces a novel adaptive central-upwind scheme designed for simulating compressible flows with discontinuities in the flow field. The proposed approach offers significant improvements in computational efficiency over the central gradient-based reconstruction approach presented in <cit.> (Hoffmann, Chamarthi and Frankel,  JCP 2024). By leveraging a combination of conservative and characteristic variable reconstruction, the proposed approach demonstrates oscillation-free results while effectively reducing computational costs and improving the results. Furthermore, the adaptive central-upwind algorithm is generalized to be compatible not only with the gradient-based reconstruction as in <cit.> but with other existing methods. In this regard, with the proposed algorithm, the standard fifth/sixth-order reconstruction scheme has also been shown to outperform existing schemes with a 20-30% reduction in computational expense with improved results. Notably, the proposed approach has successfully prevented the generation of spurious vortices in the double shear layer test cases, even with linear schemes, showcasing its robustness and effectiveness.



Low dissipation, Gradient-Based Reconstruction, Monotonocity Preserving, Central scheme, Spurious vortices





§ INTRODUCTION


Adaptive central-upwind schemes are required to simulate compressible flows when both turbulence and discontinuities, like shock waves or contact discontinuities, are present. The challenge in these simulations is to capture discontinuities (where upwind schemes with dissipative properties are required) while maintaining low dissipation (an attribute of central schemes) in turbulent flow regions, which is a difficult balance to achieve. This is particularly crucial in large eddy simulations (LES) on coarse meshes, where resolving a wide range of turbulence scales with minimal dissipation is essential. Adaptive central-upwind schemes address this by combining the strengths of upwind and central schemes, allowing them to adapt to the complex demands of such simulations <cit.>.

The adaptive central-upwind schemes are minimally dissipative, leading to sharp shocks and well-resolved broadband turbulence, but rely on an appropriate shock sensor. The shock or discontinuity sensor is the main challenge for the hybrid central/upwind method, as mentioned in <cit.>. Various researchers have proposed different approaches to address this challenge. For example, Kim and Kwon suggested a hybrid scheme that combines a central scheme with the numerical dissipation of the upwind weighted essentially non-oscillatory (WENO) scheme using a weighing function <cit.>. Hu, Wang, and Adams proposed the adaptive central-upwind WENO scheme, known as WENO-CU6, which leverages a low dissipation central scheme in smooth regions while maintaining strong shock-capturing capabilities through the WENO approach <cit.>. Subsequently, Fu et al. developed the Targeted ENO family of schemes, which further improves upon the capabilities of WENO-CU6 <cit.>. Other researchers, such as Liu et al. <cit.>, Wong and Lele <cit.>, and Subramaniam et al. <cit.>, have also made contributions by developing hybrid weighted interpolation schemes and improved versions of WENO-CU6 with localized dissipative interpolation. Subramaniam et al. <cit.> proposed an explicit-compact interpolation along with compact finite differences that provide higher resolution and more localized dissipation compared to that of Wong and Lele <cit.>. Furthermore, Chamarthi and Frankel <cit.> introduced an adaptive central-upwind scheme, which employs a sixth-order linear-compact scheme and the fifth-order Monotonicity Preserving scheme (MP) using the Boundary Variation Diminishing (BVD) algorithm. Hoffmann, Chamarthi, and Frankel <cit.> have proposed a wave-appropriate central-upwind scheme by taking advantage of the wave structure of the Euler equations <cit.>. The proposed scheme successfully simulates the hypersonic transitional flows on coarse grids. While it outperformed existing methods for benchmark and hypersonic flow test cases, the wave-appropriate centralization approach is still expensive as it requires the reconstruction of the cell-interface values to be carried out for characteristic variables. It is well known that transforming conservative variables to characteristic variables, while necessary for shock-capturing, is expensive and is not required for regions with no discontinuities. On a different note, Feng et al. have employed optimization techniques to develop an adaptive central-upwind scheme that balances dissipation and shock-capturing <cit.>. Their approach involves using multi-objective Bayesian optimization to construct a WENO-based subgrid-scale model for implicit Large Eddy Simulations of compressible flows. Feng et al. further improved the approach by using deep reinforcement learning <cit.>  and Bayesian optimization <cit.> techniques. In their papers, Feng et al. simulated the double periodic shear layer test case and observed that the proposed optimized schemes are free of spurious vortices.

 The objective of the current paper is to propose a generalized and efficient adaptive central-upwind scheme based on the wave-appropriate scheme proposed by Hoffmann, Chamarthi, and Frankel <cit.>. The proposed algorithm aims to reduce computational expenses by reconstructing conservative variables in regions without discontinuities and to prevent spurious vortices in the double periodic shear layer test case.

The rest of the manuscript is organized as follows: In Section <ref>, the governing equations are presented. Section <ref> presents the numerical methods, including the novel adaptive central-upwind approach with conservative and characteristic variable reconstruction. Section <ref> consists of the numerical results with the proposed algorithm and the conclusions are presented in Section <ref>.



§ GOVERNING EQUATIONS
 

In this study, the three-dimensional compressible Navier-Stokes equations are solved in Cartesian coordinates:

    ∂𝐔/∂ t+∂𝐅^𝐜/∂ x+∂𝐆^𝐜/∂ y+∂𝐇^𝐜/∂ z=∂𝐅^𝐯/∂ x+∂𝐆^𝐯/∂ y+∂𝐇^𝐯/∂ z,

where  U is the conservative variable vector, 𝐅^𝐜, 𝐆^𝐜,  𝐇^𝐜 and 𝐅^𝐯, 𝐆^𝐯, 𝐇^𝐯, are the convective (superscript c) and viscous (superscript v) flux vectors in each coordinate direction, respectively.  The conservative variable, convective, and viscous flux vectors are given as:


    
    𝐔 = [   ρ; ρ u; ρ v; ρ w; ρ E ],
              𝐅^c = [       ρ u; ρ u u + p; ρ u v + p; ρ u w + p;     ρ u H ],
              𝐆^c = [       ρ v; ρ v u + p; ρ v v + p; ρ v w + p;     ρ v H ],
              𝐇^c = [       ρ w; ρ w u + p; ρ w v + p; ρ w w + p;     ρ w H ],
             a– d





    [ 𝐅^𝐯=[0, τ_x x, τ_x y, τ_x z, u τ_x x+v τ_x y+w τ_x z-q_x]^T,; 𝐆^𝐯=[0, τ_x y, τ_y y, τ_y z, u τ_y x+v τ_y y+w τ_y z-q_y]^T,; 𝐇^𝐯=[0, τ_x z, τ_y z, τ_z z, u τ_z x+v τ_z y+w τ_z z-q_z]^T, ]

where ρ is density, u, v, and w are the velocities in the x, y, and z directions, respectively, p is the pressure, E = e + (u^2 + v^2 + w^2 )/2 is the specific total energy, and H = E + p/ρ is the specific total enthalpy. The equation of state is for a calorically perfect gas so that e = p/ [ ρ (γ-1) ]^-1 is the internal energy, where γ is the ratio of specific heats. The components of the viscous stress tensor τ and the heat flux q are defined in tensor notations as:

    τ_i j=μ/Re(∂ u_i/∂ x_j+∂ u_j/∂ x_i-2/3∂ u_k/∂ x_kδ_i j),


    q_i=-μ/Re Pr Ma(γ-1)∂ T/∂ x_i,    T= Ma^2γp/ρ,

where μ is the dynamic viscosity, T is the temperature, Ma and Re are the Mach number and Reynolds number, and Pr is the Prandtl number. 






§ NUMERICAL METHODS




Using a conservative numerical method, the governing equations cast in semi-discrete form for a Cartesian cell I_i,j,k = [ x_i-1/2, x_i+1/2] ×[ y_i-1/2, y_i+1/2] ×[ z_i-1/2, z_i+1/2] can be expressed via the following ordinary differential equation: 


    d/d t𝐔̌_i,j,k = 𝐑𝐞𝐬_i,j,k =    - . d𝐅̌^c/d x|_i,j,k - . d𝐆̌^c/d y|_i,j,k - . d𝐇̌^c/d z|_i,j,k
       + . d𝐅̌^v/d x|_i,j,k + . d𝐆̌^v/d y|_i,j,k + . d𝐇̌^v/d z|_i,j,k,


where the check accent, (̌·̌)̌, indicates a numerical approximation of physical quantities, 𝐑𝐞𝐬_i,j,k is the residual function, and the remaining terms are cell centre numerical flux derivatives of the physical fluxes in Equation (<ref>). For brevity, we continue with only the x-direction as it is straightforward to extend to all three dimensions in a dimension-by-dimension manner. The indices  j and k are also dropped for simplicity.





 §.§ Viscous Flux Spatial Discretization Scheme


Viscous fluxes, 𝐅̌^̌v̌, are computed using the fourth-order α-damping scheme of Chamarthi <cit.>, which is based on the α-damping approach of Nishikawa <cit.>. In the one-dimensional scenario, the cell center numerical viscous flux derivative is:


    . d𝐅^v/d x|_i = 1/Δ x( 𝐅^v_i+1/2 - 𝐅^v_i-1/2).


The cell interface numerical viscous flux is computed as,


    𝐅^v_i+1/2 = 
        [                           0;                    -τ_i+1/2; -τ_i+1/2 u_i+1/2 + q_i+1/2,;                             ],         τ_i+1/2 = 4/3μ_i+1/2. ∂ u/∂ x|_i+1/2,
              
            q_i+1/2 = -κ_i+1/2. ∂ T/∂ x|_i+1/2.




For an arbitrary variable, ϕ, the α-damping approach computes cell interface gradients as:


    . ∂ϕ/∂ x|_i+ 1/2 = 1/2( . ∂ϕ/∂ x|_i + . ∂ϕ/∂ x|_i+1) + α/2 Δ x( ϕ_R - ϕ_L ),         ϕ_L = ϕ_i + . ∂ϕ/∂ x|_iΔ x/2, ϕ_R = ϕ_i+1 - . ∂ϕ/∂ x|_i+1Δ x/2,



where, in this work, α = 4. The gradients at cell centers are computed using the fourth order compact finite-difference scheme, as in <cit.>, which is as follows:

    5/14. ∂ϕ/∂ x|_i-1 + . ∂ϕ/∂ x|_i + 5/14. ∂ϕ/∂ x|_i+1 = 1/28 Δ x( ϕ_i+2 - ϕ_i-2) + 11/14 Δ x( ϕ_i+1 - ϕ_i-1).





 §.§ Convective Flux Spatial Discretization Scheme

Similar to the viscous flux discretization, the cell centre numerical convective flux derivative is expressed as:

    . d𝐅̌^c/d x|_i = 1/Δ x( 𝐅̌^c_i+1/2 - 𝐅̌^c_i-1/2),


where i ±1/2 indicates right and left cell interface values, respectively. 𝐅̌^c_i ±1/2 are computed using an approximate Riemann solver. This work uses the componentwise local Lax-Friedrichs (cLLF) <cit.> approximate Riemann solver unless otherwise explicitly stated. The numerical fluxes at cell boundaries computed using a Riemann solver can be expressed in the following standard form:


    𝐅̌^c_i ±1/2 = 1/2[ 𝐅̌^c ( 𝐔^L_i ±1/2) + 𝐅^c ( 𝐔^R_i ±1/2) ] - 1/2| 𝐀_i ±1/2| ( 𝐔^R_i ±1/2 - 𝐔^L_i ±1/2),


where the L and R superscripts denote the left- and right-biased states, respectively, and | 𝐀_i ±1/2| denotes the convective flux Jacobian. The objective is to obtain the left- and right-biased states, 𝐔^L_i ±1/2 and 𝐔^R_i ±1/2. The procedure to obtain these interface values is elucidated in the following sections.




  §.§.§ Linear Schemes:


In this subsection, we provide the details of the calculations of candidate polynomials that can be used to approximate the values of 𝐔^L_i ±1/2 and 𝐔^R_i ±1/2. While the adaptive central-upwind approach to be presented later is a general approach that might potentially be used with most of the schemes, we consider the standard fifth-order upwind scheme <cit.> and the corresponding sixth-order scheme, and the Gradient-Based Reconstruction (GBR) schemes proposed in <cit.>.

Fifth-order upwind scheme: The well-known fifth-order upwind schemes for obtaining the values of the left and right interfaces are as follows:

    ϕ_i+1 / 2^L5, Linear=1/30ϕ_i-2-13/60ϕ_i-1+47/60ϕ_i+0+9/20ϕ_i+1-1/20ϕ_i+2, 
       ϕ_i+1 / 2^R5, Linear=1/30ϕ_i+3-13/60ϕ_i+2+47/60ϕ_i+1+9/20ϕ_i+0-1/20ϕ_i-1,

where ϕ is an arbitrary variable, either conservative (𝐔) or characteristic (𝐂) variables are used in the present paper. The superscripts L5 and R5 represent left and right-biased fifth-order formulas, and the Linear superscript represents the scheme is linear. Linear schemes cannot be used for flows with discontinuities as they will lead to oscillations. The details of the non-linear scheme will be presented later.

Upwind-biased GBR scheme: The GBR schemes differ slightly from the fifth-order mentioned above scheme. As such, the interface values depend not only on the cell centre values but also on the gradients at the cell centres. The advantage of having gradients is that they can be used in the parts of the subroutine, such as viscous flux discretization, shock capturing, post-processing quantities, and so on <cit.>. GBR method employs the first two moments of the Legendre polynomial evaluated on x_i-1 / 2≤ x ≤ x_i+1 / 2 and written for a general variable, ϕ, as:


    ϕ(x) = ϕ_i + ϕ'_i(x-x_i) + 3 ϕ”_i/2κ[(x-x_i)^2-Δ x_i^2/12],


where ϕ'_i and ϕ”_i respectively represent the first and second derivatives of ϕ_i. If x = x_i ±Δ x/2 and κ = 1/3, the following equations for the left- and right-biased states are obtained:
    

    
    ϕ^L,GBR, Linear_i+1/2   = ϕ_i+0 + Δ x/2ϕ'_i+0 + Δ x^2/12ϕ”_i+0,
      ϕ^R,GBR, Linear_i+1/2   = ϕ_i+1 - Δ x/2ϕ'_i+1 + Δ x^2/12ϕ”_i+1.

    

Like the fifth-order upwind scheme, the superscripts L, GBR and R, GBR represent left and right-biased gradient-based reconstruction formulas. In this work, ϕ'_i was computed using eighth order explicit central differences <cit.>:
    

    ϕ'_i = 1/Δ x( 1/280ϕ_i-4 - 4/105ϕ_i-3 + 1/5ϕ_i-2 - 4/5ϕ_i-1 + 4/5ϕ_i+1 - 1/5ϕ_i+2 +4/105ϕ_i+3 - 1/280ϕ_i+4).


ϕ”_i was computed by the following formula <cit.>:


    ϕ”_i = 2/Δ x( ϕ_i+1 - 2 ϕ_i + ϕ_i-1) - 1/2 Δ x( ϕ'_i+1 - ϕ'_i-1).




Central schemes: The abovementioned schemes are left and right-biased upwind schemes. They provide necessary dissipation in some regions but are unsuitable for turbulent simulations as they require a low dissipation central scheme. A central scheme can be obtained by averaging the left- and right-biased upwind reconstructions (they can also be derived in other possible ways as well).

    
    ϕ^L_i+1/2 = ϕ^C_i+1/2 = ( 1 - η) ϕ^L, Linear_i+1/2 + ηϕ^R, Linear_i+1/2,
             a
    ϕ^R_i+1/2 = ϕ^C_i+1/2 = ηϕ^L, Linear_i+1/2 + ( 1 - η) ϕ^R, Linear_i+1/2,
             b


where η = 0.5 and ( ·)^C denotes the centralized reconstruction. The well-known sixth-order central scheme can be obtained by averaging the left and right-biased reconstruction formulas given by Equation (<ref>), and is as follows:

    ϕ^C6,Linear_i+1/2=1/2(ϕ_i+1 / 2^L5, Linear + ϕ_i+1 / 2^R5, Linear)=1/60(ϕ_i-2-8 ϕ_i-1+37 ϕ_i+37 ϕ_i+1-8 ϕ_i+2+ϕ_i+3).

Similarly, the centralized GBR scheme can be obtained by averaging the Equations (<ref>), and is as follows:


    ϕ^C,GBR,Linear_i+1/2=1/2(ϕ_i+1 / 2^L,GBR, Linear + ϕ_i+1 / 2^R,GBR, Linear).






  §.§.§ Monotonicity Preserving approach


The linear reconstruction schemes represented by Equations (<ref>) and (<ref>) are susceptible to oscillations when there are discontinuities present. To address this issue, we utilized MP limiting. The following details the MP limiting procedure specifically for the left-biased state, although the procedure is identical for the right-biased state. The limiting procedure is typically carried out for the characteristic variables, 𝐂, for oscillation-free results, as explained in <cit.>. The transformation from conservative variables to characteristic variables will be explained later. The MP limiting criterion for the GBR scheme <cit.> is:


    𝐂^L_i+1/2 = 
        𝐂^L,Linear_i+1/2   if ( 𝐂^L,Linear_i+1/2 - 𝐂_i ) ( 𝐂^L,Linear_i+1/2 - 𝐂^L,MP_i+1/2) ≤ 10^-40, 
    𝐂^L,Non-Linear_i+1/2   otherwise,


where 𝐂^L,Linear_i+1/2 corresponds to Eqn. <ref> and the remaining terms are:


    𝐂^L,Non-Linear_i+1/2   = 𝐂^L,Linear_i+1/2 + minmod( 𝐂^L,MIN_i+1/2 - 𝐂^L,Linear_i+1/2, 𝐂^L,MAX_i+1/2 - 𝐂^L,Linear_i+1/2),
    𝐂^L,MP_i+1/2   = 𝐂^L,Linear_i+1/2 + minmod[ 𝐂_i+1-𝐂_i, 𝒜( 𝐂_i-𝐂_i-1) ], 
    𝐂^L,MIN_i+1/2   = max[ min( 𝐂_i, 𝐂_i+1, 𝐂^L,MD_i+1/2), min( 𝐂_i, 𝐂^L,UL_i+1/2, 𝐂^L,LC_i+1/2) ],
    𝐂^L,MAX_i+1/2   = min[ max( 𝐂_i, 𝐂_i+1, 𝐂^L,MD_i+1/2), max( 𝐂_i, 𝐂^L,UL_i+1/2, 𝐂^L,LC_i+1/2) ], 
    𝐂^L,MD_i+1/2   = 1/2( 𝐂_i + 𝐂_i+1) - 1/2 d^L,M_i+1/2, 𝐂^L,UL_i+1/2 = 𝐂_i + 4 ( 𝐂_i - 𝐂_i-1), 
    𝐂^L,LC_i+1/2   = 1/2( 3 𝐂_i - 𝐂_i-1) + 4/3 d^L,M_i-1/2,          d^L,M_i+1/2 = minmod( 0.5(d_i + d_i+1),2d_i, 2d_i+1), 
    
            d_i    = 2 ( 𝐂_i+1 - 2𝐂_i + 𝐂_i-1) - Δ x/2( 𝐂'_i+1 - 𝐂'_i-1),


where 𝒜 = 4 and minmod( a,b ) = 1/2[ sgn(a) + sgn(b) ] min( | a |, | b | ). The GBR method that employs explicit finite differences and the above MP limiter is called MEG (Monotonicity Preserving Explicit Gradient) as in <cit.>. As it uses the eighth-order explicit gradients, it is denoted as MEG8. For the fifth-order upwind scheme, the computations of the d_i+1 / 2^M and d_i differ from that of the GBR scheme and are as follows <cit.>:


    d_i+1 / 2^M     =minmod(4 𝐂_i-𝐂_i+1, 4 𝐂_i+1-𝐂_i, 𝐂_i, 𝐂_i+1),
    
        d_i     =𝐂_i-1-2 𝐂_i+𝐂_i+1.

The fifth-order upwind method with MP limiting is known as MP5, as in <cit.>, in this paper. Only the adaptive central-upwind scheme is considered in this paper and is denoted as MP6.




  §.§.§ Novel adaptive central-upwind and conservative-characteristic variable reconstruction


This section presents the novel central-upwind scheme with adaptive conservative-characteristic variable reconstruction. Before proceeding to the novel algorithm, a brief background of the approach is presented. 


	
  * In Ref. <cit.> authors took advantage of the wave structure of the Euler equations. Once the variables are transformed from physical to characteristic space, the characteristic variables have specific properties. The first and last variable is known as the acoustic waves, the second variable is the entropy wave, and the rest are known as shear waves. Together, the entropy and shear waves are known as linearly degenerate waves. The density varies across the entropy wave in characteristic space, and the rest of the variables remain unchanged. Ref. <cit.> took advantage of this and reduced the frequent activation of the MP criterion, Equation (<ref>). Shockwaves are detected by the Ducros sensor, which is used for waves (characteristic variables) other than entropy waves. However, all the waves are computed using the upwind, inherently dissipative schemes.
	
  * The approach is improved in <cit.> that all the linearly degenerate waves are computed using a central scheme if the necessary discontinuity detection criteria are met. The acoustic waves are still computed using an upwind scheme for stability, and the proposed approach significantly improved benchmark test cases and was able to predict hypersonic transitional flows.


One disadvantage of Ref. <cit.> and <cit.> is that the variables are transformed from physical (conservative) to characteristic variables, which is computationally expensive. It would be efficient to convert to characteristic variables only in regions where there are discontinuities, which are, strictly speaking, localized to some regions of the flow, and the rest of the regions are computed using conservative variables. While the Ducros sensor detects shockwaves, the contact discontinuities must be detected by a different approach and done in physical space. In physical space, across the contact discontinuity, only density varies, and the rest of the variables remain continuous; therefore, if a contact sensor coupled with the Ducros sensor can identify the discontinuities, then the regions without discontinuities can be computed using conservative variables thereby avoiding the variable transformation. The following steps describe the novel algorithm that avoids this expensive transformation procedure:



  Step 1. Compute the linear left and right-biased reconstructions of the variable density either by the fifth-order upwind schemes or the GBR schemes (depending on the linear scheme of choice the user has decided upon, either MEG8 and MP). 

  Step 2. Evaluate the MP limiting criterion (Equation (<ref>))  for the variable density (ρ = 𝐔_1) using the following equations

    L_ρ   =( 𝐔^L,Linear_i+1/2,1 - 𝐔_i+0,1) ( 𝐔^L,Linear_i+1/2,1 - 𝐔^L,MP_i+1/2,1),
     
    
    R_ρ   = ( 𝐔^R,Linear_i+1/2,1 - 𝐔_i+1,1) ( 𝐔^R,Linear_i+1/2,1 - 𝐔^R,MP_i+1/2,1).


The above equations can be considered as contact discontinuity detectors in physical space. Also evaluate the Durcos sensor, Ω_d, which is computed as follows: 

    Ω_d = | -p_i-2 + 16 p_i-1 - 30 p_i + 16 p_i+1 - p_i+2|/| +p_i-2 + 16 p_i-1 + 30 p_i + 16 p_i+1 + p_i+2|( ∇·𝐮)^2/( ∇·𝐮)^2 + | ∇×𝐮|^2,

where 𝐮 is the velocity vector, and the derivatives of velocities are computed by the Equation (<ref>) as in <cit.>.  We modify Ω_d by using its maximum value in a three-cell neighborhood:


    Ω_d = max( Ω_i+m),   for  m = -1,0,1.



  Step 3. If both the Ducros sensor and the contact discontinuity detectors identify that there are no discontinuities in the region as per the following conditions:


    if (Ω_d > 0.01)  & (L_ρ≤ 10^-40)  & (R_ρ≤ 10^-40).


One can reconstruct the conservative variables according to the following algorithm:

In all directions:


    𝐔^L,R_i+1/2,b = 
        {[                                              if  b = 1,5: 𝐔^C,Linear_i+1/2,b  i.e. ρ & ρ E  (η =0.5, in Eqn. <ref>) ].


In x-direction:


    𝐔^L,R_i+1/2,b = 
        {[                                                   if  b = 2:        𝐔^L,R,Linear_i+1/2,b  i.e. ρ u  (η =1, in Eqn. <ref>);                                                 if  b = 3,4: 𝐔^C,Linear_i+1/2,b  i.e. ρ v  & ρ w  (η =0.5, in Eqn. <ref>) ].


In y-direction:


    𝐔^L,R_i+1/2,b = 
        {[                                                   if  b = 3:        𝐔^L,R,Linear_i+1/2,b  i.e. ρ v  (η =1, in Eqn. <ref>);                                                 if  b = 2,4: 𝐔^C,Linear_i+1/2,b  i.e. ρ u  & ρ w  (η =0.5, in Eqn. <ref>) ].


In z-direction:


    𝐔^L_i+1/2,b = 
        {[                                                   if  b = 4:        𝐔^L,R,Linear_i+1/2,b  i.e. ρ w  (η =1, in Eqn. <ref>);                                                 if  b = 2,3: 𝐔^C,Linear_i+1/2,b  i.e. ρ u  & ρ v  (η =0.5, in Eqn. <ref>) ].


Only linear schemes are required in the above procedure, as there are no discontinuities in the computational region.



  Step 4. If the criterion in Equation (<ref>) is not met, then one can proceed with the characteristic variable reconstruction as in Ref. <cit.>. Compute Roe-averaged variables at the interface to construct the left, 𝐋_n, and right, 𝐑_n, eigenvectors of the normal convective flux Jacobian. 

    

  Step 5. For the GBR scheme transform 𝐔_i, 𝐔'_i, and 𝐔”_i to characteristic space by multiplying them by 𝐋_n:

    
        
    𝐂_i+m,b = 𝐋_n,i+1/2𝐔_i+m, 
                  𝐂'_i+m,b = 𝐋_n,i+1/2𝐔'_i+m,
                  𝐂”_i+m,b = 𝐋_n,i+1/2𝐔”_i+m.
                 a– c

    
    
 For the fifth-order upwind scheme transform only 𝐔_i to characteristic space by multiplying them by 𝐋_n
    
        
        
    𝐂_i+m,b = 𝐋_n,i+1/2𝐔_i+m,

    

    for m = -2,-1,0,1,2,3 and b = 1,2,3,4,5, representing the vector of characteristic variables which are defined as follows in the current implementation: 

    
        
 
  Step 6. Using Equations (<ref>), obtain the unlimited reconstruction to cell interfaces in characteristic space via:


    
    𝐂^L, GBR, Linear_i+1/2,b   = 𝐂_i+0,b + Δ x/2𝐂'_i+0,b + Δ x^2/12𝐂”_i+0,b,
      𝐂^R, GBR, Linear_i+1/2,b   = 𝐂_i+1,b - Δ x/2𝐂'_i+1,b + Δ x^2/12𝐂”_i+1,b.

    


For the fifth-order reconstruction schemes Equations (<ref>) are used for the interface values  and are as follows:

    
    𝐂_i+1/2,b^L5, Linear   =1/30𝐂_i-2,b-13/60𝐂_i-1,b+47/60𝐂_i+0,b+9/20𝐂_i+1,b-1/20𝐂_i+2,b,
      𝐂_i+1/2,b^R5, Linear   =1/30𝐂_i+3,b-13/60𝐂_i+2,b+47/60𝐂_i+1,b+9/20𝐂_i+0,b-1/20𝐂_i-1,b.

    


The left-biased reconstruction is then treated by the following algorithm:


    𝐂^L_i+1/2,b = 
        {[                                                                                                                                      if  b = 1,5:                                                             𝐂^L,Non-Linear_i+1/2,b   if Ω_d > 0.01, 
                
    𝐂^L,Linear_i+1/2,b   otherwise,;                                                                                                                                        if  b = 2: 𝐂^L,Non-Linear_i+1/2,b   if ( 𝐂^L,Linear_i+1/2 - 𝐂_i ) ( 𝐂^L,Linear_i+1/2 - 𝐂^L,MP_i+1/2) ≥ 10^-40, 
                
    𝐂^C,Linear_i+1/2,b   otherwise,;                                                                                                                                      if  b = 3,4:                                                             𝐂^L,Non-Linear_i+1/2,b   if Ω_d > 0.01, 
                
    𝐂^C,Linear_i+1/2,b   otherwise. ].


A similar procedure is carried out for the right-biased reconstruction.



  Step 7. After obtaining 𝐂^L,R_i+1/2,b, the variables are transformed back to physical fields:


    𝐔^L,R_i+1/2 = 𝐑_n,i+1/2𝐂^L,R_i+1/2.





The nonlinear scheme with the above procedure along with the GBR schemes (Equations (<ref>)) as linear schemes is denoted as MEG8-CC (CC for conservative-characteristic variables). Similarly, if fifth-order schemes (Equations (<ref>)) are used as linear schemes, then the nonlinear scheme is denoted as MP6-CC in this paper. The approach is generic, so even higher-order schemes like the seventh-order upwind schemes or implicit gradient schemes as in <cit.> can also be used. For simplicity, only the MEG8-CC and MP6-CC schemes are presented in this paper. Steps 4-7 are the same as that of the wave-appropriate centralization approach presented in <cit.>, and such scheme is denoted as MEG8-C as in <cit.>. The novel approach here is steps 1-3, where the conservative variables are reconstructed in a unique central-upwind approach if the shock and contact discontinuity sensors are satisfied. 

 
  
 It is essential to note that a complete central scheme is not used even though the criteria in Equation (<ref>) are satisfied. In Equations (<ref>), (<ref>), (<ref>), and (<ref>) not all variables are centralized. While ρ and  ρ E are centralized in all directions, the variables ρ u, ρ v and ρ w are upwind or centralized depending on the direction. The characteristic variable centralization explained in steps 4-7 and the conservative variable reconstruction in steps 1-3 above have some similar characteristics and will be explained below. A complete central scheme would still give spurious vortices for the double shear layer test case, and the proposed algorithm would avoid spurious vortices.
 
 
 
 Relation between centralization in characteristic and conservative variable space:

In <cit.>, it was found that using a centralized reconstruction for all but the characteristic acoustic waves, i.e. all the linearly degenerate waves; entropy wave and the shear waves, was a robust and superior solution. The similarity between the central-upwind scheme in conservative variable space and the central-upwind scheme in characteristic space is explained below. In order to obtain charactertic variables, 𝐂, the conservative variables,  𝐔, are multiplied by the left eigenvectors, 𝐋_n. For simplicity we consider the two-dimensional scenario:
    
        
    𝐂_b = 𝐋_n𝐔_b,




    ( [ 𝐂_1; 𝐂_2; 𝐂_3; 𝐂_4;     ]) 
            = ([ 1/2(γ-1/2 c^2 q^2+q_n/c)    -1/2(γ-1/c^2 u+n_x/c)    -1/2(γ-1/c^2 v+n_y/c)                γ-1/2 c^2;                         ;          1-γ-1/2 c^2 q^2                γ-1/c^2 u                γ-1/c^2 v                 -γ-1/c^2;                         ;                     -q_ℓ                      ℓ_x                      ℓ_y                        0;                         ; 1/2(γ-1/2 c^2 q^2-q_n/c)    -1/2(γ-1/c^2 u-n_x/c)    -1/2(γ-1/c^2 v-n_y/c)                γ-1/2 c^2 ])
            ( [ 𝐔_1; 𝐔_2; 𝐔_3; 𝐔_4 ]),

    where n = [n_x  n_y]^t and [l_x  l_y]^t is a tangent vector (perpendicular to n) such as [l_x  l_y]^t = [-n_y  n_x]^t. By taking n = [1, 0]^t and [0, 1]^t we obtain the corresponding eigenvectors in x and y directions (in the x-direction, n_x = 1, whereas n_y=0. The y-direction is analogous.). Furthermore,  q^2 = u^2 + v^2, q_l = u l_x + v l_y, and q_n = u n_x + v n_y. In the matrix (<ref>), the 𝐂_1 and 𝐂_4 are acoustic waves, 𝐂_2 is the entropy wave, and the 𝐂_3 is the shear wave. The characteristic variable 𝐂_3 in the characteristic space is as follows:
    
    𝐂_3 = -q_l 𝐔_1 + l_x𝐔_2 + l_y𝐔_3+0𝐔_4.

    
    In x- direction, [l_x  l_y]^t = [-n_y  n_x]^t= [0  1]  which means q_l = u l_x + v l_y = v. Therefore,
    
        
    𝐂_3 = -v 𝐔_1 + 1𝐔_3 = -ρ v + ρ v ≈ 0  (a small value in computations).

    
        In y- direction, [l_x  l_y]^t = [-n_y  n_x]^t= [-1  0]  which means q_l = u l_x + v l_y = -u. Therefore,
    
        
    𝐂_3 = u 𝐔_1 + (-1)𝐔_2 = ρ u - ρ u ≈ 0  (a small value in computations).

    
   The critical aspect is that in x- direction ρ v is being reconstructed in the characteristic space and centralized if the Ducros sensor criterion is satisfied. Likewise, in the y- direction, ρ u is reconstructed in the characteristic space and centralized if the Ducros sensor criterion is satisfied. Similarly, the variable ρ v is centralized in conservative variable space in x- direction (Equations (<ref>)), and  ρ u is centralized in conservative variable space in y- direction (Equations (<ref>)), respectively. The analysis straightforwardly extends to the three-dimensional scenario.
       
    Similarly, the characteristic variable 𝐂_2 represents the contact discontinuity in characteristic space and across a contact discontinuity density changes in physical space. The corresponding variables are centralized in their respective computational spaces if there is no contact discontinuity in the computation region. It explains the relation between the centralization in characteristic and conservative variable space.
    


§ RESULTS AND DISCUSSION


This section tests the proposed numerical schemes for various benchmark test cases. For time integration, we use the explicit third-order TVD Runge-Kutta method <cit.>. Time integration is performed with a CFL = 0.4 for all the problems. The advantages of the proposed algorithm are compared with the TENO5 scheme <cit.>. All the computations are carried out on the authors' Mac mini with an M1 processor.

Double periodic Shear layer

In the first test case, we demonstrate that the proposed algorithm will prevent spurious vortices for the double periodic shear layer test case <cit.>. The test involves two initially parallel shear layers that develop into two significant vortices at t = 1. The non-dimensional parameters for this test case are presented in Table <ref>.


The initial conditions were:

    
    p    = 1/γMa^2, ρ=1,  u= 
            tanh[ θ (y-0.25) ],     if  (y ≤ 0.5), 
    tanh[ θ (0.75-y) ],     if  (y > 0.5),
            
    
            v    = 0.05 sin[ 2 π(x+0.25) ] for θ  = 80, and v = 0.05 sin[ 2 π(x) ] for θ  = 120,


The parameter θ defines the initial shear layer width. Two configurations are considered with θ =80 and θ =120, as shown in the above initial conditions. The reference solutions for θ =80 and θ =120, shown in Fig. <ref>, are computed with MEG8-C scheme on a grid sizes of of 800 × 800  and 512 × 512, respectively. For this test case, if the grid is under-resolved, unphysical braid vortices and oscillations can occur on the shear layers according to the literature <cit.>. Observing the well-resolved reference solution, there are no braid vortices or oscillations that form on the shear layers. 


Fig. <ref> displays the z-vorticity computed for MEG8-CC, MEG8-C, MEG8-central, MP6-CC, and TENO5 schemes on a grid size of 160 × 160 (except for Fig. <ref>). From Fig. <ref>, it can be seen that MEG8-CC best approximates the reference solution, as there are no unphysical braid vortices and oscillations. MEG8-C, Fig. <ref> is also relatively free of spurious vortices. MEG8-CC and MEG8-C use the unique adaptive central-upwind algorithm proposed, i.e. not all variables are reconstructed using the central scheme. Fig. <ref> shows the results obtained by the fully central scheme for all the variables if both Ducros and contact discontinuity sensors criteria are satisfied (Equation <ref>). The fully central scheme produced unphysical vortical structures unlike the adaptive central-upwind scheme (Equations (<ref>), (<ref>), (<ref>), and (<ref>) in conservative variable computational space). Likewise, the MP6-CC scheme, vorticity contours shown in Fig. <ref>, also has no spurious vortices. The stencils of the MEG8-CC and MP6-CC schemes are significantly different, yet both schemes resulted in results free of spurious vortices.

Figs <ref> and <ref> show the vorticity contours computed by the TENO5 scheme on grid sizes of 160 × 160 and 320 × 320, respectively. Even on a grid size of 320 × 320, the TENO5 scheme produced spurious results in the primary vortex. It is important to note that the MP6-CC and TENO5 schemes have the same stencil for reconstructing cell-interface values. It indicates that the proposed algorithm produces results free of oscillations and spurious vortices for this test case. Feng et al. <cit.> simulated this test case on a grid size of 320 × 320 using TENO5 and TENO8. Even on that fine grid resolution, the TENO8 scheme produced spurious vortices along the shear layer (readers can refer to Fig. 13 of <cit.>). These findings suggest that the proposed algorithm with a unique central-upwind scheme can prevent spurious vortices for this test case, which has been the subject of many papers in the literature. 


Table <ref> shows the computational cost of the proposed algorithms, MEG8-CC and MP6-CC, and MEG8-C and TENO5. The MP6-CC and MEG8-CC are almost similar in computational cost, whereas MEG8-C is twice as expensive due to the reconstruction of characteristic variables only. The TENO5 scheme will require four times more grid points and is 13 times more costly than the proposed schemes.



This test case does not involve shocks or contact discontinuities; therefore, the Ducros sensor or contact sensor was unnecessary. Therefore, a linear combination of the central-upwind scheme algorithm given by Equations (<ref>), (<ref>), (<ref>), and (<ref>)  alone was sufficient. Central reconstruction scheme for all the variables except for ρ u in x- direction and ρ v in y- direction will suffice. In this regard, simulation is also carried out using the third-order upwind and fourth-order central scheme combination to show the uniqueness of the proposed algorithm. The third-order upwind scheme for left and right interfaces is as follows:


    ϕ_i+1 / 2^L3, Linear   = 1/6(-ϕ_i-1 + 5ϕ_i + 2ϕ_i+1), 
    ϕ_i+1 / 2^R3, Linear   = 1/6(2ϕ_i + 5ϕ_i+1 -ϕ_i+2),

and the fourth-order central scheme is:


    ϕ^C4,Linear_i+1/2=1/2(ϕ_i+1 / 2^L3, Linear + ϕ_i+1 / 2^R3, Linear)=1/12(- ϕ_i-1+7 ϕ_i+7 ϕ_i+1- ϕ_i+2).


Figs <ref> and <ref> show the vorticity contours computed by the linear third-fourth order scheme and linear fifth-sixth-order scheme using the proposed central-upwind algorithm on grid size of 160 × 160, respectively. Even the linear schemes are free of spurious vortices for this test case using the proposed algorithm. The only difference between the third-fourth-order scheme and the fifth-sixth-order scheme is the thickness of the shear layer. The third-fourth-order scheme is thicker than the fifth-sixth-order scheme, which can be attributed to the difference in spectral properties and order of accuracy. With the proposed algorithm, even a third-fourth-order linear scheme prevented spurious vortices, unlike the TENO5 scheme. Furthermore, for the simulations using the linear schemes (it is applicable also for the MEG8-CC and MP6-CC schemes), the fourth-order alpha-damping scheme of Nishikawa <cit.> was used instead of the compact viscous scheme proposed in <cit.> (described in Section <ref>) and yet there are no spurious vortices. It indicates that the viscous flux discretization had little impact on this test case as long as the discretization has good spectral properties <cit.>.


Finally, Fig. <ref> displays the z-vorticity computed for MEG8-CC, linear third-fourth scheme and linear fifth-sixth-order scheme on a grid size of 96 × 96  for θ=80. It can be observed that all the schemes are free of spurious vortices on such a coarse grid. These results show the advantages and uniqueness of the proposed algorithm. The algorithm also produced significantly better results than the kinetic energy preserving schemes considered in <cit.> on much coarser grids.


Finally, simulation is also carried out for inviscid, Re=∞, scenario  with θ=80 for v = 0.05 sin[ 2 π(x) ] in the initial conditions as in <cit.>. Fig. <ref> displays the z-vorticity computed for MEG8-CC, MP6-CC and the linear fifth-sixth-order scheme on a grid size of 320 × 320. In Ref. <cit.>, authors have computed the simulation on a grid size of 512 × 512 for these initial conditions. With the present algorithm, even the inviscid case is free of spurious vortices on a grid size that is two and a half times smaller than <cit.>. The fifth-sixth-order scheme has a slightly thicker shear layer width than MEG8-CC. The proposed algorithm performed better than the optimized TENO8 scheme considered in <cit.> (see Fig. 13 of the concerned reference). All the schemes, both convective and viscous flux discretization, are only second-order accurate <cit.>, as fluxes are not reconstructed directly, making these results even more significant.





Possible reasons for lack of spurious vortices:


	
  * Contact discontinuities are characterized by the condition that no mass flows through the discontinuity (normal velocities are zero) and continuous pressure. This results in a jump in density and tangential velocity (page 137 of <cit.>). The tangential velocities are v in the x- direction and u in the y- direction. In the absence of contact discontinuity, the tangential velocities will be continuous. However, while tangential velocities are discontinuous for inviscid flows, they are not discontinuous for viscous flow simulations. The driving force behind interfacial instabilities is shear. In the interface between two different fluids, the conservation of momentum requires that the tangential viscous stresses be equal, leading to the continuity of the tangential velocity. The fluids would slip in the actual inviscid limit, resulting in discontinuous tangential velocity. However, in the presence of viscosity (even artificial viscosity due to upwind reconstructions), the tangential velocity remains continuous across an interface according to Batchelor <cit.>. In the present algorithm, the tangential velocities are deemed continuous if the MP criterion for density, and Ducros sensor for shocks are met and, therefore, reconstructed using a central scheme in conservative variable space - Equations (<ref>), (<ref>), and (<ref>) or wave appropriate centraliztion in characteristic space.
	
  * The shear layer test case does not involve shocks, and the density is almost constant and therefore does not require any shock or contact discontinuity sensor and can be simulated using linear schemes itself, i.e. directly using the Equations (<ref>), (<ref>), (<ref>), and (<ref>). The flow is viscous, and the density is almost constant; therefore, the tangential velocities are also continuous, which is what the proposed algorithm is incidentally doing: centralizing the ρv in x- direction and ρu in y- direction.
	
  * Potential improvement of optimization of the numerical schemes: In <cit.> authors have optimized the parameter η in the following (a similar one) equation:


    ϕ^C_i+1/2 = ( 1 - η) ϕ^L, Linear_i+1/2 + ηϕ^R, Linear_i+1/2


The parameter η is constant, once the optimum value is found, for all the variables. Despite the optimization procedure, the optimized schemes, including the Ducros sensor, lead to spurious vortices even on much larger grid sizes than the current paper (see Figure 13 of Ref. <cit.> and the corresponding discussion). In this paper, the value of η is either 0.5 or 1 and is different for different variables in various directions. If one has to optimize, then η has to be an independent value for each variable in each direction. It implies that instead of optimizing one parameter η, there will be η_1 to η_8, eight parameters to optimize in a two-dimensional scenario (four variables in two directions)
	
  * It is only a partial explanation because there is no clear understanding of why some variables are reconstructed using upwind reconstruction for this test case (it is necessary for test cases with discontinuities as in <cit.>). Nevertheless, there are no spurious vortices even with the linear schemes using the current central-upwind algorithm, which is a significant outcome.




The shear layer test case does not have shock waves or contact discontinuities. The algorithm is tested for such cases in the following test cases, and its computational benefits and robustness are shown in the subsequent test cases.
Kelvin Helmholtz instability


First, the Kelvin-Helmholtz instability case is considered. This instability arises due to an unstable velocity gradient at the interface between the two fluids, leading to vortices and mixing patterns. The test case has the following initial conditions over a periodic domain of [0, 1] × [0, 1],


    [                        p=2.5,   ρ(x, y)={[ 2,  If  0.25<y ≤ 0.75,;              1,  else, ].;                         u(x, y)={[ 0.5,    If  0.25<y ≤ 0.75,;              -0.5,  else , ].; v(x, y)=0.1 sin (4 π x){exp[-(y-0.75)^2/2 σ^2]+exp[-(y-0.25)^2/2 σ^2]}, where σ = 0.05/√(2). ]

The computational domain is discretized with 512 cells in each direction, and the final time is t = 0.8. The numerical solutions are computed using various schemes and are depicted in Fig. <ref>.  We can observe that all the schemes capture complex structures and small-scale vortices. Both MP6-CC and MEG8-CC are free of oscillations despite being reconstructed with the adaptive conservative-characteristic variable reconstruction. Despite sharing the same stencil, MP6-CC shows more vortical structures than the TENO5 scheme. 


Table <ref> shows the schemes' computational times. MP6-CC is also 15% faster than the TENO5 scheme. For computational efficiency, the MEG8-CC scheme takes the same time as the TENO5 scheme with improved results. MEG8-CC is also 36 % faster than the MEG8-C presented in <cit.>, indicating the proposed approach's computational benefit.





Rayleigh-Taylor instability


The Rayleigh-Taylor instability is a hydrodynamic phenomenon that occurs when a dense fluid is above a less dense fluid in a gravitational field. This instability arises from an unstable density gradient, causing the denser fluid to accelerate downward and the lighter fluid to accelerate upward. As a result, complex and turbulent mixing patterns emerge between the two fluids. According to <cit.>, the initial conditions for the Rayleigh-Taylor instability are as follows:


    (ρ,u,v,p)=
       (2.0, 0, -0.025√(5p/3ρcos(8π x)), 2y+1.0),   0≤ y< 0.5,
       (1.0, 0, -0.025√(5p/3ρcos(8π x)), 1y+1.5),     0≤ y≤ 0.5,


In this test case, simulations were conducted using a uniform mesh with a 128 × 512 resolution until t=1.95. The computational domain for the simulations was [0, 1/4]× [0,1], and the adiabatic constant γ was set to 5/3. The flow conditions were specified as ρ=1, p=2.5, and u=v=0 on the top boundary, and ρ=2, p=1.0, and u=v=0 on the bottom boundary. Source terms S=(0,0,ρ, ρ v) were added to the Euler equations. The density distribution of the Rayleigh-Taylor instability problem is shown in Figs. <ref>. The proposed approaches, MP6-CC and MEG8-CC, are free of oscillations as the MP criterion with density as the variable could detect contact discontinuity and prevent oscillations. The MP6-CC is again significantly faster, as shown in Table <ref>, than the TENO5 scheme with the thinner material interface and more vortical structures. MEG8-CC is 35 % faster than the MEG8-C approach as it avoids the costly characteristic variable transformation and produces slightly more vortical structures.







Double Mach Reflection



In this test example, we will examine the proposed algorithm for the Double Mach Reflection (DMR) case originally presented by Woodward and Collela <cit.>. This scenario involves a Mach 10 unsteady planar shock wave interacting with a 30-degree inclined surface, forming intricate flow patterns. The computational domain for this particular case spans [x,y] = [0, 3] × [0, 1], and our simulations will be conducted on a uniform grid of 768 × 256. The specified initial conditions are as follows:

    (ρ, u, v, p)={[               (1.4,0,0,1),      if  y<1.732(x-0.1667); (8,7.145,-4.125,116.8333),                  otherwise ].


For this test case, we utilized the shock-stable cLLF Riemann solver from the reference <cit.> to avoid the “carbuncle” phenomenon. The simulations were conducted up to a final time of t=0.3. Reflecting wall conditions were applied to the bottom boundary for 0.1667 < x ≤ 4.0, and post-shock conditions for 0.0 ≤ x ≤ 0.1667. The top boundary was set using the exact solution of the time-dependent oblique shock. Density contours from all the considered schemes are depicted in Figs. <ref>. Among the results, the TENO5 and MP6-CC schemes exhibited the highest dissipation and failed to resolve small-scale vortical structures. Due to their superior dispersion properties, the MEG schemes captured the shocks sharply compared to the TENO and MP6-CC schemes. Furthermore, within the MEG family of schemes, the MEG8-CC outperformed the MEG8-C scheme in capturing vortical structures and near-wall jets. Similar to the previous Example <ref>, it can be concluded that the MP detection criterion with the density as the variable and the Ducros sensor effectively identify discontinuities, aiding in the reconstruction of conservative variables. 

Once again, from Table <ref>, the computational cost of the MEG8-CC scheme is cheaper than the MEG8-C and TENO5 schemes. MP6-CC is the most efficient of all the schemes, and MEG8-CC is a close second.




Riemann Problem



This example considers the two-dimensional Riemann problem of configuration 3 <cit.> for testing the proposed schemes. The initial conditions for the configuration considered are as follows:


    (ρ, u, v, p)={[             (1.5,0,0,1.5),           if  x>0.5, y>0.5;      (0.5323,1.206,0,0.3),           if  x<0.5, y>0.5; (0.138,1.206,1.206,0.029),           if  x<0.5, y<0.5;      (0.5323,0,1.206,0.3),           if  x>0.5, y<0.5 ].


The initial conditions described above result in four shocks at the interfaces of the four quadrants, and the Kelvin-Helmholtz instability along the slip lines gives rise to small-scale structures, often used as a benchmark for evaluating the numerical dissipation of a particular scheme. In this test case, simulations are conducted within a computational domain of [x,y] = [0,1]× [0,1], using a uniform grid of 1024 × 1024, until a final time, t=0.3. 


The density contours obtained by the considered schemes are presented in Figs. <ref>. The small-scale structures produced by the TENO scheme, Fig. <ref>, are dissipative compared to that of MP6-CC, Fig. <ref>, and they both share the same stencil for reconstruction. These results indicate that a properly designed scheme with limiters can also outperform the ENO-type schemes (considered standard and most popular schemes in the literature for compressible flow simulations). Furthermore, the small-scale vortices obtained by the MEG8-CC scheme are significantly better than the MEG8-C and the other schemes. 

It is important to note that this test case does not involve any contact discontinuities, and there is little difference between the MEG8-C and MEG8-CC schemes regarding computational cost, shown in Table <ref>. MP6-CC is the cheapest of all the schemes, about 36% faster than the TENO5.










Shock-Bubble interaction

In the sixth example, we focus on the interaction between a shock wave and a bubble at Mach 6, as described in Hu (2010). The shock wave and the helium bubble are assumed to behave as ideal gases. The helium bubble is positioned at coordinates x = 0.25, y = 0 within a domain spanning [0, 1] × [-0.5, 0.5]. The initial conditions for the problem are as follows:

 
    (ρ,u,v,p)=
       (1.0, -3, 0, 1),                                pre-shocked  air,
       (216/41, (1645/286)-3, 0, 251/6)        post-shocked  air,
       (0.138, -3, 0, 1)                               helium  bubble.


The initial radius of the bubble is 0.15, while the initial position of the shock front is at x = 0.05. Inflow and outflow conditions are enforced at the left and right boundaries. Neumann conditions with zero gradients for all quantities are applied at the remaining boundaries. The simulation resolution is set to 800 × 800. Simulations are also carried out with shock-stable cLLF Riemann solver <cit.> to prevent the “carbuncle” phenomenon. Once again, from Fig. <ref>, one can observe that the MEG8-CC scheme resolves more small-scale vortices than TENO5 and MP6-CC schemes. Results obtained by the MEG8-C and MEG8-CC are very similar, but the MEG8-CC scheme is 12% faster than the MEG8-C scheme.













Inviscid Taylor-Green Vortex


The Taylor-Green vortex problem is a common benchmark for evaluating the performance of the adaptive central-upwind schemes. It is a three-dimensional inviscid test case with the following initial conditions for the simulation:


    [ ρ; u; v; w; p;   ]
    =
    [                                                 1;                                      sinxcosycosz;                                     -cosxsinycosz;                                                 0; 100 + ( cos(2z) + 2 ) ( cos(2x) + cos(2y)) - 2/16 ].


 The simulations are conducted until time t=10 on a grid size of 64^3, with a specific heat ratio of γ=5/3 on a periodic domain of size x,y,z ∈ [0,2π). The flow is incompressible since the mean pressure is significantly considerable. The study aims to evaluate the ability of the proposed schemes to preserve kinetic energy and the growth of enstrophy over time. Enstrophy is the integral of the square of the vorticity and is used as a measure of the scheme's ability to preserve vortical structures. The fourth-order compact scheme, Equation (<ref>), computes the velocity derivatives required for the enstrophy computation for all the schemes, and the HLLC Riemann solver is used for this test case.
 

The kinetic energy evolution of the numerical schemes MP6-CC, MP6-C, TENO5 and the linear fifth-order upwind scheme (U5) is shown in the Fig. <ref>, indicating that the MP6-CC scheme better preserves kinetic energy than the TENO5 scheme and even the linear U5 scheme. Furthermore, the MEG8-C and MEG8-CC schemes preserve the kinetic energy significantly better than all the other schemes, as shown in Fig. <ref>. There is little to no difference between the MEG8-C and MEG8-CC schemes, as this test case has no contact discontinuity or shockwaves. Figures <ref> and <ref> depict the computed enstrophy plots of the MP6-CC, MP6-C, TENO5, U5, MEG8-C and MEG8-CC schemes. The schemes with the proposed algorithm MP6-CC and MEG8-CC perform better than TENO and U5 schemes. There is a small difference between the MEG8-C and MEG-CC schemes for enstrophy values at late times, with MEG8-CC marginally better than MEG8-C. From Table <ref>, one can infer that the MP6-CC scheme is significantly cheaper than the TENO5 scheme and better preserves kinetic energy. While the computational costs of the TENO5 and MEG8-CC are almost the same, the MEG8-CC scheme outperforms the TENO5 scheme significantly in terms of both kinetic energy preservation and enstrophy values.














§ CONCLUSIONS

This work proposed a general adaptive central-upwind scheme for the simulations of compressible flows with discontinuities in the flow field. The approach significantly improves over the centralized gradient-based reconstruction proposed by Hoffmann, Chamarthi and Frankel <cit.> with reduced computational cost where the reconstruction is carried out for characteristic variables. In the current approach, conservative variables are reconstructed if no discontinuities are detected by the Durcos shock sensor and the MP criterion-based contact discontinuity sensor; otherwise, the reconstruction is the same as that of <cit.>. While it is well known in the literature that reconstruction of the characteristic variables will result in cleaner and oscillation-free results, the conversion of the conservative variables to characteristic variables will incur significant computational expense. The current approach overcomes this by reconstructing conservative variables in the regions with no discontinuities and characteristic variables in the areas of discontinuities, significantly reducing the computational cost and improving the results in some instances. The generalized approach also works with a standard fifth/sixth-order reconstruction scheme, MP6-CC. It has been shown that the MP6-CC scheme is 20-30% cheaper than the TENO5 scheme (which has the same stencil) and yet produces significantly superior results (judging by the vortical structures produced, which is the norm in the literature for assessing a low-dissipation scheme). The computational expense of the MEG8-CC scheme is the same as that of the TENO5 scheme, yet it produces superior results despite additional computations for the derivatives and Ducros sensor. 

The proposed approach did not produce spurious vortices for the periodic double shear layer test case with the proposed algorithm for all the schemes on coarse grids, unlike that of the TENO5 and the other existing schemes in the literature. It has been shown that the current algorithm (a combination of upwind and central schemes for appropriate variables in different directions )  can also prevent spurious vortices even with linear schemes (even the third/fourth scheme has prevented spurious vortices).




§ INSPIRATION AND ACKNOWLEDGEMENTS:


The work is inspired by the paper written by Van Leer <cit.>, where the author wrote - “Presumably, the advent of high-performance computing and promise of massively parallel computing has quelled any drive toward a systematic modernization of CFD algorithms. I personally believe, though, that the next round of gains in CFD will not come from hardware improvement but, once again, from method development”, coupled with the lack of resources. A.S. thanks his wife and son for checking the code and the manuscript.



elsarticle-num




