基于k-ω湍流模型的翼尖渦演化過(guò)程數(shù)值模擬

錢(qián) 宇,蔣 皓

(中國(guó)民用航空飛行學(xué)院飛行技術(shù)學(xué)院,廣漢 618307)

當(dāng)機(jī)翼產(chǎn)生正升力時(shí),由于氣體流過(guò)上下翼面產(chǎn)生不同的壓強(qiáng),從而導(dǎo)致氣流的縱向偏移,并在翼尖脫落時(shí)形成翼尖渦[1],翼尖渦向后發(fā)展繼而形成尾流跡。由于大氣環(huán)境的特性,翼尖渦所產(chǎn)生的尾流會(huì)持續(xù)很長(zhǎng)一段時(shí)間,其擾動(dòng)的不穩(wěn)定空氣具有較大的動(dòng)能,為后機(jī)飛行安全造成很大的事故隱患[2]。在中國(guó)經(jīng)濟(jì)高速發(fā)展的現(xiàn)在,有限的空域資源逐漸不能滿足民航運(yùn)輸?shù)男枨?而尾流將直接影響到終端區(qū)的容量,影響飛機(jī)的起降架次,造成航班的延誤。在旋翼的運(yùn)動(dòng)過(guò)程中,翼尖渦的存在將造成機(jī)翼的顫振,從而產(chǎn)生嚴(yán)重的噪聲污染[3]。為解決上述問(wèn)題,對(duì)于翼尖渦演化過(guò)程的研究顯得尤為重要。

目前,對(duì)于翼尖渦的研究方法主要有風(fēng)洞實(shí)驗(yàn)、激光雷達(dá)以及計(jì)算流體力學(xué)等一系列手段[4]。其中計(jì)算流體力學(xué)較其他方法成本低、周期短、計(jì)算精度高,且能較好地捕捉到翼尖渦的變化過(guò)程。對(duì)于翼尖渦演化過(guò)程的仿真,不同的模型方法運(yùn)用于不同的計(jì)算環(huán)境,由此中外學(xué)者進(jìn)行了大量的研究。楊祥生等[5]采用雙方程SSTk-ω模型對(duì)風(fēng)力機(jī)尾流流場(chǎng)進(jìn)行了數(shù)值分析,以期證明該模型在翼尖渦預(yù)測(cè)方面具有較好的應(yīng)用效果。Anusonti-Inthra等[6]利用雷諾平均和基于粒子的渦量輸運(yùn)方法,用來(lái)預(yù)測(cè)低速流動(dòng)條件下孤立機(jī)翼的性能和尾流參數(shù)。Uhl等[7]采用DLR TAU流體求解器與雷諾應(yīng)力模型結(jié)合,研究了不同翼尖形狀對(duì)翼尖渦的演化過(guò)程的影響,并在此基礎(chǔ)上預(yù)測(cè)了渦核軸向速度的變化規(guī)律。中國(guó)對(duì)于翼尖渦流跡的研究現(xiàn)大多數(shù)還處于試驗(yàn)或驗(yàn)證階段[8],在仿真計(jì)算方面也多為研究單一翼型。

研究基于SSTk-ω湍流模型,對(duì)整機(jī)模型進(jìn)行模擬仿真與數(shù)值計(jì)算,充分考慮機(jī)身、水平安定面、垂直安定面對(duì)于機(jī)翼近場(chǎng)翼尖渦的影響,得到機(jī)翼上兩個(gè)共轉(zhuǎn)融合渦的形成過(guò)程及其在近翼流場(chǎng)轉(zhuǎn)化為尾流跡的連續(xù)演化過(guò)程,從氣動(dòng)力性能方面為近場(chǎng)及遠(yuǎn)場(chǎng)翼尖渦流跡的研究提供相應(yīng)的理論基礎(chǔ),也能為民航制定相應(yīng)的規(guī)章制度提供一定的技術(shù)參考。

1 模型與方法

1.1 計(jì)算模型

數(shù)值計(jì)算前需對(duì)計(jì)算模型進(jìn)行三維建模,為保證結(jié)果的準(zhǔn)確性,在采用SolidWorks對(duì)計(jì)算模型進(jìn)行整體建模,所得計(jì)算模型如圖1(a)所示,建模參數(shù)及運(yùn)行實(shí)況如表1所示?？紤]到計(jì)算流體動(dòng)力學(xué)( computational fluid dynamics，CFD)方法對(duì)計(jì)算資源的要求較大,且主要研究近翼流場(chǎng)的翼尖渦,為提高計(jì)算效率,縮短計(jì)算時(shí)間,因此簡(jiǎn)化計(jì)算模型,即采用半模,保留機(jī)身、機(jī)翼、水平安定面及垂直安定面,如圖1(b)所示。

