龔志斌, 李杰, 單繼祥, 張恒

(1.中國工程物理研究院 總體工程研究所, 四川 綿陽 621900; 2.西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安 710072)

發(fā)動(dòng)機(jī)噴流是動(dòng)力推進(jìn)、氣動(dòng)噪聲、燃燒以及低可探測技術(shù)的關(guān)鍵部分,對其流場結(jié)構(gòu)更好的理解和更為準(zhǔn)確的預(yù)測手段對于發(fā)展高效、靜音和環(huán)保飛行器具有重要的意義。對短距起降(STOL)飛行器而言,動(dòng)力噴流和增升裝置之間的氣動(dòng)干擾尤為強(qiáng)烈,流動(dòng)現(xiàn)象更為復(fù)雜,噪聲問題更加突出,對動(dòng)力噴流流場進(jìn)行準(zhǔn)確模擬,可為STOL構(gòu)型全機(jī)氣動(dòng)設(shè)計(jì)和噪聲預(yù)測奠定良好基礎(chǔ)。

當(dāng)前航空領(lǐng)域?qū)τ趪娏髁鲌龅念A(yù)測主要采用RANS方法,但經(jīng)過幾十年發(fā)展,RANS方法仍然不能夠?yàn)閲娏髁鲌鎏峁┛煽康念A(yù)測,標(biāo)準(zhǔn)湍流模型仍然無法準(zhǔn)確預(yù)測剪切層增長速率和不受剪切層影響的噴流無黏速勢核心長度。高溫條件下,噴流流動(dòng)中湍流Prantl數(shù)從近場到遠(yuǎn)場變化劇烈,噴流附近區(qū)域Prantl數(shù)低,摻混現(xiàn)象明顯,模擬起來更加困難。

國外不少學(xué)者采用LES方法對噴流湍流流場進(jìn)行數(shù)值模擬研究。Pope[1]指出LES方法可以很好地模擬自由剪切流動(dòng)中的動(dòng)量、熱量和質(zhì)量交換過程。噴流發(fā)動(dòng)機(jī)噴口壁面所產(chǎn)生的湍流影響通常被忽略或者假定很小,但Brunet[2]指出它們對于噴流自由剪切層的影響非常顯著,當(dāng)發(fā)動(dòng)機(jī)軸向長度減小或者涵道比增大時(shí)噴口壁面的湍流影響還會(huì)增大。在應(yīng)用 LES方法求解包含發(fā)動(dòng)機(jī)壁面的湍流問題時(shí),為了模擬壁面附近很小尺度的運(yùn)動(dòng),邊界層內(nèi)所需的網(wǎng)格量幾乎與DNS方法相當(dāng),準(zhǔn)確模擬邊界層流動(dòng)所需的計(jì)算資源也很龐大。因此,LES方法主要用于簡單的噴流流動(dòng)和低雷諾數(shù)情形,還不能應(yīng)用于工程實(shí)際中的復(fù)雜發(fā)動(dòng)機(jī)噴流流動(dòng)。此外,LES方法近壁面及可壓縮亞格子應(yīng)力模式的發(fā)展還不夠完善,對大多數(shù)包含壁面的高雷諾數(shù)噴流流動(dòng)模擬能力不足。同時(shí)LES方法噴流預(yù)測結(jié)果的調(diào)查統(tǒng)計(jì)[3]表明,噴流中心線上的速度和湍流強(qiáng)度等信息的預(yù)測非常不理想,其主要原因在于計(jì)算網(wǎng)格質(zhì)量不高且求解方法的格式離散精度不足。

