袁湘江，李 筠，田俊武，涂國華

(1.中國航天空氣動力技術(shù)研究院，北京 100074;2.中國空氣動力研究與發(fā)展中心，四川 綿陽 621000)

0 引言

傳統(tǒng)的線性穩(wěn)定性理論基于平行流假設，其適用范圍必然受到一定限制。對于零壓力梯度邊界層，只有當考查范圍較大時，非平行流的影響才比較明顯。但在壓力梯度較大時，即使在不大的范圍內(nèi)，流動的非平行性也會對流動的穩(wěn)定性起重要影響。所以，對于復雜外形，邊界層非平行性對穩(wěn)定性的影響不能忽略。在非平行流的穩(wěn)定性研究中，近十多年來發(fā)展起來的拋物化穩(wěn)定性分析方法直接以發(fā)展的剪切層流動為研究對象，在揭示流動失穩(wěn)機理以及轉(zhuǎn)捩機理方面顯示出了其更寬的適用性和更高的準確性［1－4］。已有的研究成果表明:它能像DNS一樣刻畫擾動的傳播、失穩(wěn)及非線性作用直到層流的崩潰階段(以壁面摩擦力的突然增加來判斷)，并與實驗結(jié)果相一致。

美國、歐洲和日本等都開展了拋物化穩(wěn)定性研究。在國內(nèi)，南京航空航天大學的王偉志和唐登斌等人利用PSE方法對不可壓縮流動的失穩(wěn)問題進行了探討，天津大學和中國空氣動力研究與發(fā)展中心分別利用PSE方法對可壓縮平板邊界層的失穩(wěn)與轉(zhuǎn)捩問題進行了研究，中科院力學所則對拋物化穩(wěn)定性的理論方面進行了討論。但總的說來距離國際先進水平還有一定差距，未達到工程應用的要求。

為了推廣PSE方法在航空航天方面的運用，本文開展了軸對稱流動的拋物化穩(wěn)定性研究，推導出了軸對稱拋物化穩(wěn)定性方程，以具有較強工程背景的高超音速球鈍錐邊界層為研究對象，通過線性穩(wěn)定性方法得到起始位置處的擾動波形作為初始條件，利用PSE取得了有別于線性穩(wěn)定性分析的穩(wěn)定性初步分析結(jié)果。

1 軸對稱拋物化穩(wěn)定性方程

如圖1所示，設x0－y0－z0為笛卡兒直角坐標系，x－y－z為曲線坐標系，R為子午線的曲率半徑，rb為子午線上某點到對稱軸的距離，那么拉梅系數(shù)可以表示為:

其中θ為x軸與對稱軸的夾角。另外，對于帶曲率的平板，可取H3=1。

圖1 軸對稱體以及坐標系Fig.1 Axisymmetric body and coordinate

在通常情況下，y的值是邊界層厚度的量級，而R遠遠大于邊界層厚度，同時認為R沿x方向的變化遠小于R本身，于是有

采用無窮遠來流速度U∞，密度ρ∞，溫度T∞，長度l0等為參考值，將N-S方程無量綱化，并定義雷諾數(shù)R=ρ∞U∞l0/μ∞。設流動可用基本流加上擾動量構(gòu)成，如:

對于軸對稱流動，有:

其中“*”表示任一基本流場量。

把(1)式代入可壓N-S方程，考慮到(2)式，并省略非線性項后，可以得到線化的擾動方程:

其中 φ 由擾動量所組成，φ =［ρ，u，v，w，T］T，系數(shù) Γ、A、B、C、D、Vxx等只與基本流場有關。

假設擾動量可表示為:

其中

2 程序驗證

