周恒

(天津大學 力學系，天津 300072)

0 引言

對飛行器的設計來說，在初步選定外形后，準確算出其氣動力，包括升力、阻力、力矩，以及壁面熱流等，是進一步設計的基礎。早年這主要靠半經(jīng)驗的辦法解決，要做大量的研究性實驗和模型實驗。但模型實驗有很大的局限性，受制于在地面無法完全實現(xiàn)真實飛行時遇到的條件，無法做全尺寸實驗，也無法同時滿足所有的模型律。近年來，CFD(computational fluid dynamics)在研制各種飛行器中起到了越來越大的作用，已在很大程度上取代實驗的作用。但在計算飛行器所受的摩阻及熱載荷時，其結(jié)果的正確性有賴于對邊界層轉(zhuǎn)捩的預測能力和湍流計算能力。而如果飛行器因攻角較大而產(chǎn)生分離，則轉(zhuǎn)捩預測不準還將影響到升力的預測。

邊界層轉(zhuǎn)捩預測及湍流計算這2個問題都是已歷經(jīng)100年之久還未徹底解決的問題。在2006年的Annual Review of Fluid Mechanics 的文章“Critical Hypersonic Aerothermodynamic Phenomena”[1]中引用NASA的Bushnell在1997年的話說：“歷史上人類在預測所有高超聲速(甚至超聲速)飛行器的轉(zhuǎn)捩時幾乎從來沒有成功過”。并指出，制約CFD在高超聲速飛行器設計中應用的有四大因素，其中就有“模式化轉(zhuǎn)捩及湍流的能力”。

1 轉(zhuǎn)捩預測問題

1.1 轉(zhuǎn)捩的基本概念及關鍵問題

以平板邊界層為例。一般在前緣后有一段為層流，但其中有小的擾動。如果小擾動是增長的，則到下游某處，會經(jīng)傳捩過程而變?yōu)橥牧鳎鐖D1所示。

圖1 邊界層轉(zhuǎn)捩過程示意圖Fig.1 Sketch of the process of boundary layer transition

圖2所示為轉(zhuǎn)捩前后壁面摩擦系數(shù)和傳熱系數(shù)的變化?？梢娡牧鞫伪葘恿鞫味家髷?shù)倍。如果轉(zhuǎn)捩位置預測不準，顯然會帶來誤差。實際上，轉(zhuǎn)捩預測不準，對總摩阻的影響也可能不會太大，因為終究只占總值的一小部分。但對熱防護設計來說則影響很大。因為，如實際已是湍流的地方仍按層流設計，則很可能會被燒壞，而實際仍是層流的地方按湍流設計，則會增加不必要的防熱層的質(zhì)量。此外，對有攻角的彈體，如果轉(zhuǎn)捩線預測不準，則不僅會影響總阻力的計算，還會影響氣動力矩的計算。

圖2 轉(zhuǎn)捩前后壁面摩擦系數(shù)Cf及傳熱系數(shù)qw的變化Fig.2 Change of friction Cf and heat transfer coefficients qw across the transition

轉(zhuǎn)捩有2種。一種是如上所述由小擾動增長導致的轉(zhuǎn)捩；一種是由大的擾動引起的轉(zhuǎn)捩。前者要經(jīng)歷較長的小擾動演化過程，稱自然轉(zhuǎn)捩，而后者則在有了擾動后很快發(fā)生轉(zhuǎn)捩，稱by-pass轉(zhuǎn)捩。高空飛行時的邊界層一般是自然轉(zhuǎn)捩。

自然轉(zhuǎn)捩研究需解決的關鍵問題：

(1) 邊界層中的小擾動是如何產(chǎn)生的——感受性問題

這個問題約從20世紀80年代開始，研究還不夠充分。

對不可壓流動，邊界層外的擾動傳播速度就是自由流速度，而邊界層中的以不穩(wěn)定波形式出現(xiàn)的小擾動，其傳播速度小于自由流速度，二者不匹配。因而邊界層外緣引起的邊界層內(nèi)的擾動不能直接轉(zhuǎn)化為邊界層內(nèi)的不穩(wěn)定波，而要通過與邊界層內(nèi)平均流在局部(尺度較不穩(wěn)定波的波長小很多)發(fā)生較大變化處的相互作用才能在邊界層內(nèi)激發(fā)不穩(wěn)定波。一般注意的是前緣處，邊界層壁面的幾何性質(zhì)發(fā)生突變處，如2段不光滑(可以是一階導數(shù)不光滑)相接處，或有粗糙度處。

