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

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

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

當前航空領域對于噴流流場的預測主要采用RANS方法,但經(jīng)過幾十年發(fā)展,RANS方法仍然不能夠為噴流流場提供可靠的預測,標準湍流模型仍然無法準確預測剪切層增長速率和不受剪切層影響的噴流無黏速勢核心長度。高溫條件下,噴流流動中湍流Prantl數(shù)從近場到遠場變化劇烈,噴流附近區(qū)域Prantl數(shù)低,摻混現(xiàn)象明顯,模擬起來更加困難。

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

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

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

1 數(shù)值模擬方法

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

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

Dk=ρk1.5/lIDDES

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

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

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

2 自由發(fā)展噴管流動計算驗證分析

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

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

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

圖1 噴管網(wǎng)格拓撲結構

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

圖2 時均流場噴管下游中心線上速度計算結果

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

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

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

3 ARN2噴管流動計算驗證分析

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

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

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

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

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

圖7 瞬態(tài)Q等值面噴流流場結構圖

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

圖8 噴流中心線上的時均速度與試驗結果對比

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

圖9 時均流場流向速度分布與試驗結果對比

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

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

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

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

圖10 剪切層厚度分布與試驗結果對比 圖11 噴流中心線上脈動速度與試驗結果對比

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

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

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

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

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

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

5 結 論

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

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

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