為了檢驗程序的正確性，分別用不同方法計算了平板邊界層擾動波的增長因子和考查了流向曲率對穩(wěn)定性的影響。當曲率半徑R趨向于無窮大時，曲面演變成平面。本文考慮Ma=4.5的平板問題，在當?shù)乩字Z數(shù)Rx=1500(參考長度為的地方引入T-S波，分別利用LST、PSE和DNS求解該T-S波的增長因子N沿流向的發(fā)展。圖2(a)給出了在不同流向位置處不同方法獲得的N值比較，圖2(b)給出了在不同流向位置處DNS與PSE方法獲得的波形比較，結(jié)果是吻合的。從而驗證了程序的正確性。

圖2 不同方法獲得的N與波形在不同流向位置的比較Fig.2 Comparing N-factor and disturbance shapes acquired by different methods on different streamwise positions

為檢驗本文PSE方法研究曲率對穩(wěn)定性影響的正確性，圖3給出了不同流向曲率半徑下擾動增長因子沿流向的增長情況。其中正的曲率半徑表示凸外形，負的曲率半徑表示凹外形。從圖中可以明顯看出，隨著曲率減小(R的絕對值增大)，擾動增長情況向無曲率的平板靠近;正曲率具有增強穩(wěn)定性的作用，曲率越大擾動增長越慢;負曲率具有失穩(wěn)作用，曲率的絕對值越大擾動增長越快。這與袁湘江等人［5］采用線性穩(wěn)定性(LST)的分析結(jié)果一致，與Christian［9］等人采用其他分析方法所得的結(jié)果也一致，還與Malik＆Balakumar［10］的不可壓LST結(jié)果一致。文獻［11］采用Malik＆Poll方法研究了曲率對穩(wěn)定性的影響，也得出正曲率具有穩(wěn)定性作用，負曲率具有失穩(wěn)作用。

圖3 增長因子隨Rx的變化Fig.3 Increasing factor N as a function of Rx

3 球鈍錐穩(wěn)定性分析

由于根據(jù)線性穩(wěn)定性理論分析，鈍體頭部流場中的邊界層總是穩(wěn)定的，利用這種方法對鈍體進行的轉(zhuǎn)捩預測，其轉(zhuǎn)捩位置都在遠離頭部的錐身上。但Widhopf和Hall的地面實驗表明，在某些條件下，轉(zhuǎn)捩位置是在鈍錐的頭部，并不是在錐身。這一點也可以從美國的NRV、回收的Gemini返回艙以及我國神州一號返回艙的燒蝕外形上得到了證實。近年來，人們經(jīng)過長期研究逐漸認識到，鈍體端頭邊界層的失穩(wěn)機制可能與邊界層的非平行性以及流向速度變化有關，也可能其轉(zhuǎn)捩機理根本不是自然轉(zhuǎn)捩，即躍過線性發(fā)展階段直接轉(zhuǎn)捩，即:“Bypass”轉(zhuǎn)捩。鑒于這些考慮，本文利用上述軸對稱拋物化穩(wěn)定性方程，對鈍錐端頭邊界層進行穩(wěn)定性分析，并與線性穩(wěn)定性理論的分析結(jié)果進行比較。初始流場的狀態(tài)參數(shù)為:M=6.0，Re=2 ×106，T∞=100K，TW=300K。計算外形是5°鈍錐，基本流場的計算區(qū)域和網(wǎng)格系統(tǒng)見圖4。為提高效率，在計算時，采用軸對稱NS方程和激波裝配技術(shù)。計算網(wǎng)格點為241×251。

圖4 基本流場的計算區(qū)域和網(wǎng)格系統(tǒng)Fig.4 Computational zone of base flow and mesh

在穩(wěn)定性分析中，以球鈍錐的駐點作為坐標原點，球鈍錐的表面弧長作為流向坐標。為了方便對兩種穩(wěn)定性理論的分析結(jié)果進行比對，取與體軸夾角為18°處，為擾動波的入口，入口處的擾動波由LST確定。圖5、圖6分別表示根據(jù)LST得到的放大因子與根據(jù)PSE得到的放大因子。由圖可知，根據(jù)LST，擾動始終衰減;而根據(jù)PSE，擾動波在傳播一段距離后，出現(xiàn)增長。可見對于球鈍錐邊界層非平行性以及流向速度變化對擾動增長與衰減特性有較大影響。圖7表示根據(jù)LST確定的在流向各站位的規(guī)一化波形比較，可見各站位波形峰值點沿流動方向逐漸向邊界層外緣移動，除此之外，波形幾乎沒有變化。圖8、圖9分別表示根據(jù)PSE確定的在流向前半程與后半程各站位的規(guī)一化波形比較，可見前半程各站位波形變化較大，與根據(jù)LST確定的前半程波形變化不同，根據(jù)LST確定的前半程波形一開始波形峰值點沿流動方向向壁面方向移動，壁面附近出現(xiàn)第二峰值，并且第二峰值逐漸增長到超過第一峰值后，第二峰值點逐漸向邊界層外緣移動。而后半程波形沿流動方向變化較小，只看到第二峰值點逐漸向邊界層外緣移動。圖10到圖15為波傳播前半程幾個站位由LST確定的規(guī)一化波形與PSE確定的對應站位規(guī)一化波形的比較。由圖可知，在擾動增長特性變化前，LST確定的波形與PSE確定的波形基本相似，而在擾動增長特性變化點附近PSE確定的波形變化較大，開始出現(xiàn)多峰值。在擾動增長特性變化后，擾動波模態(tài)改變，波形相對穩(wěn)定。

圖5 LST確定的增長因子Fig.5 Increasing factor N decided by LST

4 結(jié)論

本文推導了平板/軸對稱拋物化穩(wěn)定性方程，并通過與LST和DNS的結(jié)果比較驗證了程序的正確性，并系統(tǒng)研究了超音速鈍體邊界層的失穩(wěn)機制，結(jié)果表明:1)隨著擾動波向下游的傳播，曲率對擾動波的影響增大，正曲率起失穩(wěn)的作用，負曲率起穩(wěn)定的作用;2)流動的非平行性以及流向速度變化對擾動增長與衰減特性有較大影響，在一定條件下，可以改變擾動的增減性質(zhì);3)流動的非平行性以及流向速度變化對波形的有較大影響，尤其擾動波在增減性質(zhì)發(fā)生變化時，波形變化較大;4)球鈍錐邊界層中的失穩(wěn)擾動波是TS波通過演化形成的。線性理論不能反映TS波的演化，因此，無論在什么地方都找不到失穩(wěn)的解;5)在球鈍錐邊界層中疊加TS波時，應有足夠大的幅值，否則將在開始階段被完全衰減。

