Pω增強(qiáng)型k－ω湍流模型在三角翼旋渦流動(dòng)的應(yīng)用

張冬云，李喜樂，楊 永，張 強(qiáng)

（1.中國(guó)商用飛機(jī)有限責(zé)任公司上海飛機(jī)設(shè)計(jì)研究院，上海 201210；2.中國(guó)航天空氣動(dòng)力技術(shù)研究院，北京 100074；3.西北工業(yè)大學(xué)翼型、葉柵空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西西安 710072）

Pω增強(qiáng)型k－ω湍流模型在三角翼旋渦流動(dòng)的應(yīng)用

張冬云1，李喜樂2，＊，楊 永3，張 強(qiáng)3

（1.中國(guó)商用飛機(jī)有限責(zé)任公司上海飛機(jī)設(shè)計(jì)研究院，上海 201210；2.中國(guó)航天空氣動(dòng)力技術(shù)研究院，北京 100074；3.西北工業(yè)大學(xué)翼型、葉柵空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西西安 710072）

在三角翼旋渦繞流數(shù)值模擬中，標(biāo)準(zhǔn)Wilcoxk－ω湍流模型生成項(xiàng)未考慮旋度的影響而導(dǎo)致預(yù)測(cè)的旋渦強(qiáng)度較弱。通過引入探測(cè)因子區(qū)分剪切層和渦核，在旋渦流動(dòng)的高旋度區(qū)域增加ω方程生成項(xiàng)的方法，基于結(jié)構(gòu)化網(wǎng)格上的RANS求解器，加入了Pω增強(qiáng)型k－ω湍流模型，對(duì)繞尖前緣三角翼亞聲速和跨聲速旋渦流場(chǎng)進(jìn)行了數(shù)值模擬。計(jì)算結(jié)果與NASA的NTF風(fēng)洞和DLR的DNW－TWG風(fēng)洞試驗(yàn)數(shù)據(jù)進(jìn)行了對(duì)比分析，結(jié)果表明：不論在亞聲速還是跨聲速自由來流條件下，Pω增強(qiáng)型k－ω湍流模型計(jì)算的壓力分布、渦破裂位置均與試驗(yàn)數(shù)據(jù)吻合良好，準(zhǔn)確地預(yù)測(cè)出了三角翼上翼面的主渦、二次渦結(jié)構(gòu)，特別是跨聲速條件下激波干擾導(dǎo)致的渦破裂的臨界迎角及渦破裂位置，表明Pω增強(qiáng)型k－ω湍流模型在繞三角翼旋渦流動(dòng)數(shù)值模擬中具有良好的適用性。

k－ω湍流模型；生成項(xiàng)；三角翼；跨聲速；臨界迎角；渦破裂

0 引 言

現(xiàn)代先進(jìn)戰(zhàn)斗機(jī)常采用三角翼布局，繞三角翼大迎角旋渦流動(dòng)包含復(fù)雜的流動(dòng)現(xiàn)象［1］，它們對(duì)戰(zhàn)斗機(jī)的空氣動(dòng)力特性、操縱性與穩(wěn)定性、氣動(dòng)彈性和結(jié)構(gòu)動(dòng)力學(xué)特性均會(huì)產(chǎn)生重要的影響［2］。因此，近些年科研人員對(duì)繞三角翼旋渦流場(chǎng)進(jìn)行了大量的實(shí)驗(yàn)［3－8］和數(shù)值模擬［9－14］研究。研究表明在大迎角下，尖前緣三角翼上翼面的旋渦流場(chǎng)包含主分離渦、二次渦，甚至三次渦［15］。隨著迎角的增加，分離渦會(huì)發(fā)生渦破裂現(xiàn)象。亞聲速情況下，渦破裂位置隨迎角增加逐漸從三角翼后緣向上游移動(dòng)［16－17］；與亞聲速渦破裂現(xiàn)象不同的是，跨聲速條件下，由于激波對(duì)主渦的干擾作用，在某一臨界迎角下，三角翼上會(huì)發(fā)生渦破裂位置突然前移的現(xiàn)象，使得繞三角翼的跨聲速旋渦流動(dòng)變得更為復(fù)雜［11－12，14］。

