90.dvi - PDF 免费下载

Commun. Comput. Phys. do: 10.408/ccp.041109.160910a Vol. 10, No. 1, pp. 90-119 July 011 Performance of Low-Dsspaton Euler Fluxes and Precondtoned LU-SGS at Low Speeds Kech Ktamura 1,, E Shma 1, Kechro Fumoto 1, and Z. J. Wang 3 1 JAXA s Engneerng Dgtal Innovaton (JEDI) Center, Japan Aerospace Exploraton Agency, 3-1-1 Yoshnoda, Chuuou, Sagamhara, Kanagawa, 5-510, Japan. JAXA s Engneerng Dgtal Innovaton (JEDI) Center, Japan Aerospace Exploraton Agency, -1-1 Sengen, Tsukuba, Ibarak, 305-8505, Japan. 3 Department of Aerospace Engneerng, Iowa State Unversty, 71 Howe Hall, Ames, IA 50011, USA. Receved 4 November 009; Accepted (n revsed verson) 16 September 010 Avalable onlne 7 March 011 Abstract. In low speed flow computatons, compressble fnte-volume solvers are known to a) fal to converge n acceptable tme and b) reach unphyscal solutons. These problems are known to be cured by A) precondtonng on the tme-dervatve term, and B) control of numercal dsspaton, respectvely. There have been several methods of A) and B) proposed separately. However, t s unclear whch combnaton s the most accurate, robust, and effcent for low speed flows. We carred out a comparatve study of several well-known or recently-developed low-dsspaton Euler fluxes coupled wth a precondtoned LU-SGS (Lower-Upper Symmetrc Gauss-Sedel) mplct tme ntegraton scheme to compute steady flows. Through a seres of numercal experments, accurate, effcent, and robust methods are suggested for low speed flow computatons. AMS subect classfcatons: 65Z05, 76M1, 76N99 Key words: All-speed scheme, low-dsspaton, precondtonng, LU-SGS. 1 Introducton In recent years, compressble fnte-volume methods (FVMs) have been used n a wde spectrum of flow regmes, ncludng low speed flows n whch compressblty plays no Correspondng author. Emal addresses: ktamura.kech@axa.p (K. Ktamura), shma.e@axa.p (E. Shma),fumoto.kechro@axa.p (K. Fumoto),zw@astate.edu (Z. J. Wang) http://www.global-sc.com/ 90 c 011 Global-Scence Press

K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 91 sgnfcant role. The applcaton of compressble flow solvers to low speeds has been motvated by the fact that users need only slght modfcatons to the exstng (compressble) codes for computatons of such low speed flows, and that ths extenson has the followng potental applcatons of engneerng nterests: 1. Analyss of flows nvolvng both low speeds (M < 0.1) and hgh speeds (M 10 or even 100), e.g., a cavtatng flow n a rocket engne [1,].. Aeroacoustc analyss n low speed flows [3]. When appled to low speed flow computatons, however, compressble solvers are known to a) fal to converge n acceptable tme (stffness problem), and b) reach unphyscal solutons. These problems are known to be cured by A) precondtonng on the tmedervatve term so that acoustc wave speed s properly scaled, and B) control of dsspaton n numercal fluxes, respectvely. There have been several methods of A) [4, 5] and B) [6 9] proposed separately. However, t s unclear whch combnaton s the most accurate, robust, and effcent n low speed flows. It s dffcult to prove ths mathematcally because, for nstance, the amount of dsspaton added to the computaton s dependent not only on the adopted methods, but on the computatonal grd, flow condtons, and so forth. If a combnaton of methods A) and B) has nsuffcent dsspaton for the gven condtons, the calculaton wll suffer from numercal oscllaton/nstablty, and may eventually dverge. If the method s too dsspatve, on the other hand, ts accuracy s sgnfcantly lost. Therefore, n the present paper, we pursue an expermental approach by performng a comparatve study of dfferent methods of A) along wth B) for dfferent grds and dfferent flow condtons of low speeds. We wll pay partcular attenton to several well-known or recently-developed low-dsspaton Euler fluxes coupled wth a precondtoned LU- SGS (Lower-Upper Symmetrc Gauss-Sedel) mplct scheme [10, 11] n the framework of steady flows. Smlar comparsons have already been conducted by others (n [1], for example), but ther dscussons were lmted to only a few methods/cases and lacked concrete conclusons. In ths study, through an extensve seres of numercal experments, accurate, effcent, and robust methods among 14 dfferent approaches wll be suggested for low speed flow computatons. The paper s organzed as follows: n Secton, numercal methods and flow condtons adopted here wll be descrbed. Then, n Secton 3, numercal results and dscussons wll be presented from a vscous, moderate speed case (Case 1: M =0.5, Re =5,000 n 3.1), nvscd, low speeds cases (Cases A-C: M = 0.1 0.001 n 3.), and a vscous, low speed case (Case 3: M =0.01, Re =,000 n 3.3). CFL effects wll also be dscussed n 3.4. These computatons wll be conducted wth global tme steppng so that dscussons theren could be appled (or at least referenced) to unsteady flow computatons wth the use of dual-tme steppng n whch temporal convergence s attaned n each tme step [4] (not actually covered n ths work, though). On the other hand, t s natural to use local tme steppng technque f one s nterested n steady solutons. Thus, we wll address the local tme steppng ssue as a separate nvestgaton n 3.5. Features of each method wll be summarzed n 3.6, and Secton 4 wll conclude the present artcle.

9 K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 Numercal methods and flow condtons.1 Governng equatons The governng equatons are the compressble Naver-Stokes equatons as follows, ncludng the precondtonng matrx Γ of Wess and Smth [4]. Γ s gven n the Appendx, and t s smply elmnated n the non-precondtoned form. In general, conservatve varables are commonly used n compressble flow solvers and occasonally t s not straghtforward to change varables used n the whole code (our code s not an excepton). Thus, we kept usng conservatve varables Q as dependent varables [8, 13], nstead of usng prmtve ones q, although some drawbacks of usng conservatve varables were ponted out n [14]. Γ Q t + F k = Fv k, x k x k ρ ρu k Q= ρu l, F k = ρu l u k + pδ lk ρe ρu k H ( ul τ lk = µ + u ) k x k x l 3 µ u n δ lk, x n, Fv k = 0 τ lk u m τ mk +κ T x k (.1a), (.1b) (.1c) where ρ s the densty, u velocty components n Cartesan coordnates, E total energy, p pressure, H total enthalpy,.e., H = E+ p/ρ, and T temperature. The workng gas s ar approxmated by the calorcally perfect gas model wth the specfc heat rato γ=1.4. The Prandtl number s Pr = 0.7. The molecular vscosty µ and thermal conductvty κ are related as κ= c p µ/pr, where c p s specfc heat at constant pressure. Fgure 1: Schematc of cell geometrc propertes. Eq. (.1) s solved wth a fnte-volume code, and can be wrtten n the delta form as: V t Q +Γ 1 Σ (F, Fv, )S, =0, (.) where V stands for the volume of the cell, t the (local) tme step, Q change of conservatve varables n tme, F, and Fv, the nvscd (Euler) and vscous fluxes through