圖1 A320飛機(jī)等比例縮放與計(jì)算模型Fig.1 Proportional scaling model and computational model of A320 airplane

表1 模型尺寸大小及其運(yùn)行實(shí)況Table 1 Model size and operating conditions

1.2 網(wǎng)格的劃分

網(wǎng)格的劃分將計(jì)算模型離散為若干網(wǎng)格,通過(guò)對(duì)每一個(gè)網(wǎng)格進(jìn)行計(jì)算及分析,整合后可得計(jì)算模型的整體氣動(dòng)力性能。利用ICEM-CFD對(duì)計(jì)算模型進(jìn)行結(jié)構(gòu)化網(wǎng)格劃分。結(jié)構(gòu)化網(wǎng)格較非結(jié)構(gòu)化網(wǎng)格而言更為精細(xì),對(duì)于區(qū)域邊界的擬合能更好地實(shí)現(xiàn),且在網(wǎng)格生成的過(guò)程中速度快、質(zhì)量好、表面光滑、更加貼近實(shí)際模型,能很好地適用于流體和表面應(yīng)力的計(jì)算。結(jié)構(gòu)化網(wǎng)格的網(wǎng)格質(zhì)量較好,較好的網(wǎng)格質(zhì)量能充分地提高數(shù)值計(jì)算過(guò)程中各項(xiàng)參數(shù)的收斂性,使得計(jì)算所得的結(jié)果更為準(zhǔn)確、合理。劃分后的計(jì)算網(wǎng)格如圖2所示。

圖2 計(jì)算流場(chǎng)及計(jì)算模型局部放大圖Fig.2 Calculation flow field and local enlarged drawing of calculation model

在劃分網(wǎng)格前,為保證模擬真實(shí)著陸過(guò)程,屏蔽其他因素的干擾。因此設(shè)置流場(chǎng)入口距離模型4c(c為翼展長(zhǎng)度),流場(chǎng)出口離模型20c,上壁面離模型5c,下表面離模型1c。為了使計(jì)算結(jié)果更為精確,則機(jī)翼處需生成更加精細(xì)的網(wǎng)格。考慮到單純?cè)黾泳W(wǎng)格節(jié)點(diǎn)數(shù)或者減小網(wǎng)格增長(zhǎng)率會(huì)大大增加網(wǎng)格數(shù)量,為數(shù)值計(jì)算帶來(lái)相應(yīng)的難度。為了解決這個(gè)問(wèn)題,采用兩次O-Block方法對(duì)機(jī)翼附面層網(wǎng)格進(jìn)行相應(yīng)的加密處理。該方法不僅能對(duì)局部進(jìn)行加密,而且可以較好地避免復(fù)雜形狀的Block頂點(diǎn)處發(fā)生的網(wǎng)格扭曲,從而在機(jī)翼表面得到所需的邊界加密層。O-Block方法劃分的機(jī)翼網(wǎng)格如圖3所示,參數(shù)如表2所示。

圖3 O-Block劃分機(jī)翼網(wǎng)格Fig.3 Generating wing mesh of O-Block

表2 網(wǎng)格劃分參數(shù)Table 2 Meshing parameters

最后一定要保證附面層網(wǎng)格的高度,這將直接影響計(jì)算的收斂性,附面層網(wǎng)格高度y的計(jì)算公式[9]為

(1)

式(1)中:y+取1；μ為流體的動(dòng)力黏度；Uτ為流體的估計(jì)速度；ρ為流體密度。

通過(guò)式(1)的計(jì)算,可以得到附面層網(wǎng)格高度為1.979×10-3mm,因此取y=1.5×10-3mm，滿足計(jì)算要求。最終生成的結(jié)構(gòu)化網(wǎng)格數(shù)總量為7.057×106。