近年來RANS/LES混合方法逐步發(fā)展并應(yīng)用于動(dòng)力噴流流場的數(shù)值模擬,其基本思想是在噴管壁面附近采用高效的RANS方法來對邊界層內(nèi)的小尺度湍流進(jìn)行模擬,下游噴流中的大尺度湍流結(jié)構(gòu)則采用LES方法進(jìn)行模擬,基準(zhǔn)RANS湍流模式和LES區(qū)域亞格子應(yīng)力模型在很大程度上決定了RANS/LES混合方法特性。結(jié)合RANS和LES方法的各自優(yōu)點(diǎn)發(fā)展而來的混合方法是兼顧模擬精度和計(jì)算效率的最優(yōu)選擇,它可以認(rèn)為是LES方法處理高雷諾數(shù)噴流流動(dòng)的可行方法。Shur和Spalart等人[4]針對復(fù)雜噴管噴流流動(dòng)開展了RANS/LES混合方法數(shù)值模擬研究,Cai等[5]基于該方法結(jié)合高階格式討論了網(wǎng)格分布、邊界條件處理等對噴流流動(dòng)數(shù)值模擬的影響。Brunet等[6]基于Zonal-DES方法,針對較為實(shí)際的機(jī)翼/掛架/短艙構(gòu)型,開展了跨聲速情況下的動(dòng)力噴流數(shù)值模擬研究,基于簡單的SA、SST湍流模型和復(fù)雜的DRSM湍流模型的混合方法都取得了與試驗(yàn)較為吻合的結(jié)果。國外還有一些其他研究人員,如Grinstein和Fureby等[7]認(rèn)為不使用亞格子模式的隱式LES方法(ILES方法,又稱NLES方法)能夠較好地捕捉復(fù)雜噴流湍流結(jié)構(gòu),Eastwood等[8-9]采用RANS/NLES方法對簡單軸對稱噴管和發(fā)動(dòng)機(jī)復(fù)雜噴流進(jìn)行模擬并得到良好效果。從公開的文獻(xiàn)資料來看,國內(nèi)采用RANS/LES混合方法對動(dòng)力噴流流場數(shù)值模擬研究成果仍相對較少。

本文擬開展動(dòng)力噴流效應(yīng)混合方法分析及應(yīng)用研究。基于自由發(fā)展噴管流動(dòng)和ARN2噴管流動(dòng)算例,分別討論網(wǎng)格密度和空間離散格式精度對混合方法數(shù)值模擬結(jié)果的影響。在此基礎(chǔ)上,對較為實(shí)際的渦扇發(fā)動(dòng)機(jī)噴流流場進(jìn)行混合方法模擬分析。

1 數(shù)值模擬方法

在有限體積法基礎(chǔ)上,對三維可壓縮非定常N-S方程進(jìn)行求解。對無黏項(xiàng)采用三階MUSCL格式,黏性項(xiàng)采用二階中心差分格式進(jìn)行空間離散;時(shí)間推進(jìn)采用隱式LU-SGS格式。為了提高空間離散格式精度,采用Jiang等[10]提出的有限體積五階WENO-JS格式進(jìn)行左右狀態(tài)變量重構(gòu)來實(shí)現(xiàn)。

采用Shur等[11]的方法,在k-ωSST湍流模型的基礎(chǔ)上,結(jié)合DDES方法和LES壁面模型的特點(diǎn),重新定義亞格子尺度,解決邊界層附近對數(shù)邊界層不匹配問題,加快RANS方法到LES方法的轉(zhuǎn)換,實(shí)現(xiàn)了SST-IDDES方法。其中湍動(dòng)能輸運(yùn)方程中引入長度尺度lIDDES,耗散項(xiàng)為

Dk=ρk1.5/lIDDES

Δ=min[max(CwΔmax,Cwd,Δmin),Δmax)]

Cw=0.15,Δmax和Δmin分別為網(wǎng)格單元三向最大和最小尺度,d為網(wǎng)格單元與壁面距離。

發(fā)動(dòng)機(jī)進(jìn)排氣邊界條件有不同的構(gòu)造方法。其中一種將進(jìn)氣口作為流場出口,給定質(zhì)量流量,或者是靜壓和總溫,排氣口作為流場入口,給定質(zhì)量流量和總溫條件。但是這種方法在整個(gè)網(wǎng)格邊界面上需要指定均勻的質(zhì)量流量,會(huì)影響計(jì)算過程收斂性和計(jì)算結(jié)果精度。本文根據(jù)自由來流條件、進(jìn)氣口質(zhì)量流量和噴口總溫總壓確定發(fā)動(dòng)機(jī)工作狀態(tài),以總壓比(Po/P∞)和總溫比(To/T∞)統(tǒng)一定義進(jìn)排氣動(dòng)力邊界條件,其中Po和To為當(dāng)?shù)乜倝汉涂倻?P∞和T∞為自由來流靜壓和靜溫。在風(fēng)扇入口指定質(zhì)量流量比,在已知自由來流條件下,由能量守恒關(guān)系式和等熵關(guān)系式,可迭代求出風(fēng)扇入口馬赫數(shù),進(jìn)而得到風(fēng)扇入口總壓比和總溫比,內(nèi)外涵道尾噴出口的總壓比和總溫比在已知涵道比的情況下可以類似得到,詳細(xì)推導(dǎo)參見文獻(xiàn)[12-13]。