K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 93 the cell-nterface S, (whch separates the cell and ts neghbor cell ), respectvely (see Fg. 1). Detals are explaned below.. Numercal methods The computatonal code employed here s LS-FLOW : JAXA s n-house, unstructured, compressble Naver-Stokes solver for arbtrary polygons. LS-FLOW has many optons for spatal reconstructon and temporal evoluton. Included n Table 1 are only the methods adopted for the present study. The second order of spatal accuracy s guaranteed [18]. The Euler fluxes and the mplct schemes are summarzed n Table. Vscous fluxes are computed by usng Wang s second-order method [17]. Formulaton of each Euler flux s brefly ntroduced below, followed by that of tme evoluton methods (for detals, see the orgnal lterature). Governng Equatons Spatal Dscretzaton Table 1: Numercal methods. Compressble Euler/Naver-Stokes Equatons Cell-centered FVM Green-Gauss Method [15, 16] (wthout slope lmter) Gradents Spatal Reconstructon Invscd Term see Euler Fluxes n Table Vscous Term Wang [17] Temporal Evoluton see Implct Schemes n Table Table : Euler fluxes and mplct schemes. Baselne Low-Dsspaton/ Precondtoned (Specfed reference values) Roe [19] P-Roe [4] (M co ), A-Roe [9] (ρ u ) Euler Fluxes AUSM + [0] AUSM + -up [6] (M ) SHUS [1] SLAU [8] (No reference values requred) Implct Schemes LU-SGS [10] plu-sgs [4] (M co )..1 Euler fluxes Invscd numercal fluxes at cell-nterfaces F 1/ are calculated by one of the followng Euler fluxes. 1) Roe [19], Precondtoned Roe (P-Roe) [4], and All-Speed-Roe (A-Roe) [9]: Usng the dfference of varables ()=() R () L, and the Roe-averaged [19] values (ˆ), the Roe flux s expressed n the followng form of Lu and Vnokur []: where F 1 = 1 ( FL +F R ˆR ˆΛ ˆL Q ), (.3a) ˆR ˆΛ ˆL Q= ˆ λ 1 Q+δ 1 ˆQ +δ N, (.3b) ˆQ =(1 u v w H) T, N=(0 n x n y n z V n ) T, Γ=dag(λ 1, λ,, λ 5 ), (.4a)

94 K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 δ 1 = λ+ ( p/ĉ)+λ ˆρ V n, δ = λ + ˆρ V n +λ ( p ), ĉ ĉ (.4b) λ + = λˆ 1 + λˆ + λˆ 3, λ = λˆ λˆ 3, (.4c) λ 1,4,5 =V n, λ =V n +c, λ 3 =V n c, (.4d) V n = un x +vn y +wn z. The mass flux (frst row of F 1/ ), for example, s wrtten as follows: ˆm 1 = 1 ( ρ L V nl +ρ R V nr ˆV n ρ ˆM+1 + ˆM 1 ˆM p ĉ (.4e) ˆM+1 + ˆM 1 ˆρ V n ), (.5) The Roe s Remann solver s one of the most wdely-used numercal fluxes, but ths flux s known to suffer from the carbuncle phenomenon [3, 4] or an expanson shock at hgh speeds, and as demonstrated later, unphyscal oscllatons at low speeds. To overcome the defect n low speed flow computatons, Wess and Smth derved a verson of Roe flux for a precondtoned system, called Precondtoned Roe (P-Roe) [4], by multplyng nverse precondtonng matrx Γ 1 to the numercal dsspaton term as follows: where wth F 1 = 1 ( FL +F R Γ 1 ˆR ˆΛ ˆL Q ), (.6a) Γ 1 ˆR ˆΛ ˆL Q= ˆλ 1 Q+δ 1 ˆ Q +δ N, (.6b) ˆ Q =(1 u v w H) H, N=(0 n x n y n z V n ) T, Γ 1 Λ=dag(λ 1, λ, λ 3, λ 4, λ 5 ), (.7a) [ δ 1= λ ˆ + ˆ [ δ = ˆ λ 3 λ + λ ˆ 3 ε V n 1 ε ˆ V n + 1 ε ˆ λ 1,4,5 =V n, V,3=V n±c = 1 λ λ ˆ 3 Vn ] p ĉ εĉ + λ ˆ λ ˆ 3 λ λ ˆ 3 Vn ] ĉ ˆρ V n + λ ˆ λ ˆ 3 p ĉ, { (1+ε) V n ± ˆρ V n ĉ, (.7b) (.7c) } (ε 1) V n +4εc, (.7d) ε=mn ( 1, max(km,m co) ), (.8) where K s taken as 1.0, and M co s cutoff Mach number whch s set as freestream Mach number M n the present study, leadng to ε = mn(1,max(m,m )). The above expresson s borrowed from a precondtonng matrx Γ (see Eq. (A.) n Appendx) so that