三角翼上翼面渦破裂的發(fā)生與分離渦的許多參數(shù)有關(guān)（例如基于臨界態(tài)理論和旋渦穩(wěn)定性理論的羅斯比數(shù)［18－19］，旋轉(zhuǎn)速度、軸向速度和逆壓梯度等內(nèi)部參數(shù)的關(guān)系等）。因此，在數(shù)值模擬研究中，為了準(zhǔn)確預(yù)測(cè)繞三角翼旋渦流動(dòng)中的渦破裂現(xiàn)象，必須準(zhǔn)確預(yù)測(cè)出分離渦結(jié)構(gòu)。在基于求解RANS方程的三角翼旋渦繞流數(shù)值模擬研究中，湍流模型在對(duì)渦結(jié)構(gòu)的捕捉中發(fā)揮著至關(guān)重要的作用，進(jìn)而影響著對(duì)渦破裂的預(yù)測(cè)。標(biāo)準(zhǔn)的一方程和兩方程湍流模型都存在在渦核區(qū)域預(yù)測(cè)的湍流黏性水平過高的問題，導(dǎo)致分離渦較弱，耗散過快，難以準(zhǔn)確預(yù)測(cè)出分離渦結(jié)構(gòu)和渦破裂現(xiàn)象。本文采用Brandsma等［2］提出的方法，在流動(dòng)的高旋度區(qū)域增加ω（耗散率）方程生成項(xiàng)，在基于結(jié)構(gòu)化網(wǎng)格的RANS方程求解器上，加入Pω增強(qiáng)型k－ω湍流模型，對(duì)65°后掠角尖前緣三角翼的亞聲速和跨聲速旋渦流場(chǎng)進(jìn)行了數(shù)值模擬。計(jì)算了不同狀態(tài)下三角翼上翼面的物面壓力分布、渦破裂位置，分析了流場(chǎng)結(jié)構(gòu)，通過與試驗(yàn)數(shù)據(jù)和標(biāo)準(zhǔn)Wilcoxk－ω湍流模型［20］計(jì)算結(jié)果對(duì)比分析，考核了Pω增強(qiáng)型k－ω湍流模型對(duì)繞三角翼旋渦流動(dòng)的數(shù)值模擬能力。

1 數(shù)值模擬方法

1.1流體力學(xué)控制方程

繞三角翼旋渦流場(chǎng)的流動(dòng)控制方程采用RANS方程，其積分形式為：

其中，V為控制體體積，Q為守恒變量矢量，σ為控制體表面，F(xiàn)為通過表面σ的凈通量矢量，包含黏性項(xiàng)和無黏項(xiàng)，n為表面σ的單位外法向矢量。

流動(dòng)控制方程空間離散采用二階有限體積方法。無黏通量項(xiàng)的離散采用迎風(fēng)型通量差分分裂格式（Flux－Difference Splitting）；而對(duì)于黏性項(xiàng)，采用二階中心格式進(jìn)行離散。定常計(jì)算的時(shí)間推進(jìn)采用隱式近似因子化法（Approximate Factorization Method），在每個(gè)空間方向上獨(dú)立進(jìn)行隱式求解運(yùn)算。

1.2 標(biāo)準(zhǔn)Wilcoxk－ω湍流模型及其改進(jìn)

首先，應(yīng)變率張量和旋轉(zhuǎn)張量分別定義為：

一般情況下，在剪切層內(nèi)，速度梯度主要為物面法向方向上的梯度，因此，應(yīng)變率張量和旋轉(zhuǎn)張量大致相等。但是，在接近渦核的區(qū)域，流動(dòng)基本處于旋轉(zhuǎn)狀態(tài)，旋轉(zhuǎn)張量將會(huì)增大。

式（4）從數(shù)學(xué)上表示了能量從平均流動(dòng)到脈動(dòng)速度流場(chǎng)的傳遞關(guān)系，這是由于在渦團(tuán)拉伸過程中，平均速度梯度和雷諾應(yīng)力相互作用引起的。其中，只有速度梯度的對(duì)稱分量（即應(yīng)變率張量珚Sij）和雷諾應(yīng)力張量的各向異性分量aij對(duì)湍動(dòng)能有貢獻(xiàn)。因此，式（4）可以寫為：

可以看出，湍動(dòng)能的產(chǎn)生與應(yīng)變率張量是成比例的。對(duì)于三角翼旋渦流動(dòng)，式（5）表明湍流大多出現(xiàn)在剪切層和其周圍的流場(chǎng)，不會(huì)出現(xiàn)在高速旋轉(zhuǎn)的渦核區(qū)域。因此，可以認(rèn)為渦核區(qū)域的湍流度較低，甚至可以認(rèn)為是層流。

標(biāo)準(zhǔn)Wilcoxk－ω湍流模型是Wilcox［20］基于Boussinesq假設(shè)提出的兩方程湍流模型，該模型根據(jù)湍動(dòng)能k和湍動(dòng)能的比耗散率ω來計(jì)算渦黏系數(shù)，以封閉RANS方程。k－ω模型的湍流渦黏系數(shù)定義如下：