2 自由發(fā)展噴管流動(dòng)計(jì)算驗(yàn)證分析

Panda等[14]在NASA的Glenn研究中心對自由噴管流動(dòng)進(jìn)行了試驗(yàn)測量,文獻(xiàn)[15]給出了相應(yīng)的噴管直徑、試驗(yàn)環(huán)境和噴管出口流動(dòng)條件。計(jì)算過程中除噴管出口給定排氣邊界外,噴管其他物面邊界給定滑移壁面邊界。試驗(yàn)條件下出口馬赫數(shù)為1.4,對應(yīng)噴管出口邊界總壓比Poj/P∞=3.192、總溫比Toj/T∞=1.012 4,其中Poj和Toj為噴管出口當(dāng)?shù)乜倝汉涂倻?。自由來流馬赫數(shù)為0,為了使得計(jì)算穩(wěn)定,給定來流馬赫數(shù)為0.01。計(jì)算過程中給定雷諾數(shù)1×106。

為準(zhǔn)確計(jì)算噴口下游的動(dòng)力噴流,對噴口下游15D(D為噴管出口直徑)范圍內(nèi)的流場網(wǎng)格進(jìn)行了適當(dāng)加密,生成了粗、細(xì)2套計(jì)算網(wǎng)格,噴流區(qū)域內(nèi)粗網(wǎng)格3個(gè)方向上最大尺寸為0.08D,網(wǎng)格規(guī)模約為700萬;細(xì)網(wǎng)格流向最大尺寸為0.04D,網(wǎng)格規(guī)模約為2 100萬,網(wǎng)格空間截面如圖1所示。

在非定常RANS方法計(jì)算結(jié)果的基礎(chǔ)上采用結(jié)合了三階MUSCL插值格式的SST-IDDES方法進(jìn)行續(xù)算,無量綱時(shí)間步長Δt取為0.000 5。

圖1 噴管網(wǎng)格拓?fù)浣Y(jié)構(gòu)

圖2給出了噴管下游中心線上的時(shí)均速度與試驗(yàn)結(jié)果的對比。X/D(噴管出口下游位置X與噴管直徑D之比)小于8的范圍內(nèi),時(shí)均速度接近于波動(dòng)范圍內(nèi)的平均值。X/D大于8時(shí),時(shí)均速度結(jié)果與試驗(yàn)存在偏差,網(wǎng)格密度增加有利于結(jié)果改善。從圖3不同位置速度分布結(jié)果來看,粗網(wǎng)格情況噴口下游X/D=2和X/D=4處速度峰值與試驗(yàn)接近,但是計(jì)算得到的噴流直徑比試驗(yàn)情況小,當(dāng)?shù)丶羟袑铀俣忍荻群艽?。網(wǎng)格加密后,X/D=2處速度分布沒有變化,X/D=4和X/D=6位置速度分布結(jié)果得到較大改善。

圖2 時(shí)均流場噴管下游中心線上速度計(jì)算結(jié)果

圖4給出了噴流流向截面內(nèi)的瞬態(tài)馬赫數(shù)云圖。從粗網(wǎng)格結(jié)果來看,在噴口下游保持了較長一段相對穩(wěn)定的流動(dòng)狀態(tài),然后高速噴流與周圍靜止氣流在黏性的作用下,產(chǎn)生了剪切層失穩(wěn)。當(dāng)網(wǎng)格加密后,剪切層失穩(wěn)更早出現(xiàn)。結(jié)合圖4結(jié)果來看,噴口下游出現(xiàn)的較長流動(dòng)穩(wěn)定狀態(tài)與物理情況不符,網(wǎng)格加密在一定程度上能夠更加逼近實(shí)際情況。圖6給出了以馬赫數(shù)渲染的流場瞬態(tài)Q等值面圖,網(wǎng)格密度的影響體現(xiàn)得更為明顯,在細(xì)網(wǎng)格的基礎(chǔ)上可獲得噴流流場更多、更豐富的三維小尺度湍流結(jié)構(gòu)。