對可壓縮流，特別是高超聲速流，邊界層外的擾動除了有以自由流速度傳播的外，還可以以快聲波或慢聲波的形式出現(xiàn)，其傳播速度有可能和邊界層內(nèi)不穩(wěn)定波傳播速度相同或相近，可以直接或間接激發(fā)邊界層內(nèi)不穩(wěn)定波。

(2) 小擾動在邊界層中如何演化——擾動演化的線性及非線性理論

這個問題的研究比較充分，已有近100年歷史。

已知邊界層中的不穩(wěn)定波，可以用線性或非線性穩(wěn)定性理論，包括拋物化穩(wěn)定性方程，研究其演化，特別是其放大過程。這時要求解穩(wěn)定性方程。對不可壓流，這類方程的邊界條件是不難確定的，一般在壁面是不可滑移條件，邊界層外為零邊界條件。但對高超聲速邊界層，在壁面處多了一個溫度條件，在邊界層外會由于存在激波而較難處理(主要是馬赫數(shù)非常高的高超聲速邊界層)。特別是不穩(wěn)定波的演化對壁面溫度條件較敏感，而壁面溫度條件又往往很難確定，一般既不是定溫，也不是絕熱。

(3) 擾動演化至什么程度將觸發(fā)轉(zhuǎn)捩——轉(zhuǎn)捩判據(jù)

邊界層內(nèi)擾動演化至何時會觸發(fā)轉(zhuǎn)捩，是轉(zhuǎn)捩判據(jù)問題。由于轉(zhuǎn)捩不僅依賴于擾動大小，而且還與擾動頻率、形狀等有關，增加了提出轉(zhuǎn)捩判據(jù)的難度。國外目前還主要依靠經(jīng)驗或半經(jīng)驗方法。

高超聲速邊界層內(nèi)溫度可以很高，導致需考慮真實氣體效應。好在需要較長時間在大氣層中飛行的飛行器，其速度一般不會太高。如果Ma<7，則真實氣體效應主要表現(xiàn)為氣體的內(nèi)能必須考慮振動能，比熱將依賴于溫度，但還不需要考慮分子的離解及化學反應。相對來說處理還不太難。

1.2 轉(zhuǎn)捩預測方法

目前已有的轉(zhuǎn)捩預測方法有2類：

(1) 基于線性穩(wěn)定性理論的預測方法

基于線性穩(wěn)定性理論的預測方法具有一定的合理性，如圖1所示，在轉(zhuǎn)捩發(fā)生前，要經(jīng)歷一段小擾動增長的過程。這一過程，占了從開始有小擾動到實際發(fā)生轉(zhuǎn)捩的整個過程的大部分。但按原來的做法，要靠經(jīng)驗確定作為轉(zhuǎn)捩判據(jù)的一個參數(shù)。因此，這是一種半經(jīng)驗法，國外稱為eN法，N就是那個參數(shù)。而N的值，要針對某一類特定的流動，由經(jīng)驗得到。對eN方法的內(nèi)容和應用可參考文獻[2-3]。

在認為這一方法可靠的條件下，國外有人針對一類問題，將由其所得結(jié)果和邊界層的某些參數(shù)，如邊界層外沿的馬赫數(shù)及邊界層厚度等做關聯(lián)，使得更便于納入CFD的計算中，更便于工程技術人員的應用。

近年來，對基于線性穩(wěn)定性理論的轉(zhuǎn)捩預測方法做了較大改進，大大減少了對經(jīng)驗的依靠。

(2) 將轉(zhuǎn)捩前的擾動演化植入湍流模式計算中，使得轉(zhuǎn)捩和湍流合并計算。但擾動植入湍流計算這一點不太能令人信服，而且它有若干參數(shù)要靠經(jīng)驗確定，且對不同流動沒有通用性。

美國和歐洲的航空航天界基本都采用eN方法。

1.3 eN方法及其改進簡介

以超聲速錐體邊界層轉(zhuǎn)捩為例(如圖3所示)。