［1］HERBERT T.Parabolized stability equations［J］.Annual Review of Fluid Mechanics，1997，29:245-283.

［2］HU S，ZHONG X.Nonparallel stability analysis of compressible boundary layer using 3-D PSE［R］.AIAA 99 －0813，1999.

［3］CHANG C，MALIK M R.Non-parallel stability of compressible boundary layers［R］.AIAA 93 －2912，1993.

［4］LI F，MALIK M R.Spectral analysis of parabolized stability equations［J］.Computers ＆ Fluids，1997，26(3):279-297.

［5］袁湘江，周恒.考慮流向曲率和壓力梯度的可壓縮邊界層穩(wěn)定性分析［J］.空氣動力學學報，1998，16(3):276－281.

［6］CHRISTIAN M，STEPHAN B，MARC L，et al.Curvature and nonparallel effects in 3-D compressible transition analysis［R］.AIAA97－0824，1997.

［7］MALIK M R，BALAKUMAR P.Linear stability of three-dimensional boundary layers-effects of curvature and non-parallelism［R］.AIAA 93－0079，1993.

［8］LI M J，GAO Z.Analysis and application of ellipticity of stability equations on fluid mechanics［J］.Applied Mathematics and Mechanics(English Edition)，2003，24(11):1334－1341.

［9］CEBECI T，CHEN H H，KAUPS K.Further consideration of the effect of curvature on the stability of three-dimensional flows［J］.Computers＆ Fluids，1992，21(4):491 －502.

［10］涂國華，袁湘江.超聲速平板邊界層流動的二維拋物化穩(wěn)定性研究［A］.中國第一屆近代空氣動力學與氣動熱力學會議［C］.四川綿陽:國防工業(yè)出版社，2006.809－813.