An improved low diffusion E-CUSP upwind scheme

Computers & Fluids 48 (2011) 214–220 Contents lists available at ScienceDirect Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d Technical note An improved low diffusion E-CUSP upwind scheme Ge-Cheng Zha ⇑, Yiqing Shen, Baoyuan Wang Dept. of Mechanical and Aerospace Engineering, University of Miami, Coral Gables, FL 33124, United States a r t i c l e i n f o Article history: Received 24 July 2010 Received in revised form 23 December 2010 Accepted 21 March 2011 Available online 2 April 2011 a b s t r a c t An improved low diffusion E-CUSP (LDE) scheme is presented. The E-CUSP scheme can capture crisp shock profile and exact contact surface. Several numerical cases are presented to demonstrate the accuracy and robustness of the new scheme. Ó 2011 Elsevier Ltd. All rights reserved. Keywords: Riemann solver Low diffusion 1. Introduction With the application of computational fluid dynamics becoming more and more popular, the demand on high accuracy and high efficiency CFD solutions also becomes stronger to satisfy the needs of the broad engineering problems. To achieve efficiency, accuracy and simplicity, many efforts have been made to develop upwind schemes only using scalar dissipation instead of matrix dissipation such as that of the Roe’s flux difference splitting (FDS) scheme [1]. The examples include AUSM family schemes of Liou [2–6], the Van Leer–Hänel scheme [7], Edwards’s LDFSS schemes [8,9], Jameson’s CUSP schemes and limiters [10–12], and the E-CUSP schemes developed by Zha et al. [13–16], etc. Zha and Hu suggested an E-CUSP schemes, which has low diffusion and can capture crisp shock wave profiles and exact contact discontinuities [16]. The scheme is consistent with the characteristic directions due to the nature of E-CUSP scheme. The scheme shows the highest stability for two shock tube tests problems compared with several other popularly used upwind schemes for the explicit Euler time marching scheme. The scheme also works well when extended to multi-dimensions [16]. However, the E-CUSP scheme of Zha–Hu may generate temperature oscillation near the computation boundary, in particular when the mesh is skewed. Zha was able to remove the temperature oscillation by introducing the total enthalpy in the smooth factor for the energy equation [17,18]. However, the scheme loses the capability to capture the exact contact surface due to the modification. The purpose of this paper is to present an improved low diffusion and efficient E-CUSP upwind scheme that is able to capture crisp shock profile and exact contact surface, and is smooth for ⇑ Corresponding author. Tel.: +1 305 284 3328; fax: +1 305 284 2580. E-mail address: gzha@miami.edu (G.-C. Zha). 0045-7930/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2011.03.012 multi-dimensional flow calculations. This paper modifies the Zha–Hu E-CUSP scheme [16] by using the Mach number splitting of Edwards’s LDFSS schemes [8] for the convective flux. The scheme is shown to be accurate, robust and efficient by the cases tested in this paper. 2. The numerical scheme 2.1. Governing equations To describe the new scheme, we will begin with the quasi-1D Euler equations in Cartesian coordinates for inviscid flow: @t U þ @xE  H ¼ 0 ð1Þ q where U ¼ SQ ; Q ¼ @ qu A; E ¼ SF, qe 0 0 qu 1 C B F ¼ @ qu2 þ p A; ðqe þ pÞu 1 0 1 0 dS B C H ¼ @pA dx 0 ð2Þ In above equations, q is the density, u is the velocity, p is the static pressure, e is the total energy per unit mass and S is the cross sectional area of the 1D duct. The following equation is also employed to relate pressure with the conservative variables:   1 p ¼ ðc  1Þ qe  qu2 2 ð3Þ where c is the specific heat ratio with the value of 1.4 for ideal gas. The finite volume method with the explicit Euler temporal integration is used to discretize the governing equations. It yields the following formulation at cell i: 215 G.-C. Zha et al. / Computers & Fluids 48 (2011) 214–220 DQ nþ1 i    H n i ¼ Dt C Eiþ1  Ei1 þ 2 2 Si ð4Þ where C = 1/(Dx Si), n is the time level index. A numerical scheme is needed to evaluate the interface flux: Eiþ1 ¼ SFiþ1 2 ð5Þ 2 1 1 ðML þ 1Þ2  ðM L  1Þ2 1 4 4 1 1 C  ðMÞ ¼  ðMR  1Þ2 þ ðM L  1Þ2 2 4 4 C þ ðMÞ ¼ 1 1 0 aR þ aL / aR þ aL aL þ aR =/ ¼ aR þ aL 1 ¼ ð14Þ 2 ð15Þ where, 0 qu C B C B F ¼ Fc þ Fp ¼ @ qu2 A þ @ p A pu qeu 0 ð13Þ It can be seen, if 1 = 2 = 1, the Mach number splitting is the Van Leer’s Mach number splitting [19]. The 1 and 2 are evaluated as [8]: 2.2. The low diffusion E-CUSP (LDE) scheme In [13,14,16], the characteristic analysis is given as the foundation to construct the E-CUSP scheme. The basic idea is to split the flux F to the convective flux Fc and the pressure flux Fp. That is: ð12Þ ð6Þ /¼ ðqa2 ÞR ðqa2 ÞL ð16Þ c The Jacobian matrix of vector F has the eigenvalues of the velocity (u, u, u), which represent the convective terms. The Jacobian matrix of vector Fp has the eigenvalues of the speed of sound (a, 0, a), which represents the acoustic waves propagating in each direction at subsonic regime. Based on the above separation of convective and wave terms in subsonic regime, Zha et al. [13–16] suggested to treat the convective term Fc in an upwind manner and to average the wave term Fp in both upwind and downwind direction with the weight of u ± a. Such a construction of the upwind scheme is consistent with the characteristic directions. In this paper, the modified E-CUSP scheme will have the same pressure splitting as that in [16,18]. The convective flux will be based on the Mach number splitting instead of the mass flux splitting. It can be written as: q C B c Fc ¼ u@ qu A ¼ aM qf qe 1 0 ð7Þ where a is the speed of sound, M is the Mach number, and 0 1 1 B C c f ¼ @uA e ð17Þ ð18Þ The pressure splitting coefficient borrowed from Van Leer [19] is used: P L;R ¼ 1 ðM L;R  1Þ2 ð2  M L;R Þ 4 ð19Þ For supersonic flows, if u > a; F1 ¼ FL ; If u < a, F1 ¼ FR 2 2 Above formulations can be written as a general form considering the Mach number range 1 < M < 1, which is given in the following section describing the extension of the scheme to 3D. 2.3. Extension to multi-dimensions  c c F1 ¼ a1 qL C þ ðMÞf L þ qR C  ðMÞf R 2 0 B þB @ 1 0 C Pp C  A 1 p u  a1 2 2 L R 1 ðaL þ aR Þ 2 uL uR ML ¼ ; MR ¼ a1 a1 2 In subsonic regime, if the flow is in the direction, 0 6 u 6 a, ð20Þ where 3 2 3 2 3 Other fluxes F, G, S, and T can be obtained following the symmetric rule. The stress s and heat flux q are, þ The interface speed of sound a1 , Mach number, uþ L , and uR are eval2 uated as: 2   @Q @E @F @G 1 @R @S @T þ þ þ ¼ þ þ @t @x @y @z Re @x @y @z 0 q qu 7 7 6 6 6 2 sxx 7 7 6 6 qu 7 6 qu þ p 7 7 7 7 6 6 6 7 7 6 6 6 Q ¼ 6 qv 7; E ¼ 6 quv 7; R ¼ 6 sxy 7 7 7 7 6 6 6 sxz 7 5 4 4 qw 5 4 quw 5 qe uk sxk  qx ðqe þ pÞu 0 B C þ P B þ @  p C A 1 p u þ a1 2 a1 ¼ 2.3.1. Governing equations The normalized Navier–Stokes equations governing compressible viscous flows in three-dimensions can be written in the Cartesian coordinate as: 2 1 2 0 2 1 1 ðML þ 1Þ2  ðM R þ 1Þ2 1 4 4 1 1 2  C ðMÞ ¼  ðMR  1Þ þ ðM R þ 1Þ2 2 4 4 C þ ðMÞ ¼ ð8Þ Similar to the treatment in [5,14,8], the interface speed of sound will be used in order to capture the contact discontinuities. The present E-CUSP scheme will adopt the Mach number splitting of Edwards [8], which is derived from the van Leer’s Mach number splitting formulation [19] and is made to capture the contact discontinuities. The interface flux of the present E-CUSP scheme is given as the following in the subsonic regime: 2 For most of the cases, 1 and 2 are on the order of 1 except at the shock discontinuities. The convective flux is hence primarily evaluated by the upwind flux. In subsonic regime, if the flow is in the direction, a 6 u 6 0, the Mach number splitting is constructed following the symmetric rule, which is: ð10Þ ð11Þ   @ui @uk 2 @uj  dik þ @xk @xi 3 @xj 1 l @T qj ¼ ðc  1ÞM 21 Pr @xj sik ¼ l  The repeated index k stands for the Einstein summation over x, y and z. The following equation is also employed to relate pressure with the conservative variables: 216 G.-C. Zha et al. / Computers & Fluids 48 (2011) 214–220 qe ¼ p 1 þ qðu2 þ v 2 þ w2 Þ The following relations borrowed from Edwards LDFFS scheme [8,9] to express the formulations from 1 < M < 1 are used, c1 2 In the above equations, q is the density, u, v, and w are the Cartesian velocity components in x, y and z directions, p is the static pressure, and e is the total energy per unit mass, l is the molecular viscosity, J is the transformation Jacobian, c, Re, M1, and Pr are the ratio of specific heat, Reynolds number, freestream Mach number, and Prandtl number, respectively. The Navier–Stokes equations Eq. (20) can be solved either by finite volume or finite difference method in the generalized coordinates (n, g, f). In this paper, when the 5th order weighted essentially non-oscillatory (WENO) scheme is used to discretize the inviscid fluxes, a finite difference method described in [20] is employed to save CPU time. When the first order scheme is used for the inviscid fluxes, a standard finite volume method is utilized. 2.3.2. The LDE Scheme in 3-dimensions[21] The extension of the LED scheme to 3D is straightforward and is basically the same as the 1D scheme. In generalized coordinate system, the flux E can be split as the following 0 qU B C B C B quU C B lx p C C B C B c p C B C E¼E þE ¼B B qv U C þ B ly p C C B C B @ qwU A @ lz p A qeU pU 0 1 0 1 ð21Þ where U is the contravariant velocities in n direction expressed as U ¼ lt þ lx u þ ly v þ lz w ð22Þ where l is the normal vector on n surfaces with its magnitude equal to the elemental surface area and pointing to the direction of increasing n. l¼ rn J dgdf ð23Þ nt dgdf J ð24Þ The velocity U is defined as U ¼ U  lt ¼ lx u þ ly v þ lz w ð25Þ In multi-dimensions, the formulations given in Section 2.2 can be written in a general form from subsonic to supersonic with 1 < M < 1 as the following. The convective term, Ec is evaluated as, 0 1 1 B C BuC B C c C Ec ¼ qU B B v C ¼ qUf ; B C @wA e 0 1 1 B C BuC B C C fc ¼ B Bv C B C @wA e ð26Þ let the contravariant speed of sound C be expressed as:  1 2 2 2 2 C ¼ a lx þ ly þ lz ð27Þ Ec1 ¼ C 1 qL C þ fLc þ qR C  fRc ð28Þ pffiffiffiffiffiffiffiffiffi where a ¼ cRT is the speed of sound. Then the convective flux at interface 12 is evaluated as: 2 2 where the subscripts L and R represent the left and right hand sides of the interface. The interface speed of sound is C1 ¼ 2 1 ðC L þ C R Þ 2 2 C  ¼ aR ð1 þ bR ÞM R  bR M R þ M1 2 UL UR ML ¼ ; MR ¼ C1 C1 2 2 1 aL;R ¼ ½1  signðML;R Þ 2 bL;R ¼ max½0; 1  intðjM L;R jÞ   qC 2 CR þ CLU C L þ C R U1  R M1 ¼ M1 ; M1 ¼ M1 ;U ¼ 2 C þ C 2 2 2 CR þ CL R L ðqC 2 ÞL ð30Þ þ M 1 ¼ bL dþ M L  bR d M þR 2 1 M L;R ¼  ðML;R  1Þ2 4   1 1 d ¼ 1  sign ðML þ MR Þ 2 2 The pressure flux, Ep is evaluated as the following 1 1 0 0 0  þ  C B C B B plx C B DL pL þ DR pR lx C C B C B þ  p B D p þ D p l B C E1 ¼ B ply C ¼ B  L L R R yC C 2  C B C B þ @ plz A @ DL pL þ DR pR lz A  þ  C 1 S L pL þ S R pR pU ð31Þ DL;R ¼ ½að1 þ bÞ  bP  L;R ð32Þ 0 2 where, The pressure splitting coefficient is: P L;R ¼ 1 ðM L;R  1Þ2 ð2  ML;R Þ 4 ð33Þ For the pressure term in the energy equation, the contravariant speed of sound C is consistent with U and is calculated as: lt stands for the grid moving velocity and is defined as lt ¼ C þ ¼ aþL ð1 þ bL ÞM L  bL M þL  Mþ1 ð29Þ C ¼ C  lt ð34Þ     ð1 þ bÞM S L;R ¼ ½a  bM L;R ð35Þ where M¼ U C ð36Þ  are evaluated based on M using the formulations  and b and the a given in Eq. (30). The use of U; C and M in the pressure term for the energy equation is to take into account of the grid speed so that the flux will transit from subsonic to supersonic smoothly. When the grid is stationary, lt ¼ 0; C ¼ C; U ¼ U. In this paper, the 5th order weighted essentially non-oscillatory scheme of Jiang and Shu [22] is used to achieve high order accuracy for the inviscid flux. The third-order TVD Runge–Kutta scheme [23] is used for time marching of the 1D unsteady problems. For multidimensional steady state problems, the implicit time marching scheme is used with the unfactored Gauss–Seidel line relaxation [24,25]. 3. Results and discussion According to Godunov [26], when there are discontinuities in the solutions, monotone behavior of a solution can not be assured with higher than first order scheme. Hence, for an upwind scheme to be used as a Riemann solver, it is essential to examine the performance of the scheme using first order accuracy. For the follow- 217 G.-C. Zha et al. / Computers & Fluids 48 (2011) 214–220 ing test cases, both the 1st order scheme and the 5th WENO scheme are used. For all the 2D flows calculated in this paper, the strategy of using increased value of e = 102 to minimize the numerical dissipation and to improve convergence [20] is employed. A comparison of the computational cost of the improved LDE scheme and the previous scheme [17] shows no distinguishable difference. 1.4 Present Scheme 1.3 CFL=0.9 CFL=0.09 Theory 1.2 T/TL 1.1 3.1. Shock tubes 3.1.1. The sod problem The Sod shock tube problem has the initial solution set to be at rest with a diaphragm located in the middle of the shock tube. The pressure on the left side of the diaphragm is 10 times higher than the pressure on the right side. At time level t = 0, the diaphragm breaks. A shock wave propagates to the right side of the tube. A contact surface follows the shock tube traveling toward the right side at a lower speed. An expansion wave propagates to the left side of the tube. Since the computation stops before the waves reach either end of the shock tube, the first order extrapolation boundary conditions are used at both ends of the shock tube. Fig. 1 is the first order accuracy results of the Sod shock tube [27] temperature profiles computed by the present scheme, Roe scheme and the analytical results. The maximum CFL numbers without generating oscillations for this shock tube computation is 1.0 for the present scheme and 0.95 for the Roe scheme. A uniform mesh of 201 points is used. Fig. 1 shows that the present scheme captures the shock profile using 3 grid points while the Roe scheme uses 4 grid points. The contact surface profiles resolved by both the schemes with 1st order accuracy are severely smeared and are virtually identical. Fig. 2 is the temperature profiles calculated by the 5th order WENO and 4th order Runge–Kutta method with the present ECUSP scheme at two different CFL numbers varied by 10 times. Fig. 2 indicates that the two profiles are basically identical, which means that the 5th order WENO scheme with 4th order Runge– Kutta method preserves the time accuracy well and the solutions are insensitive to the time steps. Obviously, the overall profiles are much more accurate than the 1st order scheme. The expansion front and tail are accurately resolved and the contact discontinuity is sharply captured. The shock profile is about the same as the 1st order result and is primarily determined by the Riemann solver. The slight oscillation at the contact surface may be removed by using the characteristic variables as shown in [28]. 1 0.9 0.8 0.7 0.6 0 5 10 X/L Fig. 2. Computed temperature distributions of the Sod 1D shock tube using 5th order WENO schemes with different time steps. 3.1.2. Slowly moving contact surface This is a shock tube case used in [3] to demonstrate the improvement of the present scheme over the previous scheme [17] to capture the contact surface. The initial conditions are [q, u, p]L = [0.125, 0.112, 1.0], [q, u, p]R = [10.0, 0.112, 1.0]. All the results are first order accuracy. A uniform mesh of 201 points is used. Fig. 3 shows that both the present scheme and the Roe scheme can resolve the exact contact surface as they are designed. The results of those schemes are at time level 0.01. The velocity is constant and the density discontinuity is monotone. The results calculated using the previous E-CUSP scheme [17] show that it can not capture the constant velocity of the contact surface with large velocity oscillation. The new E-CUSP scheme also has far higher CFL number than the Roe schemes with the CFL value of 1.0. The Roe scheme has the maximum CFL = 0.3 without generating oscillations. Fig. 3 indicates that as a Riemann solver, both the present scheme and the Roe scheme can capture the exact contact discontinuity. Unfortunately, for the 5th order WENO scheme with the 4th order Runge–Kutta method, it can not achieve the constant velocity and pressure across the contact discontinuity with no oscillation for both the new E-CUSP and Roe scheme. 1.4 0.2 10 present CFL=1.0 Roe CFL =0.95 Theory 9 0.175 CFL=1.0 Previous 8 1.2 Density CFL=0.3 Roe Scheme 7 1.1 0.125 Density T/TL 0.15 CFL=0.3 1 0.9 6 0.1 5 4 velocity Velocity 1.3 Present Scheme 0.075 3 0.8 0.05 2 1 0.7 0.025 0 0.6 10 5 X/L Fig. 1. Computed temperature distributions of the Sod 1D shock tube using 1st order schemes. 5 10 0 X/L Fig. 3. Computed density and velocity profiles of the slowly moving contact discontinuity using different schemes. 218 G.-C. Zha et al. / Computers & Fluids 48 (2011) 214–220 3.1.3. Shock interaction with an entropy wave This test case is taken from Ref. [29], which represents a Mach 3 shock wave interacting with a sine entropy wave. The initial condition is 3 Present Theory Roe Scheme 2.5 ð1 þ esinð5xÞ; 0; 1Þ x P 4 where e = 0.2 is used. Fig. 4 is the density distribution calculated by the new scheme and the Roe scheme with the 5th-WENO scheme and the 4th order Runge–Kutta method based on the mesh size of 400 points. Both schemes achieve virtually the same profile. Compared with the accurate solutions calculated using 8000 mesh points, the shock profile and the front of the entropy wave are resolved well, but some of the wave amplitudes are under-predicted due to the numerical dissipation. Mach Number ð3:857143; 2:629369; 10:3333Þ x < 4 ðq; u; pÞ ¼ 2 1.5 1 0.5 0 2 4 6 8 X/h 3.2. Entropy condition This case is to test if a scheme violates the entropy condition by allowing expansion shocks. The test case is a simple quasi-1D converging–diverging transonic nozzle [14,15]. The correct solution should be a smooth flow from subsonic to supersonic with no shock. However, for an upwind scheme which does not satisfy the entropy condition, an expansion shock may be produced. For the subsonic boundary conditions at the entrance, the velocity is extrapolated from the inner domain and the other variables are determined by the total temperature and total pressure. For supersonic exit boundary conditions, all the variables are extrapolated from inside of the nozzle. The analytical solution was used as the initial flow field. Explicit Euler time marching scheme was used to seek the steady state solutions. All the schemes use first order differencing. Fig. 5 is the comparison of the analytical and computed Mach number distributions with 201 mesh points using the present scheme and the Roe scheme. The analytical solution is smooth throughout the nozzle and reaches the sonic speed at the throat (the minimum area of the nozzle, located at X/h = 4.22). It is seen that the Roe scheme generates a strong expansion shock at the nozzle throat. Both schemes can converge to machine zero (12 order of magnitude) with CFL = 0.95 even with the expansion shock waves. Similar to the original Zha CUSP scheme [16,17], the present scheme does not have an expansion shock wave at the sonic point, but is not smooth and has a little glitch. As reported in Fig. 5. Computed Mach number distributions for the quasi-1D converging–diverging nozzle using 1st order schemes. [14], when the mesh is refined, the amplitude of the expansion shock waves is reduced. When the 2nd order scheme is used, the expansion shock waves disappear. 3.3. Supersonic Laminar wall boundary layer To examine the numerical dissipation of the new scheme, a laminar supersonic boundary layer on an adiabatic flat plate is calculated using first order accuracy. The incoming Mach number is 2.0. The Reynolds number based on the length of the flat plate is 40,000. Adiabatic no-slip wall boundary conditions are used on the wall. The Prandtl number of 1.0 is used in order to compare the numerical solutions with the analytical solution of compressible Blasius solution [30]. The baseline mesh size is 41  31 in the direction along the plate and normal to the plate respectively. The height of the computational domain is 0.82 times of the flat plate length. There are 13 points within the boundary layer. Fig. 6 is the comparison between the computed velocity profiles and the Blasius solution. The solutions of the present scheme and the Roe scheme all agree precisely with the analytical solution. Fig. 7 is the comparison between the computed temperature profiles and the Blasius solution. Again, the present scheme and 10 4.5 Present Blasius Roe Scheme 8 4 3.5 eta Density 6 3 present scheme Roe scheme accurate 2.5 4 2 2 1.5 1 0 -4 -2 0 2 4 X/L Fig. 4. Computed density profiles of the shock interaction with an entropy wave using 5th order WENO schemes. 0 0.2 0.4 0.6 0.8 1 U/Uf Fig. 6. Computed velocity profiles of the laminar boundary layer using 1st order schemes. 219 G.-C. Zha et al. / Computers & Fluids 48 (2011) 214–220 10 -3 present Blasius Roe Scheme 9 8 -4 -5 -6 7 -7 Cp eta 6 5 -8 Roe’s scheme LDE -9 4 -10 3 -11 2 -12 1 -13 0 1 1.2 1.4 1.6 5 1.8 Fig. 7. Computed temperature profiles of the laminar boundary layer using 1st order schemes. Roe scheme accurately predict the temperature profiles and the computed solutions basically go through the analytical solution. The 1st order accuracy results presented here indicate that the present E-CUSP scheme has low diffusion as the Roe scheme. For the upwind schemes with high diffusion such as the Van Leer scheme, they cannot resolve the boundary layer profiles accurately [16,18]. 10 x T/Tf Fig. 9. Computed nozzle surface pressure coefficient using the 5th order WENO schemes. 175  50. The grid is clustered near the wall. The inlet Mach number is 0.22. The CFL number of 10 is used. Fig. 8 is the computed pressure contours showing the oblique shock and its reflection on the wall. Fig. 9 is the computed nozzle surface pressure coefficient using the 5th order WENO schemes with the LDE scheme and the Roe scheme, which shows that the results computed by the two schemes are indistinguishable. 4. Conclusions 3.4. Transonic converging–diverging nozzle To examine the performance of the new scheme with 5th order WENO scheme in two-dimensional flow and the capability to capture the shock waves which do not align with the mesh lines, a transonic converging–diverging nozzle is calculated as inviscid flow. The nozzle was designed and tested at NASA and was named as Nozzle A1 [31]. The solution is steady state and the time marching is achieved by the unfactored implicit Gauss–Seidel iteration. The nozzle is symmetric about the centerline. Hence only upper half of the nozzle is calculated. The upper boundary uses the slip wall boundary conditions and the lower boundary of the center line uses the symmetric boundary conditions. The mesh size is The improved LDE scheme can capture crisp shock profiles and exact contact surfaces. The computational cost of the improved scheme is basically the same as the previous E-CUSP scheme. The new scheme is shown to have low diffusion and is able to resolve wall boundary layer accurately. For a 1D transonic nozzle, the new LDE scheme does not generates any expansion shock wave as the Roe scheme does. Several 1D shock tube flows and 2D transonic flows are calculated using the 5th order WENO scheme with the LDE scheme. The results of the LDE scheme are in good agreement with experiment and theoretical results. Acknowledgment This work is partially supported by AFOSR Grant FA9550-06-10198 monitored by Dr. Fariba Fahroo, by ARO Grant 56827-RT-ISP monitored by and John D. Schmisseur and Lt Col. R. Jefferies at AFOSR and Peggy A. Lacewell at ARO, and by Miami WindTR at University of Miami. 10 9 8 7 References y 6 5 4 3 2 1 0 0 5 10 x Fig. 8. Computed nozzle pressure contours using the 5th WENO and the new E-CUSP scheme. [1] Roe P. Approximate riemann solvers, parameter vectors, and difference schemes. J Comput Phys 1981;43:357–72. [2] Liou M-S, Steffen CJ. A new flux splitting scheme. J Comput Phys 1993;107:1–23. [3] Wada Y, Liou MS. An accurate and robust splitting scheme for shock and contact discontinuities. AIAA Paper 94-0083; 1994. [4] Liou MS. Progress towards an improved CFD methods: AUSM+.’’ AIAA Paper 951701-CP; June, 1995. [5] Liou M-S. A sequel to AUSM: AUSM+. J Comput Phys 1996;129:364–82. [6] Liou M-S. Ten years in the making-AUSM-family. AIAA 2001:2001–521. [7] Hänel D, Schwane R, Seider G. On the accuracy of upwind schemes for the solution of the Navier–Stokes eqautions. AIAA paper 87-1105 CP; 1987. [8] Edwards JR. A low-diffusion flux-splitting scheme for Navier–Stokes calculations. AIAA Paper 95-1703-CP; June, 1995. [9] Edwards JR. A low-diffusion flux-splitting scheme for Navier–Stokes calculations. Comput Fluids 1997;6:635–59. 220 G.-C. Zha et al. / Computers & Fluids 48 (2011) 214–220 [10] Jameson A. Analysis and design of numerical schemes for gas dynamics I: artificial diffusion, upwind biasing, limiters and their effect on accuracy and multigrid convergence in transonic and hypersonic flow. AIAA Paper 93-3359; July, 1993. [11] Jameson A. Analysis and design of numerical schemes for gas dynamics I: artificial diffusion, upwind biasing, limiters and their effect on accuracy and multigrid convergence in transonic and hypersonic flow. J Comput Fluid Dyn 1995;4:171–218. [12] Jameson A. Analysis and design of numerical schemes for gas dynamics II: artificial diffusion and discrete shock structure. J Comput Fluid Dyn 1995;5:1–38. [13] Zha G-C, Bilgen E. Numerical solutions of euler equations by using a new flux vector splitting scheme. Int J Numer Methods Fluids 1993;17:115–44. [14] Zha G-C. Numerical tests of upwind scheme performance for entropy condition. AIAA J 1999;37:1005–7. [15] Zha GC. Comparative study of upwind scheme performance for entropy condition and discontinuities. AIAA Paper 99-CP-3348; June 28–July 1, 1999. [16] Zha G-C, Hu Z-J. Calculation of transonic internal flows using an efficient high resolution upwind scheme. AIAA J 2004;42(2):205–14. [17] Zha G-C. Low diffusion efficient upwind scheme. AIAA J 2005;43. [18] Zha GC. A low diffusion e-cusp upwind scheme for transonic flows. AIAA Paper 2004-2707, 34th AIAA Fluid Dynamics Conference; June 28–July 1 2004. [19] Van Leer B. Flux–vector splitting for the euler equations. Lect Note Phys 1982;170:507–12. [20] Shen Y-Q, Zha G-C, Wang B-Y. Improvement of stability and accuracy of implicit WENO scheme. AIAA J 2009;47(2):331–44. [21] Zha GC, Shen Y, Wang B. Calculation of transonic flows using WENO method with a low diffusion E-CUSP upwind scheme. AIAA Paper 2008-0745, 46th AIAA Aerospace Sciences Meeting, Reno, NV, January 2008. [22] Jiang GS, Shu CW. Efficient implementation of weighted ENO schemes. J Comput Phys 1996;126:202–28. [23] Shu CW. Essentially Non-Oscillatory and Weighted Essentially Schemes for Hyperbolic Conservation Laws. [24] Hu Z-J, Zha G-C. Calculations of 3D compressible using an efficient low diffusion upwind scheme. Int J Numer Methods Fluids 2004;47:253–69. [25] Shen YQ, Zha GC. A comparison study of Gauss–Seidel iteration methods for internal and external flows. AIAA Paper 2007-4332; 2007. [26] Godunov S. Finite-difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics. Mat Sb 1959;47:271–306. [27] Sod G. A survey of several finite difference methods for systems of nonlinear hyperbolic conversation laws. J Comput Phys 1978;27:1–31. [28] Shen YQ, Zha GZ. Generalized finite compact difference scheme for shock/ complex flowfield interaction. J Comput Phys. doi:10.1016/j.jcp.2011.01.039. [29] Shu C-W, Osher O. Efficient implementation of essentially non-oscillatory shock capturing schemes. J Comput Phys 1988;77:439–71. [30] White FM. Viscous fluid flow. 2nd ed. McGraw-Hill Inc; 1991. [31] Mason ML, Putnam LE. The effect of throat contouring on two-dimensional converging–diverging nozzles at static conditions. NASA Technical Paper 1704; 1980.