圖3 小攻角錐體邊界層轉(zhuǎn)捩問題坐標圖Fig.3 Coordinates used for the transition prediction of the boundary layer of a cone with small angle of attack

預測轉(zhuǎn)捩要研究從每一個子午面出發(fā)的各種小擾動可能導致轉(zhuǎn)捩的地方。

對每一個子午面，找出所謂的ZARF曲線(如圖4所示)。其上的每一個點，都代表一個頻率為F，在當?shù)?對應橫坐的標x值)放大率為0但卻是局部最大(指頻率相同但波數(shù)不同條件下)的擾動。

圖4 某子午面中的ZARF曲線Fig.4 ZARF curves in a certain meridian plane

對ZARF曲線上每一點代表的擾動波，計算它沿其傳播方向的放大倍數(shù)eN，N顯然為x的函數(shù)。當N達到某一由經(jīng)驗確定的預設的值，則對應的x值就是一個可能的轉(zhuǎn)捩位置。對所有ZARF上的點，可以找到一個可能轉(zhuǎn)捩位置的點集合。其中對應于最小的x的點就決定了由該子午面出發(fā)的擾動導致的轉(zhuǎn)捩點。由所有的子午面確定的轉(zhuǎn)捩點就組成轉(zhuǎn)捩線。

但是，在將這一方法用于超聲速小攻角圓錐的轉(zhuǎn)捩預測時，得到的轉(zhuǎn)捩線如圖3所示。其形狀甚至在定性上也和一般實驗所得結(jié)果不符。

圖5 某一Ma=6，攻角為1°的錐體邊界層， 用傳統(tǒng)eN法所得轉(zhuǎn)捩線Fig.5 Transition line obtained by conventional eN method for the boundary layer of a cone with angle of attack 1° and Ma=6

一個合理的轉(zhuǎn)捩預測方法，應該考慮前述與轉(zhuǎn)捩有關的3個關鍵因素，而原來的方法只是考慮了其中的一個因素，即小擾動的增長過程。

針對原來方法的不足，文獻[4]提出了2方面的改進意見。

(1) 考慮小擾動的實際產(chǎn)生機理

不從ZARF開始，而假設在離錐頭不遠的某一x處，所有子午面中的各種頻率擾動的幅值均相同。對每一子午面，求所有從該x點出發(fā)的擾動波的幅值沿x的變化。當其幅值達到邊界層外緣速度的0.015倍處，就認為該擾動會觸發(fā)轉(zhuǎn)捩。

在某一給定的初始擾動幅值下，得到了圖6a)中的轉(zhuǎn)捩線(Ma=6，攻角為1°錐體邊界層)。而文獻[5]在同樣基本流條件下，采用在錐頭下游不遠處用壁面吹吸方法引入復雜擾動波系，得到了圖6b)中的轉(zhuǎn)捩線。二者十分相似。

文獻[6]還用改進后的eN方法對其他能有實驗結(jié)果的情況做了驗算，結(jié)果也都滿意。圖7是其中的一個驗證例子。

1992年，King做了一個超聲速小攻角圓錐的轉(zhuǎn)捩實驗[7]。他所用的風洞有一個裝置，在該裝置不用的時候，風洞中的背景擾動主要是洞壁的湍流邊界層發(fā)出的噪聲擾動。而在該裝置使用時，實驗段前洞壁邊界層的噪聲擾動被抑制，此時的風洞稱為靜風洞，其主要背景擾動就是一般流動中的擾動。圖8中帶有實驗點的2條線是對應于這2種情況的2種轉(zhuǎn)捩線。靠上的是對應于靜風洞的結(jié)果，靠下的是非靜風洞的結(jié)果。對應于靜風洞的轉(zhuǎn)捩線更靠下游，這是因為這時的擾動更小。用上述改進后的eN方法，可以在適當選擇擾動初值的情況下得到和實驗很接近的結(jié)果。但對非靜風洞情況，無論是原來的eN方法，還是改進了的eN方法，都不能得到甚至只是定性上相符的結(jié)果，見圖8a)中最下面三角形點和圖8b)中最下面的長方形點。文獻[8]認為，這是因為對以聲擾動為主的情況，應考慮邊界層對聲擾動的感受性。這時擾動將不是在錐頭部被感受，而是在慢聲波的波速和邊界層中同頻率的小擾動波的波速相同的地方被感受。采用這一因素后，得到的結(jié)果如圖8c)中空心圓圈所示，和實驗基本符合。要指出的是，完全符合是不可能的，因為風洞中聲擾動在不同位置的幅值是不同的，而實驗中并未給出這么全面的數(shù)據(jù)。計算中不得不做一定的假設。