圖4 瞬態(tài)速度云圖計(jì)算結(jié)果

圖5 流場瞬態(tài)Q等值面圖

3 ARN2噴管流動(dòng)計(jì)算驗(yàn)證分析

NASA的Glenn研究中心設(shè)計(jì)了一系列亞聲速流動(dòng)噴管ARN(acoustic reference nozzles),并在此基礎(chǔ)上進(jìn)行了大量的噴流及噪聲試驗(yàn),為CFD數(shù)值驗(yàn)證提供了有用的數(shù)據(jù),是驗(yàn)證噴流湍流流動(dòng)的標(biāo)準(zhǔn)算例之一[16]。ARN系列收縮噴管為軸對稱外形,其中 ARN2出口直徑為5.08 cm。

針對ARN2噴管外形與流動(dòng)特點(diǎn),建立與上一節(jié)類似的三維軸對稱柱形計(jì)算域,計(jì)算域大小為50D×80D(D為噴管出口直徑)。噴管位置與試驗(yàn)保持一致,噴流入口位于X=0處,上游邊界距離噴管出口30D,噴管外壁面向上游延伸直到上游邊界,下游邊界距離噴管出口的距離為50D,周向邊界到噴管中心軸線的距離為25D。采用點(diǎn)對接多塊結(jié)構(gòu)化網(wǎng)格技術(shù)生成流場結(jié)構(gòu)網(wǎng)格,在噴管出口下游20D范圍內(nèi)對網(wǎng)格適當(dāng)加密,壁面附近網(wǎng)格尺度為1×10-4D,唇口厚度方向布置33個(gè)網(wǎng)格點(diǎn),在出口下游20D范圍的關(guān)注區(qū)域內(nèi)網(wǎng)格最大尺度不超過0.04D,網(wǎng)格規(guī)模為2 200萬。

在非定常RANS方法計(jì)算結(jié)果的基礎(chǔ)上采用結(jié)合了三階MUSCL格式和五階WENO格式的SST-IDDES方法進(jìn)行續(xù)算。在噴管壁面給定無滑移邊界條件,延伸的外壁面給定滑移邊界條件,在周向和上、下游邊界給定黎曼遠(yuǎn)場邊界條件。不考慮溫度效應(yīng),在噴流入口給定總溫比To/T∞=1.0,總壓比Po/P∞=1.8,對應(yīng)噴管出口的馬赫數(shù)為0.9。給定來流馬赫數(shù)為0.01,試驗(yàn)雷諾數(shù)為2.0×106。

圖6 瞬態(tài)馬赫數(shù)云圖

圖7給出了計(jì)算得到的瞬態(tài)馬赫數(shù)云圖。噴流入口到噴管出口之間的區(qū)域內(nèi),流動(dòng)逐步加速到馬赫數(shù)0.9,渦量幾乎為零,流動(dòng)是穩(wěn)定的。采用三階MUSCL格式計(jì)算時(shí),噴管出口下游保持了較長的流動(dòng)穩(wěn)定狀態(tài),然后才由于黏性和速度差導(dǎo)致剪切層失穩(wěn),湍流混合作用增強(qiáng)。采用五階WENO格式計(jì)算時(shí),在唇口下游立即出現(xiàn)較大尺度的流場結(jié)構(gòu),這些流場結(jié)構(gòu)的大小與唇口厚度處于一個(gè)量級,噴流流場中的渦結(jié)構(gòu)更加細(xì)致,渦量耗散要小得多。圖9給出了瞬態(tài)Q等值面結(jié)果,相對于三階MUSCL格式而言,采用五階WENO格式計(jì)算時(shí)整個(gè)噴流流場小尺度渦結(jié)構(gòu)要更加豐富細(xì)致,在唇口附近表現(xiàn)尤為明顯。