為了理解k－ω模型在三角翼旋渦流動(dòng)中的作用，下面對(duì)k－ω模型的生成項(xiàng)進(jìn)行分析。湍動(dòng)能的生成項(xiàng)見式（4），與之相關(guān)的比耗散率方程生成項(xiàng)為：

由于該模型采用Boussinesq假設(shè)計(jì)算雷諾應(yīng)力，湍動(dòng)能的生成項(xiàng)變?yōu)椋?/p>

從k和ω方程的生成項(xiàng)定義中，可以明顯地看出，它只依賴于平均流動(dòng)的應(yīng)變率而沒有考慮旋度的影響。然而渦核區(qū)域?yàn)楦咝葏^(qū)域，這種過度的簡(jiǎn)化導(dǎo)致渦核區(qū)域黏性過大，使得預(yù)測(cè)的旋渦耗散偏大，導(dǎo)致旋渦較弱，并且不能持續(xù)，迅速耗散。

為了使標(biāo)準(zhǔn)Wilcoxk－ω模型可以準(zhǔn)確預(yù)測(cè)渦核區(qū)域的流動(dòng)結(jié)構(gòu)，Brandsma等［2］提出了兩種不同的旋轉(zhuǎn)修正方法。這些方法建議控制湍動(dòng)能的生成項(xiàng)，進(jìn)而控制渦核區(qū)域的湍流渦黏系數(shù)。第一種方法是直接把耗散項(xiàng)作為限制器來限制k方程的生成項(xiàng)，第二種方法是在流動(dòng)的高旋度區(qū)域增加ω方程的生成項(xiàng)。本文采用第二種旋轉(zhuǎn)修正方法，為了在合理的區(qū)域應(yīng)用旋轉(zhuǎn)修正，引入一個(gè)探測(cè)因子來區(qū)分剪切層和渦核。這個(gè)因子定義為平均應(yīng)變率張量和平均旋轉(zhuǎn)張量大小的比值：

在剪切層內(nèi)，應(yīng)變率張量和旋轉(zhuǎn)張量的大小幾乎相等，r≈1，在渦核區(qū)域，流動(dòng)接近旋轉(zhuǎn)狀態(tài)，r1。修正后的生成項(xiàng)定義為：

Pωnew等于ω方程的生成項(xiàng)除以min （r2，1），即：

采用這一修正后，加強(qiáng)了ω方程的生成項(xiàng)，使ω增加，也就更快的耗散k，減小了渦核區(qū)域的湍流渦黏性，渦核就不會(huì)被過快耗散，改進(jìn)之后的模型稱之為Pω增強(qiáng)型k－ω湍流模型。

2 計(jì)算模型

幾何模型采用VFE－2實(shí)驗(yàn)中Chu和Luckring［4］在美國(guó)國(guó)家航空航天局（NASA）國(guó)家跨聲速實(shí)驗(yàn)室（National Transonic Facility，NTF）風(fēng)洞中使用的尖前緣三角翼模型。模型幾何尺寸如圖1所示，三角翼的根弦長(zhǎng)Cr＝0.980m，翼展b＝0.914m，前緣后掠角φ＝65°，機(jī)翼面積S＝0.448m2，展弦比AR＝1.86。

3 數(shù)值模擬結(jié)果分析

首先研究計(jì)算網(wǎng)格對(duì)計(jì)算結(jié)果的影響，然后選用合適的計(jì)算網(wǎng)格，與標(biāo)準(zhǔn)Wilcoxk－ω湍流模型計(jì)算結(jié)果進(jìn)行對(duì)比分析，驗(yàn)證Pω增強(qiáng)型k－ω湍流模型的適用性。本文對(duì)VFE－2實(shí)驗(yàn)中亞聲速和跨聲速條件下不同的流動(dòng)狀態(tài)進(jìn)行了數(shù)值模擬。為方便敘述，各計(jì)算狀態(tài)用表1中對(duì)應(yīng)編號(hào)表示。

圖1 65°后掠三角翼模型［12］Fig.1 Geometry of the 65°swept leading edge delta wing

表1 計(jì)算狀態(tài)對(duì)比Table 1 CFD cases chosen for summary comparison

3.1 網(wǎng)格影響分析

本文網(wǎng)格采用了H－O型網(wǎng)格拓?fù)?。為了分析網(wǎng)格影響，相對(duì)粗網(wǎng)格，對(duì)三角翼前、后緣，支架頂點(diǎn)和物面法向區(qū)域網(wǎng)格進(jìn)行加密處理。圖2為密網(wǎng)格的拓?fù)浣Y(jié)構(gòu)、物面網(wǎng)格和空間網(wǎng)格切片，部分網(wǎng)格分布參數(shù)見表2。

圖2 三角翼拓?fù)浣Y(jié)構(gòu)及計(jì)算網(wǎng)格示意圖（fine grid）Fig.2 Mesh and grid topology around the delta wing