圖6 不同方法得到的轉(zhuǎn)捩曲線Fig.6 Transition curves obtained with different methods

圖7 另一情況與實驗結(jié)果的比較Fig.7 Comparison with experimental observation for another case

圖8 實驗與預測結(jié)果比較
Fig.8 Comparisons between experiments and theoretical predictions

(2) 給出基本不依賴經(jīng)驗的轉(zhuǎn)捩判據(jù)

eN方法雖基于線性穩(wěn)定性理論，但對從小擾動開始的自然轉(zhuǎn)捩，在擾動小于0.01時，線性理論可以很好地描述擾動增長過程。而不少轉(zhuǎn)捩的DNS結(jié)果顯示，轉(zhuǎn)捩開始時，一般擾動的大小也就在0.01～0.02之間。因此，建議以擾動幅值達到邊界層外沿流速的0.015為轉(zhuǎn)捩開始的判據(jù)。而不再像原來的eN方法中要針對不同流動靠實驗或經(jīng)驗去確定作為轉(zhuǎn)捩判據(jù)的N值。文獻[9]對這一判據(jù)是否可靠做了進一步的驗證，結(jié)果是肯定的。

2 湍流計算問題

2.1 目前湍流模式用于超聲速/高超聲速湍流邊界層時的問題

對高超聲速湍流邊界層的湍流計算，目前還要靠湍流模式。最常用的有：①代數(shù)模式，如Baldwin-Lomax (B-L)模式，Cebeci-Smith模式；②一方程模式，如Spalart-Allmaras(S-A)模式；③ 兩方程模式，如k-ε模式，k-ω模式及其改進型，Menter SST模式等。這些的共同特點是，模式最終都提供一個渦粘系數(shù)。代數(shù)模式最容易用，其次是一方程模式，再次是兩方程模式。 C. J. Roy & F. G. Blotther[10]在2006年， P. A. Gnoffo等[11]以及J. L. Brown[12]在2013年分別將各種模式計算結(jié)果與實驗結(jié)果進行對比，結(jié)果發(fā)現(xiàn)，如果邊界層沒有由于和激波相互作用而發(fā)生分離，則對壁面壓力和壁面剪應力，各種模式都能給出還算合理的結(jié)果，與實驗結(jié)果相比的誤差在10%以內(nèi)。但對壁面熱流，則各種模式和實驗結(jié)果的差別都較大，在某些局部地方，甚至可達200%。出乎意料的是，代數(shù)模式在熱流計算上，比一方程和兩方程模式的結(jié)果更好。但如果邊界層由于激波干擾而出現(xiàn)分離，則代數(shù)模式就會失效。

但是，只以實驗結(jié)果為根據(jù)，很難分析上述各種模式不足的原因。近年來由于計算機及計算技術的發(fā)展，已可以對湍流邊界層做湍流的直接數(shù)值模擬(DNS)。而DNS的結(jié)果可以提供詳盡的數(shù)據(jù)。而且，對高超聲速邊界層做實驗，壁面的熱邊界條件很難確切給出。一般它既不是絕熱，也不是定溫。即使其他實驗條件都相同，壁面的材料傳熱系數(shù)及壁面內(nèi)部溫度的不同，也可導致不同的壁面熱流條件。而DNS則可確切規(guī)定壁面溫度或熱流條件。因此，用DNS結(jié)果和用模式計算結(jié)果做對比，則可避免由于壁面熱邊界條件不完全對應帶來的不確定性。

目前對湍流做直接數(shù)值模擬，還只能做雷諾數(shù)比較小的。但是，有2種因素使得這一限制影響減小。一是對飛行在高空，例如30 km以上的情況，單位長度雷諾數(shù)實際不大。因為雷諾數(shù)Reunit=v/ν=vρ/μ，其中v是速度，μ是粘性系數(shù)，ρ是密度。對高超聲速而言，v雖然很大，比亞聲速可以大一個量級，但在30 km以上的高空，密度約是地面空氣密度的1/67，因此單位雷諾數(shù)比在地表附近做亞聲速飛行的要小得多。二是戰(zhàn)術導彈的尺寸一般不太大，因此實際雷諾數(shù)也不會很大。

