張浩南，李國輝，徐 宇

(空軍航空大學(xué) 航空作戰(zhàn)勤務(wù)學(xué)院，長春 130000)

細長體繞流在大攻角下會出現(xiàn)一種現(xiàn)象，即旋渦對稱流動發(fā)展為非對稱流動。當(dāng)旋渦表現(xiàn)為非對稱流動時，細長體會產(chǎn)生一個側(cè)向力，進而產(chǎn)生偏航力矩。Zilliac等[1]進行的實驗研究顯示了四種不同的流動模式：α<30°時的穩(wěn)定對稱流動、30°<α<50°時的不穩(wěn)定對稱流動、50°<α<65°時的雙穩(wěn)態(tài)模式以及α>65°時的類卡門渦街(或隨機尾跡)流態(tài)。這種流動狀態(tài)轉(zhuǎn)換的另一種表現(xiàn)為側(cè)向力的轉(zhuǎn)換[1-3]，了解這種現(xiàn)象的原理將對于航空航天應(yīng)用具有重要價值,例如對于導(dǎo)彈和高機動性飛機而言，可以通過流動控制進而達到減阻的目的。大攻角下非對稱流發(fā)展的研究是一項極具挑戰(zhàn)性的工作，由于缺乏對流動機理的認識，在實驗和數(shù)值模擬時常常會遇到很多困難，包括模型缺陷、來流條件的選擇、和測量技術(shù)的局限性等等。為研究這種現(xiàn)象，眾多學(xué)者對影響流動非對稱性的因素進行了研究分析，例如馬赫數(shù)[4]、雷諾數(shù)[5]、鈍度[6]和細長比[7]等等。由于在傳統(tǒng)風(fēng)洞實驗研究中，實驗數(shù)據(jù)重現(xiàn)性很差,隨著計算機技術(shù)的發(fā)展，數(shù)值模擬方法開始逐漸顯示出它的優(yōu)勢。大多數(shù)數(shù)值模擬使用RANS方法進行研究，其中引入了人工擾動[8]或幾何缺陷[9]是產(chǎn)生流動非對稱所必需的；同樣，部分學(xué)者使用LES方法研究[10]也必須添加人工擾動，而也有學(xué)者則認為，在沒有任何擾動以及缺陷的情況下也可以得到非對稱的旋渦[11]。

本文的目的是使用RANS方法研究細長體在大攻角下尾渦的發(fā)展，而既不引入幾何缺陷，也不引入人為擾動。在五個攻角α=30°，α=40°，α=50°，α=55°和α=60°上進行研究，使用壓力分布分析流場沿細長體的截面法向力和側(cè)面力，此外，提出一種驗證細長旋成體截面繞流拓撲結(jié)構(gòu)的方法，結(jié)果證明，數(shù)值仿真結(jié)果與理論分析保持一致[12]。

本文的目的是使用RANS方法研究細長體在大攻角下尾渦的發(fā)展，既不引入幾何缺陷，也不引入人工擾動。在5個攻角α=30°，α=40°，α=50°，α=55°和α=60°下進行研究，使用壓力分布分析流場沿細長體的截面法向力和側(cè)面力，此外，提出一種驗證細長旋成體截面繞流拓撲結(jié)構(gòu)的方法。

1 數(shù)值模擬方法

控制方程為三維粘性N-S方程，其守恒通量形式為

(1)

式中：U為非定常項；E、F、G是為無粘項；Eμ、Fμ、Gμ為粘性項，本文對N-S方程進行時間平均進而模擬湍流流動。通過Reynolds平均法對流動控制方程式(1)進行時均處理得到的雷諾平均N-S方程(RANS)為封閉RANS方程，選用SSTk-ω湍流模型，該模型對網(wǎng)格適應(yīng)性較好[13]，其湍流動能k和比耗散率ω的輸運方程為

(2)

Gω-Yω+Dω+Sω

(3)

式中：Gk表示湍流的動能；Gω為ω的生成項；Γk和Γω分別代表k與ω的有效擴散項；Yk和Yω分別代表k與ω的發(fā)散項；Dω代表正交發(fā)散項，有效擴散方程為

(4)

(5)

式中：σk和σω是方程的湍流能量普朗特數(shù)；μt為湍流黏度，計算公式如下：

(6)

其中S為旋率，系數(shù)α*使得湍流粘度產(chǎn)生低雷諾數(shù)修正，a1=0.31。該模型在近壁區(qū)等價于k-ω模型，在遠離壁面的區(qū)域自動轉(zhuǎn)換為標準的k-ε模型。對控制方程采用有限體積法進行求解，對于控制方程的空間離散方式，擴散項采用二階中心差分格式，對流項采用三階MUSCL格式，速度與壓力的耦合方式選用PISO算法。

2 網(wǎng)格及邊界條件