K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 95 the resultant flux formula has all the egenvalues scaled from the order of the speed of sound to that of the local flud velocty at low speeds. Note that at supersonc speeds (max(m,m ) > 1.0), ε s unty and orgnal Roe flux s recovered; otherwse, even at a moderate speed (e.g., M = 0.5), the egenvalues are no longer the same as those n Eq. (.4d) and the resultant flux s expected to behave n a dfferent manner as unprecondtoned Roe. All-Speed-Roe (A-Roe), whch was developed recently by L and Gu [9], modfed the Roe flux by ntroducng the swtchng functon f(m) for all speeds as follows: F 1 = 1 ( ) f(m)(f L +F R )+[1 f(m)]fc press ˆR ˆΛ A Roe ˆL Q, (.9a) ( 4+(1 M ) M ) 1+M,1, (.9b) f(m)=mn Fc press = U c(φ L +Φ R ) + P L+P R, (.9c) U c = (V nl+v nr ) c (ρ u ) (p R p L ), (.9d) Φ=(ρ, ρu, ρv, ρw, ρh) T, P=(0, pn x, pn y, pn z, 0) T, (.9e) λ A Roe 1,4,5 =V n, λ A Roe =V n + f(m)c, λ A Roe 3 =V n f(m)c, (.9f) where Fc press s a pressure stablzaton term wth c = 0.05 and ρ u = ρ u (accordng to the orgnal paper [9], these should be the maxmum values n the whole computatonal doman; however, they are smply set to be freestream values here, snce contnuty equaton ustfes ρ u = ρ u = const n ths work). Ths scheme does not rely on cutoff Mach number M co, whch s typcally borrowed from precondtonng matrx Γ (see Eq. (A.) n Appendx) and s ncluded n some other all speed schemes (such as Eq. (.8) n P-Roe [4], as explaned above), though reference values ρ u should be specfed. Furthermore, as P-Roe, ths flux s supposed to behave dfferently compared wth orgnal Roe even at a moderate speed. ) AUSM + [0] and AUSM + -up [6]: AUSM-famly schemes [6 8,0,1] are another set of wdely-used fluxes featurng smplcty and relatve robustness aganst shock-related anomales (e.g., carbuncle phenomenon) [3, 4]. Among AUSM-famly, we ntroduce two representatve methods,.e., AUSM + and ts all-speed extenson, AUSM + -up. Formulaton of AUSM + s gven as: F 1 = ṁ+ ṁ Ψ + + ṁ ṁ Ψ + pn, (.10a) Ψ=(1, u, v, w, H) T, N=(0, n x, n y, n z, 0) T, (.10b) { ρl, f M 1 >0, ṁ 1 = M 1 c 1 (.10c) ρ R, otherwse, M 1 = f + M,L β= 1 + f 8 M,R β= 1, p= f + 8 p,l α= 3 p L + f 16 p,r α= 3 p R, (.10d) 16

96 K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 where { f ± 1 (M± M ), f M 1, M β = ± 1 4 (M±1) ±β(m 1), otherwse, { f p ± 1 (1±sgn(M)), f M 1, α = 1 4 (M±1) ( M)±αM(M 1), otherwse, c 1 =mn( c L, c R ), c= c max(c, u ), (.10e) (.10f) c = (γ 1) H. (.11) (γ+1) Ths scheme was extended for all speeds as AUSM + -up, wth ntroducton of addtonal user-specfed parameters: where F 1 = ṁ+ ṁ Ψ + + ṁ ṁ Ψ + pn, (.1a) Ψ=(1, u, v, w, H) T, N=(0, n x, n y, n z, 0) T, (.1b) { ρl, f M 1 >0, ṁ 1 = M 1 c 1 (.1c) ρ R, otherwse, M 1 = f + M,L β= 1 + f 8 M,R β= 1 +M p, p= f + 8 p,l α p L + f p,r α p R + p u, (.1d) c 1 =mn( c L, c R ), c L = c max(c,v nl ), c R = c max(c, V nr ), (.13a) M p = K p max(1 σ M,0) p R p L f α ρ 1 c, ρ 1 = ρ L+ρ R, (.13b) 1 wth p u = K u f + pl f pr (ρ L+ρ R )( f α c 1)(V nr V nl ), (.13c) K p =0.5, K u =0.75, σ=1.0, (.14a) M = V nl +V nr c 1, α= 3 16 ( 4+5 f α), (.14b) f α (M o )= M o ( M 0 ), M o =mn ( 1,max( M,M ) ). (.14c) Ths scheme also excludes cutoff Mach number M co, though freestream Mach number M s requred. In ths paper, as n [6, 9], cutoff Mach number M co and freestream Mach number M are used n dfferent meanngs, because M co s not necessarly equal to M, but arbtrary chosen from the values of the order of M.

K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 97 3) SHUS [1] and SLAU [8]: SHUS (Smple Hgh-resoluton Upwnd Scheme) s one of AUSM-famly schemes, whch replaced mass flux of AUSM + wth that of Roe (Eq. (.5) wth the use of arthmetc averaged values rather than Roe averaged ones. Ths scheme acheved accuracy of Roe flux whle keepng the robustness of AUSM + aganst shock anomales. = ṁ+ ṁ Ψ + + ṁ ṁ Ψ + pn, (.15a) F 1 Ψ=(1, u, v, w, H) T, N=(0, n x, n y, n z, 0) T. (.15b) The mass flux functon of SHUS s gven as: ṁ= 1 { (ρv n ) L +(ρv n ) R V nl+v nr ρ M+1 M 1 ρ V n M+1 + M 1 M } p, (.16) c where The pressure term s: ()=() R () L, ( )= () L+() R. (.17) p= f + pl α=0 p L + f pr α=0 p R, (.18) SHUS was further developed to gve more relable solutons both at low and hgh speeds. The latest verson s named SLAU (Smple Low-dsspaton AUSM): F 1 = ṁ+ ṁ Ψ + + ṁ ṁ Ψ + pn, (.19a) Ψ=(1, u, v, w, H) T, N=(0, n x, n y, n z, 0) T. (.19b) The mass flux functon of SLAU s: ṁ= 1 { ρ L (V nl + V n + )+ρ L (V nl V n ) χ c }, p V n = ρ L V nl +ρ R V nr, (.0a) ρ L +ρ R V n + =(1 g) V n +g V nl, V n =(1 g) V n +g V nr, (.0b) g= max [ mn(m L,0), 1 ] mn [ max(m R,0),1 ] [0,1], (.0c) and the pressure flux s p= p L+p R + f + pl α=0 f pr α=0 (p L p R )+(1 χ) ( f + pl α=0 + f pr α=0 1 ) p L + p R, (.1a) χ=(1 M), (.1b) ( M=mn 1.0, 1 c u ) L +v L +w L +u R +v R +w R, (.1c) M= V n c = un x+vn y +wn z. (.1d) c

98 K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 SLAU needs no cutoff Mach number M co or freestream Mach number M. To the best of the authors knowledge, ths flux s the only method among all speed schemes whch s totally free from restrctons of specfyng reference values. Ths property s desrable for computatons of flows nvolvng no unform flow, such as turbopump nternal flows [1]... Tme evoluton methods 1) LU-SGS and plu-sgs (precondtoned LU-SGS) Implct Schemes: Tme ntegraton s conducted by usng LU-SGS mplct method or ts precondtoned verson, precondtoned LU-SGS [4], whch s referred to as plu-sgs for brevty here. Its formulaton starts from Eq. (.), expressed wth tme step ndex n ncluded: where and Res s the rght-hand sde resdual, V Q n t +Γ 1 ( F n, Fvn, )S,= Γ 1 Res n, (.) () n =() n+1 () n, (.3) Res n = (F n, Fvn, )S,. (.4) Agan, n the case wthout precondtonng, Γ 1 s smply dropped. Then, Eq. (.) s rewrtten n the form of Gauss-Sedel (GS) teratve method by decomposton of new (updated) and old (non-updated) values where [ Q n+1 =D n 1 +Γ 1 Γ 1 U pper Lower S, A +, S, A + new/old, Q new/old new/old Q new/old ] Γ 1 Res n, (.5) A, = F, Q Fv, Q, (.6) s flux Jacoban from cell to cell through the cell-nterface S,. A + has only the postve components of the egenvectors. The dagonal matrx D s gven as D = V I+Γ 1 t S, A +,. (.7)

K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 99 Specfcally, n LU-SGS (Lower Upper Symmetrc Gauss-Sedel), Eq. (.5) s further rewrtten as [ ] Q =D 1 Γ 1 S, A +, Q Γ 1 Res n, (.8a) [ Q n+1 =D 1 Lower Γ 1 = Q +D 1 Lower Γ 1 S, A +, Q +Γ 1 U pper U pper S, A +, n+1 Q n+1 ] Γ 1 Res n S, A + n+1, Q n+1. (.8b) Then, A + s approxmated as the followng as n Jameson and Turkel s LU-SGS [10] (ths verson s commonly referred to as LU-SGS ). where A +, A,+σ, I, (.9) σ, = σ, (A, )= V n, +c + (µ +µ ) (ρ +ρ ) h,, (.30) where h, s dstance between cell-centers of and. In plu-sgs, the spectral radus σ, s scaled as σ,, thus, Γ 1 A +, Γ 1 A, +σ, I, (.31) where σ, = σ (Γ 1 A, )= 1 { (1+ε) V n, + (ε 1) Vn, +4εc }+ (µ +µ ). (.3) (ρ +ρ ) h, The precondtonng coeffcent ε, whch should be of the order of M, appears n the precondtonng matrx as gven by ε = mn(1,max(km,m co)) (Eq. (A.3) n Appendx). Varables n ths equaton can be arbtrary chosen, and here, K s taken as unty and M co = M, leadng to ε=mn ( 1,max(M,M ) ). (.33) Wth the above approxmaton, the dagonal block D (n Eq. (.7))) s transformed nto merely a scalar D = V t I+ S, A, =0. ( S, Γ 1 A +, = V t + σ ), S, I, (.34a) (.34b)

100 K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 Substtutng Eqs. (.6), (.31), and (.34a), Eqs. (.8a) and (.8b) becomes ( Q = V t + ( V = ( V = t + t + ( Q n+1 = Q + V σ ), 1 [ S, σ ), 1 [ S, Γ 1 Lower Lower σ ) [, 1 S, 1 Lower t + Γ 1 S, σ ), S, 1 Γ 1 S, A +, ] Q Γ 1 Res n A, +σ, I ] Q Γ 1 Res n ( S, Γ 1( F, Fv, )+σ, Q Upper S, A + n+1, Q n+1 ( = Q + V σ t + ), 1 ( S, 1 S, Γ 1 ( F n+1 Upper, Fv n+1, ) ] Γ 1 Res n, (.35a) )+σ, ) Q. (.35b) Note that the computatonal cost for the mplementaton of Γ 1 F s trval accordng to Turkel [5], by usng the followng form. Γ 1 F= F (1 ε)dp c (1 u v w H) T, (.36a) ( u dp=(γ 1) +v +w F 1 u F v F 3 w F 4 + F 5 ), (.36b) where F l stands for the l-th row of F (e.g., F 1 = (ρv n )). Wth the above expresson one can easly obtan precondtoned varables by avodng actual matrx operaton of Γ 1 F. ) Tme step: Tme step t s gven by the followng formula: t = CFL V [ ], (.37) σ, S, where CFL s Courant-Fredrchs-Lewy (CFL) number, and the spectral radus σ can be replaced by σ for precondtoned systems. The use of Eq. (.37)) s called local tme steppng (used n 3.5), whereas the global tme steppng (used n 3.1-3.4) takes the form of t= CFL max, (σ, / h, ). (.38) 3) Sub-teraton procedure: The sub-teraton (sometmes called Newton-teraton) s used to enhance the convergence rate outsde the LU-SGS loop. Eq. (.), V Q n t +Γ 1 ( F n, Fvn, )S,= Γ 1 Res n, (.39)