上面列舉的兩方程模式，其最終所給出的渦粘系數(shù)，都是基于渦粘系數(shù)μT正比于k2/ε的假設，即μT=Ck2/ε，其中k是湍動能，ε是湍動能的耗散速度，C為常數(shù)。文獻[13]對Ma為6，壁面溫度tW為來流溫度5.5倍的平板邊界層分別用DNS和幾種模式計算k和ε，結(jié)果發(fā)現(xiàn)k的分布還相差不多，只是量的差別，而ε的分布則定性上都差別很大，見圖9a),9b)。而用不正確的ε組合出正確的渦粘系數(shù)，從道理上是說不通的。文獻[5]也針對Ma為6的平板湍流邊界層，把由DNS得到的k與ε代入k-ε模式的公式所算出的渦粘系數(shù)(圖9c)中的實線)，與直接由DNS結(jié)果推出的渦粘系數(shù)(圖9c)中的空心圓連線)進行比較，發(fā)現(xiàn)無論是定量還是定性都相差很多。可見原來的基本假定有問題。近年來，有不少研究發(fā)現(xiàn)，對平板邊界層、槽道流及圓管流動這類有剪切的壁湍流，沿壁面方向可分為若干層，每一層的湍流動力學特征不同。因此μT=Ck2/ε這一假設并不正確。

圖9 湍動能k、湍動能耗散率ε及渦粘系數(shù)μT的比較Fig.9 Comparisons for the turbulent kinetic energy k, dissipation rate ε of the turbulent kinetic energy and the eddy viscosity μT

2.2 對B-L模式所做的改進

由于航天用飛行器外形相對航空飛行器簡單，很多部位由直線、平面、錐體或接近于他們的幾何形狀組成。在沒有分離的地方，邊界層湍流接近于平衡或充分發(fā)展湍流，B-L模式可用，而B-L模式是最便于應用的湍流模式。因此，B-L模式迄今在航天系統(tǒng)還常被采用。據(jù)此，文獻[13]決定先針對B-L模式，看是否有改進的可能。他們在對高超聲速平板湍流邊界層做DNS的結(jié)果分析的基礎上發(fā)現(xiàn)，傳統(tǒng)的B-L模式有3個方面要做改進。

(2) 湍流Prandtl數(shù)不應是常數(shù)。層流狀態(tài)下的粘性和傳熱，其起因都是分子間的碰撞引起的，機理相同，因此其系數(shù)之比為常數(shù)的假定合理。但渦粘性和渦傳熱是流體微團間的動量和熱量傳遞。雖然最直接的原因都是由接觸處的分子碰撞導致，但在微團內(nèi)動量傳遞比熱量傳遞效率高。例如設想一個一維的情況如圖10 所示。有2個微團，各有速度v1，v2及溫度T1，T2，且v1>v2。當微團1追上微團

2，動量及熱量傳遞開始。如果不考慮流體的壓縮性，則動量傳遞瞬間完成，兩微團的速度立即變得相等。但熱量傳遞不可能瞬間完成，而仍要通過熱擴散完成。三維情況二者的差別要小一些。湍流強度越大的地方，二者差別越大。根據(jù)這一考慮，他們以邊界層中湍動能的分布(各種參數(shù)下不完全相同，但差別不大)為依據(jù)，給出了一個修正的表達式。

圖10 2流體微團相碰情況Fig.10 Sketch for the collision of two fluid lumps

(3)

對冷卻壁，要對B-L模式中的混合長公式在壁面附近做一定的修正。經(jīng)過這3個修正之后，B-L模式不僅能給出更精確的壁面摩擦系數(shù)和壁面?zhèn)鳠嵯禂?shù)，而且也能給出更精確的整個平均流和平均溫度的剖面。圖11是一個例子。