研究中使用的幾何模型由細長比為3的尖拱狀前身和細長比為10的圓柱狀后身組成。模型的總長度為26 cm，直徑D=2 cm。分別在尖拱狀前身的11個截面(OG段)和圓柱狀后身的9個截面(CY段)收集數(shù)據(jù)以便于檢查分析[如圖1(a)所示]，截面位置與文獻[14]所提供的位置相匹配，尖拱段部分截面用以監(jiān)視初始旋渦的發(fā)展。在選擇計算域時，細長體前端延長至6D處，后端與徑向均延長至15D處[圖1(b)]。整體網(wǎng)格選用結(jié)構(gòu)網(wǎng)格，第一層網(wǎng)格高度為 0.000 1D，對模型頭部網(wǎng)格進行了細化，總網(wǎng)格數(shù)為150萬(圖2)。計算采用Fluent軟件，邊界條件為速度入口與壓力出口，壁面選用無滑移絕熱壁面，此外，自由來流速度固定為0.2馬赫，湍流度固定為0.2%。

圖1 截面選取以及計算域示意圖

3 模擬結(jié)果與分析

在本節(jié)中，通過繪制x渦量等值線圖，定性地評價了計算結(jié)果；通過截面?zhèn)认蛄头ㄏ蛄ο禂?shù)定量地分析了計算結(jié)果；并且根據(jù)截面上奇點的數(shù)量與位置，對細長體的截面拓撲結(jié)構(gòu)進行了驗證。

圖3顯示了所有計算攻角下的x渦量等值線?？梢钥闯觯悍菍ΨQ性沿軸向發(fā)展，并隨著攻角的增加而增加，在30°≤α<50°時旋渦基本保持對稱渦流型，非對稱性的發(fā)展在α=50°時開始，并且起初較為微弱，在尖拱段時相對較低，到圓柱段時加強。在α=60°時旋渦已發(fā)展為高度非對稱的，非對稱性在尖拱段末端已經(jīng)開始顯現(xiàn)。

圖3 平均x渦量等值線

側(cè)向力和法向力系數(shù)(分別為Cy和Cz)相對于攻角的變化如圖4所示?？梢钥闯?，Cz隨著攻角的增加而增加，這是由于隨著攻角的增加，背風(fēng)側(cè)的吸力會越來越大，因此Cz會有所增加。而且，隨著攻角的增加，Cz最大值的位置沿軸向向下游移動。在α=30°時，Cz最大值出現(xiàn)在OG(6)截面；在α=40°到50°時，Cz最大值的位置移到OG(7)截面；在α=60°時，Cz最大值出現(xiàn)在OG(8)截面，并且這也是所有計算攻角中的最大值。在所有計算攻角下，當(dāng)達Cz到最大幅值后均沿軸向減小。另一方面，在30°≤α≤40°時，整個模型Cy值保持在較低水平，這表明尾流旋渦仍然是對稱的。但是，Cy在α=50°時沿軸向顯著增加，在OG(11)截面達到最大值，然后，側(cè)向力沿軸向在負向和正向之間波動。當(dāng)攻角增加到α=60°時，Cy達到所有攻角下的最大幅值。與Cz相反，Cy最大值的位置隨著攻角增加沿軸向向前移動，在α=60°時Cy最大值與Cz最大值發(fā)生在相同位置，即OG(8)截面。α=55°和α=60°的曲線趨勢類似于α=50°，即Cy方向不斷在正向與負向之間改變，這表明尾流旋渦具有強烈的非對稱性，非對稱性隨著攻角的增加而增強，并沿軸向從圓柱段移動到尖拱段。

圖4 側(cè)向力和法向力系數(shù)曲線

圖5為攻角α=30°時OG(3)截面、OG(6)截面以及CY(1)截面的速度矢量分布示意圖。從圖中可以看出，在OG(3)截面,流動還未產(chǎn)生分離,截面流型為附著結(jié)構(gòu)(圖5(a))；沿軸向往下游發(fā)展，這種附著結(jié)構(gòu)開始發(fā)生變化，細長旋成體背渦出現(xiàn)了流動分離現(xiàn)象，到OG(6)截面，在細長旋成體左右兩側(cè)產(chǎn)生一對對稱旋渦，由于此時對稱旋渦剛開始出現(xiàn)，流動的分離還較為微弱，卷入旋渦中的渦量較少，使得產(chǎn)生的旋渦為體積小且強度弱的一對對稱旋渦，故而在卷成旋渦之后很快就在背風(fēng)側(cè)的對稱面附近發(fā)生再附(圖5(b))，當(dāng)旋渦處于對稱流態(tài)時，兩個頭渦之間的誘導(dǎo)作用是相互平衡的，并且他們之間相距較遠，相互擠壓的作用微乎其微，這就保證了這種對稱渦流態(tài)是一種穩(wěn)定的平衡狀態(tài)；隨著流動沿軸向繼續(xù)往下游發(fā)展，分離現(xiàn)象進一步增強，分離剪切層不斷向兩個頭渦輸入渦量，使得旋渦強度增強體積增大，并且相互靠近進而伴隨著擠壓現(xiàn)象的發(fā)生(圖5(c))，此時旋渦仍保持對稱渦流態(tài)，兩個頭渦之間的誘導(dǎo)作用依舊是相互平衡的，但是由于他們不斷靠近，已經(jīng)產(chǎn)生了較為強烈的擠壓作用，使得此時的對稱渦流態(tài)成為一種不穩(wěn)定的平衡狀態(tài)。

