进气道1

American Institute of Aeronautics and Astronautics 1 A New Parabolized Navier-Stokes Algorithm and its Applications in Some Hypersonic Propulsion Aerodynamic Problems Chen Bing*, Xu Xu†, and Cai Guobiao‡ Faculty 403, School of Astronautics, Beijing University of Aeronautics and Astronautics, Beijing, 100083, PR China A new implicit finite-volume Single-Sweep Parabolized Navier-Stokes (SSPNS) algorithm is developed. Theoretical analysis is focused on the mathematic properties about the parabolized Navier-Stokes (PNS) Equations, especially on the treatment of streamwise pressure gradient. Then the original implicit time iterative Lower-Upper Symmetric Gauss- Seidel (LU-SGS) method is successfully extended to integrate the PNS Equations in the streamwise direction. The hybrid upwind schemes, including Advection Upstream Splitting Method (AUSM) family schemes and Low-Diffusion Flux-Splitting (LDFSS) schemes, are used to compute the crossflow inviscid fluxes, while central schemes for the viscous fluxes. Three typical flows, i.e., supersonic flat plate flow, 15° ramp hypersonic flow, and cone flows with different angles of attack, are calculated with the SSPNS codes. Numerical results agree well with those obtained from NASA’s UPS PNS codes and experimental results by Tracy or Holden et al. Furthermore, several scramjet component flowfields, including 3 hypersonic inlet flows and 2 Single-Expansion Ramp Nozzle (SERN) flows, are also obtained with the SSPNS codes. Results of inlets, such as flow structures, wall pressure distributions, and heat transfer coefficients, show good agreement with those of NASA UPS codes, IMPNS codes, SCRAMIN NS codes, and experimental data by Holland et al. SSPNS results of the 2D and 3D SERN flowfields also agree well with those of NASA’s experiments. By comparison with the traditional time-iterative Full Navier-Stokes (FNS) flow solvers, the SSPNS codes show 1~2 order of magnitude of computational speed faster and at least 1 order of magnitude of storage saving in the 3D hypersonic inlet flow field simulation. All the numerical results indicate that SSPNS is a highly efficient, highly accurate, and also highly robust algorithm for steady supersonic/hypersonic flows without any large streamwise separation, and it is appropriate to be used in the aerodynamic optimization design of scramjet inlet and nozzle. Nomenclature A, B, C = flux Jacobian matrices c = speed of sound, m/s e, h = total energy and enthalpy per unit mass, J/kg Ma = Mach number p = pressure, N/m2 T = temperature, K ρ = density, kg/m3 u, v, w = velocity components in the x, y and z directions, m/s Cp, Ch = wall pressure coefficient and heat transfer coefficient γ = specific heat ratio κ = thermal conductivity, W/(m·K) μ = viscosity coefficient, Pa·s Pr = Prandtl number * Postdoctoral Research Assistant, School of Astronautics, BUAA, E-mail: chenbing@sa.buaa.edu.cn. † Professor, School of Astronautics, BUAA. ‡ Professor, School of Astronautics, BUAA. 42nd AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit 9-12 July 2006, Sacramento, California AIAA 2006-4353 Copyright © 2006 by the author(s). Published by the American Institute of Aeronautics and Astronautics, Inc., with permission. American Institute of Aeronautics and Astronautics 2 I. Introduction -43A had achieved two flights in 2004, and it is believed that a full-scale hypersonic vehicle will be developed in the near future1. At a pace with this research, much effort has been put into computational fluid dynamics (CFD) to unveil the physics of hypersonic flows in the hypersonic propulsion systems1~4. However, the conventional time iterative Full Navier-Stokes (FNS) methods do suffer from a very slow convergence rate, which leads to high computational cost (in terms of CPU time and memory storage) and inefficient numerical analysis for steady supersonic/hypersonic flows. Some of the reasons are the severe viscous dissipation in the boundary layers, strong shock waves, shock-boundary interactions, and broad spectrum of the timescale in the flowfields. However, in the steady state supersonic/hypersonic viscous flow cases, the computational efficiency can be greatly promoted by using the approximate form of the FNS equations, known as Parabolized Navier-Stokes (PNS) equations5~42. The PNS equations can be integrated in a space-marching manner because of their hyperbolic-parabolic characteristics in streamwise direction. The parabolized Navier-Stokes algorithm has been proved to be accurate enough in supersonic and hypersonic steady flows in the absence of streamwise separation, and has been used successfully in many cases, including perfect gas flows5~25, chemically reaction flows26~32, and magnetohydro- dynamic (MHD) flows33~34, etc. Most of the space-marching PNS solutions are almost the same accurate as those of the corresponding time iterative FNS algorithms. And in some circumstances, the former ones are even more accurate due to much finer grids (especially in streamwise direction) employed in the space-marching algorithms. The space-marching PNS algorithms have gained much progress in the latter half of last century5~42. However, the instability problem is still one of the most challengeable problems in the algorithms now. This involves several factors; among which, the streamwise integration methods, the crossflow numerical flux calculation schemes, and streamwise pressure gradient treatment may be the most serious ones. In this paper a new Single-Sweep Parabolized Navier-Stokes (SSPNS) algorithm is developed. The original implicit time iterative Lower-Upper Symmetric Gauss-Seidel (LU-SGS) method by Jameson et al.43~44 is successfully extended to integrate the PNS equations in the streamwise direction. The hybrid upwind schemes, including AUSM- family (including AUSM+45, AUSMPW46, and AUSMPW+1~2, etc.) schemes and LFDSS schemes47~48, are used to discrete the crossflow inviscid fluxes, while the conventional central scheme for the viscous fluxes. The present paper is organized as follows. A brief description on the governing equations of gasdynamic is given in Section II with special emphasis on the treatment of the streamwise pressure gradient. In Section III, LU-SGS integration method and the AUSM-family, LFDSS schemes used in the SSPNS algorithm, are introduced in detail. Numerical test cases from simple hypersonic flows to complex hypersonic inlet and nozzle flowfields are presented to verify the properties (i.e., accuracy and robustness, etc.) of the SSPNS algorithm in Section IV. Finally, some conclusions based on the results of the previous sections are drawn in Section V. II. Governing Equations A. Three-Dimensional Full Navier-Stokes Equations The three-dimensional compressible full Navier-Stokes equations can be written in a general nonorthogonal co- ordinate systems (ξ, η, ζ) as 0) ˆˆ()ˆˆ()ˆˆ(ˆ = ∂ −∂ + ∂ −∂ + ∂ −∂ + ∂ ∂ ζηξ vivivi t GGFFEEQ , (1) where t denotes time, and the subscripts i and v indicate inviscid and viscous fluxes, respectively. The flow and flux vectors are J/ˆ QQ = , Jiziyixi /)(ˆ GFEE ξξξ ++= , Jiziyixi /)(ˆ GFEF ηηη ++= , Jiziyixi /)(ˆ GFEG ζζζ ++= , Jvzvyvxv /)(ˆ GFEE ξξξ ++= , Jvzvyvxv /)(ˆ GFEF ηηη ++= , Jvzvyvxv /)(ˆ GFEG ζζζ ++= , with the Cartesian flow and flux vectors Tewvu ],,,,[ ρρρρρ=Q , X American Institute of Aeronautics and Astronautics 3 Ti uhuwuvpuuu ],,,,[ ρρρρρ +=E , Ti vhvwpvvvuv ],,,,[ ρρρρρ +=F , Ti whpwwwvwuw ],,,,[ ρρρρρ +=G , Txxzxyxxxzxyxxv qwvu ],,,,0[ −++= ττττττE , Tyyzyyxyyzyyxyv qwvu ],,,,0[ −++= ττττττF , Tzzzyzxzzzyzxzv qwvu ],,,,0[ −++= ττττττG , and ρ is density; u, v, w are velocity components corresponding to Cartesian coordinates x, y, z; e and h are total energy and enthalpy per unit mass, respectively; p is pressure; τ is viscous stress; q is heat flux; and J is coordinate transform Jacobian. The variables have been nondimensionalized with the free stream variables as the reference values using the following formulas ∞ = UL tt / * , L xx * = , L yy * = , L zz * = , ∞ = U uu * , ∞ = U vv * , ∞ = U ww * , ∞ = ρ ρρ * , (2) 2 * ∞∞ = U pp ρ , ∞ = T TT * , 2 * ∞ = U ee , ∞ = μ μμ * , where T is temperature, μ is viscosity, U∞ is free stream velocity magnitude, L∞ is a characteristic length, superscript * denotes dimensional quantity, and subscript ∞ represents dimensional free stream conditions. Then the non- dimensional form of viscous stress τ and heat flux q are expressed as )2( Re3 2 zyxxx wvu −−= ∞ μ τ , )( Re xyxy vu += ∞ μ τ , )2( Re3 2 zxyyy wuv −−= ∞ μ τ , )( Re xzxz wu += ∞ μ τ , )2( Re3 2 yxzzz vuw −−= ∞ μ τ , )( Re zyyz vw += ∞ μ τ , (3) xx Tq ∞ −= Re κ , yy Tq ∞ −= Re κ , zz Tq ∞ −= Re κ , and Re∞ is Reynolds number, κ is the thermal conductivity. The subscripts x, y, and z represent partial differentiation with respect to Cartesian coordinates. The total energy e, enthalpy h, and speed of sound c are, respectively, written as )( 2 1 )1( 222 wvupe +++ − = ργ , ρ/peh += , ∞ = MaTc / , (4) where γ is specific heat ratio and Ma is Mach number. The viscosity μ can be divided into laminar (subscript l) and turbulence (subscript t) components. The turbulence part μt can be determined by turbulence models35, such as B-L model49 and k-ε model50, etc. The other can be calculated with Sutherland’s law51 ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ + + = ∞ ∞ TT TTl /4.110 /4.11012/3μ , (5) and the thermal conductivity is obtained according to the relation ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ + − = ∞ t t l l Ma PrPr)1( 1 2 μμ γ κ , (6) where Prl and Prt are Prandtl number for laminar and turbulence flows, respectively. Finally, the perfect gas state equation is used to complete the set of equations ρ γ pMaT 2 ∞ = . (7) American Institute of Aeronautics and Astronautics 4 B. Three-Dimensional Parabolized Navier-Stokes Equations According to the outstanding work of Rudman and Rubin39, Davis40, Lubard and Helliwell41, Gao42, et al., in the high Reynolds number flows the magnitude of viscous and heat diffusions are O(1) and O(ReL1/2) in streamwise and crossflow directions, respectively. Therefore, the velocity and temperature partial differentiations in the streamwise direction can be omitted. That is, the three-dimensional parabolized Navier-Stokes equations are derived from the steady-state full Navier-Stokes equation by neglecting any derivative in the streamwise direction contained in the stress tensor, and all viscous and heat fluxes in the streamwise direction. Assume ξ is the streamwise direction, the three-dimensional parabolized Navier-Stokes equations are reduced to the following form 0) ˆˆ()ˆˆ(ˆ ** = ∂ −∂ + ∂ −∂ + ∂ ∂ ζηξ vivii GGFFE , (8) where the superscript * of crossflow viscous fluxes denote viscous and heat diffusions having been neglected, this is always named as diffusion parabolizing assumption17, 42. Based on the differential chain rules, the derivatives in the stress tensor τ* and heat fluxes q* in Eq. (3) can be expressed as xxx ζϕηϕϕ ζη ⋅+⋅= , yyy ζϕηϕϕ ζη ⋅+⋅= , (9) zzz ζϕηϕϕ ζη ⋅+⋅= , where ϕ represents u, v, w or T, and subscripts ξ, η, ζ indicate partial differentiation with respect to the generalized coordinates. The PNS equations, i.e. Eq. (8), is a mixed set of hyperbolic-parabolic equations with respect to the streamwise direction, provided that 1) the inviscid core flow is supersonic; 2) the streamwise velocity components is everywhere greater than zero; 3) and the streamwise pressure gradient is properly handled18. This set of PNS equations are widely used in the space-marching codes all over the world, such as AFWAL PNS7, NASA UPS16~18, 24~25, 27~34, IMPNS13~15, TORPEDO10, 19, 35, etc. C. Streamwise pressure gradient treatments In the subsonic boundary layer, information can propagate upstream because of the presence of the pressure gradient term in the streamwise moment equations. As a consequence, a space-marching method is not well-posed and will come front with the instability problem. So it is necessary to do some special treatment to eliminate this difficulty. Vigneron et al.37 had proposed a successful method. This approach involves separating the streamwise inviscid flux vector into two parts PEE += *ˆˆ ii , (10) where Jiziyixi /)(ˆ **** GFEE ξξξ ++= , Tzyx JJJp ]0,/,/,/,0[)1( ξξξω−=P , Ti uhuwuvpuuu ],,,,[ * ρρρωρρ +=E , Ti vhvwpvvvuv ],,,,[ * ρρωρρρ +=F , Ti whpwwwvwuw ],,,,[ * ρωρρρρ +=G , and the streamwise pressure-gradient coefficient ω is determined through a linear eigenvalue analysis with the constraint that the hyperbolic-parabolic nature of the governing equations is preserved with respect to the stream- wise direction. This coefficient ω is expressed as37 ⎪⎭ ⎪⎬⎫⎪⎩ ⎪⎨⎧ −+ = 2 2 0 )1(1 ,1min ξ ξ γ γσ ω Ma Ma v , (11) where Maξ is the Mach number in the streamwise direction; σ0 is a safety factor (usually ranging between 0.75 and 0.9) included to account for nonlinearities that are not considered in the eigenvalue analysis; and subscript v is an abbreviation of “Vigneron”. It should be noticed that a value of ωv, a little smaller than 1, implies that small fractions of the streamwise pressure gradient are neglected in corresponding with Maξ slightly greater than 1. Let σv be the derivative of ωv with respect to Maξ2, that is ⎪⎩ ⎪⎨ ⎧ < −+ ≥ = ∂ ∂ = ,, ])1(1[ ,,0 )( 2 , 2 22 0 2 , 2 2 vsw vsw v v MaMa Ma MaMa Ma ξ ξ ξ ξ γ γσωσ (12) American Institute of Aeronautics and Astronautics 5 and Masw,v performs as a switching streamwise Mach number to determine pressure gradient to be split or not. Thivet38 pointed out that the streamwise pressure-gradient coefficient obtained from Eq. (11) is not smoothing with respect to Maξ2. This would introduce a sacrifice in the robustness of the algorithms, especially for high order schemes. Thivet proposed another streamwise pressure-gradient coefficient ωt, which is a C2 function of Maξ2 and specified according to the relation38 ⎪⎩ ⎪⎨ ⎧ < =≤≤Ω = ,,1 , 2 0),,( )( 22 , 0 2 , 220 2 ξ ξξ ξ κ π κ ω MaMa MaMaMa Ma tsw t tswtt t (13) where the trigonometric function, Ωt, is defined as [ ]1)cos( 12/ 1),( 202020 −+ − =Ω ξξξ κκπ κ MaMaMa tttt , and the constant κt0 is σ0αtγ(π/2-1), where σ0 is the safety factor as in Vigneron’s proposal, and αt=0.999. Then its derivative with respect to Maξ2 is ⎪⎩ ⎪⎨ ⎧ > =≤≤ ∂ Ω∂ = ,,0 , 2 0, )( ),( )( 2 , 2 0 2 , 2 2 20 2 tsw t tsw tt t MaMa MaMa Ma Ma Ma ξ ξ ξ ξ ξ κ πκ σ (14) and [ ])sin( 12/ 1 )( ),( 2000 2 20 ξ ξ ξ κκκ π κ Ma Ma Ma ttt tt − − = ∂ Ω∂ . Figure 1 displays the pressure gradient coefficients, the derivatives, and the switching Mach number of the two methods. It’s obvious that result of Thivet’s approach obeys the Vigneron’s limit. Moreover, the streamwise pressure gradient coefficient of Thivet’s approach is two-order of smooth with respect to Maξ2. 2 ξMa ωv ωt ω σ σv σt σ0 = 1.0 1.2 1.0 0.8 0.6 0.4 0.2 0.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 0.0 0.5 1.0 1.5 2.0 2 ξMa ωv ωt ω σ σv σt σ0 = 0.7 1.2 1.0 0.8 0.6 0.4 0.2 0.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Mav Mat σ0 1.0 1.5 2.0 2.5 3.0 3.5 4.0 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Vigneron Thivet M a s w (a) Safety factor σ0=1.0 (b) Safety factor σ0=1.0 (c) Switching Mach number Figure 1. Streamwise pressure gradient coefficient obtained with different methods38. Except for the pressure gradient splitting approaches, Chang et al.52 noticed that flux-vector splitting (FVS) schemes based on characteristic splitting, such as Steger-Warming53 scheme, also can be used to separate the elliptic information carrying part from the streamwise flux vector. The elliptic part in Eq. (8) is usually neglected or treated as a source term to preserve the PNS equations’ hyperbolic-parabolic nature. III. Numerical Method A. Discretization and Implicit Integration In present single-sweep space-marching algorithm, the elliptic part of the streamwise flux is simply omitted, i.e. substitute Eq. (10) into Eq. (9), and neglect vector P, then the PNS equations are of the form 0) ˆˆ()ˆˆ(ˆ *** = ∂ −∂ + ∂ −∂ + ∂ ∂ ζηξ vivii GGFFE . (15) 1. Discretization and Linearization of the PNS Equations The PNS equations, Eq. (15), are to be integrated on the finite-volume abcdefgh as Lawrence et al.16~18 did in NASA’s UPS algorithm, see Figure 2. The flow field is discretized by successively adding slabs of thickness Δξ as the solution proceeds. The indices in ξ, η, and ζ are n, j, and k, respectively. A typical control volume cell is constructed by two successive slabs, e.g. the nth and (n+1)th slabs. In the space-marching algorithm the area-averaged American Institute of Aeronautics and Astronautics 6 flow field properties are reserved on the second grid points54 (labeled as “●”), while in the traditional time iterative finite-volume methods on cell vertices (primary grid points, labeled as “□”) or cell center (third grid points, labeled as “∆”). ζ, k ξ, n η, j Primary Grid Point Second Grid Point Third Grid Point a b c d e f g h (n, j, k) (n+1, j, k) (n+1/2, j-1/2, k) (n+1/2, j, k-1/2) (n+1/2, j, k+1/2) (n+1/2, j+1/2, k) Figure 2. A typical finite-volume cell for the discretization of the PNS equations. On the structured grid shown in Figure 2, the discretization of the PNS equations, Eq. (15), can be written as ,0])ˆˆ()ˆˆ[( ])ˆˆ()ˆˆ[( ])ˆ()ˆ[( 2/1 2/1, *2/1 2/1, * 2/1 ,2/1 *2/1 ,2/1 * , *1 , * =−−−+ −−−+ − + − + + + − + + + n kjvi n kjvi n kjvi n kjvi n kji n kji GGGG FFFF EE (16) where n+1/2, j±1/2 and k±1/2 are used to label the other four crossflow surfaces of the control volume. It is convenient to represent the streamwise flux vector at a given ξ station using the functional notation given by ),(ˆ)ˆ( 1, 1 , *1 , * +++ = n kj n kji n kji dSQEE , (17) and n+1 denotes the spatial index in the streamwise direction where the solution is currently being calculated. Eq. (17) means that *ˆ iE is function of the grid metrics (i.e. geometry) designated by (dS) and the conservative flow variables (Q). It is also advisable to change the dependent variable from *ˆ iE to Q to avoid the complex mathematic calculation in extracting the flow properties from *ˆ iE . According to Eq. (17) the streamwise flux vector can be linearized spatially in the following manner , ),(ˆ ),(ˆ ),(ˆ)ˆ( , , 1 ,, * 1 ,, * 1 , 1 , *1 , * n kjn kj n kj n kjin kj n kji n kj n kji n kji d d d Q Q SQE SQE SQEE Δ ∂ ∂ += = + + +++ where the increment n kj n ji n ji , 1 ,, QQQ −=Δ + , and the flux Jacobians are n kj n kj n kjin kj n kj d d , 1 ,, * 1 ,, * ),(ˆ),(ˆ Q SQE SQA ∂ ∂ = + + , n kj n kj n kjin kj n kj n kj d d , ,, * ,, * , * ),(ˆ),(ˆˆ Q SQE SQAA ∂ ∂ == , and it is obvious that we have the identities n kj n kj n kj n kj n kj dd , 1 ,, *1 ,, * ),(ˆ),(ˆ QSQASQE ++ = , n kj n kj n kji ,, * , * ˆ)ˆ( QAE = . After substituting the above linearizations and flux Jacobians, the streamwise flux gradient becomes n kj n kj n kj n kj n kj n kj n kj n kji n kji dd , 1 ,, * ,, *1 ,, * , *1 , * ),(ˆ]ˆ),(ˆ[)ˆ()ˆ( QSQAQASQAEE Δ+−=− +++ . (18) General