超聲速邊界層中的模態轉換及壁溫影響效應

蘇彩虹, 宋明真

(天津大學 機械學院 高速空氣動力學研究室, 天津 300072)

0 引 言

邊界層轉捩預測是飛行器設計中的重要問題，因為湍流狀態下壁面摩阻比層流下大得多，而高速飛行時，湍流狀態下的壁面熱流也比層流大好幾倍。因此，只有準確預測轉捩位置，才能合理預估摩阻和表面熱流[1-3]。從小擾動開始的自然轉捩具有以下幾個典型的特征階段：1)感受性，即自由流擾動通過一定機制轉換為邊界層內不穩定波的過程；2)邊界層內不穩定波的線性增長；3)幅值增長到一定程度的非線性作用及轉捩至湍流。因此，想要從理論的角度做完全基于物理機制的轉捩預測，必須要回答以下三個問題[4]：1)自由流擾動如何激發邊界層內不穩定波；2)邊界層內的不穩定波如何演化；3)不穩定波發展到何種程度轉捩便會發生。

目前為止，以上問題中，只有第二個問題的線性演化階段有成熟的理論進行描述，即線性穩定性理論。線性穩定性理論基于局部平行性假設[5]，給定來流條件和擾動，求解Orr-Sommerfeld方程的特征值問題，可以計算邊界層內不穩定模態的增長過程。目前，基于線性穩定性理論的eN方法被認為是最有理論依據的轉捩預測方法，已在低速流的轉捩預測中取得了豐碩成果。然而，對于高速流的轉捩預測，遠沒有低速流的成功。主要有兩方面的原因：一是，在高速流中用以確定轉捩判據N值的實驗和經驗遠沒有低速流中的多；二是，在高速流中感受性的機制與低速流中顯著不同，導致eN方法所隱含的感受性的考慮，即所有擾動波在中性點處都具有相同的初始幅值，不再成立。Su & Zhou曾在對小攻角高超聲速圓錐進行轉捩預測時，對eN方法提出過兩個改進意見[6]，改進后的eN方法證實可以得到與實驗相比更合理的結果[7-8]。

不少研究表明，正確考慮感受性機制成為高速邊界層轉捩預測成功的關鍵[9]。與不可壓邊界層中的“尺度轉換”機制不同[10-12]，高速邊界層中不穩定波的激發是通過“同步”或“共振”機制實現的[13-18]。同步是指，頻率相同的自由流擾動與邊界層模態具有相同或相近的相速度，后者就會由前者激發。高超聲速圓錐邊界層的感受性研究表明[19]，自由流中的快、慢聲波會分別激發邊界層內的快、慢模態。在向下游的演化過程中，快模態的相速度逐漸減小，慢模態的相速度逐漸增大，直到與第二模態波相同(或相近)，從而將后者激發。由于高空中不存在聲源，因而并不存在聲波，主要的擾動源可能是由湍流脈動及密度或溫度不均引起的渦波和熵波。數值研究表明[20-21]，渦波和熵波可以通過與頭部激波作用，產生激波后的聲波，激波后的聲波再通過同步機制激發出快模態，快模態向下游演化，通過模態轉換激發不穩定的第二模態。

對于馬赫數大于4的超聲速邊界層，Mack第二模態是邊界層中起主導作用的不穩定模態。Fedorov[13]曾采用多尺度法研究了快模態到第二模態的轉換機制，表明在模態轉換區域附近，非平行性扮演著非常重要的作用。因此，基于局部平行假設的線性穩定性理論無法描述快模態到第二模態的轉換過程。若能發展從快模態出發、考慮模態轉換的轉捩預測方法將比原有的半經驗方法更有物理依據。

本文以超聲速平板邊界層為研究對象，采用數值方法研究了快模態與第二模態間的模態轉換過程，構建了模態轉換系數和模態轉換區域與擾動頻率之間的模型公式。基于此，發展了從快模態出發、考慮模態轉換的擾動幅值的計算方法。該方法可以準確描述多個壁面溫度條件下包含模態轉換過程的擾動演化。

1 計算參數和數值方法

1.1 計算工況

1.2 直接數值模擬方法

本文分別采用直接數值模擬方法(DNS)和線性拋物化穩定性方程(PSE)的方法計算快模態向下游的演化。

DNS采用的是二維N-S方程的擾動形式，即:

(1)

(2)

其中，C0=110.4 K。

采用有限差分方法離散控制方程。對流項采用五階迎風差分格式，黏性項采用六階中心差分格式，時間推進采用四階Runge-Kutta格式。壁面分別采用無滑移和等溫壁面邊界條件，上邊界為無反射特征邊界條件。計算域入口的擾動由線性穩定性理論(LST)分析給出，即:

(3)

1.3 拋物化穩定性方程

線性的拋物化穩定性方程(PSE)具有如下形式：

(4)

2 超聲速邊界層中的模態轉換過程

