高超聲速尖錐邊界層轉捩數(shù)值模擬

宋 博 李椿萱

(北京航空航天大學 航空科學與工程學院，北京100191)

邊界層轉捩對高超聲速飛行器的氣動力、熱產生重要的影響，因此轉捩的精確預估對飛行器氣動和熱防護系統(tǒng)的設計極為重要［1－2］.然而，邊界層的轉捩受當?shù)豈ach數(shù)、單位Reynolds數(shù)、壁面冷卻效應、頭部鈍度 (包括entropy swallowing效應)、橫流或三維流動效應、壁面粗糙度、吹吸等流動和幾何參數(shù)的影響［3－4］，且這些因素相互干擾、耦合，難以精確模擬，以致迄今仍未能建立起完整的轉捩理論，特別是高超聲速邊界層轉捩的預測仍是尚未解決的流體力學難題.

基于層流脈動能 (laminar kinetic energy)的轉捩位置判斷方法在近年來受到了較多的關注.該方法的基本思想是轉捩的發(fā)生是由層流中存在的不穩(wěn)定波動引起的，從而通過模擬層流區(qū)域的脈動能發(fā)展過程來構造轉捩的預測模型.這種基于不同物理機理引發(fā)的轉捩模擬模式是對以往強烈依賴于經驗的模式的一個重大改進.文獻［5－7］最先應用穩(wěn)定性理論導出了不同擾動模態(tài)的渦粘性.文獻［8－10］通過一系列數(shù)值實驗驗證了高超聲速條件下該模式是相當有效的.

然而，文獻［8－10］的所有算例均假設轉捩過程為自然轉捩，但所參考的實驗結果均出自常規(guī)風洞，其壁面湍流邊界層所產生的噪聲要比飛行實驗中觀測到的高出一個量級，造成小擾動誘發(fā)自然轉捩的實驗呈現(xiàn)出逾越型 (bypass)轉捩機制，導致轉捩位置較飛行實驗的結果大大提前，同時機制的不同可令轉捩參數(shù)的影響趨勢隨之改變［11］.靜風洞的實驗段近壁流動為層流，其環(huán)境噪聲接近于高空的真實情況.本文作者曾在Favré(密度加權)平均的框架下推導了可壓縮流動的層流脈動能輸運方程，在一定的假設下結合尺度分析針對方程中出現(xiàn)的顯式可壓縮項進行了封閉，并采用靜風洞轉捩實驗的數(shù)據(jù)對所建立模型的模型系數(shù)進行了標定［12－13］.本文將采用上述已標定系數(shù)的模型開展高超聲速尖錐邊界層轉捩位置的數(shù)值預測研究.

1 基于Favré平均的層流脈動能方程

Favré平均N-S方程的具體形式參見文獻［4］.通過推導Favré平均可壓縮Reynolds應力方程，再對方程中的Reynolds應力張量進行收縮可導出笛卡兒坐標系下基于Favré平均的層流脈動能輸運方程:

其中，上標“～”與“'”分別代表流場變量的Favré平均量及其脈動;上標“－”與“″”分別代表流場變量的時間平均量及其脈動.

方程 (1)右端項的簡化和封閉可參閱文獻［12］，在此僅給出其最終形式:

其中

其中MaTL為脈動馬赫數(shù)，定義為

μtL為層流脈動的渦粘性，與誘導轉捩發(fā)生的不穩(wěn)定波動的特性有關，可寫成［7］

Cμ為常數(shù)，取值0.09;τtL為層流脈動的時間尺度，可表達為

ω1和ω2分別表示最不穩(wěn)定的第一、第二模態(tài)波動的頻率，可表示為

δ為速度邊界層厚度;Ue和Up分別為邊界層的邊緣速度和相速度，有Up≈0.94Ue;λ為流場的特征長度，取λ≈2δ.

在上述模型中，模型系數(shù)的取值為［12］

其中Tu表示來流脈動強度，定義為

轉捩的預測可通過對式 (2)與Favré平均Navier-Stokes方程組的聯(lián)立求解給出.其轉捩準則設定為:流場中滿足

