超聲速鈍錐三維分離條件的數(shù)值分析

李 沁，張涵信

（1.中國空氣動力研究與發(fā)展中心，四川綿陽 621000；2.北京航空航天大學國家計算流體力學實驗室，北京 100083）

0 引 言

分離和旋渦是常見的流動現(xiàn)象。對于二維分離，流場通常會出現(xiàn)逆壓梯度，并產生分離旋渦等流動結構；三維分離則復雜得多，首先會存在不同的分離類型，如開式分離與閉式分離；與此相聯(lián)系，分離起始會出現(xiàn)不同的奇點形態(tài)，如正常點起始、鞍點起始以及鞍結點組合起始；旋渦的渦軸既可以與來流垂直（與二維情況接近），也可能與來流平行。在實際問題中，在一定迎角下的機翼分離、二維拐角／臺階分離都屬于二維情況下的分離，逆壓梯度引起的壓力增加對流體力學特性尤其是氣動性能產生重要影響；而三維分離產生的旋渦，既可能增加升力——如流向渦產生的渦吸力，又可能由于分離流態(tài)的變化，引起氣動性能的降低，如渦軸抬升以及渦破裂。因此，準確預測分離，開展分離發(fā)生、發(fā)展的理論研究，具有重要的實際意義。其中分離理論的研究成果，可以直接用于計算結果的檢驗和分析；同時當前數(shù)值計算所采用的邊界條件，由于在分離點附近理論上不嚴格、數(shù)值上存在誤差，往往造成計算得到的分離與實驗存在較大差距，從而影響氣動特性的準確計算；分離理論研究將有助于這些難題的解決。

不可壓二維分離的理論由Prandtl［1］給出，其標志是Prandtl分離條件，該條件既被用在定常也被用在非定常的表面分離；對于更為復雜的物面外分離，Moor、Rott和Sears提出 MRS準則［2］。張涵信院士對此也作了研究并得到更為一般的形式［3］。三維定常分離的研究學者主要有Hesieh、K Wang、Tobak和Peake等，其中關于分離線性質的研究曾經存在長時間的爭論，最后由張涵信院士給出的結果平息了學術爭論［3－6］，更為重要的是研究還給出了分離判則。在文獻［7］中，張涵信院士給出了關于流動分離詳細和完整的理論。

最近張涵信院士［8］通過將非定常流動中密度影響吸收轉換為一個類似速度梯度項，從而利用定常分離的研究成果來闡述非定常分離。比較非定常分離判則和定常分離判則，可以發(fā)現(xiàn)非定常分離的前三個條件實際上與定常分離是一致的，所以深入開展定常分離的數(shù)值分析仍然具有必要性。

關于如何開展三維定常分離理論的數(shù)值分析，目前有兩個技術性問題：第一個是三維分離判則給出的結論基于分離坐標系，三個坐標軸分別是分離線、分離線法向和物面法向，而實際計算曲線坐標系除個別簡單情況外，坐標系與分離坐標系之間缺乏直接聯(lián)系，對分離判則進行計算分析與驗證較為困難；第二個是密度的影響問題，密度在高馬赫數(shù)下的分離判則中所起作用需要深入開展研究。

本文的研究正是基于這樣的考慮，開展了15°半錐角鈍錐超聲速層流分離條件的數(shù)值研究。來流馬赫數(shù)為1.8和10.6，迎角從10°到40°。研究一方面探索分離條件的數(shù)值分析方法，給出分析的結果并與理論進行比較；另一方面也研究不同馬赫數(shù)對分離形態(tài)的影響，并給出不同情況下流動的具體分離特征。

在第1節(jié)，首先介紹數(shù)值模擬使用的方法并給出驗證算例；在第2節(jié)給出三維數(shù)值分離判則的數(shù)值分析方法；應用數(shù)值分析方法，在第3節(jié)中對超聲速鈍錐計算結果進行了分析，對各種條件下分離特征進行了研究；最后在第4節(jié)給出研究結論。

1 數(shù)值方法

1.1 計算方法

數(shù)值模擬控制方程為守恒型的積分Navier－Stokes方程：

其中U是守恒變量向量，F(xiàn)是控制體邊界上的通量向量。V和S是控制體的體積和邊界面面積是S的法向單位向量。采用通?；诹骟w的有限體積方法進行離散，控制方程可表示為：