表2 計(jì)算網(wǎng)格量對(duì)比Table 2 Comparison of computational grid size

進(jìn)行網(wǎng)格影響分析選用的計(jì)算狀態(tài)為：Ma∞＝0.85，α＝23.0°，Re∞＝6×106。在文獻(xiàn)［12］中，NASA NTF風(fēng)洞試驗(yàn)數(shù)據(jù)顯示三角翼上翼面并未發(fā)生主渦破裂現(xiàn)象，該迎角與發(fā)生渦破裂的臨界迎角接近，而文獻(xiàn)［12］的CFD方法預(yù)測(cè)出了渦破裂現(xiàn)象。

圖3為選用Pω增強(qiáng)型k－ω湍流模型，粗、密網(wǎng)格的CFD計(jì)算結(jié)果對(duì)比，表明粗、密網(wǎng)格均預(yù)測(cè)出了三角翼上翼面的主渦峰值和二次渦峰值，前三個(gè)站位處，密網(wǎng)格預(yù)測(cè)的主渦和二次渦峰值均高于粗網(wǎng)格，更接近于試驗(yàn)數(shù)據(jù)。x／cr＝0.8站位處，粗網(wǎng)格預(yù)測(cè)的主渦峰值及站位與試驗(yàn)結(jié)果吻合良好，但沿展向外側(cè)壓力分布均低于試驗(yàn)數(shù)據(jù)，二次渦峰值不明顯；密網(wǎng)格預(yù)測(cè)的主渦和二次渦峰值均高于試驗(yàn)數(shù)據(jù)。x／cr＝0.95站位處，粗網(wǎng)格預(yù)測(cè)的主渦峰值與試驗(yàn)數(shù)據(jù)吻合較好，但主渦形態(tài)較為平滑，密網(wǎng)格預(yù)測(cè)的主渦峰值偏高，但主渦形態(tài)與試驗(yàn)數(shù)據(jù)接近（即沿展向方向，在主渦外側(cè)的壓力系數(shù)曲線更為陡峭）。

上述對(duì)比表明，兩套網(wǎng)格均較為準(zhǔn)確地預(yù)測(cè)出了主渦和二次渦結(jié)構(gòu)，網(wǎng)格密度達(dá)到研究要求。為了更為準(zhǔn)確地預(yù)測(cè)繞三角翼旋渦流動(dòng)中的分離渦、剪切層、激波與渦的干擾現(xiàn)象等，本文后續(xù)的研究均采用密網(wǎng)格。

圖3 不同網(wǎng)格密度在不同展向站位物面壓力系數(shù)曲線對(duì)比Fig.3 Comparison of surface pressure distribution between coarse and fine grid

3.2 亞聲速旋渦流動(dòng)

3.2.1 渦破裂前的狀態(tài)——狀態(tài)1

狀態(tài)1的自由來流馬赫數(shù)為Ma∞＝0.4，迎角α＝18.5°，為中度迎角，該狀態(tài)下繞三角翼上翼面旋渦流場(chǎng)未發(fā)生渦破裂。圖3給出了該狀態(tài)的數(shù)值模擬結(jié)果，其中圖4（a）左側(cè)為采用標(biāo)準(zhǔn)k－ω湍流模型計(jì)算的物面壓力系數(shù)分布，右側(cè)為Pω增強(qiáng)型k－ω湍流模型的計(jì)算結(jié)果，并給出了沿三角翼5個(gè)流向站位（x／cr＝0.2，0.4，0.6，0.8，0.95）的物面壓力系數(shù)分布與NASA NTF風(fēng)洞試驗(yàn)數(shù)據(jù)的比較。通過對(duì)比可以看出，標(biāo)準(zhǔn)k－ω模型和Pω增強(qiáng)型k－ω模型均預(yù)測(cè)出了三角翼上翼面的主分離渦結(jié)構(gòu)，但k－ω模型僅預(yù)測(cè)出了主渦的位置，不同站位處的主渦吸力峰峰值均低于試驗(yàn)數(shù)據(jù)，且主渦寬度較寬，沒有預(yù)測(cè)出二次渦。Pω增強(qiáng)型k－ω湍流模型預(yù)測(cè)的主渦峰值和主渦寬度均與試驗(yàn)值吻合良好，且預(yù)測(cè)出了二次渦吸力峰，但吸力峰值低于試驗(yàn)數(shù)據(jù)，二次渦吸力峰的展向站位也比試驗(yàn)數(shù)據(jù)靠近機(jī)翼內(nèi)側(cè)，這是由于在計(jì)算中沒有考慮橫流轉(zhuǎn)捩導(dǎo)致的。