的第1個點對應的流向坐標xt即被認為是轉捩的起始位置.

上述基于層流脈動能輸運的邊界層轉捩預測模型包含了一定的物理機理.可以看到，式(2)右端各項均依賴于層流脈動渦粘性系數(shù)μtL.其計算式 (3)中顯含了流場的特征時間尺度，該時間尺度是流場中最不穩(wěn)定波動的頻率的函數(shù)，而最不穩(wěn)定波的頻率與引發(fā)轉捩的物理機制密切相關.更多關于不穩(wěn)定波動頻率的討論可參閱文獻 ［14］.

2 數(shù)值模擬參數(shù)研究

本文采用隱式方法計算定常層流流動，采用顯式方法進行轉捩位置預測計算.對控制方程中無粘和粘性通量項的差分離散分別采用了Harten-Yee TVD格式和二階中心格式.在獲取層流平均流場后，即可將Favré平均層流脈動能方程與N-S方程組進行聯(lián)立求解以給出轉捩的位置.

2.1 單位Reynolds數(shù)效應

2.1.1 算例說明

文獻 ［15］進行了半角為5°尖錐的高超聲速邊界層轉捩實驗研究.實驗模型的頭部半徑為2.54×10－6m，來流 Mach數(shù)為6.0，攻角為0°，總壓為3.34×105Pa.8組實驗 (H1～H8)的來流單位長度Reynolds數(shù)Re∞的變化范圍為3.6×106m－1～25.6×106m－1.實驗分別對絕熱 (Tw/T0∞=0.86)和等溫 (Tw/T0∞=0.59)壁面條件下的流場進行了測量.

為考察文獻 ［12］所建立的高超聲速可壓縮轉捩預測模型對Reynolds數(shù)效應的捕捉效果，本節(jié)將以與上述轉捩流動實驗相同條件的尖錐繞流為算例開展數(shù)值模擬參數(shù)研究.算例采用等溫壁條件.所采用的流向、法向和周向網(wǎng)格點數(shù)分別為201，101和37.流向網(wǎng)格采用均勻分布，法向網(wǎng)格在壁面附近進行了加密.

2.1.2 邊界層轉捩位置預測結果

在轉捩預測計算中，對邊界層厚度的估算在此提出采用基于當?shù)乜傡史逯档臏蕜t.數(shù)值實驗表明，采用該準則可給出更為準確的邊界層厚度.在此，僅以算例H1為例，圖1給出了模型子午面上 x=0.1，0.3，0.5和0.8 m處的無量綱切向速度和以來流總焓無量綱化的無量綱總焓剖面.可以看到，在速度邊界層外緣附近，總焓峰值位置接近速度邊界層外緣，并在流動下游具有更好的吻合度，從而證實了基于總焓峰值的邊界層厚度判斷準則亦適用于等溫壁面條件的流動.

考慮到本算例所引的實驗是在生產型常規(guī)風洞中進行的，模型參數(shù)b的取值在此沿用文獻［8－10］在數(shù)值模擬高噪聲環(huán)境下轉捩實驗的做法，即令b=a.然而文獻 ［8－10］在其工作中假設了第二模態(tài)與橫流模態(tài)的影響可合并表達為ω2=Up/λ'，其中 λ'≈3.5δ.在上述假設下，文獻 ［9］數(shù)值模擬了Horvath實驗，并標定了實驗所用NASA LaRC 20”Mach 6風洞的湍流度Tu為0.67%.然而，采用該Tu值在本文的轉捩預測模型框架下進行的數(shù)值計算發(fā)現(xiàn)，轉捩位置與實驗結果相比大大提前.作者認為，這緣于Papp轉捩預測模型［8－10］中 ω2的表示形式并不適用于本節(jié)所研究的軸對稱尖錐流動，因為在零攻角狀態(tài)下橫流模態(tài)并不起作用;而考慮橫流模態(tài)后令λ'≈3.5δ的處理方式相當于減小了第二模態(tài)的頻率，即增大了其特征時間尺度.由式(3)可知該時間尺度的增大將加速層流脈動渦粘性的增長，從而導致轉捩過早發(fā)生.