上述網(wǎng)格數(shù)量是在滿足附面層網(wǎng)格要求的基礎(chǔ)上,通過(guò)網(wǎng)格增長(zhǎng)率計(jì)算得到的網(wǎng)格數(shù)量。為驗(yàn)證網(wǎng)格的無(wú)關(guān)性,適量增加節(jié)點(diǎn)數(shù)目、減小網(wǎng)格增長(zhǎng)率,以此增大網(wǎng)格數(shù)量,計(jì)算增大后的758×104網(wǎng)格,并與705×104網(wǎng)格所得計(jì)算結(jié)果做比較。結(jié)果顯示758×104網(wǎng)格與705×104網(wǎng)格對(duì)于近翼流場(chǎng)渦量的相對(duì)誤差小于5%,如圖4所示,說(shuō)明網(wǎng)格數(shù)量對(duì)于結(jié)果的影響已不敏感,故可采用705×104網(wǎng)格進(jìn)行數(shù)值計(jì)算。

圖4 不同網(wǎng)格數(shù)量近翼流場(chǎng)渦量對(duì)比圖Fig.4 Vorticity diagram of near wing flow field with different mesh numbers

1.3 控制方程

控制方程采用常黏度條件下不可壓縮流體的N-S(Navier-Stokes)方程,該矢量形式[10]為

(2)

1.4 湍流模型

對(duì)于湍流的計(jì)算方法,目前主要有雷諾時(shí)均(RANS)、尺度解析以及直接數(shù)值模擬(DNS)。選取雷諾時(shí)均,該方法滿足計(jì)算要求且具有較為高效的計(jì)算速度。基于雷諾時(shí)均方法下常用的湍流模型有如下幾種:Spalart-Allmaras模型、SSTk-ω模型、標(biāo)準(zhǔn)k-ω模型、Reynolds Stress等。

在數(shù)值計(jì)算的過(guò)程中,不同的模型具有不同的適應(yīng)環(huán)境。Spalart-Allmaras模型廣泛用于航空領(lǐng)域?qū)諝鈩?dòng)力學(xué)的計(jì)算,對(duì)于有壁面邊界空氣動(dòng)力學(xué)流動(dòng)有較好的計(jì)算結(jié)果,但是該模型為單一運(yùn)輸方程模型,有較大能量耗散,可以運(yùn)用于粗網(wǎng)格的計(jì)算,所得計(jì)算精度不高。標(biāo)準(zhǔn)k-ε模型計(jì)算量適中,具有較高的計(jì)算精度,但在模擬旋流和擾流時(shí)有缺陷,整個(gè)計(jì)算過(guò)程中伴隨較大的能量耗散,且在計(jì)算壓力梯度較大的復(fù)雜流場(chǎng)時(shí)易失準(zhǔn),導(dǎo)致前端渦旋的過(guò)度生成至整片機(jī)翼上,從而影響計(jì)算結(jié)果的準(zhǔn)確性。Reynolds Stress多用于強(qiáng)渦流,通過(guò)直接使用運(yùn)輸方程來(lái)求解雷諾應(yīng)力,避免了其他模型的黏性假設(shè)。RANS模型最復(fù)雜,收斂性較差,占用過(guò)多的CPU、運(yùn)算時(shí)間及運(yùn)算內(nèi)存,對(duì)于較為復(fù)雜的三維強(qiáng)流動(dòng)較為適用。

對(duì)于SSTk-ε模型,該模型使用混合函數(shù)將標(biāo)準(zhǔn)的k-ε和標(biāo)準(zhǔn)的k-ε相結(jié)合,在近翼面區(qū)域計(jì)算時(shí)保留了標(biāo)準(zhǔn)k-ε湍流模型的特性,即適合于存在逆壓梯度情況的邊界層流動(dòng)；在遠(yuǎn)場(chǎng)區(qū)域則按照標(biāo)準(zhǔn)k-ε湍流模型進(jìn)行計(jì)算。由此該模型在計(jì)算流動(dòng)分離和逆壓方面有較好的表現(xiàn),適用于高逆壓梯度及剪切流的計(jì)算[11]。其數(shù)學(xué)表達(dá)式為

Sω+Dω

(4)

式中：Γ為有效擴(kuò)散率；G、Y、S、D分別為生成項(xiàng)、擴(kuò)散項(xiàng)、用戶自定義項(xiàng)以及交叉擴(kuò)散項(xiàng)；ρ為流體密度；k為紊流動(dòng)能；ω為耗散率，其表達(dá)式為

(9)

σk=σω=2

(10)