圖4（b）為不同站位切片上速度分量U等值線分布圖，可以看到清晰的主渦結(jié)構(gòu)，主渦沒有發(fā)生渦破裂。二次渦結(jié)構(gòu)靠近機(jī)翼表面，不夠清晰，而圖4（a）中x／cr＝0.95站位二次渦位置附近壓力峰值消失，分布較為平坦，表明在三角翼后緣附近發(fā)生了渦破裂。

3.2.2 渦破裂后的狀態(tài)——狀態(tài)2

該亞聲速狀態(tài)發(fā)生了渦破裂現(xiàn)象，迎角為α＝23°，和狀態(tài)1一樣給出了標(biāo)準(zhǔn)k－ω模型和Pω增強(qiáng)型k－ω模型的計(jì)算結(jié)果，如圖5（a）所示。通過對(duì)比同樣可以明顯地看出，標(biāo)準(zhǔn)k－ω模型預(yù)測(cè)的主渦吸力峰值偏低，寬度偏大，主渦強(qiáng)度較弱；Pω增強(qiáng)型k－ω模型預(yù)測(cè)的主渦吸力峰值與試驗(yàn)數(shù)據(jù)吻合良好，僅x／cr＝0.80站位處偏高，二次渦峰值高于試驗(yàn)數(shù)據(jù)。從圖4（b）后兩個(gè)站位切片的速度分量U等值線圖及沿渦軸流線可以明顯地看出，在三角翼后緣附近，主渦和二次渦均發(fā)生了渦破裂。

圖5 狀態(tài)2計(jì)算結(jié)果Fig.5 Computational results for case 2

3.3 跨聲速旋渦流動(dòng)

3.3.1 渦破裂前的狀態(tài)——狀態(tài)3

圖6是跨聲速狀態(tài)下渦破裂前的計(jì)算狀態(tài)，α＝18.5°。該狀態(tài)下，標(biāo)準(zhǔn)k－ω模型沒有模擬出明顯的主渦結(jié)構(gòu)，主渦吸力峰也不明顯，x／cr＝0.20和0.4站位壓力系數(shù)峰值遠(yuǎn)低于試驗(yàn)結(jié)果；而Pω增強(qiáng)型kω模型預(yù)測(cè)出了清晰的主渦和二次渦結(jié)構(gòu)，x／cr＝0.20和0.4兩個(gè)站位處主渦峰值稍低于試驗(yàn)數(shù)據(jù)，且沿展向方向稍靠近機(jī)翼內(nèi)側(cè)，x／cr＝0.60處預(yù)測(cè)峰值稍高，沿展向方向稍靠近機(jī)翼外側(cè)，最后兩個(gè)站位與試驗(yàn)數(shù)據(jù)吻合良好。此外，各個(gè)站位處的二次渦峰值均高于試驗(yàn)值。

圖6 狀態(tài)3計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)對(duì)比Fig.6 The comparison between the computational results and the experimental data for case 3

3.3.2 渦破裂后的狀態(tài)——狀態(tài)4、5

由文獻(xiàn)［11－12，14］可知，跨聲速下，隨著迎角的增加，由于正激波的干擾作用，三角翼上主渦會(huì)發(fā)生突然的渦破裂現(xiàn)象。當(dāng)Ma∞＝0.85時(shí)，VFE－2尖前緣三角翼發(fā)生渦破裂的臨界迎角為24.6°，在該迎角下，支架前端正激波的干擾導(dǎo)致主渦在x／cr＝0.60附近發(fā)生了渦破裂。

圖7為采用標(biāo)準(zhǔn)k－ω模型預(yù)測(cè)α＝23°時(shí)的物面壓力分布，從x／cr＝0.60站位處主渦峰值減弱，最后兩個(gè)站位處主渦峰值消失，變得較為平坦可以明顯地看出，當(dāng)α＝23°時(shí)已經(jīng)發(fā)生了渦破裂，表明標(biāo)準(zhǔn)k－ω模型預(yù)測(cè)的臨界迎角偏小。此外，與試驗(yàn)數(shù)據(jù)對(duì)比可知，標(biāo)準(zhǔn)k－ω模型并未預(yù)測(cè)出明顯的二次渦結(jié)構(gòu)，壓力系數(shù)曲線與α＝24.6°的試驗(yàn)數(shù)據(jù)更為接近。

圖7 標(biāo)準(zhǔn)k－ω模型計(jì)算的狀態(tài)4結(jié)果Fig.7 The comparison between computational results of standard Wilcoxk－ωmodel and experimental data for case 4