圖1 尖錐子午面上無量綱切向速度和無量綱總焓剖面

基于上述討論，在此采用本文建立的可壓縮轉捩預測模型對算例H4進行了一系列不同來流湍流度下的數(shù)值模擬實驗，并由此將來流湍流度重新標定為0.435%.應用該標定值進行數(shù)值模擬所得的算例H1～H8的轉捩起始位置見表1所列.

表1 Reynolds數(shù)效應算例轉捩起始位置計算結果 m

可以看到，計算結果與實驗值吻合較好，特別是算例H2～H4，二者的誤差小于2%.隨著來流單位Reynolds數(shù)的增大，轉捩位置逐漸向頭部靠近.然而當來流單位Reynolds數(shù)較大時(算例H5～H8)，計算轉捩位置與實驗結果相比存在一定的滯后.一個可能的原因是此時層流脈動能在模型頭部下游快速增長，而數(shù)值計算采用的網(wǎng)格在頭部附近較稀，其流向分辯率不足，從而令脈動能的耗散過大，造成了轉捩的延遲.

2.2 攻角效應

三維效應對高超聲速飛行器邊界層內不穩(wěn)定波動的發(fā)展規(guī)律和轉捩位置的影響一直是研究的重點.在高超聲速條件下，三維邊界層轉捩的主導模態(tài)是橫流模態(tài).有攻角圓錐邊界層是典型的三維邊界層.考察高超聲速尖錐 (半錐角小于10°)攻角對邊界層轉捩影響的風洞實驗可以得到較為一致的結論:小攻角 (小于5°)條件下，與零攻角時相比，背風子午線上的轉捩位置前移，迎風子午線上的轉捩位置后移，且該趨勢隨攻角的增大而增大.然而，由于三維高超聲速邊界層轉捩過程受實驗環(huán)境噪聲、來流Reynolds數(shù)、Mach數(shù)、攻角以及錐角等多重因素的影響，實驗數(shù)據(jù)存在相當大的離散度.靜風洞的環(huán)境噪聲水平接近于真實飛行環(huán)境的噪聲，其實驗數(shù)據(jù)應具有更高的參考價值，但作者未搜集到研究高超聲速尖錐攻角效應的靜風洞實驗的相關文獻.

2.2.1 算例說明

本節(jié)將選用5°半角的尖錐模型，對其在攻角 α =0°，0.5°，1°，2°下的繞流流場進行數(shù)值模擬參數(shù)研究，分析小攻角狀態(tài)下的流場特征并采用本文所建立的可壓縮轉捩預測模型預估不同攻角下的轉捩位置.鑒于實施相關實驗的常規(guī)風洞的噪聲水平可對邊界層轉捩位置造成顯著影響，而目前對其尚無準確可信的標定值，且為便于對比分析攻角對轉捩位置的影響，在此計算的來流條件采用與靜風洞尖裙錐繞流實驗相同的條件［12］;轉捩預測計算時來流湍流度取值為Tu=0.1%，參數(shù)b取值0.54;壁面采用絕熱壁邊界條件.

2.2.2 邊界層橫流轉捩預測的數(shù)值模擬結果

應該指出，為了顯式表征高超聲速三維邊界層轉捩的橫流模態(tài)，在三維邊界層轉捩預測計算時，層流脈動的特征時間尺度應修正為

其中τcf表示橫流模態(tài)的特征時間尺度，以特征頻率表達有

其橫流模態(tài)的特征頻率ωcf=Up/3.5δ.

采用本文作者所提出的可壓縮轉捩預測模型［12］計算得到的迎風子午面和背風子午面上的轉捩起始位置見表2.可以發(fā)現(xiàn)，隨著攻角的增大，其迎風子午面上計算所得到的轉捩起始位置相比于零攻角呈單調增大，相當于迎風面上轉捩延遲發(fā)生;在背風子午面上，轉捩起始位置則單調減小，相當于背風面上轉捩提前發(fā)生.但轉捩延遲和提前的程度并不大，主要是由于計算所采用的網(wǎng)格在展向的細密度不足所致.即使如此，本項小攻角范圍的參數(shù)研究所顯示的攻角對轉捩位置的影響趨勢與高超聲速邊界層轉捩風洞實驗的結論是一致的.