式中：ε為紊流耗散項(xiàng)。將N-S方程與SSTk-ω湍流模型相結(jié)合,保證數(shù)值計(jì)算過(guò)程的整體封閉,即可對(duì)生成的結(jié)構(gòu)化網(wǎng)格進(jìn)行數(shù)值計(jì)算。

2 仿真計(jì)算與分析

2.1 方法可行性驗(yàn)證

采用Fluent對(duì)結(jié)構(gòu)化網(wǎng)格進(jìn)行數(shù)值求解,選用基于渦度的轉(zhuǎn)捩修正的SSTk-ω湍流模型；流場(chǎng)采用理想氣體,黏度選擇sutherland；邊界條件采用遠(yuǎn)場(chǎng)壓力入口,設(shè)置標(biāo)準(zhǔn)大氣壓為101.325 kPa,速度為130 n mile/h,迎角為5.5°,邊界出口選用壓力出口。求解方法采用基于耦合的隱式求解,又由于一階迎風(fēng)夸大了擴(kuò)散的影響,容易偏離真正的場(chǎng)分布,因此以一階迎風(fēng)迭代400步作為初始計(jì)算值,再轉(zhuǎn)為二階迎風(fēng)迭代500步,總迭代步長(zhǎng)為900步。最后對(duì)網(wǎng)格標(biāo)準(zhǔn)初始化后開(kāi)始計(jì)算。

在數(shù)值計(jì)算的過(guò)程中,為了證明數(shù)值方法的可行性,運(yùn)用上述方法計(jì)算來(lái)流速度124 n mile/h,雷諾數(shù)4.6×106,迎角10°工況下NACA0012翼型的靜壓力系數(shù)CP,并繪制翼展方向1/3處曲線圖,與Chow等[12]實(shí)驗(yàn)進(jìn)行對(duì)比驗(yàn)證,如圖5所示。

圖5 翼型靜壓力系數(shù)曲線圖Fig.5 Diagram of static pressure coefficient on airfoil

計(jì)算結(jié)果顯示,仿真值與實(shí)驗(yàn)值的平均誤差小于0.05,所得結(jié)果置信度高,誤差在接受范圍之類,因此該數(shù)值計(jì)算方法可行。

2.2 機(jī)翼表面翼尖渦的演化

圖6所示為機(jī)翼側(cè)表面不同截面處渦量分布。通過(guò)該圖可知翼尖渦在機(jī)翼上的演化大體可以分為4個(gè)階段。

圖6 機(jī)翼不同截面處渦量分布Fig.6 Vorticity distribution at different sections of airfoil

第1階段,次渦旋于機(jī)翼前緣形成后,渦量沿前緣到后緣有不斷增加的趨勢(shì),當(dāng)渦旋位于x/b=0.20(b為翼弦長(zhǎng)度)時(shí),次渦旋開(kāi)始從機(jī)翼表面脫落且向上翼面不斷移動(dòng),渦量表現(xiàn)出減小的趨勢(shì)。第2階段,主渦旋在x/b=0.20時(shí)形成,且有一定的增加趨勢(shì),但渦旋的當(dāng)量值總體來(lái)說(shuō)小于次渦旋。第3階段,當(dāng)次渦旋持續(xù)后移,在位于x/b=0.6時(shí)次渦旋開(kāi)始和主渦旋相互糾纏融合,逐漸形成同一個(gè)渦旋,即第1共轉(zhuǎn)融合渦。第4階段為第2共轉(zhuǎn)融合渦的形成,在位于x/b=0.90時(shí)機(jī)翼表面次渦旋和主渦再一次形成,并于x/b=0.95時(shí)再一次糾纏融合,形成第2共轉(zhuǎn)融合渦。

圖7所示為機(jī)翼表面渦核靜壓力系數(shù)變化曲線,從圖7可以看出,次渦旋與次級(jí)尾跡渦的渦核靜壓力系數(shù)有一個(gè)當(dāng)量先減小后增加的過(guò)程,原因是因?yàn)闄C(jī)翼前緣附近有較大的逆壓梯度差,在位于x/b=0.05和x/b=0.10時(shí)存在一對(duì)旋向相反的渦旋。而當(dāng)位于x/b=0.15時(shí),正向的渦旋消失,從而導(dǎo)致負(fù)向渦旋增加,至此次渦旋的渦核靜壓力系數(shù)開(kāi)始不斷上升。

圖7 機(jī)翼表面渦核靜壓力系數(shù)變化曲線Fig.7 Variation curve of vortex core static pressure coefficient on wing surface