選取壁面溫度為140 K，擾動頻率為1，圖1(a)、圖1(b)中分別給出了LST分析得到的邊界層內快、慢模態的相速度和增長率沿流向的變化。可見，該工況下，快、慢模態的相速度并沒有嚴格的相交點，快模態與第二模態在x=100左右相速度最為接近，可認為是同步點。事實上，從圖1(b)可知，該位置與擾動開始放大的中性點非常接近。快模態在下游的演化過程中始終是衰減的，而慢模態，即Mack第一模態，其相速度和增長率曲線在下游與增長的第二模態相連。

(a)相速度

圖2 入口處快模態的特征函數

圖3 快模態向下游的演化

(a)x=165

3 考慮模態轉換的擾動幅值計算方法

基于線性穩定性理論(LST)的eN方法，通過累計擾動的增長率來計算擾動的演化，即:

(5)

-αi表示擾動的增長率，xN表示中性點的位置，An表示中性點處擾動的初始幅值。由于線性穩定性理論忽略了非平行性，而在模態轉換區域，非平行性的影響不再是可以忽略的高階量[13]，因此，基于線性穩定性理論的eN方法無法從快模態開始，描述擾動向下游第二模態的轉換過程。而PSE方法實質上求解的是一個演化問題，而非特征值問題，在實際的應用過程中遠不如LST方便[26]。并且，若擴展到三維問題，可能需要采用面推進的方法，計算量很大，無法被工程實際采用[27]。為此，我們仍以LST為基礎，發展可以考慮模態轉換過程的擾動幅值計算方法。

考慮模態轉換對第二模態波演化的影響應包含兩個方面。一是，考慮模態轉換發生的區間。模態轉換發生在一定的流向區間內，轉換完成后第二模態波才開始放大；二是，考慮第二模態開始放大的初始幅值。由于模態轉換發生前，快模態經歷一個衰減的過程，不同頻率的快模態衰減的程度不同，因此，不同頻率的第二模態開始放大的初始幅值也不同。因此，考慮模態轉換過程的擾動幅值計算方法可寫為:

(6)

其中，A0表示某一位置x0處快模態的初始幅值，xN為中性點的位置，Γ為模態轉換系數，表示在中性點處通過模態轉換生成的第二模態的幅值與入口快模態的幅值之比，其大小與頻率有關。由于模態轉換并非準確地在中性點處完成，模態轉換系數中生成的第二模態的幅值實質上是一種等效幅值。Δx代表模態轉換區間，也與頻率有關。假設在x0處激發的不同頻率的快模態具有相同的初始幅值，取A0=3×10-4。本文僅考慮從快模態開始的演化過程，并不包含其激發過程。若能通過感受性研究估計出快模態擾動幅值與頻率的關系，只要將式(6)中的A0替換為A0(w)即可。

根據PSE方法得到的擾動速度幅值的演化結果，以及采用LST從中性點處積分第二模態增長率的結果，如圖5所示，可計算Γ和Δx。具體步驟為：1)根據PSE得到的幅值最大值Amax，計算對應的LST在中性點處的幅值，即An=Amax/eN,N為LST從中性點開始積分擾動增長率得到的最大N值，由式(5)給出；2)計算模態轉換系數Γ，Γ=An/A0；3)由于LST沒有考慮模態轉換的過程，給出的第二模態增長比PSE的結果靠前，將其整體向下游平移Δx，兩者結果一致，Δx即為模態轉換區間。經過2)和3)兩步得到的結果，也就是利用Γ、Δx采用式(6)得到的第二模態增長的結果，如圖5中虛線所示。

圖5 模態轉換系數Γ和模態轉換區間Δx的示意圖

分別考慮壁面溫度為140 K、160 K、180 K和200 K的情況，采用LST得到的中性曲線如圖6所示。由于低頻擾動相對高頻擾動而言具有更大的增長區域，我們主要考察頻率較低的快模態。顯然，快模態在模態轉換前經歷一個顯著的衰減過程，由于我們主要的關注點是靠近中性曲線下支界區域的模態轉換過程，為了減小計算量，將計算域入口移至x=500處。在該處引入快模態擾動，分別采用PSE方法計算不同頻率的快模態向下游的演化，以及式(5)計算相應頻率的第二模態波在不穩定區域內的增長，按照上述方法得到模態轉換系數和轉換區間，結果如圖7和圖8所示。圖7中豎線標出了各壁面溫度下對應的中性頻率(中性點對應的頻率)。由于實際中，快模態總是在不穩定區上游被激發，因此，我們只考慮頻率小于中性頻率的快模態。由圖7可見，不同壁面溫度下的模態轉換系數隨頻率的變化具有相似的規律，即隨著頻率增加，模態轉換效率增大，直至某一頻率(略小于中性頻率)。這是因為，頻率越高的擾動，入口的位置距離中性點越近，快模態的衰減過程越短，從而對應的轉換效率越高。而對于非常接近中性頻率的擾動而言，模態轉換需要在一定區域內完成，反而占據了一部分不穩定區域，因而降低了轉換效率。圖8給出的結果表明，模態轉換需要的流向距離隨著頻率的增加而減小。

圖6 不同壁面溫度條件下的中性曲線

圖7 模態轉換系數Γ隨頻率ω的變化

圖8 模態轉換區間Δx隨頻率ω的變化

