尖拱細(xì)長(zhǎng)旋成體低速大攻角繞流數(shù)值模擬

張浩南，李國(guó)輝，徐 宇

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

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

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

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

1 數(shù)值模擬方法

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

(1)

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

(2)

Gω-Yω+Dω+Sω

(3)

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

(4)

(5)

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

(6)

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

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

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

圖1 截面選取以及計(jì)算域示意圖

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

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

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

圖3 平均x渦量等值線

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

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

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

圖5 速度矢量分布示意圖

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

圖6 細(xì)長(zhǎng)體繞流拓?fù)浣Y(jié)構(gòu)示意圖

圖7 切向速度曲線

圖8 空間鞍點(diǎn)的縱向位置與橫向位置曲線

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

圖9 鞍點(diǎn)位置示意圖

圖10 臨界狀態(tài)拓?fù)浣Y(jié)構(gòu)示意圖

4 結(jié)論

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

2) 確定截面繞流中奇點(diǎn)的精確數(shù)量及位置對(duì)截面繞流拓?fù)浣Y(jié)構(gòu)的驗(yàn)證，結(jié)果與理論分析一致。

3) 在數(shù)值模擬中通過(guò)截面繞流中特殊奇點(diǎn)的位置變化可觀察到對(duì)稱(chēng)渦流型的臨界流動(dòng)狀態(tài)。