圖5 速度矢量分布示意圖

為探究細長旋成體對稱渦流態(tài)的演變，基于Lowson和Ponton[15]給出的細長旋成體繞流截面流型拓撲結(jié)構(gòu)(圖6)，為獲得繞流截面流型拓撲結(jié)構(gòu)中鞍點的精確位置，圖7給出了α=30°攻角下OG(3)截面、OG(6)截面以及CY(1)截面上壁面切向速度隨方位角θ的變化曲線，進而確定近壁面奇點位置，除此之外，通過尋求速度矢量與截面法向量夾角的0點來確定空間鞍點的精確位置，圖8(a)與圖8(b)分別為該攻角下空間奇點沿軸向不同截面縱向與橫向位置變化曲線，由于細長旋成體頭部為尖拱型，不同截面對應(yīng)的模型半徑不同，為準確表示空間奇點與壁面的距離，在圖8(a)中選用空間奇點縱向坐標與該截面模型半徑的差值作為y軸坐標值，從圖8可以看出，隨著截面沿軸向往下游發(fā)展，空間奇點逐漸遠離細長體壁面，線性程度較高，且一直處于背風(fēng)側(cè)對稱面附近很小的范圍之內(nèi)。通過上述分析可得，對于OG(3)截面，切向速度存在兩個0點，分別存在于θ=0°和θ=180°處，而由于壁面是無法穿過的，所以細長旋成體壁面的徑向速度始終為0，進而可以得出，在θ=0°和θ=180°處存在兩個奇點，并且通過壁面切向速度的流動方向可以判斷θ=0°處為再附型半鞍點Rh1，θ=180°處為分離型半鞍點Sh1，空間中無奇點；同理，對于OG(6)截面，除去θ=0°和θ=180°處的兩個奇點，在背風(fēng)側(cè)對稱面兩側(cè)附近存在兩個位置對稱的再附型半鞍點Rh2和Rh3(θ=140°與θ=220°)，而在距離背風(fēng)側(cè)對稱面兩側(cè)更遠一些的位置存在兩個位置對稱的分離型半鞍Sh2和Sh3(θ=90°與θ=270°)，空間中無奇點，在整個細長旋成體壁面上，對稱面兩側(cè)附近的兩個再附型半鞍點在旋成體壁面上代表兩條再附線；對于CY(1)截面，從圖8可以看出，相比于OG(6)截面，θ=180°處的奇點由分離型半鞍點轉(zhuǎn)變?yōu)樵俑叫桶氚包cRh2，并且細長體背風(fēng)側(cè)對稱面左右兩側(cè)的一對再附型半鞍點也不復(fù)存在，這代表細長旋成體壁面上將不再有再附線，此時空間中存在一個奇點S1。不難得出，OG(3)截面、OG(6)截面以及CY(1)截面的對稱渦流型分別對應(yīng)圖6所示的3種流動拓撲機構(gòu)。

圖6 細長體繞流拓撲結(jié)構(gòu)示意圖

圖7 切向速度曲線

圖8 空間鞍點的縱向位置與橫向位置曲線

在上述分析中，驗證了兩渦對稱流型的兩種拓撲結(jié)構(gòu)，它們的主要區(qū)別在于兩個再附型半鞍點的消失以及θ=180°處奇點類型的轉(zhuǎn)變，為此對這3個奇點的演變進行進一步研究，由于對稱渦流型的轉(zhuǎn)變是一種漸變的過程，對OG(6)截面與CY(1)截面中間數(shù)個截面的壁面切向速度進行了分析，圖9為其中部分截面上3個奇點的位置示意圖，從圖中可以看出細長體背風(fēng)側(cè)對稱面左右兩側(cè)的一對再附型半鞍點Rh2和Rh3會不斷的靠近θ=180°處的分離型半鞍點Sh1，直至3個奇點重合，轉(zhuǎn)變?yōu)楦唠A奇點O，此時在θ=180°附近已再無任何奇點存在，稱這個狀態(tài)為臨界流動狀態(tài)，從上述分析可以得出，臨界流動狀態(tài)的高階奇點O由圖6(b)中Rh2、Rh3以及Sh1重合而成，此時壁面存在兩個分離型半鞍Sh1和Sh2，一個再附型半鞍點Rh1和一個高階奇點O，且空間中無奇點存在，拓撲結(jié)構(gòu)如圖10所示，與文獻[12]提出的臨界狀態(tài)保持一致。

圖9 鞍點位置示意圖

圖10 臨界狀態(tài)拓撲結(jié)構(gòu)示意圖

4 結(jié)論

1) 在無任何擾動與不規(guī)則性的前提下，可以得到非對稱的旋渦，且非對稱性隨著攻角的增加而增強，并沿軸向逐漸增強。

2) 確定截面繞流中奇點的精確數(shù)量及位置對截面繞流拓撲結(jié)構(gòu)的驗證，結(jié)果與理論分析一致。

3) 在數(shù)值模擬中通過截面繞流中特殊奇點的位置變化可觀察到對稱渦流型的臨界流動狀態(tài)。