式中，V代表網格單位的體積，F(xiàn)、Fv、S分別代表網格邊界面上的無粘、粘性通量和面積向量；由于所采用方法為標準有限體積算法，文中將不對這些變量作詳細解釋。進一步，有限體積求解采用格心法，即流動參數(shù)控制點取在控制體中心，單元邊界面上的流動參量需要插值格式得到。對于無粘通量，其計算方法采用MUSCL格式，具體形式為：

其中，κ和b是格式調節(jié)因子，1≤b≤，κ為時格式為三階精度；通量F（UL，UR）采用了耗散小的有限差分裂Roe格式。粘性通量采用中心差分進行計算。時間算法部分，由于求解問題為定常問題，采用LU－SGS隱式進行迭代求解。

物面速度邊界條件為無滑移條件，在M＝1.8情況下，溫度邊界條件為絕熱壁條件；在M＝10.6情況下，為等溫壁條件，壁面溫度為294.44K。外邊界為無反射邊界條件。

1.2 數(shù)值驗證

由于本文所選取的鈍錐外形，國外曾經在M＝10.6、迎角為20°時開展過驗證實驗（Cleary［9］），本文也針對該狀態(tài)進行驗證計算，其目的是通過考察程序對CFD的難點——熱流的計算準確性來開展程序驗證。該算例Re數(shù)為1.1×105，流動為層流，計算條件詳見3.1節(jié)。

圖1給出了用駐點熱流歸一化后的熱流分布及其與實驗結果比較，其中φ＝0°子午線為背風區(qū)物面對稱線。比較顯示計算與實驗符合良好。

圖1 計算熱流與Cleary實驗［9］比較Fig.1 The comparison of heat flux between the computation and Cleary’s experiment［9］

驗證結果表明，計算具有較高的數(shù)值模擬精度，可以用來開展進一步的鈍錐分離研究。

2 三維分離判則的數(shù)值分析

雖然在文獻［3－7］中，張涵信院士早就給出了三維定常分離的判則，但由于分析是在分離坐標系下進行，所以除了一些十分規(guī)則的幾何外形外（如方腔），分離判則的直接數(shù)值分析存在一定困難。本節(jié)擬通過對分離判則開展分析，得到不依賴于具體坐標系的、易于數(shù)值驗證的形式。

圖2給出了理論分析所用的曲線坐標系，其中y軸沿分離線，軸正方向為分離線流向，x軸為與y軸垂直的物面曲線，z軸為與x、y正交的曲線坐標；速度u、v、w分別與x、y、z軸對應。值得注意的是，z軸通常并不位于分離流面。再附理論分析所用坐標系與之類似。

圖2 分離坐標系示意圖Fig.2 The schematic of the coordinate system of the separation

文獻［3－7］給出的三維定常分離條件為：流動在物面分離線上滿足

其中下標“SL”表示分離線。

與之對應，三維定常再附條件為：流動在物面再附線上滿足

其中下標“AL”表示再附線。通常計算中觀察到的是第一類再附（（?2u／?x?z）AL＞0）。

從式（4）～式（9）可以看出，分離和再附條件的關鍵變量是：?u／?z、?2u／?x?z和?2w／?z2。下面分別對這些變量與式（4）～（9）中的基本關系進行分析，討論分離判則在一般計算坐標系下的數(shù)值分析方法。

2.1 ?u／?z＝0

為了能夠脫離具體的分離坐標系，考慮到分離坐標系下速度矢量的正交性，首先補充v、w的速度梯度組成正交的空間矢量：（?u／?n，?v／?n，?w／?n）＝?V／?n。如果去掉矢量的法向部分，則可得到?V／?n在物面切平面上的投影。“?u／?z＝0”表示投影矢量在分離線上將與該線相切。為此考慮位于物面切平面的新矢量：?＝?V／?n－n（n·?V／?n），分離條件（4）意味在分離線上?與之相切；新矢量脫離具體的坐標表達，可以方便在計算坐標系下計算。

2.2 ?2u／?x?z＞0與?2u／?x?z＜0

考慮式（4）或式（7），如果面向分離坐標系y軸的正方向，?2u／?x?z＜0意味著在y軸兩側?有指向y軸的分量。這樣考慮2.1節(jié)的分析，如果在分離線附近以?為“速度場”畫表面流線，則流線將向分離線匯聚，所畫出的極限流線或包絡線應該與分離線重合。這個結論可以用來數(shù)值驗證前兩個分離條件的正確性。實際上?對應表面摩擦力確定的矢量場［7］，所以如果分離線由速度場確定的極限流線給出，分離條件式（4）、式（5）也表述為：極限流線與摩擦力線給出的分離線相同。類似地，對于與?2u／?x?z＜0條件對應的再附線，可知用矢量場?畫出的表面流線將由分離線向兩側發(fā)散，再附線是?矢量場的不穩(wěn)定軌線。