圖11中，Ma=6，壁面溫度Tw=0.61Taw，Taw為絕熱壁溫度。DNS為直接數(shù)值模擬結(jié)果，BL為傳統(tǒng)BL模式結(jié)果，modify1為僅第一種改進后的結(jié)果，modify2為僅有1，2兩種改進的結(jié)果，modify3為3種改進都采用后的結(jié)果。

3 結(jié)束語

需要指出的是，迄今所有湍流模式，都不能同時給出精確的壁面摩擦系數(shù)、壁面?zhèn)鳠嵯禂?shù)和整個平均流和平均溫度的剖面。而對有分離的問題，如果分離不是由入射激波或壓縮拐角等引起的，而是由于如攻角較大等因素引起的，則分離位置與分離前的邊界層厚度等有關。如果分離前的平均流剖面不準確，則計算出的分離位置一般也不準確。

如前所述，對有分離的情況，代數(shù)模式不能用。是否可能對其再做改進，以使之能用，還是必需用帶方程的模式，還需要研究。而如必需用帶方程的模式，則它們也需要有實質(zhì)性的改進。

注：圖中橫坐標分別為：a)無量綱流向平均速度；b)無量綱溫度；c)無量綱流向坐標；d) 無量綱流向坐標圖11 平板湍流邊界層的速度、溫度、壁面摩擦系數(shù)及壁面?zhèn)鳠嵯禂?shù)的分布Fig.11 Distribution of the mean flow velocity, temperature and heat transfer coefficient at the wall

參考文獻：

[1] Bertin John J, Cummings Russell M. Critical Hypersonic Aerothermodynamicc Phenomena[J]. Annul Review of Fluid Mechanics, 2006, 38：129-157.

[2] CEBECI T, STEWARTSON K. On Stability and Transition in Three-Dimensional Flows[J]. AIAA, 1980, 18(4): 398-405.

[3] CEBECI T, SHAO J P. CHEN H H, et al. The Preferred Approach for Calculating Transition by Stability Theory[C]∥ Proceeding of International Conference on Boundary and Interior Layers, ONERA, Toulouse France,2004.

[4] SU Cai-hong, ZHOU Heng. Transition Prediction of a Hypersonic Boundary Layer Over a Cone at Small Angle of Attck—with the Improvement ofeNMethod[J]. Science in China Series, 2009,52(1):115-123.

[5] LI X L, FU D X, MA Y W. Direct Numerical Simulation of a Spatially Evolving Supersonic Turbulent Boundary Layer at Ma=6[J]. Chinese Physics Letters, 2006,23(6):1519-1522.

[6] SU Cai-hong, ZHOU Heng.Transition Prediction for Supersonic and Hypersonic Boundary Layers on a Cone with Angle of Attack[J]. Science China Physics, Mechanics & Astronomy, 2009, 52 (8):1223-1232.

[7] KING R A. Three-Dimensional Boundary-Layer Transition on a Cone at Mach 3.5[J]. Experiments in Fluids, 1992(13):305-314.

[8] SU Cai-hong, ZHOU Heng. Transition Prediction of the Supersonic Boundary Layer on a Cone Under the Consideration of Receptivity to Slow Acoustic Waves[J]. Science China, Physics, Mechanics & Astronomy, 2011, 54(10):1875-1882.

[9] SU Cai-hong. The Reliability of the Improved eNMethod for the Transition Prediction of Boundary layer on a flat Plate[J]. Science China Physics, Mechanics & Astronomy, 2012, 55 (5)：837-843.

[10] ROY C J,BLOTTNER F G. Review and Assessment of Turbulence Models for Hypersonic Flows[J]. Progress in Aerospce Sciences, 2006,42:469-530.

[11] GNOFFO P A, BERRY S A, Van Norman J W. Uncertainty Assessment of Hypersonic Shock Wave-Boudary-Layer Interaction at Compression Corner[J]. Journal of Spacecraft and Rockets, 2013, 50(1):69-95.

[12] BROWN J L. Hypersonic Shock Wave Impingement on Turbulent Boundary Layer: Computational Analysis and Uncertainty[J].Journal of Spacecraft and Rockets, 2013, 50(1):96-123.

[13] DONG Ming, ZHOU Heng.The Improvement of Turbulence Modeling For the Aerothermal Computation of Hypersonic Turbulent Boundary Layers[J]. Science China Physics, Mechanics & Astronomy, 2010, 53(2)：369-379.