K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 101 wth Q n Q m, and applyng three-pont backward dfference (subscrpts and are omtted for clarty), leads to Thus, 3Q m+1 4Q n +Q n 1 Q m = 3 [ 3Q m 4Q n +Q n 1 = 3 Qm +3Q m 4Q n +Q n 1 (.40) = t ( F n Fv n )S }. V (.41) + t V Γ 1 {Res m ( F n Fv n )S }], (.4) where m s the number of sub-teratons, and when m reaches the specfed maxmum teraton number or the Q m reduced to the threshold value, the sub-teraton process s termnated as Q m Q n. Note that ths procedure acheves second-order temporal accuracy f t s frozen throughout the computaton. In addton, wth precondtonng matrx Γ, dual tme steppng s usually adopted for unsteady calculatons [4]. However, we dd not take ths strategy but used sub-teratons only to accelerate and stablze computatons of steady flows. The same dea s shared by developers of OVERFLOW code (see [5] and Prvate Communcaton wth Dr. Chrs Nelson, Mar. 010). In the subsequent sectons, CFL s chosen as CFL=0 n consderaton of both stablty and effcency, and no sub-teratons (= one sub-teraton) or three sub-teratons are employed, f not mentoned otherwse. The global tme steppng technque s usually used (unless stated otherwse). Based on the flow condtons explaned below (n Table 3), no slope lmters or turbulence models are used. Table 3: Test cases and condtons. Cases Condtons Comments 1) Vscous M =0.5, Re =5,000 [6, 7] Moderate Speed (Valdaton) ) Invscd M =0.001 0.1 Low Speeds 3) Vscous M =0.01, Re =,000 [3, 11] Low Speed.3 Flow condtons Computatons are conducted for a subsonc or a low-mach-number flow over NACA001 arfol, under the condtons gven n Table 3. The arfol has no angle-of-attack throughout the present study. In order to hghlght dfferences among methods, we generated the followng two grds that are not too fne or too coarse (Fg. ): We found that three nner-teratons were enough to accelerate convergence and/or stablze computatons for steady state computatons wth sgnfcant reducton ( orders drop) of RHS resdual at the thrd nner loop. Further teratons reduced resdual slghtly more, but dd not mprove accuracy or convergence dramatcally.

10 K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 Fgure : Computatonal grds for (a) nvscd (01 31 ponts) and (b) vscous (01 51 ponts) smulatons. Two-dmensonal, O-type, structured grds. 01 ponts n the crcumferental drecton, and 31 ponts (nvscd) or 51 ponts (vscous) n the radal (wall-normal) drecton, respectvely. The mnmum spacng near the wall for vscous cases s δ = 1.0e 3, based on the chord length of 1. Ths spacng acheves suffcent resoluton for boundary-layers consdered here. Far feld boundary s 50 tmes chord length away from the wall. 3 Results and dscussons The results are summarzed n Tables 4 and 5, n whch the followng notatons are used: S (Successful): The L -norm of densty resdual dropped at least four orders wth physcally correct soluton. U (Unsatsfactory): The soluton reached an unphyscal one wth poor qualty, and/or the resdual remaned sgnfcant. F (Falure): The calculaton dverged. In the Case 1, separatons ponts (desgnated as % chord length from the leadng edge) are also ncluded n the table. The result of each case wll be dscussed n the subsequent sub-sectons. 3.1 Case 1: M =0.5 vscous (lamnar) flow for code valdaton Ths test case has been wdely used as a benchmark [6, 7]. The computatons were conducted for 10,000 tmesteps. Typcal computed flow feld s dsplayed n Fg. 3. Fg. 4 shows hstores of drag coeffcent C D and L -norm of densty resduals for the successful cases. For all the successful cases, the computed flows were almost dentcal to each other, wth slghtly dfferent separaton ponts [6, 7] near the tralng edge. These locatons,

K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 103 (a) (b) Fgure 3: Computed flow feld by LU-SGS/SLAU, no sub-teratons, Case 1 (M = 0.5, vscous): (a) Iso- Mach-contours (0< M <0.59), (b) u-velocty contours; blow-up vew of separaton regon near the tralng-edge ( 0.01< u<0). resdual CD CPU tme [s] (a) CPU tme [s] (b) resdual CD tme steps (c) tme steps Fgure 4: Resdual and drag coeffcent hstores, Case 1 (M = 0.5, vscous): (a) Resdual vs. CPU tme, (b) C D vs. CPU tme, (c) Resdual vs. tme steps, and (d) C D vs. tme steps. (d) ncluded n Tables 4 and 5, are n good agreement wth reference separaton ponts of 80%-89% chord length (comparsons to other methods n lterature were made n a separate work [18]).