圖8為采用Pω增強(qiáng)型k－ω模型計(jì)算的狀態(tài)4、5狀態(tài)結(jié)果。圖8左側(cè)為α＝23°迎角下的計(jì)算結(jié)果，該狀態(tài)未發(fā)生渦破裂現(xiàn)象，通過物面壓力系數(shù)云圖可以看出清晰的主渦和二次渦結(jié)構(gòu)，與試驗(yàn)數(shù)據(jù)對(duì)比可知，各個(gè)站位計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)吻合良好；圖8右側(cè)為α＝24.6°迎角下的計(jì)算結(jié)果，由物面壓力系數(shù)云圖可以看出，該狀態(tài)在x／cr＝0.60位置附近發(fā)生了渦破裂現(xiàn)象。各個(gè)站位壓力系數(shù)曲線與試驗(yàn)值吻合良好，前兩個(gè)站位處主渦和二次渦吸力峰均比較清晰，x／cr＝0.60站位處主渦峰值降低，這是由于激波干擾導(dǎo)致主渦強(qiáng)度減弱引起的，后兩個(gè)站位處由于發(fā)生了渦破裂現(xiàn)象，主分離渦結(jié)構(gòu)消失，壓力分布變得較為平坦。Pω增強(qiáng)型k－ω模型準(zhǔn)確地模擬出了這一復(fù)雜現(xiàn)象，表明Pω增強(qiáng)型k－ω模型適用于繞三角翼跨聲速流場(chǎng)的數(shù)值模擬。

圖8Pω增強(qiáng)型k－ω模型計(jì)算的狀態(tài)4、5結(jié)果與試驗(yàn)數(shù)據(jù)對(duì)比Fig.8 The comparison between the computational results ofk－ωwithPωenhancer model and the experimental data for case 4and case5

3.3.3 軸向流動(dòng)預(yù)測(cè)的驗(yàn)證——狀態(tài)6

為了進(jìn)一步考核Pω增強(qiáng)型k－ω模型對(duì)渦破裂前旋渦流場(chǎng)軸向流動(dòng)的預(yù)測(cè)能力，與konrath等［21］在DLR的PIV試驗(yàn)結(jié)果進(jìn)行了對(duì)比。選擇對(duì)比的實(shí)驗(yàn)狀態(tài)為Ma＝0.8，α＝25.9°，Re＝3×106。CFD計(jì)算選用迎角為α＝26°，其他狀態(tài)參數(shù)與試驗(yàn)相同。

圖9（a）～（f）為x／cr＝0.5、0.55、0.6、0.7、0.75和0.8站位處采用PIV測(cè)量出的速度分量u′（u′＝U／V∞）切片云圖，試驗(yàn)中，該狀態(tài)的渦破裂位置在x／cr＝0.6～0.7之間，圖10（a）～（f）給出了采用Pω增強(qiáng)型k－ω模型計(jì)算出的與試驗(yàn)相同站位的無量綱速度分量u′切片云圖。通過對(duì)比可以看出，在x／cr＝0.5站位，數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果一樣，主渦均呈扁平形狀，試驗(yàn)的主渦形狀更為扁平，渦核的展向和法向位置基本一致。在x／cr＝0.55處，主渦特性與前一站位基本一致，CFD和試驗(yàn)預(yù)測(cè)的二次渦均出現(xiàn)了渦破裂。x／cr＝0.6處，兩者出現(xiàn)了較大差異，從圖9（c）可以看出，主渦結(jié)構(gòu)仍然存在，但與前兩個(gè)站位相比，由于接近渦破裂位置，主渦強(qiáng)度已經(jīng)明顯減弱，變得不穩(wěn)定。CFD方法預(yù)測(cè)的渦破裂位置在x／cr＝0.59處（從圖10可以看出），因此，CFD方法預(yù)測(cè)結(jié)果中，x／cr＝0. 6已處于渦破裂之后。后三個(gè)站位均處于渦破裂之后，PIV和CFD的速度云圖基本一致。

圖9 PIV顯示的速度分量u′切片云圖［21］Fig.9 Velocity componentu′contours for PIV

圖10 狀態(tài)6計(jì)算所得的速度分量u′切片云圖Fig.10 Velocity componentu′contours for computational results of case 6