第1共轉(zhuǎn)融合渦在形成后其渦核的靜壓力系數(shù)的數(shù)值變化較小,且渦核的靜壓力系數(shù)要小于第2共轉(zhuǎn)融合渦,主要是因?yàn)榈?共轉(zhuǎn)融合渦形成的時(shí)間及距離較第2共轉(zhuǎn)融合渦更長(zhǎng),因此有一定的能量耗散,且在形成第2共轉(zhuǎn)融合渦時(shí),主渦所提供的能量要大于第1次融合時(shí)所提供的能量。

綜上,由主渦和次級(jí)渦兩次糾纏融合形成的兩個(gè)共轉(zhuǎn)融合渦旋與次級(jí)尾跡渦將共同影響近翼流場(chǎng)渦旋的形成。

2.3 近場(chǎng)翼尖渦的演化過(guò)程

通過(guò)上述渦旋的相互作用產(chǎn)生的負(fù)向渦旋流動(dòng),最終在x/c=0.10處形成一個(gè)完整的翼尖渦,如圖8所示。在模擬的過(guò)程中也可以發(fā)現(xiàn),來(lái)流對(duì)近場(chǎng)渦旋移動(dòng)有較大的影響,無(wú)風(fēng)情況下,渦旋應(yīng)有下降的趨勢(shì),而5.5°的迎風(fēng)來(lái)流使得翼尖渦有一定程度的向上偏移,因此翼尖渦在近翼流場(chǎng)的演化也很大程度上受環(huán)境因素的影響。

圖8 近翼流場(chǎng)翼尖渦的渦量分布Fig.8 Vorticity distribution of wingtip vortex at near wing flow field

由圖4、圖9可知,翼尖渦在近場(chǎng)的形成及傳遞會(huì)有一定的能量轉(zhuǎn)換。當(dāng)渦旋從翼面脫離后,在近翼流場(chǎng)形成渦旋的過(guò)程中,近場(chǎng)渦核的靜壓力系數(shù)呈二次曲線減少。在x/c=0.40的位置時(shí),渦核的靜壓力系數(shù)才逐漸趨于平穩(wěn)值,直到在x/L=0.5(L為機(jī)身長(zhǎng)度)的位置時(shí),其數(shù)值變化較小呈現(xiàn)穩(wěn)定狀態(tài)。而翼尖渦也隨著渦核靜壓力系數(shù)的減小,其體積在不斷膨脹增加,繼而逐步生成一個(gè)穩(wěn)定的翼尖渦。

圖9 不同截面處?kù)o壓力系數(shù)分布云圖Fig.9 Static pressure coefficient distribution at different sections

綜上,在機(jī)翼表面形成的融合渦旋受次渦旋、主渦旋共同作用。且融合后渦旋直至離開(kāi)機(jī)翼表面無(wú)明顯的能量損耗。在負(fù)向渦流脫離機(jī)翼后,翼尖渦在近翼流場(chǎng)的形成與運(yùn)動(dòng)過(guò)程受環(huán)境因素影響較大。翼尖渦的生成過(guò)程中,伴隨著明顯的能量轉(zhuǎn)換。

3 結(jié)論

通過(guò)仿真計(jì)算,獲得了機(jī)翼表面及近翼流場(chǎng)翼尖渦的演化過(guò)程,得到了如下結(jié)論。

(1)O-Block方法能得到較好的近壁面網(wǎng)格,有利于近翼流場(chǎng)的計(jì)算；SSTk-ω湍流模型將標(biāo)準(zhǔn)k-ε模型和標(biāo)準(zhǔn)k-ω模型相結(jié)合,能較好地適用于渦旋演化過(guò)程的數(shù)值計(jì)算。

(2)在近翼流場(chǎng)的數(shù)值計(jì)算過(guò)程中,翼尖渦演化至x/L=0.50時(shí)才有較為穩(wěn)定的渦核靜壓力,該處的計(jì)算參數(shù)則能成為遠(yuǎn)渦流場(chǎng)的初始計(jì)算值。在著陸階段,主要的影響因素是大氣的質(zhì)量與空氣的流動(dòng)。對(duì)于遠(yuǎn)場(chǎng)尾跡的計(jì)算則需要考慮地效等因素,這也是今后研究遠(yuǎn)場(chǎng)尾渦演化過(guò)程的重點(diǎn)。