圖7 瞬態(tài)Q等值面噴流流場結(jié)構(gòu)圖

圖8給出了時(shí)均流場噴流中心線上的速度分布。在X/D小于5的速勢中心區(qū)域內(nèi),不同離散格式的結(jié)果基本相同。X/D大于5時(shí),相對于三階MUSCL格式而言,采用五階WENO格式計(jì)算得到的速勢中心區(qū)域要更長,中心線上的速度衰減速率要小得多,量值上與試驗(yàn)結(jié)果更為貼近。

圖8 噴流中心線上的時(shí)均速度與試驗(yàn)結(jié)果對比

圖9給出了時(shí)均流場噴口下游X/D為1,4,10的流向速度沿軸向的分布。速度分布結(jié)果與試驗(yàn)較為接近,但不同離散格式的結(jié)果還是存在一定差異。在X/D=1位置,流向速度峰值與試驗(yàn)幾乎一致,三階MUSCL格式結(jié)果在y/D=±0.5位置處速度存在階躍,y/D=±0.5時(shí)速度很快減小并接近為零,這意味著當(dāng)?shù)丶羟袑雍穸缺葘?shí)際情況要小得多,五階WENO格式結(jié)果則不存在階躍,與試驗(yàn)結(jié)果較為接近。在X/D=4位置,2種格式均能夠獲得較好的結(jié)果。在X/D=10位置,速度峰值相對于試驗(yàn)結(jié)果偏低而速勢中心區(qū)以外的速度偏高,尤其是三階MUSCL格式的結(jié)果差異增大。

圖9 時(shí)均流場流向速度分布與試驗(yàn)結(jié)果對比

為了定量的表示噴流與周圍氣流的摻混,定義剪切層厚度參數(shù)b為

b=r|u=0.1Uj-r|u=0.9Uj

式中，r為與噴流中心線的軸向間距。圖10給出了X/D小于8范圍內(nèi)的剪切層厚度分布,速勢中心區(qū)域基本在此范圍。試驗(yàn)結(jié)果顯示剪切層厚度基本上呈線性增長。三階MUSCL格式得到的剪切層厚度在唇口附近接近于零,呈類似拋物型增長,在X/D大于4以后呈線性增長且增長速率要明顯大于試驗(yàn)結(jié)果。五階WENO格式得到的剪切層厚度比試驗(yàn)結(jié)果略高,總體增長趨勢基本一致。

圖11給出了噴流中心線上的脈動(dòng)速度(即湍流強(qiáng)度)分布。中心線上的脈動(dòng)速度在唇口附近最小,往下游逐漸增大,在達(dá)到一定峰值后慢慢下降。三階MUSCL格式和五階WENO格式計(jì)算得到的中心線上脈動(dòng)速度分布與試驗(yàn)結(jié)果趨勢上基本一致。三階MUSCL格式得到的脈動(dòng)速度峰值出現(xiàn)較早,X/D在5到10的范圍內(nèi)要高于試驗(yàn)結(jié)果,當(dāng)X/D大于10時(shí)要明顯低于試驗(yàn)結(jié)果。五階WENO格式計(jì)算得到的中心線上脈動(dòng)速度與試驗(yàn)更為接近。

圖10 剪切層厚度分布與試驗(yàn)結(jié)果對比 圖11 噴流中心線上脈動(dòng)速度與試驗(yàn)結(jié)果對比

4 渦扇發(fā)動(dòng)機(jī)動(dòng)力噴流數(shù)值模擬分析

下面針對較為實(shí)際的渦扇發(fā)動(dòng)機(jī)模型,采用SST-IDDES混合方法,對動(dòng)力噴流效應(yīng)進(jìn)行數(shù)值模擬分析。計(jì)算網(wǎng)格在噴流區(qū)域內(nèi)適當(dāng)加密,在發(fā)動(dòng)機(jī)下游3倍發(fā)動(dòng)機(jī)長度范圍內(nèi)網(wǎng)格尺度不超過發(fā)動(dòng)機(jī)最大直徑的2%,網(wǎng)格總數(shù)約為3 500萬。取海平面標(biāo)準(zhǔn)大氣條件,來流狀態(tài)為:Ma=0.2,Re=4.66×106,發(fā)動(dòng)機(jī)軸線與來流平行。發(fā)動(dòng)機(jī)動(dòng)力邊界條件為:質(zhì)量流量比RMF=1.6,涵道比RB=6.0,風(fēng)扇入口總壓比Po/P∞=1.028 3,總溫比To/T∞=1.008;風(fēng)扇噴流出口總壓比Po/P∞=1.8,總溫比To/T∞=1.2;渦輪噴流出口總壓比Po/P∞=1.6,總溫比To/T∞=2.8。在非定常RANS方法計(jì)算結(jié)果的基礎(chǔ)上采用結(jié)合了三階MUSCL格式和五階WENO格式的SST-IDDES方法進(jìn)行續(xù)算。