圖11給出了本文CFD方法和Glassgow大學(xué)的CFD方法［21］（使用PMB（Parallel Multi－Block）流場(chǎng)求解器，基于多塊結(jié)構(gòu)化網(wǎng)格求解非定常RANS方程，控制方程的空間離散采用格心格式的有限體積法。其中，對(duì)流通量項(xiàng)的離散采用Osher和Roe迎風(fēng)型數(shù)值通量，采用MUSCL變量插值方法獲得二階空間離散精度，Van Albada限制器消除激波附近的非物理振蕩；粘性項(xiàng)離散采用二階中心差分格式；定常計(jì)算采用隱式時(shí)間推進(jìn)格式；湍流模型的處理與本文類似，采用對(duì)耗散率方程進(jìn)行旋轉(zhuǎn)修正的k－ω模型）計(jì)算的渦軸處無量綱速度分量u′分布與DLR試驗(yàn)數(shù)據(jù)的對(duì)比，試驗(yàn)給出了x／cr＝0.5、0.55和0.6站位處渦軸上三個(gè)點(diǎn)的速度u′，分別為：1.962（x／cr＝0.5處）、1.870（x／cr＝0.55處）、1.522（x／cr＝0.60處）。Glassgow大學(xué)和本文CFD方法預(yù)測(cè)的主渦軸線上u′最大值分別為1.864（x／cr＝0.378處）和1.879（x／cr＝0.55處），可以看出，本文的預(yù)測(cè)值更接近試驗(yàn)結(jié)果。CFD方法預(yù)測(cè)的渦核軸向速度最大值均低于試驗(yàn)數(shù)據(jù)，這可能是CFD方法預(yù)測(cè)的渦破裂位置靠前的原因。

圖11 渦核軸線上速度分量u′與試驗(yàn)PIV數(shù)據(jù)對(duì)比Fig.11 Comparison of velocity componentu′through vortex core between computational results and experimental PIV data for case 6

4 結(jié) 論

采用Pω增強(qiáng)型k－ω湍流模型，通過對(duì)尖前緣三角翼亞聲速和跨聲速繞流的數(shù)值模擬，可以得出以下結(jié)論：

1）Pω增強(qiáng)型k－ω湍流模型能夠增加旋渦流動(dòng)高旋度區(qū)域ω方程生成項(xiàng)，減小旋渦區(qū)域的湍流渦黏系數(shù)，而其他區(qū)域保持標(biāo)準(zhǔn)Wilcoxk－ω湍流模型基本特性，適用于繞三角翼的亞聲速和跨聲速旋渦流動(dòng)數(shù)值模擬。

2）Pω增強(qiáng)型k－ω湍流模型能夠準(zhǔn)確預(yù)測(cè)出主渦吸力峰值和主渦的展向位置，但預(yù)測(cè)的二次渦峰值偏高，這與采用全湍流計(jì)算，未考慮轉(zhuǎn)捩影響有關(guān)。

3）Pω增強(qiáng)型k－ω湍流模型能夠準(zhǔn)確地預(yù)測(cè)出跨聲速下激波干擾導(dǎo)致主分離渦破裂的臨界迎角及渦破裂位置。

［1］Nelson R C，Pelletier A.The unsteady aerodynamics of slender wings and aircraft undergoing large amplitude maneuvers［J］.Progress in Aerospace Sciences，2003，39（2－3）：185－248.

［2］Brandsma F J，Kok J C，Dol H S，et al.Leading edge vortex flow computations and comparison with DNW－HST wind tunnel data［C］／／Proceeds of the RTO／AVT Symposium on Vortex Flows and High Angle of Attack.NATO RTO／AVT，2001.

［3］Earnshaw P B.An experimental investigation of the structure of a leading edge vortex［R］.RAE Technical Note No.Aero.2740，Royal Aeronatical Establisment，March 1961.

［4］Chu J，Luckring J M.Experimental surface pressure data obtained on a 65°delta wing across Reynolds number and Mach number ranges：Volume 1－sharp leading edge［R］.NASA Technical Memorandum 4645，NASA Langley Research center，F(xiàn)ebruary 1996.

［5］Klute S M，Vlachos P P，Telionis D P.High speed digital particle image velocimtery study of vortex breakdown［J］.AIAA Journal，March 2005，43（3）：642－650.

［6］Jobe C E.Vortex breakdown location over 65°delta wing empiricism and experiment［R］.AIAA 1998－2526.

［7］Konrath R，Schrder A，Kompenhans J.Analysis of PIV results obtained for the VFE－2 65°delta wing configuration at sub－and transonic speeds［R］.AIAA Paper 2006－3003.

［8］Schrder A，Agocs J，F(xiàn)rahnert H，et al.Application of stereo PIV to the VFE－2 65°delta wing configuration at sub－and transonic speeds［R］.AIAA 2006－3486.

［9］Morton S A，F(xiàn)orsythe J，Mitchell A M，et al.DES and RANS simulations of delta wing vertical flows［R］.AIAA 2002－0587.

［10］Mitchell A M，Morton S A，F(xiàn)orsythe J R，et al.Analysis of delta wing vortical substructures using detached eddy simulation［J］.AIAA Journal，2006，44（5）：964－972.