2.3 ?2 w／?z2

類似地，考慮矢量（?2V／?n2）w＝（?2u／?n2，?2v／?n2，?2w／?n2）w，其中?／?n表示物面法向的方向導數(shù)、“w”表示物面，則分離坐標系下的?2w／?z2可表示為：n·?2V／?n2，該式脫離具體坐標表達，可由數(shù)值計算結果方便地得到。具體地，根據(jù)n·?2V／?n2大于零或小于零將物面進行“二值化”著色，通過考察分離線和再附線所在區(qū)域開展與理論的對比分析。

3 鈍錐三維分離條件的數(shù)值分析

3.1 計算條件與網格生成

計算外形為15°鈍錐，頭部半徑為1，球心到底部長15。網格為含奇性軸的旋成體網格。具體計算條件見表1。

表1 計算條件與網格Table 1 The condition of computation and grids

其中“＊”表示該狀態(tài)同時用于1.2中的驗證計算。

3.2 鈍錐M＝1.8、Re＝1×105下系列迎角分離拓撲

以下將分三方面對該條件下的鈍錐分離進行分析。其中在圖例中，“S”或“SP”表示鞍點，“N”表示結點，“SL”表示分離線，“AL”表示再附線。

3.2.1 基本分離特性

對于旋成體有迎角表面分離形態(tài)，可以根據(jù)兩個基本特征進行分類：分離類型和多次分離的級數(shù)。圖3給出了不同迎角下極限流線表示的流動表面分離。結果顯示，鈍錐分離由兩個基本組成部分構成：頭部分離和錐體分離。

（1）頭部分離。在球錐頭部，起始分離類型為閉式分離；分離由對稱子午線上的鞍點開始，結束于結點（螺旋點、正?；蛲嘶Y點）。其中10°迎角為背風區(qū)分離剛剛開始的情況（圖3a），頭部分離由位于對稱子午線上的結點開始，通過鞍、結分叉，發(fā)展出“鞍點→分離線→結點”結構的典型閉式分離形態(tài)。隨著迎角的增大，頭部分離發(fā)展出更加復雜的形態(tài)（將在下面討論），起始分離鞍點不在對稱子午線上，流動不對稱，但閉式分離的基本特征仍然保持。

（2）錐體分離。由于錐體分離隨迎角的增大會出現(xiàn)多次分離，所以錐體分離僅討論最外側的主分離。錐體主分離類型在15°、20°和30°為正常點起始的開式分離；當迎角增大到40°時，錐體主分離發(fā)生鞍結點分叉，分離線分裂成兩部分，第一部分一側仍為由正常點起始的開式分離，另一側則由鞍點起始，分離在退化結點結束；第二部分由同一鞍點起始，分離線移植延伸到鈍錐地面，如圖3（d）所示。

圖3 系列迎角下速度極限流線給出的表面分離拓撲Fig.3 The surface separation topology by limiting streamlines at series of angles of attack（α）

從表面流線看，錐體的多次分離由10°迎角的主分離出現(xiàn)開始，逐漸發(fā)展出兩個二次以上分離（20°下的SL2和SL3），進一步發(fā)展出三個二次以上分離（20°下的SL2、SL3和SL5），復雜程度不斷增加，但是新出現(xiàn)的多次分離都呈近似直線分布，體現(xiàn)類錐形流特征。

3.2.2 分離線的奇點分布規(guī)律

文獻［7］給出分離線上奇點分布規(guī)律：分離線上若存在奇點，則為交替分布的鞍點和結點；若分離由奇點開始，則奇點為鞍點；若分離線進入奇點結束，則奇點為結點。下面以鈍錐頭部分離線上奇點分布為例，對分布規(guī)律開展數(shù)值對比分析。

圖4顯示的10°迎角背風區(qū)表面極限流線給出了十分少見的、分離剛開始的結構。頭部流線拓撲由鞍點開始，結束于對稱子午線上的穩(wěn)定結點；鞍點上方、對稱子午線上存在不穩(wěn)定的結點；鞍點向上的奇軸構成一段很小的分離線。鞍點和結點通過鞍點的奇軸連接。

圖4 10°迎角頭部表面分離的奇點分布Fig.4 The singular－point distribution of the surface separation atα＝10°