此外，與零攻角尖裙錐絕熱壁條件下的繞流計算結果［12］相對比可以發(fā)現(xiàn)，零攻角尖錐上的轉捩位置存在顯著延遲.該結果再次證實了文獻［12］所得出的結論:逆壓梯度可增強流動不穩(wěn)定性，使轉捩提前發(fā)生.

表2 攻角效應算例轉捩起始位置計算結果

3 結論

本文通過對不同Reynolds數(shù)下繞尖錐的高超聲速流動的系列數(shù)值實驗研究了單位Reynolds數(shù)對軸對稱邊界層轉捩位置的影響，結果表明增大單位Reynolds數(shù)可使轉捩提前發(fā)生.隨后對攻角為0.5°，1°和2°的三維尖錐邊界層展開了數(shù)值模擬研究，計算并分析了轉捩位置在圓錐迎、背風面的不同變化趨勢，得到了與實驗相一致的結論:迎風面上轉捩延遲，背風面上則提前.

References)

［1］Lau K Y.Hypersonic boundary layer transition-application to high speed vehicle design［R］.AIAA-2007-310，2007

［2］Papp J L，Dash S M.Hypersonic transitional modeling for scramjet and missile applications［R］.AIAA-2002-0155，2002

［3］Morkovin M V.Transition at hypersonic speeds［R］.NASACR-178315，1987

［4］Wilcox D C.Turbulence modeling for CFD ［M］.2nd ed.California:DCW Inc，2000:174－179

［5］Warren E S，Hassan H A.Transition closure model for predicting transition onset［J］.Journal of Aircraft，1998，35(5):769－775

［6］Xiao X，Edwards J R，Hassan H A.Transitional flow over an elliptic cone at Mach 8［R］.AIAA-2001-0276，2001

［7］McDaniel R D，Nance R P，Hassan H A.Transition onset prediction for high speed flow［R］.AIAA-99-3792，1999

［8］Papp J L，Dash S M.Rapid engineering approach to modeling hypersonic laminar to turbulent transitional flows［J］.Journal of Spacecraft and Rockets，2005，42(3):467－475

［9］Papp J L，Kenzakowski D C，Dash S M.Extensions of a rapid engineering approach to modeling hypersonic laminar to turbulent transitional flows［R］.AIAA-2005-0892，2005

［10］Papp J L，Dash S M.A rapid engineering approach to modeling hypersonic laminar to turbulent transitional flows for 2D and 3D geometries［R］.AIAA-2008-2600，2008

［11］Schneider S P.Effects of high-speed tunnel noise on laminarturbulent transition ［J］.Journal of Spacecraft and Rockets，2001，38:323－333

［12］Song B，Lee C H.A favré averaged transition prediction model for hypersonic flows［J］.Science China，Technological Sciences，2010，53(8):2049－2056

［13］宋博，李椿萱.其于Favré平均的高超聲速可壓縮轉捩預測模型［J］.中國科學:技術科學，2010，40(8):879－885

Song B，Lee C H.A Favré averaged transition prediction model for hypersonic flows［J］.Science China:Technological Sciences，2010，40(8):879－885(in Chinese)

［14］宋博，李椿萱.基于非湍流脈動動能方程的高超聲速轉捩預測［J］.北京航空航天大學學報，2010，33(10):1162－1165

Song Bo，Lee Chunhian.A transition prediction model for hypersonic flows based on laminar fluctuation energy transport equation ［J］.Journal of Beijing University of Aeronautics and Astronautics，2010，33(10):1162 －1165(in Chinese)

［15］Horvath T J，Berry S A，Hollis B R，et al.Boundary layer transition on slender cones in conventional and low disturbance Mach 6 wind tunnels［R］.AIAA-2002-2743，2002