可以發現，不同壁面溫度條件下，圖7給出的曲線具有相似的性質。為了在計算擾動演化時考慮模態轉換的影響，需要建立Γ-ω的關系。考慮以下變換:

ω*=ωn-ω

(7)

ωn為中性頻率，Δω代表了擾動波的頻率范圍，即ωn-Δω為考察的最低頻率。本例中取Δω=0.3。Γm是Γ-ω曲線的峰值。假定Γm，Γ|ωn-Δω隨壁面溫度Tw變化是線性關系，利用圖7中Tw=140 K和Tw=200 K的結果，可得:

Γm=-0.1492Tw/T∞+0.7144

Γ|ωn-Δω=-0.0424Tw/T∞+0.1622

(8)

在本文的計算中，ωn通過對基本流場進行穩定性分析得到。在實際應用中，簡單起見，也可利用已知的兩個壁面溫度下的中性頻率做線性插值得到。對于本文的算例，結果差別很小，未在此列出。采用PSE和LST計算更多壁面溫度的Γ*-ω*曲線，所得結果如圖9所示。可以發現，對于很廣泛的壁面溫度范圍，Γ*-ω*曲線幾乎是重合的。可以利用數據擬合的辦法，建立Γ*-ω*的關系，即:

圖9 不同壁面溫度條件下的Γ*-ω*曲線

Γ*=1/(C1ω*2+C2ω*+C3)+C4ω*

(9)

其中，C1=296.7722，C2=-10.2683，C3=1.0874，C4=-0.1319。這樣，根據式(9)，結合式(7)～(8)可計算給定壁面溫度條件下的模態轉換系數。

類似地，由圖8給出的Tw=140 K和Tw=200 K的數據，可建立模態轉換區間Δx與頻率和壁面溫度的關系，即:

Δx=-6.248Tw/T∞+14.8/ω+8.56

(10)

Fedorov等人[13-14]采用多尺度法研究模態轉換時發現，模態轉換區間Δx滿足:

Δx=Cε2/3

(11)

其中，ε=1/ReLn，ReLn表示以中性點距平板前緣的距離Ln定義的雷諾數，C為常數。由于不同壁面溫度下不同頻率對應的中性點位置不同，因此ReLn與頻率有關，即得到的Δx是頻率的函數。在應用式(11)時，需要針對每個壁面溫度條件，采用PSE方法計算某一個頻率下的Δx，利用式(11)確定常數C。表1中給出了不同壁面溫度下C的取值。圖10給出了分別由式(10)和(11)計算的Δx與圖8給出的數值結果的比較。可見，相比于擾動線性增長的區域(本例中為幾百到幾千個無量綱長度)，這兩種方法計算出的Δx的差別很小。

表1 不同壁面溫度下C的取值

圖10 不同方法計算的Δx隨頻率的變化曲線

4 擾動幅值計算方法的驗證

考慮模態轉換過程的擾動幅值計算方法可以描述為：給定壁面溫度Tw，利用式(7)～(9)計算不同頻率快模態的模態轉換系數Γ，再結合式(10)或(11)給出的Δx，代入到式(6)中，計算由模態轉換激發的不同頻率的第二模態波的幅值演化。為了驗證該方法的有效性，采用PSE計算快模態的演化進行比較。圖11給出了壁面溫度為150 K時，分別采用式(10)和式(11)計算得到的幅值演化曲線與PSE所得結果的比較。可見，考慮模態轉換的擾動計算方法可以得到與PSE方法一致的結果。并且，分別采用我們建立的經驗公式(10)與Fedorov給出的式(11)得到的模態轉換區間，所得的幅值演化結果幾乎沒有差別。類似地，圖12和圖13分別給出了壁面溫度為170 K和190 K的幅值演化結果。同樣，考慮模態轉換后的擾動演化結果與PSE的結果吻合很好。因此，采用該方法計算擾動的演化，結合一定的轉捩判據，可以預測從快模態開始的包含模態轉換過程的邊界層轉捩位置。

圖11 擾動幅值演化結果(Tw=150 K)

圖12 擾動幅值演化結果(Tw=170 K)

圖13 擾動幅值演化結果(Tw=190 K)

5 結 論

本文采用直接數值模擬和拋物化穩定性方程的方法，對馬赫數4.5的超聲速平板邊界層中快模態到第二模態的模態轉換過程開展了研究。研究表明，PSE方法可以得到與DNS一致的擾動演化結果。利用PSE所得的結果，定義了模態轉換系數和區間。研究發現，不同壁面溫度下，模態轉換系數、轉換區間與擾動頻率的依賴關系存在相似的規律。在此基礎上，建立了考慮壁面溫度條件的模態轉換系數和轉換區間的計算模型，發展了考慮模態轉換的擾動幅值計算方法，并采用PSE對新發展的方法進行了驗證。結果表明該方法可以在較廣泛的壁面溫度范圍內，準確地描述從快模態到第二模態轉換的擾動演化過程。

致謝:感謝天津大學吳雪松教授，多次與他進行的關于感受性問題的討論給予作者很多啟發性的思考。感謝天津大學韓宇峰博士在論文研究過程中給予的熱情幫助。