104 K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 Table 4: Computed results of no sub-teraton cases (CFL=0), S (Successful), U (Unphyscal or oscllatory), and F (Falure, dverged). Sub-teratons: 1 Case 1 Case Case 3 Implct Euler M =0.5 A) M =0.1 B) M =0.01 C) M =0.001 M =0.01 Schemes Fluxes Re =5,000[Sep.pont(%)] Re = (Invscd) Re =,000 Roe S[84.9] AUSM + S[85.0] U U SHUS S[84.9] LU-SGS P-Roe S[83.8] F F F F A-Roe S[8.7] S F F F AUSM + -up S[84.9] F F F F SLAU S[83.8] S S S S Roe S[84.9] AUSM + F F F SHUS U[89.7] LU-SGS P-Roe S[83.8] S S S S A-Roe S[8.7] F F F S AUSM + -up F F S S S SLAU S[83.8] S S S S Table 5: Computed results of three sub-teraton cases (CFL=0), S (Successful), U (Unphyscal or oscllatory), and F (Falure, dverged). Sub-teratons: 3 Case 1 Case Case 3 Implct Euler M =0.5 A) M =0.1 B) M =0.01 C) M =0.001 M =0.01 Schemes Fluxes Re =5,000[Sep.pont(%)] Re = (Invscd) Re =,000 Roe S[84.9] AUSM + S[84.9] U U SHUS S[84.9] LU-SGS P-Roe S[84.1] F F F F A-Roe S[83.0] S F F F AUSM + -up S[84.9] F F F F SLAU S[83.8] S S S S Roe S[84.9] AUSM + F F F SHUS U[85.4] LU-SGS P-Roe S[83.8] S S S S A-Roe S[8.6] S S S S AUSM + -up F F S S S SLAU S[83.8] S S S S From Tables 4 and 5 and Fgs. 3 and 4, the followng features are noteworthy: Coupled wth LU-SGS, all the fluxes yelded physcally correct solutons. In ths test case, the convergence rate was not practcally mproved by precondtonng of LU-SGS, although hstores of the drag coeffcent and resdual are slghtly affected (Fg. 4(a) and (b)). Even worse, calculatons dverged n some cases (AUSM +, SHUS, and AUSM + -up, see Tables 4 and 5). Ths would be because ) some combnatons, such as plu-sgs/ausm +, clashed and produced an nsuffcent amount of dsspaton (explaned later), or ) the scalng functon of AUSM + -up dd not work well n conuncton wth plu-sgs under the present flow condtons. Effect of Euler fluxes seemed to be mnor (Fgs. 4(c) and (d)), compared wth the above mentoned factors.

K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 105 For further code valdaton, we also conducted M =.0, α=0, nvscd flow computatons on a grd n Fg. (a). We pont out that all the 14 methods produced smlar and satsfactory results wth the help of Venkatakrshnan s lmter [8] wth Wang s correcton (ε =0.05) [9] (not shown). 3. Case : Low speed (M =0.1, 0.01, and 0.001), nvscd flow In ths secton, nvscd computatons were carred out for,000 tmesteps wth the freestream Mach number as a parameter: M = 0.1 (Case A), 0.01 (B), and 0.001 (C). Solutons and convergence rates are compared for dfferent methods. Fg. 5 shows the typcal computed flowfelds by LU-SGS/Roe, LU-SGS/SLAU and plu-sgs/slau. In Fg. 6, drag coeffcent hstores are shown for the three sub-teraton cases. Under the current flow condtons, the computed drag s regarded as an ndcator of numercal error. For example, n LU-SGS/Roe calculaton wth M = 0.01 (Case B), the drag coeffcent hstory reached a plateau at a sgnfcant value (Fg. 6(b)) wth an apparently unphyscal soluton shown n Fg. 5(a), even though the correspondng densty resdual showed three orders of reducton as well as other baselne fluxes cases (Fg. 7). The fnal values of computed drag coeffcents n all the runs are summarzed n Table 6. (a) (b) (c) Fgure 5: Computed flow felds ( 1 < C p < 1), no sub-teratons, Case B (M = 0.01, nvscd): (a) LU- SGS/Roe (Unphyscal), (b) LU-SGS/SLAU (Stable, slow convergence), and (c) plu-sgs/slau (Stable, fast convergence). From those Fgs. 5 and 6 and Tables 4 and 5, the followng general remarks are confrmed: (a) If no-precondtoned system of Eq. (.1) wthout Γ s solved, such as LU-SGS/Roe, calculatons do not dverge, but reach unphyscal solutons due to excessve numercal dsspaton of the method (Fg. 5(a)) [3,8]. (b) If only precondtonng A) (tme-dervatve precondtonng) s used, such as plu-sgs/roe, calculatons dverge (usually wthn a few tme steps), because the dsspaton n the numercal flux s not scaled properly [4,8,30]. (c) If only precondtonng B) (numercal flux precondtonng) s used, such as LU-SGS/P-Roe, calculatons dverge n most cases; the only excepton s LU-SGS/SLAU combnaton, showng slow convergence but a stable soluton (Fg. 5(b)) [3,8].

106 K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 CD CD CD tme steps (a) tme steps (c) CD tme steps (b) tme steps (d) Fgure 6: Drag coeffcent hstores, Case (nvscd), 3 sub-teratons: (a) Case A (M = 0.1), (b) Case B (M = 0.01), (c) Case C (M = 0.001), and (d) Case C (M = 0.001, wder scale for vertcal axs). resdual tme steps Fgure 7: Resdual hstores, Case B (M = 0.01, nvscd), 3 sub-teratons. (d) If both precondtonng A) and B) are used, such as plu-sgs/slau, physcally correct solutons are obtaned n most cases, wth clearly mproved convergence (Fg. 5(c)). The plu-sgs/a- Roe dverged wthout sub-teratons, but t was expectedly cured by ntroducton of sub-teratons (Recall that, as stated above, P-Roe, A-Roe and Roe are supposed to behave n dfferent manners at subsonc speeds). The plu-sgs/ausm + -up, however, faled to reach a soluton at M =0.1 as n the Case 1 (M =0.5, vscous case) even wth sub-teratons. These remarks are summarzed n Table 7.

K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 107 Table 6: Calculated drag coeffcent n nvscd flow (= numercal error) over NACA001 arfol, Case (M =0.1, 0.01, and 0.001; nvscd). Case, Sub-teratons: 1 Case, Sub-teratons: 3 Implct Euler A) M =0.1 B) M =0.01 C) M =0.001 A) M =0.1 B) M =0.01 C) M =0.001 Schemes Fluxes Re = (Invscd) Re = (Invscd) Roe 0.0137 0.070 0.680 0.0137 0.0704 0.567 AUSM + 0.0177 0.0956 0.844 0.0177 0.0915 0.7157 SHUS 0.0144 0.0758 0.661 0.0144 0.0738 0.5565 LU-SGS P-Roe (Dverged) (Dverged) (Dverged) (Dverged) (Dverged) (Dverged) A-Roe 0.0019 (Dverged) (Dverged) 0.0019 (Dverged) (Dverged) AUSM + -up (Dverged) (Dverged) (Dverged) (Dverged) (Dverged) (Dverged) SLAU 0.0038 0.0041 0.0069 0.0038 0.0037 0.0051 Roe AUSM + (Dverged) (Dverged) SHUS plu-sgs P-Roe 0.0055 0.0055 0.0055 0.0055 0.0055 0.0055 A-Roe (Dverged) (Dverged) (Dverged) 0.0019 0.0019 0.0019 AUSM + -up (Dverged) 0.0049 0.0049 (Dverged) 0.0049 0.0049 SLAU 0.0038 0.0037 0.0036 0.0038 0.0037 0.0036 Table 7: Summary of typcal computed results (expect for a few exceptons), S (Successful), U (Unphyscal or oscllatory), and F (Falure, dverged). Euler Fluxes Baselne (Roe, Low-Dsspaton (P-Roe, A-Roe, AUSM + -up, AUSM +, SHUS) SLAU) (Precondtonng B) Implct LU-SGS U S (slow convergence, SLAU) or F (other fluxes) Schemes plu-sgs (Precondtonng A) F S (fast convergence) Accordng to Tables 4 and 5, qualtatve performances of most of the methods presented here were ndependent on the Mach number. Quanttatvely, however, the drag coeffcents for non-precondtoned cases ncreased wth decreasng Mach number, from the smlar order to that of the precondtoned cases (M = 0.1) to orders larger (M = 0.001), ndcatng the recommended lower lmt of non-precondtoned methods for use be over M = 0.1. Meanwhle, those C D values for the precondtoned cases stayed almost constant (Fg. 6). The drag coeffcent, the error ndcator, showed the least value for plu-sgs/a-roe of 0.0019, followed by plu-sgs/slau (0.0037), plu-sgs/ausm + -up (0.0049), and plu-sgs/p-roe (0.0055), n the M = 0.01, three sub-teraton cases (the sub-teraton dd not seem to affect solutons sgnfcantly) (Table 6). Therefore, at least from these data, plu-sgs/a-roe appeared to have produced the most accurate soluton f t successfully worked. In addton, t s confrmed that a plu-sgs/low-dsspaton-flux combnaton can handle even M = 0.001 flow. Specfcally, plu-sgs/p-roe and plu-sgs/slau produced more successful results than other methods n a range of M =0.001 0.5. 3.3 Case 3: Low speed (M =0.01), vscous flow Ths test case has also been used to nvestgate the effects of precondtonng [3, 11]. Here, however, we focus on the vscous effects. Agan, typcal computed flow felds