渦扇發(fā)動(dòng)機(jī)由于存在內(nèi)外涵道,在出口下游的噴流中會(huì)產(chǎn)生內(nèi)、外剪切層,內(nèi)側(cè)剪切層源于內(nèi)涵噴流與外涵噴流之間的黏性作用,外側(cè)剪切層源于外涵噴流與周圍環(huán)境氣流之間的黏性作用,從圖12瞬態(tài)密度分布和圖13瞬態(tài)渦量分布結(jié)果來看，內(nèi)側(cè)剪切層相對于外側(cè)剪切層更加穩(wěn)定,外側(cè)剪切層更早出現(xiàn)失穩(wěn)現(xiàn)象并誘導(dǎo)內(nèi)側(cè)剪切層失穩(wěn),從而導(dǎo)致內(nèi)外涵道噴流摻混。圖14給出了以馬赫數(shù)渲染的發(fā)動(dòng)機(jī)噴流瞬態(tài)Q等值面結(jié)果,圖中更好地顯示出采用五階WENO格式所解析出的發(fā)動(dòng)機(jī)出口下游豐富清晰的三維小尺度湍流結(jié)構(gòu)和噴流的影響范圍。

圖12 子午面瞬態(tài)密度分布云圖

圖13 子午面瞬態(tài)渦量云圖

圖14 瞬態(tài)Q等值面結(jié)果

5 結(jié) 論

本文從動(dòng)力噴流效應(yīng)精細(xì)化模擬需求出發(fā),開展了混合方法數(shù)值模擬分析與應(yīng)用研究。在k-ωSST湍流模型基礎(chǔ)上建立了IDDES混合方法,在此基礎(chǔ)上,針對自由發(fā)展噴管流動(dòng)和ARN2噴管流動(dòng)開展了模擬驗(yàn)證分析,討論了網(wǎng)格密度和空間離散格式精度對噴流流場計(jì)算結(jié)果的影響,并對較為實(shí)際的渦扇發(fā)動(dòng)機(jī)動(dòng)力噴流流場進(jìn)行了數(shù)值模擬分析。

自由發(fā)展噴管流動(dòng)混合方法計(jì)算結(jié)果表明:增加網(wǎng)格密度可有效改善噴流速度分布預(yù)測結(jié)果,縮短噴口下游非物理流動(dòng)穩(wěn)定狀態(tài),瞬態(tài)流場結(jié)構(gòu)更加清晰合理。ARN2噴管流動(dòng)混合方法計(jì)算結(jié)果表明:相比于三階MUSCL格式,采用五階WENO格式計(jì)算得到的噴流湍流結(jié)構(gòu)更加豐富細(xì)膩,噴口下游流場更加合理,中心線上的速度衰減速率要更小且與試驗(yàn)更為吻合,噴流區(qū)內(nèi)的速度分布與試驗(yàn)更為接近,湍流強(qiáng)度尤其是噴口下游附近的湍流強(qiáng)度預(yù)測更為準(zhǔn)確。

渦扇發(fā)動(dòng)機(jī)噴流流場混合方法數(shù)值模擬結(jié)果表明,外側(cè)剪切層更早出現(xiàn)失穩(wěn)現(xiàn)象并誘導(dǎo)內(nèi)側(cè)剪切層失穩(wěn),從而導(dǎo)致內(nèi)、外涵道噴流摻混。提高空間離散格式精度,噴流摻混更早更充分,湍流結(jié)構(gòu)更加細(xì)致合理。