隨著迎角的增大，在一定范圍頭部分離線奇點逐漸形成穩(wěn)定的拓撲結構，該結構以迎角20°狀態(tài)為代表，如圖5所示。圖中頭部分離線分別起始于對稱子午線上方兩側和下方線上的兩個鞍點，結束于螺旋點。在分離線上，鞍、結點交替分布。進一步的分析（見3.2.3節(jié)）傾向于認為對稱子午線兩側、向下延伸的收斂曲線是第二類再附線［7］，計算顯示該曲線的演化是與分離線結合。值得提出的是，上、下鞍點的奇軸并沒有直接連接，因此不違背結構穩(wěn)定性的要求。當迎角增大到40°時（圖6），頭部表面分離變得非常復雜，出現(xiàn)兩個主控的不穩(wěn)定節(jié)點，但仍然可以通過鞍點的奇軸將奇點連接起來，在連接線上，鞍、結點交替分布；同時底部的分離線仍由鞍點開始，結束于退化的結點；左右兩側分離結構相似但不對稱，尤其是，底部分離線上、原來在對稱子午線上的鞍點已移動到頭部右側。

圖5 20°迎角頭部表面分離的奇點分布Fig.5 The singular－point distribution of the surface separation atα＝20°

圖6 40°迎角頭部表面分離的奇點分布Fig.6 The singular－point distribution of the surface separation atα＝40°

通過數(shù)值分析表明，分離線上的奇點符合文獻［7］給出的拓撲規(guī)律。

3.2.3 分離、再附判則的數(shù)值分析

由分析可知，如果分離線和再附線用計算給出的速度極限流線表示（圖7中用黃色線條顯示），則式（4）、（5）、（7）和（8）表明它們應該與由矢量場?給出的分離和再附線（圖7中用黑色線條顯示）重合。在這樣的思路下，我們將后者作為背景，將速度極限流線給出的分離、再附圖像疊加在背景上來進行比較。圖7以40°迎角為代表，給出了比較結果。結果表明，矢量場?給出的分離和再附線與速度極限流線得到的結構幾乎完全重合：黃色線條完全覆蓋黑色線條，這表明計算與分離條件式（4）、（5）和再附條件式（7）、（9）相符，彼此得到相互印證。

關于條件式（5）和（9），可以利用數(shù)值結果直接計算n·?2V／?n2，并在物面用紅色顯示大于零的區(qū)域、用灰色顯示小于零的區(qū)域，得到結果以30°迎角為代表，如圖8所示。在頭部區(qū)域，閉式分離的分離線始終處于n·?2V／?n2大于零區(qū)域，對稱子午線兩側、向下延伸的收斂曲線由于在靠近頭部位置位于小于零區(qū)域，分析傾向認為該線是第二類再附線；主分離線始終位于n·?2V／?n2大于零區(qū)域，二次以上分離線也基本位于大于零區(qū)域——有時該區(qū)域僅表現(xiàn)為一條窄縫，而分離線恰好從中間穿過；錐體部分的再附線，基本位于n·?2V／?n2小于零區(qū)域，說明流場中再附線主要為第一類再附線。計算結果與理論分析基本上是一致的。

圖7 40°迎角物面速度極限流線與?流線疊加Fig.7 The superimposition of the surface limiting streamlines and that given by?atα＝40°

圖8 30°迎角物面n·?2V ／?n2 分布Fig.8 The surface distribution of n·?2V ／?n2 atα＝30°

3.3 鈍錐M＝10.6、Re＝1.1×105下系列迎角分離拓撲

分兩方面對該條件下鈍錐分離進行分析。

3.3.1 基本分離特性

圖9、圖10表明，在高超M＝10.6、Re＝1.1×105來流條件下，所計算迎角鈍錐背風區(qū)都出現(xiàn)分離?；痉蛛x形態(tài)是由正常點起始的開式分離；從20°左右迎角起，流動在背風區(qū)出現(xiàn)二次分離，隨迎角增大，主分離和二次分離起始位置向頭部移動，分離區(qū)域不斷擴大；分離線、再附線交替分布，將物面分割成不同區(qū)域。值得指出的是，當迎角為40°時，再附線AL12在錐體后部逐漸出現(xiàn)與分離線SL2融合的趨勢，而分離線SL2僅在靠對稱面一側流線向其匯聚（即滿足式（5）），另一側則顯示發(fā)散的特征（即分離坐標系下：（?2u／?x?z）＞0），這種現(xiàn)象的數(shù)學物理規(guī)律尚需進一步分析?？偟膩碚f，高超情況下鈍錐大迎角分離拓撲形態(tài)明顯比M＝1.8下簡單。