108 K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 (a) (b) (c) Fgure 8: Computed flow felds ( 1< C p <1), 3 sub-teratons, Case 3 (M =0.01, vscous): (a) LU-SGS/Roe (Unphyscal), (b) LU-SGS/SLAU (Stable, slow convergence), and (c) plu-sgs/slau (Stable, fast convergence). CD resdual tme steps (a) CPU tme [s] Fgure 9: Drag coeffcent and resdual hstores, Case 3 (M = 0.01, vscous), 3 sub-teratons: (a) C D, (b) Resdual. (b) and drag/resdual hstores are shown n Fgs. 8 and 9, respectvely. As can be seen from these fgures and Tables 4 and 5, these computatons behaved n a broadly smlar manner to ther nvscd counterparts (Case B; Fgs. 5, 6(b), and 7): Non-precondtoned cases showed excessve drag and/or slow convergence. Fully precondtoned computatons yelded satsfactory solutons wth fast convergence. SLAU s agan only one flux whch showed accepted solutons wthout tme-dervatve precondtonng (Precondtonng A), n expense of slow convergence, though. Therefore, vscous effects played a mnor role n the present cases wth only one excepton: A-Roe flux (wthout sub-teratons) yelded a satsfactory soluton only n the vscous case, probably because ts pressure stablzaton term (Eqs. (.9c)-(.9e)) n combnaton wth the vscous source term (Fv n Eq. (.1b)) had a favorable contrbuton to the soluton. We also pont out that the dfference between nvscd and vscous cases wll be enhanced wth the use of local tme steppng as shown later n 3.5.

K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 109 Table 8: Summary of CFL Effect, S (Successful), U (Unphyscal or oscllatory), and F (Falure, dverged). Sub-teratons: 1 Case B) M =0.01, Re = (Invscd) [C D ] Case 3) M =0.01, Re =,000 Implct Scheme Euler Fluxes CFL=0 00,000 0,000 CFL=0 00,000 0,000 P-Roe S [0.0055] S [0.0055] S [0.0055] S [0.0055] S S S S plu-sgs A-Roe F F F F S S S S AUSM + -up S [0.0049] S [0.0049] S [0.0049] U [0.0035] S S S U SLAU S [0.0037] F F F S F F F Table 9: Summary of CFL Effect, S (Successful), U (Unphyscal or oscllatory), and F (Falure, dverged). Sub-teratons: 3 Case B) M =0.01, Re = (Invscd) [C D ] Case 3) M =0.01, Re =,000 Implct Scheme Euler Fluxes CFL=0 00,000 0,000 CFL=0 00,000 0,000 P-Roe S [0.0055] S [0.0055] S [0.0055] S [0.0055] S S S S plu-sgs A-Roe S[0.0019] F F F S S S S AUSM + -up S [0.0049] S [0.0049] S [0.0049] F S S S F SLAU S [0.0037] F F F S F F F 3.4 Effect of CFL numbers Now that we have confrmed plu-sgs/low-dsspaton-flux combnatons are recommended for low speeds and that ther performances are almost rrelevant to Mach number at least n the range of M = 0.1 0.001, we compared convergence rates of these combnatons wth dfferent CFL numbers rangng from 0 to 0,000 for M = 0.01, both n nvscd (Case B) and vscous (Case 3) cases. The results are summarzed n Tables 8 and 9, respectvely. Accordng to the results, the larger the CFL number, the more oscllatory or unstable the computaton tends to be. Fg. 10 shows resdual hstores for the cases wth sub-teratons (dverged cases are excluded, e.g., plu-sgs/slau wth CFL=00). Judgng from ths fgure, plu-sgs/p-roe wth CFL=,000 and CFL=0,000 (wth no remarkable dfference), followed by plu-sgs/ausm + -up wth CFL =,000, gave the fastest convergence (to machne zero) wth a satsfactory soluton both n the nvscd and the vscous cases; whereas these combnatons of methods wth CFL=00 acheved faster convergence for four-order reducton of resdual (about 100-00 tme steps; Fgs. 10(a), and (b)), whch s twce as fast as that of the CFL=0 case. Thus, n general, CFL=00-,000 appeared to be optmum. CFL larger than such values led to stagnated acceleraton of convergence (e.g., CFL = 0,000 wth P-Roe), or oscllatory or unstable solutons (e.g., CFL = 0,000 for fluxes other than P-Roe). SLAU s the excepton whch showed the smallest maxmum allowable CFL = 0 among four low-dsspaton fluxes for stable computatons. In Fg. 11, only successful cases wth CFL=00 or,000 both wth and wthout subteratons are shown together for nvscd (Case B) and vscous (Case 3) computatons, agan (as resdual hstores versus CPU tme ). It s confrmed that usng sub-teratons generally yelded faster convergence. Wth the effect of number of sub-teratons taken nto account, the choce of plu-sgs/p-roe wth CFL=,000 (3 sub-teratons), followed by plu-sgs/p-roe wth CFL =,000 (no sub-teratons) (nvscd, Case B) or plu-

