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

李 沁，張涵信

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

0 引 言

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

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

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

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

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

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

1 數(shù)值方法

1.1 計算方法

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

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

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

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

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

1.2 數(shù)值驗證

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

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

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

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

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

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

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

圖2 分離坐標(biāo)系示意圖Fig.2 The schematic of the coordinate system of the separation

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

其中下標(biāo)“SL”表示分離線。

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

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

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

2.1 ?u／?z＝0

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

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

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

2.3 ?2 w／?z2

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

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

3.1 計算條件與網(wǎng)格生成

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

表1 計算條件與網(wǎng)格Table 1 The condition of computation and grids

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

3.2 鈍錐M＝1.8、Re＝1×105下系列迎角分離拓?fù)?/h3>

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

3.2.1 基本分離特性

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

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

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

圖3 系列迎角下速度極限流線給出的表面分離拓?fù)銯ig.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），復(fù)雜程度不斷增加，但是新出現(xiàn)的多次分離都呈近似直線分布，體現(xiàn)類錐形流特征。

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

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

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

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

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

圖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ù)值分析表明，分離線上的奇點符合文獻(xiàn)［7］給出的拓?fù)湟?guī)律。

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

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

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

圖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下系列迎角分離拓?fù)?/h3>

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

3.3.1 基本分離特性

圖9、圖10表明，在高超M＝10.6、Re＝1.1×105來流條件下，所計算迎角鈍錐背風(fēng)區(qū)都出現(xiàn)分離?；痉蛛x形態(tài)是由正常點起始的開式分離；從20°左右迎角起，流動在背風(fēng)區(qū)出現(xiàn)二次分離，隨迎角增大，主分離和二次分離起始位置向頭部移動，分離區(qū)域不斷擴大；分離線、再附線交替分布，將物面分割成不同區(qū)域。值得指出的是，當(dāng)迎角為40°時，再附線AL12在錐體后部逐漸出現(xiàn)與分離線SL2融合的趨勢，而分離線SL2僅在靠對稱面一側(cè)流線向其匯聚（即滿足式（5）），另一側(cè)則顯示發(fā)散的特征（即分離坐標(biāo)系下：（?2u／?x?z）＞0），這種現(xiàn)象的數(shù)學(xué)物理規(guī)律尚需進一步分析?？偟膩碚f，高超情況下鈍錐大迎角分離拓?fù)湫螒B(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），我們將由速度極限流線得到的分離線、再附線與由矢量場?給出的結(jié)構(gòu)疊加，結(jié)果顯示在圖11。結(jié)果同樣表明，矢量場?給出的分離和再附線與速度極限流線得到的結(jié)構(gòu)幾乎完全重合，從而在數(shù)值上驗證上述條件。

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

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

4 結(jié) 論

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

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

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

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

（4）對于來流M＝1.8的情況，鈍錐頭部會出現(xiàn)復(fù)雜的鞍、結(jié)點組合。雖然這些奇點拓?fù)浣Y(jié)構(gòu)復(fù)雜，但其沿分離線等的分布符合文獻(xiàn)［7］中給出的拓?fù)湟?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］.力學(xué)學(xué)報，1983，6：559－570.

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

［5］張涵信.分離流的某些進展［J］.航空學(xué)報，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］張涵信.分離流與渦運動的結(jié)構(gòu)分析［M］.國防工業(yè)出版社，2005.

［8］張涵信，張樹海，田浩，等.三維可壓縮非定常流的壁面分離判據(jù)及其分離線附近的流動形態(tài)［J］.空氣動力學(xué)學(xué)報，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.