圖9 迎角20°的驗證計算與表面速度極限流線Fig.9 The validating computation atα＝20°and the surface limiting streamlines

圖10 迎角30°、40°的表面速度極限流線Fig.10 The surface limiting streamlines atα＝30°and 40°

3.3.2 分離、再附判則的數(shù)值分析

類似地，為了考察分離判則式（4）、（5）和再附判則式（7）、（8），我們將由速度極限流線得到的分離線、再附線與由矢量場?給出的結構疊加，結果顯示在圖11。結果同樣表明，矢量場?給出的分離和再附線與速度極限流線得到的結構幾乎完全重合，從而在數(shù)值上驗證上述條件。

圖11 鈍錐物面速度極限流線與?流線疊加Fig.11 The superimposition of the surface limiting streamlines and that given by?of the blunt cone

但當利用數(shù)值結果直接計算n·?2V／?n2來驗證分離、再附條件式（6）、（9）時，卻發(fā)現(xiàn)n·?2V／?n2大于零、小于零的區(qū)域分布散亂，沒有明顯規(guī)律性，缺乏M＝1.8情況下分離線、再附線與不同區(qū)域間明顯存在的規(guī)律性。這方面的原因可能源于n·?2V／?n2與式（6）、（9）尚存在誤差，或者由于高馬赫數(shù)下密度引起的可壓縮效應，針對該問題的進一步分析正在進行中。

4 結 論

本文采用數(shù)值模擬的手段，對不同馬赫數(shù)、系列迎角下的鈍錐分離拓撲進行了研究，其目的是探討分離判則的數(shù)值分析方法并開展相應分析。計算經過驗證，結果準確可靠，可以用來提供數(shù)據(jù)以開展分析研究。經過計算、分析，有以下結論：

（1）文中通過數(shù)值模擬，詳細給出了普通超聲速和高超范圍鈍錐系列迎角下的分離特性。結果表明，在來流M＝1.8條件下，流動分離模式由頭部閉式分離和錐體開式分離組成，錐體分離隨迎角增大會進一步出現(xiàn)多次分離；在來流M＝10.6條件下，在計算迎角下分離模式僅有開式分離；

（2）文中通過分析文獻［7］中分離、再附判則，構造了不依賴于具體坐標的矢量?和n·?2V／?n2，給出了判則的等價表述。新的表述可以利用數(shù)值模擬結果進行直接計算，更加有利于數(shù)值分析；

（3）文中基于數(shù)值模擬結果和構造的矢量，對分離、再附判則進行了數(shù)值分析，結果表明當M＝1.8時，計算與理論一致；當M＝10.6時，數(shù)值分析與分離判則式（4）、（5）和再附判則式（7）、（8）一致；

（4）對于來流M＝1.8的情況，鈍錐頭部會出現(xiàn)復雜的鞍、結點組合。雖然這些奇點拓撲結構復雜，但其沿分離線等的分布符合文獻［7］中給出的拓撲規(guī)律。

［1］PRANDTL.über Flüssigkeitsbewegung bei sehr kleineir Reibung［M］／／Verh.III，Int.Math.Kongr.Heidelberg，1904：485－491.

［2］MOORE F K.Boundary layer research［M］.GORTTER H G，Edited.Berting，Springer Vertag，1958：296－310.

［3］張涵信.二維粘性不可壓縮流動的通用分離判據(jù)［J］.力學學報，1983，6：559－570.

［4］張涵信.三維定常粘性流動的分離條件及分離線附近流動的性狀［J］.空氣動力學學報，1985，1：1－12.

［5］張涵信.分離流的某些進展［J］.航空學報，1985，6（4）：301－312.

［6］ZHANG H X.Numerical simulation of three dimensional separated flow field and application of topological theory［J］.AdvancesinScienceofChina，Mechanics，1991，1：59－80.

［7］張涵信.分離流與渦運動的結構分析［M］.國防工業(yè)出版社，2005.

［8］張涵信，張樹海，田浩，等.三維可壓縮非定常流的壁面分離判據(jù)及其分離線附近的流動形態(tài)［J］.空氣動力學學報，2012，30（4）：421－430.

［9］CLEARY J W.Effects of angle of attack and bluntness on laminar heating－rate distributions of a 15°cone at a Mach number of 10.6［R］.NASA TN D－5450，1969.