110 K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 resdual resdual tme steps (a) tme steps Fgure 10: Resdual hstores wth dfferent CFL numbers for M = 0.01, 3 sub-teratons: (a) nvscd (Case B) and (b) vscous (Case 3) computatons. (b) resdual resdual CPU tme [s] (a) CPU tme [s] Fgure 11: Resdual hstores wth dfferent CFL numbers, wth and wthout sub-teratons, for M = 0.01: (a) nvscd (Case B) and (b) vscous (Case 3) computatons. (b) SGS/AUSM + -up wth CFL=,000 (3 sub-teratons) (vscous, Case 3), showed the fastest convergence rate towards machne zero; for 4-order drop of resdual, plu-sgs/p-roe wth CFL = 00 (3 sub-teratons) s the fastest n the nvscd cases (Case B), whereas n the vscous smulatons (Case 3) plu-sgs/p-roe wth CFL = 00 and,000 (3 subteratons), plu-sgs/a-roe wth CFL = 00 (3 sub-teratons), and plu-sgs/ausm + - up wth CFL =00 (3 sub-teratons) are n the fastest group. From those results, n terms of effcency, plu-sgs/p-roe appeared to be the best wth the maxmum allowable CFL numbers, followed by plu-sgs/ausm + -up. Based on ths lmted set of results, numercal dsspaton n P-Roe and AUSM + -up seemed to be compatble wth those produced by LU-SGS (or plu-sgs) for large CFL numbers, probably due to ther use of M. From Eqs. (.34) and (.37), the larger the CFL, the smaller the scalar D becomes, degradng ts dagonal domnance and hence, ntroducng more numercal dsspaton nto the system of equatons.

K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 111 3.5 Effect of local tme-steppng All the computatons above were conducted wth global tme steppng so that dscussons theren could be appled to unsteady flow computatons wth the use of dual-tme steppng, n whch temporal convergence s attaned n each tme step [4] (not actually covered n ths work, though). On the other hand, t s commonly known that the local tme-steppng technque (Eq. (.37)) accelerates convergence rate for steady flow computatons. Thus, n case one s nterested only n steady solutons, we employed local-tmesteppng technque for the test cases here. The computatons were conducted for the same cases and methods as n the above dscussons, and three sub-teratons were adopted. The results are summarzed n Tables 10 (CFL = 0) and 11 (CFL = 0 0,000 for plu-sgs/low-dsspaton-flux cases), and resdual hstores are shown n Fgs. 1 and 13. The portons of successful cases and others are roughly smlar to the global tme-steppng cases shown n Table 5. However, SLAU seemed to be destablzed by employng the local tme-steppng and showed no successful cases when coupled wth plu-sgs, although ths flux s only one whch showed relatve robustness at low speeds when used wth LU-SGS, agan. plu-sgs/p-roe s the most robust aganst ncreasng CFL whether the local tme-steppng s used or not. Table 10: Computed results of three sub-teraton cases (CFL=0 wth local tme steppng), S (Successful), U (Unphyscal or oscllatory), and F (Falure, dverged). Sub-teratons: 3 Case 1 Case Case 3 Implct Euler M =0.5 A) M =0.1 B) M =0.01 C) M =0.001 M =0.01 Schemes Fluxes Re =5,000[Sep.pont(%)] Re = (Invscd) Re =,000 Roe S[84.9] AUSM + S[84.9] U U SHUS S[84.9] LU-SGS P-Roe S[83.8] F F F F A-Roe S[8.7] S F F F AUSM + -up S[8.7] F F F F SLAU S[8.7] S S S U Roe S[8.7] AUSM + F F F SHUS U[8.9] plu-sgs P-Roe S[8.7] S S S S A-Roe S[8.7] S S S S AUSM + -up F F S S S SLAU U[83.8] F F F F Table 11: Summary of CFL Effect, S (Successful), U (Unphyscal or oscllatory), and F (Falure, dverged). Sub-teratons: 3 wth local tme steppng Case B) M =0.01, Re = (Invscd) [C D ] Case 3) M =0.01, Re =,000 Implct Scheme Euler Fluxes CFL=0 00,000 0,000 CFL=0 00,000 0,000 P-Roe S [0.0055] S [0.0055] S [0.0055] S [0.0055] S S S S plu-sgs A-Roe S[0.0019] F F F S S S S AUSM + -up S [0.0049] F F F S F F F SLAU F F F F F F F F

11 K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 resdual resdual resdual CPU tme [s] (a) resdual CPU tme [s] (b) CPU tme [s] (c) CPU tme [s] Fgure 1: Resdual hstores wth dfferent CFL numbers, wth and wthout sub-teratons: (a) vscous, M =0.5 (Case 1), (b), (c) nvscd, M = 0.01 (Case B), and (d) vscous, M = 0.01 (Case 3) computatons. (d) resdual resdual tme steps (a) CPU tme [s] Fgure 13: Resdual hstores wth dfferent CFL numbers, wth and wthout local tme steppng, for M =0.01: (a) nvscd (Case B) and (b) vscous (Case 3) computatons. (b) As shown n Fg. 1(a), the local tme-steppng clearly accelerated the convergence for vscous, moderate Mach number flow of M = 0.5 (Case 1). At ths flow speed, local tme-steppng appeared to be more effectve than precondtonng, and ths s explaned from the formulaton of

K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 113 Eqs. (.30), (.3), and (.37). The spectral radus, σ= V n +c+(vscous term), s domnated by V n and vares from one cell to another wth the order of c. Ths change s amplfed by changes of cell szes of the order of 10 or more (n ths case, about 100; Fg. (b)), sgnfcantly affectng the tme step t, compared wth precondtonng σ σ, whch leaves the order of the spectral radus remaned. At low Mach numbers, on the other hand, the local tme-steppng s less effectve than tmedervatve precondtonng (Fg. 1(b)-(d)). Agan, ths s clearly explaned from Eqs. (.30), (.3), and (.37). If precondtonng technque s used at low speeds, c s reduced to c wth the order of V n (whch s orders smaller than orgnal c), resultng n orders larger tme steps n the whole computatonal doman; on the contrary, f only the local tme-steppng s used, ths technque has lttle effect on the spectral radus σ= V n +c+(vscous term), snce V n c. Combnaton of the tme-dervatve precondtonng and the local tme steppng s consdered qute effectve, but ths set led to unstable or oscllatory solutons under some condtons (e.g., plu-sgs/ausm + -up wth CFL=00 or more). Accordng to Fg. 13, n whch resdual hstores for cases wth the tme-dervatve precondtonng and the local tme steppng both are shown, convergence acceleraton by ncreasng CFL s stagnated at CFL=00, showng comparable convergence rate wth CFL=,000 wthout local tme-steppng, n ether the nvscd (Case B) or vscous (Case 3) smulaton for M =0.01. 3.6 The Best schemes Accordng to the error estmaton above (M = 0.01, nvscd case of Case B), plu- SGS/A-Roe, followed by plu-sgs/slau, produced the most accurate soluton f t successfully worked. In terms of effcency, we counted mnmum requred CPU tme for machne zero convergence for each method wth ts maxmum allowable CFL (usually,000, except for plu-sgs/slau wth CFL=0) n M =0.01, vscous case of Case 3 (Fg. 10(b)). Accordng to ths crteron, plu-sgs/p-roe (357 sec.), followed by plu-sgs/ausm + -up (473 sec.), appears to be the best. To compare robustness, we smply counted numbers of successful cases marked n Tables 4 and 5 for global tme-steppng cases: plu-sgs/p-roe (10), plu-sgs/a-roe (7), plu-sgs/ausm + -up (6), plu-sgs/slau (10). Thus, plu-sgs/p-roe and plu- SGS/SLAU produced more successful cases than other methods n a range of M = 0.001 0.5, and these combnatons seem to be the most robust among all the methods. In local tme-steppng cases, on the other hand, plu-sgs/slau faled and plu-sgs/p- Roe and plu-sgs/a-roe are the best. Based on all the dscussons above, evaluaton of each combnaton of precondtoned LU-SGS scheme and a low-dsspaton flux s presented n Table 1. All n all, n low speed flow computatons, each method s suggested for use n the followng occasons: plu-sgs/p-roe: One seeks the fastest convergence, and M co s avalable. plu-sgs/a-roe: Obtanng the most accurate solutons s the top prorty. plu-sgs/ausm + -up: One seeks fast convergence, and M s avalable.