［11］Schiavetta L A，Boelens O J，F(xiàn)ritz W.Analysis of Transonic Flow on a slender delta wing using CFD［R］.AIAA 2006－3171.

［12］Schiavetta L A，Boelens O J，Crippa S，et al.Shock effects on delta wing vortex breakdown［R］.AIAA 2008－395.

［13］Wang GuangXue，Deng Xiaogang，Liu Huayong，et al.Application of high－order scheme（WCNS）at high angles of incidence for delta wing［J］.Acta Aerodynamica Sinica，2012，30（1）：28－33.（in Chinese）王光學(xué)，鄧小剛，劉化勇，等.高階精度格式WCNS在三角翼大迎角模擬中的應(yīng)用研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2012，30（1）：28－33.

［14］Li Xile，Yang Yong，Zhang Qiang，et al.Numerical simulation of shock／vortex interaction in transonic flow around a delta wing［J］.Acta Aeronautica et Astronautica Sinica，2013，34（4）：750－761.（in Chinese）李喜樂，楊永，張強(qiáng)，等.繞跨聲速三角翼的激波／渦干擾流場(chǎng)數(shù)值模擬［J］.航空學(xué)報(bào)，2013，34（4）：750－761.

［15］Anderson J D.Fundamentals of aerodynamics［M］.2nd edition，McGraw－Hill International Editions，1991，356－364.

［16］Deléry J M.Aspects of vortex breakdown［J］.Progress in Aerospace Sciences，1994，30（1）：1－59.

［17］Hall M G.Vortex breakdown［J］.Annual Review of Fluid Mechanics，1972，4：195－218.

［18］Spall R E，Gatski T B，Grosch C E.A criterion for vortex breakdown［J］.Physics of Fluids，1987，30（11）：3434－3440.

［19］Robinson B A，Barnett R M，Agrawal S.Simple numerical criterion for vortex breakdown［J］.AIAA Journal，1994，32（1）：116－122.

［20］Wilcox D C.Reassessment of the scale determining equation for advanced turbulence models［J］.AIAA Journal，26（11）：1299－1310.

［21］Lucy Schiavetta.Evaluation of URANS and DES predictions of vortical flows over slender delta wings［D］.Department of Aerospace Engineering University of Glassgow，F(xiàn)ebruary 2007.

Application of k－ωwith Pωenhancer turbulence model to delta wing vortical flows

Zhang Dongyun1，Li Xile2，＊，Yang Yong3，Zhang Qiang3
（1.Shanghai aircraft design ＆Research Institute，Commercial Aircraft Corporation of China Ltd，Shanghai 201210，China；2.China Academy of Aerospace Aerodynamics，Beijing 100074，China；3.National Key Laboratory of Science and Technology on Aerodynamic Design and Research，Northwestern Polytechnical University，Xi’an 710072，China）

The production ofkωwithin standard Wilcoxk－ωturbulence model is dependent on the mean strain－rate of the flow and does not take the rotation rate into account，resulting in the prediction of a weak vortex which can’t be sustained and diffuses quickly.A suitable sensor is defined to distinguish between shear layers and vortex cores，to increase the production of the dissipation rate（ω）within regions of highly rotational flow.k－ωwithPωenhancer turbulence model is applied to RANS solver for structured grid.Subsonic and transonic flows around a 65oswept leading edge delta wing are simulated numerically.The computational results are compared with detailed experimental data obtained from a series of wind tunnel tests in the NTF at NASA and DNW－TWG at DLR.The comparisons show that thek－ωwithPωenhancer turbulence model produces better results when it comes to capturing the vortex core，the critical angle of attack and the position of vortex breakdown at subsonic and transonic flows，thus is suitable for simulating the vortical flow around delta wing.

k－ωturbulence model；production；delta wing；transonic；critical angle of attack；vortex breakdown

V211.3

Adoi：10.7638／kqdlxxb－2015.0200

0258－1825（2016）04－0461－08

2015－11－25；

2016－03－01

張冬云（1985－），男，安徽合肥人，工程師，主要從事民用飛機(jī)氣動(dòng)設(shè)計(jì)與CFD研究.E－mail：zhangdongyun＠comac.cc

李喜樂（1983－），男，河南洛陽人，博士，工程師，主要從事飛行器設(shè)計(jì)與CFD研究.E－mail：lxl1027＠163.com

張冬云，李喜樂，楊永，等.Pω增強(qiáng)型k－ω湍流模型在三角翼旋渦流動(dòng)的應(yīng)用［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2016，34（4）：461－467.

10.7638／kqdlxxb－2015.0200 Zhang D Y，Li X L，Yang Y，et al.Application ofk－ωwithPωenhancer turbulence model to delta wing vortical flows［J］.Acta Aerodynamica Sinica，2016，34（4）：461－467.