114 K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 Table 1: Evaluaton of Precondtoned LU-SGS Scheme/Low-Dsspaton Euler Fluxes. Accuracy Effcency Robustness (ndcated by numercal (ndcated by mnmum (numbers of successful cases) errors, C D, n M =0.01, CPU tme for machne global tme local tme step nvscd flow smulaton) zero convergence) step (10 cases) (5 cases) P-Roe 0.0055 357 sec. 10 5 plu-sgs A-Roe 0.0019 546 sec. 7 5 AUSM + -up 0.0049 473 sec. 6 3 SLAU 0.0037 More than 4,000 sec. 10 0 plu-sgs/slau: One s not sure whether the computaton reaches stable solutons, or there s no reference (unform) flow present. In addton, f computatonal tme really does not matter, LU-SGS/SLAU would be an alternatve choce. Therefore, t s expected that a promsng flux functon can be developed f, for nstance, SLAU s mproved by ncorporatng numercal dsspaton whle ts robustness s mantaned, by usng reference flow values as n P-Roe, A-Roe or AUSM + -up only when they are avalable. 4 Conclusons We carred out a comparatve study for several well-known or recently-developed low-dsspaton Euler fluxes coupled wth precondtoned LU-SGS (plu-sgs) mplct scheme n the framework of steady flows. It s confrmed that plu-sgs along wth low-dsspaton Euler fluxes gave accurate solutons wth sgnfcant mprovement of the computatonal effcency. The system of non-precondtoned counterparts, on the other hand, suffered from unphyscal solutons (no precondtonng at all), dvergence or slow convergence (control of dsspaton n numercal flux only), or dvergence of calculatons (precondtonng of tme ntegraton only). It s also confrmed that the recommended lower lmt of non-precondtoned methods for use be over M = 0.1. All n all, n low speed flow computatons, plu-sgs/p-roe, plu-sgs/a-roe, plu-sgs/slau or plu- SGS/AUSM + -up combnaton s suggested for use n the followng occasons: plu-sgs/p-roe: One seeks the fastest convergence, and the cutoff Mach number M co are avalable. plu-sgs/a-roe: Obtanng the most accurate solutons s the top prorty. plu-sgs/ausm + : One seeks fast convergence, and the freestream Mach number M s avalable. plu-sgs/slau: One s not sure whether the computaton reaches stable solutons, or there s no reference (unform) flow present. Moreover, f computatonal tme really does not matter, LU-SGS/SLAU would be an alternatve choce. In addton, SLAU s the only all-speed scheme whch s totally free from restrctons of specfyng reference values, such as M co or M.

K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 115 Therefore, t s expected that a promsng flux functon can be developed f, for nstance, SLAU s mproved by ncorporatng numercal dsspaton whle ts robustness s mantaned, by usng reference flow values as n A-Roe or AUSM + -up only when they are avalable. Furthermore, we would lke to pont out the followng: 1. Vscous effects played a mnor role.. The local tme steppng technque was proven to be effectve to accelerate convergence, but ts effect decreased wth decreasng Mach number M. At low speeds, the effect of local tme steppng s recovered f t s coupled wth precondtonng of tme ntegraton, but ths combnaton led to unstable or oscllatory solutons under some condtons. Appendx The transformaton matrx from prmtve varables to conservatve varables s wrtten as γ c 0 0 0 ρ T γ c u ρ 0 0 ρ T u Q q = γ c v 0 ρ 0 ρ T v γ c w 0 0 ρ ρ, q = u +v +w, (A.1) T w γ c H 1 ρu ρv ρw ρ q T where Q s the conservatve state vector (ρ,ρu,ρv,ρw,ρe) T, and q s the prmtve one employng pressure [14] (ρ,u,v,w,t) T. Then, the precondtoner of Wess and Smth for conservatve varables s wrtten as follows due to Turkel [5], although ths form s not used n the actual mplementaton. Γ 1 =I (1 ε)(γ 1) c ε=mn ( 1,max(KM,M co) ), dag(1, u, v, w, H) q q q q q u v w u v w u v w u v w u v w, (A.) (A.3)

116 K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 where K s constant usually taken as 0.5-1.0, and M co s cutoff Mach number whch s as the same order as freestream Mach number M. As stated n the man text, Eq. (.36) shown below s used nstead and the computatonal cost for the mplementaton of Γ 1 F s trval accordng to Turkel [5]: Γ 1 F= F (1 ε)dp c (1 u v w H) T, (A.4) ( u +v +w dp=(γ 1) F 1 u F v F 3 w F 4 + F 5 ), (A.5) where F l stands for the l-th row of F (e.g., F 1 = (ρv n )). Wth the above expresson one can easly obtan precondtoned varables by avodng actual matrx operaton of Γ 1 F. Nomenclature C D c p C p drag coeffcent specfc heat at constant pressure pressure coeffcent c pressure stablzaton coeffcent n All-Speed-Roe, 0.05 E F,F v ε mnmum spacng of grd, 1.0e 3 n vscous cases (based on the arfol chord length of 1) total energy nvscd (Euler) and vscous flux vectors precondtonng coeffcent γ specfc heat rato, 1.4 Γ H, precondtonng matrx total enthalpy cell ndces K coeffcent n precondtonng matrx, 1.0 κ M ṁ thermal conductvty,κ= µc p /Pr, Mach number mass flux, ṁ=ρu, µ molecular vscosty n P normal vector to the cell-nterface, (n x,n y,n z ) T pressure Pr Prandtl number, 0.7 Q (conservatve) state vector

K. Ktamura, E. Shma, K. Fumoto and Z. J. Wang / Commun. Comput. Phys., 10 (011), pp. 90-119 117 Re Reynolds number (based on the arfol chord length of 1) ρ densty T temperature u,v,w velocty components n x,y,z-drectons, respectvely x,y,z Cartesan coordnates V V n volume of cell velocty component normal to the cell-nterface V n =(u,v,w), n=un x +vn y +wn z ( ) arthmetc averaged value (ˆ) Roe averaged value () precondtoned value Subscrpts co L,R cutoff left and rght runnng wave components freestream condton maxmum value n whole computatonal doman n All-Speed-Roe Superscrpts m n value at m-th sub-teraton value at n-th tmestep Acknowledgments Ths work was conducted as a ont research between JAXA and Iowa State Unversty (ISU). The authors are grateful to Prof. Kozo Fu and other staff at JAXA and ISU both for the relevant support. The authors also thank Assocate Prof. Nobuyuk Tsubo at Kyushu Insttute of Technology and Dr. Chrs Nelson at Innovatve Technology Applcatons Company for havng constructve dscussons wth us. We also thank Dr. Kazuto Kuzuu at JEDI Center, JAXA, and Mr. Junya Aono at Research Center of Computatonal Mechancs, Inc., for mprovement of our numercal code. References [1] C. C. Krs, D. Kwak, W. Chan, and J. A. Housman, Hgh-fdelty smulatons of unsteady flow through tubopumps and flowlners, Comput. Fluds., 37 (008), 536 546.