余磊,肖明雅,王潤濤,陳江鋒,周濤濤,唐志全*

1. 安徽全柴動力股份有限公司,安徽滁州 239500;2.中國科學技術大學工程科學學院,安徽合肥 230026;3.合肥工業(yè)大學汽車與交通工程學院,安徽合肥 230009

0 引言

隨著能源與環(huán)境問題日益突出,機動車污染物排放標準不斷提高,2019年開始實行的國六排放標準對燃油車排放提出更加嚴格的要求,高效與清潔技術已經成為內燃機發(fā)展的重點方向[1-2]。燃油霧化質量是影響內燃機燃燒過程的重要因素,在很大程度上決定內燃機的動力和排放特性[3]。噴霧霧化過程是極其復雜的兩相流問題,涉及到氣-液兩相間的相互耦合以及液滴的碰撞、破碎和蒸發(fā)等多個過程。常用的試驗手段只能宏觀展現燃油霧化的工作過程[4],獲取的流場信息有限,數值模擬具有經濟性和強大的信息獲取能力,可以獲得霧化過程中湍流與液滴的詳細演化特征,已經成為研究噴霧發(fā)展過程的重要手段[5-6]。

模擬湍流和液相運動是噴霧霧化過程數值仿真的重點。由于直接數值模擬和大渦模擬需要巨大的計算量,目前工程應用上仍多采用Reynolds平均法(簡稱RANS方法)處理霧化過程中的氣相湍流運動,同時耦合液滴破碎模型模擬霧化過程。萬吉安[3]基于歐拉-拉格朗日方法,湍流模型采用標準k-ε模型,液滴破碎模型采用Wave模型進行柴油機燃油噴霧過程的三維數值模擬;周乃君等[7]采用標準k-ε湍流模型耦合Wave液滴破碎模型研究了各種工況下高壓燃油噴射霧化過程中的噴霧貫穿距和錐角等宏觀特性的發(fā)展規(guī)律;邢志海等[8]采用標準k-ε湍流模型耦合Wave液滴破碎模型對高壓共軌單孔噴油器高壓燃油噴霧微觀特性數值模擬研究;王勇[9]采用標準k-ε湍流模型、KH-RT液滴破碎模型進行了超高壓燃油霧化過程仿真研究,探討了不同模型參數不同工況下數值模擬效果;王昆朋[10]采用RNGk-ε湍流模型耦合KH-RT液滴破碎模型探究了不同噴射壓力、不同噴射策略對發(fā)動機噴霧發(fā)展和混合氣形成過程的影響。雖然對噴霧霧化過程進行模擬時采用了不同的湍流模型和液滴破碎模型,但對內燃機工作過程中的高壓噴霧計算,哪種模型的精度更高并無明確結論,因此有必要開展不同湍流模型和液滴破碎模型對高壓噴霧過程模擬結果準確性的評估研究。

近年來,美國Sandia國家實驗室下的Engine Combustion Network(ECN)研究組開展了一系列燃油噴霧試驗研究,為噴霧霧化及燃燒過程的數值模擬提供了更多的標定數據,其中應用較為廣泛的是以正十二烷為燃料、代號為Spray A的試驗研究,這是由于正十二烷與柴油的碳鏈長度、沸點等性質更接近,使用正十二烷可以更準確地再現柴油噴出后的蒸發(fā)和混合過程,因此受到廣大研究人員的關注[11-16]。

本文基于歐拉-拉格朗日方法,采用不同的湍流模型和液滴破碎模型對Spray A噴霧開展數值模擬研究,對比仿真與試驗結果,分析不同模型對內燃機高壓噴霧霧化過程模擬準確性的影響,為高壓噴霧霧化過程的數值模擬研究提供參考。

1 數值模型

1.1 湍流模型

RANS方法的核心是不直接求解瞬時的Navier-Stokes (N-S) 方程,而是求解時均化的雷諾(Reynolds)方程。對流動質量、動量和能量輸運方程進行Reynolds平均,構建時均輸運方程。本文中不涉及氣體高速流動,將流體視為不可壓縮流體,故以張量表示的時均連續(xù)方程[17]和時均Reynolds方程(RANS方程)[17]為:

(1)

(2)

式(1)和(2)組成的方程組中未知量大于方程組數量,方程組不封閉,需要引入新的方程才能使方程組封閉。引入不同的方程則形成不同的模型。

1.1.1 Reynolds應力方程模型(RSM)

(3)

式中:Cij為對流項,DT,ij為湍流擴散項,DL,ij為分子黏性擴散項,Pij為剪應力產生項,Gij為浮升力產生項,φij為壓力應變項,εij為黏性耗散項,Fij為系統(tǒng)旋轉產生項。

雷諾應力模型是對湍流流動最完整的物理表示,可以直接模擬RANS方程中的流動項,適用于強旋轉和流線曲率造成的湍流穩(wěn)定性等非常復雜的流動問題。

1.1.2k-ε模型

在標準k-ε模型中,關于湍動能k的輸運方程[17]為:

(4)

式中:μt為湍動黏度,μt=ρCμk2/ε,其中Cμ為經驗常數;Prk為k對應的普朗特數;Gk為由平均速度梯度引起的k的產生項;Gb為由浮升力引起的k的產生項;YM為可壓縮湍流中的脈動擴張項;Sk為用戶定義的關于k的源項。

關于湍流耗散率ε的輸運方程[17]為:

(5)

式中:C1ε、C2ε、C3ε是由研究經驗得出的常數;Prε為ε對應的普朗特數;Sε為用戶定義的關于ε的源項。

為了彌補標準k-ε模型的不足,科研人員提出了適用不同情況的RNGk-ε模型[18]、Realizablek-ε模型[19]。RNGk-ε模型中對湍流黏度進行了修改,解決了標準k-ε模型對于旋轉流動、強曲率流動計算失真的問題。Realizablek-ε模型將湍流黏度計算中的Cμ與應變率相聯系。該模型適用于高壓射流、剪切流、大曲率流動,并且強化了分離與強壓差流動界層的性能。

1.1.3k-ω模型

標準k-ω模型中,對應的湍動能k的輸運方程[17]為:

(6)

式中:Γk為k的擴散頂,Yk為k在湍流下的耗散。

湍流耗散率ω的輸運方程為:

(7)

式中:Gω為ω產生項;Γω為ω的擴散率,Yω為ω在湍流下的耗散,Sω為源項。

相比k-ε模型,k-ω模型可以更好地模擬壁面附近的流動,而且能夠很好地模擬逆壓梯度邊界層流動和分離。但是該模型不能準確模擬自然流,對于逆壓梯度造成的剪切力預測過高。相比于標準k-ω模型,通常建議選用改進的SSTk-ω模型[20]。改進的SSTk-ω模型以壁面距離為基準,在壁面近場使用k-ω模型,在壁面遠場采用k-ε模型,更好地結合k-ε和k-ω模型的優(yōu)點。

1.2 離散相模型

噴霧過程包含液滴的噴射、破碎等復雜過程,均需要相應的模型進行?；痆21]。本文中以離散相模型 (discrete phase model,DPM) 對液滴的位置、質量、動量以及溫度進行追蹤求解。

弛豫時間[17]

(8)

式中:dp為液滴顆粒直徑,m;Re為相對雷諾數;CD為曳力系數。

笛卡爾坐標系下x方向的離散相液滴的作用力平衡方程[17]為:

(9)

式中:u為流體速度,m/s;uP為液滴的速度,m/s;gx為液滴在x方向的加速度,m/s2;ρp為液滴密度,kg/m3;Fx為液滴受到的其他力的作用,包括熱致遷移力、布朗力、薩夫曼升力等,本文中只考慮虛擬質量力即使顆粒周圍流體加速引起的附加作用力、壓力梯度力的作用。

虛擬質量力[17]

(10)

壓力梯度力[17]

(11)

考慮瞬時湍流速度脈動對顆粒軌跡的影響,對離散相顆粒的分布采用隨機軌道模型。在實際過程中,通過連續(xù)相的差得到液滴所需要的流動信息。

1.3 液滴破碎模型

由于高壓噴霧霧化過程中韋伯數較高,故選用Wave Breakup Model(簡稱Wave模型)、KH-RT Breakup Model(簡稱KH-RT模型)、Stochastic Secondary Droplet Model(簡稱SSD模型) 3種適用于高韋伯數流場的液滴破碎模型進行分析。

1.3.1 Wave模型

基于圓柱射流穩(wěn)定性的線性分析,研究人員建立了Wave破碎模型[22]。圓柱射流在環(huán)境氣體中受到的小擾動

η=η0exp(inz+wt),

(12)

式中:η0為初始振幅,n為波數,z為流動方向,w為波的增長率。

基于小擾動假設,對控制射流柱運動的N-S方程進行線性化處理,代入邊界條件可得到色散方程。擬合色散方程的解,可得表面波的最大增長速率的方程[22]為:

(13)

表面波的波長方程為:

(14)

式中:ΛKH為增長最快表面波的波長,m。

Wave破碎模型認為射流柱表面氣液間剪切力的作用使射流柱出現了不穩(wěn)定性,即KH表面波引起了射流表面的不穩(wěn)定性,并最終導致了子液滴從射流柱表面剝落。因此,Wave破碎模型也稱為KH破碎模型。

液滴破碎后的子液滴半徑[22]

(15)

式中:B0為模型常數,本文中B0=0.61。

關于母液滴的半徑表達式[22]為:

(16)

式中:tKH為KH波破碎時間,tKH=3.726B1r/ΛKHΩKH,其中,B1為與噴嘴結構、噴嘴內部流動狀態(tài)有關的模型常數,常取1.7～60.0,本文中B1=9.0。

1.3.2 KH-RT模型

KH-RT模型認為液滴的破碎不僅是因為KH表面波的擾動,還有另一種不穩(wěn)定波的擾動,即RT表面波。RT模型認為擾動是在高速射流的氣液交界面上由密度差導致液相向氣相加速而引起的。與KH模型類似,RT模型也是根據增長率最高的波的波長來決定液滴的破碎方式與破碎時間。雖然KH-RT模型中的KH波與RT波共同控制液滴破碎,但兩者之間存在一定的競爭關系。在噴霧近場,由于射流內部接觸的氣體很少,所以KH波對液滴的破碎起主要作用;在噴霧遠場,液滴與環(huán)境氣體直接接觸,首先在RT模型中,由液滴的破碎特征時間是否大于RT破碎時間決定液滴是否發(fā)生破碎,其次在KH模型中,由液滴韋伯數是否大于臨界韋伯數判斷是否發(fā)生KH破碎,臨界韋伯數設為12。

KH波最大增長率ΩKH、相應波長ΛKH及破碎時間表達式在KH破碎模型中已描述,KH-RT模型中B0=0.61,B1=60。

RT波最大增長率ΩRT、相應波長ΛRT、破碎時間tRT以及子液滴半徑rchild的表達式[22]分別為:

(17)

ΛRT=2πCRT/KRT,

(18)

tRT=Cτ/ΩRT,

(19)

rchild=πCRT/KRT,

(20)

式中:gt為液滴行進方向上的加速度,m/s2;CRT為破碎半徑常數,CRT=1;Cτ為破碎時間常數,Cτ=0.5。

1.3.3 SSD模型

SSD模型將液滴破裂視為一種離散的隨機事件導致直徑尺度在一定范圍內的分布。在SSD模型中,液滴破裂的概率與母液滴和次級液滴的大小無關。破裂模型可預測發(fā)生破裂的時間以及新液滴的數量和屬性。液滴半徑大于臨界半徑rcr時液滴發(fā)生破裂。

臨界半徑[17]

(21)

式中:Wecr為臨界韋伯數,本文中Wecr=6。

破碎時間[17]

(22)

式中:B為破碎常數,B=4。半徑大于臨界半徑的液滴的破裂時間會增大,當液滴上的破裂時間大于臨界破裂時間時,發(fā)生破裂。

當一個液滴破碎時,這個液滴破碎成數個新的包裹。包裹中顆粒的半徑通過對數分布函數隨機獲取[17]。

液滴發(fā)生分解時,創(chuàng)建足夠多的小包裹,使每個小包裹所代表的液滴數大致等于所設置的小包裹目標數Np,本文中Np=1 000。

1.4 物理模型及計算條件

噴霧過程示意如圖1所示。由圖1可知:噴霧的物理過程大致可分為噴嘴孔內流動、噴霧近場的初次霧化、噴霧遠場的二次霧化3個階段。這3個階段中包含以下幾方面的復雜過程:液柱流動、分裂形成液滴;液滴進一步破碎、碰撞、聚合、再破碎;液相蒸發(fā)、與環(huán)境氣體混合以及液相與氣相的相互作用等,使得噴霧過程的模擬計算十分復雜。

圖1 噴霧過程示意圖 圖2 平口噴嘴結構

Spray A噴霧試驗中使用單孔共軌噴油器,數值模擬選用平口霧化噴嘴模型,平口霧化噴嘴模型結構如圖2所示,噴嘴結構參數孔板長度L=0.1 mm,噴射器內直徑dj=0.09 mm,拐角曲率半徑re=0.001 mm。

為簡化計算,將噴霧過程近似看作一個空間對稱的物理過程,可進一步簡化為二維旋轉對稱的物理過程。參考Omidvar等[23]的數值計算方法,建立二維對稱計算域,計算網格如圖3所示,其中,對網格以噴口為中心采用橫向和縱向的漸進加密,網格的長和寬分別為80、40 mm,網格數為800×400個,O(0,0)處為噴嘴位置;邊界1為恒溫壁面邊界,邊界2為壓力出口邊界。本文中數值模擬基于Spray A試驗,基本參數采用Spray A的計算工況的參數,如表1所示。

表1 Spray A計算工況相關參數

a)局部放大 b)計算網格圖3 Spray A計算網格

1.5 網格無關性分析

為了保證獲得的結果不受網格質量影響,對網格進行無關性驗證,選用標準k-ε湍流模型耦合Wave破碎模型進行驗證,網格數量分別為20.6萬、32.0萬,40.5萬。以y=0中軸線上的1.5 ms時氣相流場的速度計算結果作為對比,不同網格數量下的y=0中軸線x方向不同位置處氣相流場速度曲線如圖4所示。由圖4可知:除氣相流場峰值速度外,3個網格的速度基本相同;網格數量為32.0萬時氣相流場的峰值速度與40.5萬時的峰值速度幾乎一致,但是網格數量為20.6萬時氣相流場的峰值速度明顯較小,說明網格數量為32.0萬時的網格精度可以滿足計算要求,即滿足網格無關性。因此,本文中的計算采用網格數量為32.0萬的網格。

圖4 不同網格數量時中軸線x方向不同位置處氣相流場速度對比圖

2 仿真結果及分析

2.1 湍流模型對仿真結果的影響分析

通過標準k-ε、RNGk-ε、Realizablek-ε、SSTk-ω以及RSM 5種湍流模型耦合Wave破碎模型對Spray A噴霧霧化過程進行模擬計算,分析不同湍流模型對數值結果準確性的影響,將5種湍流模型的仿真結果和試驗結果進行對比,不同湍流模型對液相貫穿距的仿真結果如圖5所示。

圖5 試驗和不同湍流模型對液相貫穿距的預測結果

由圖5可知:不同湍流模型預測的液相貫穿距發(fā)展趨勢基本與試驗結果一致,都是先快速增長,而后液相貫穿距在10 mm附近上下波動;但RNGk-ε和RSM模型仿真的液相貫穿距比試驗偏大的比例較大,這是由于模型對噴霧場中的湍動能預測偏小,流場的主流速度較大,使得液相軸向運動距離更遠;5種湍流模型中,標準k-ε模型對液相貫穿距模擬結果的穩(wěn)定性最好,隨時間的波動最小。

不同湍流模型對氣相貫穿距的預測結果和試驗結果對比如圖6所示。由圖6可知:5種模型計算的噴霧氣相貫穿距均偏小,但Realizablek-ε模型的仿真結果最接近試驗結果,其次是標準k-ε模型,RNGk-ε模型的誤差最大。

圖6 不同湍流模型對氣相貫穿距的預測結果和試驗結果對比

不同時刻噴霧形態(tài)試驗結果及不同湍流模型仿真結果如圖7～12所示,其中白色點劃線表示噴霧軸向到達的最遠位置。由圖7～12可知:仿真得到的噴霧形態(tài)與試驗圖像大致相同。仿真圖像與試驗圖像的差異主要表現在2方面:1)試驗圖像的邊界輪廓粗糙不平,仿真圖像邊界輪廓光滑平整,這是由于實際情況下液滴與連續(xù)相的湍流強相互作用形成噴霧的非對稱結構,從而導致噴霧蒸氣邊緣粗糙不平,而5種湍流模型得到的是噴霧流場的湍流平均信息,形成了噴霧結構沿中心軸線成近似軸對稱形狀;2)模擬計算的噴霧蒸氣發(fā)展一直慢于實際情況。Realizablek-ε模型計算的氣相縱向發(fā)展最接近試驗數據,與圖6中的噴霧氣相貫穿距一致;但從噴霧氣相徑向發(fā)展來看,SSTk-ω模型仿真結果更接近試驗,RNGk-ε模型得到的蒸氣形態(tài)存在一定程度的失真,主要表現在噴霧近場徑向發(fā)展過慢,噴霧遠場徑向發(fā)展過快。

a)經過64 μs b)經過106 μs c)經過298 μs d)經過510 μs e)經過808 μs圖7 氣相噴霧形態(tài)演變過程試驗結果

不同湍流模型對氣相噴霧錐角的預測結果如表2所示。由表2可知:5種湍流模型的仿真結果與試驗結果都存在一定的誤差,但是SSTk-ω模型的氣相噴霧錐角與試驗結果的誤差較小,其次是標準k-ε模型, RNGk-ε模型的誤差最大,與氣相噴霧形態(tài)演變過程的徑向發(fā)展結果相同。

綜上,Realizablek-ε模型在噴霧貫穿距、氣相發(fā)展的仿真結果與試驗結果最接近, RNGk-ε模型在所有結果的預測中誤差都相對較大。在Lu等[21]的研究中也發(fā)現,RNGk-ε模型不能準確求解大體積運動,對噴霧體內氣流速度的預測不合理,結合本文的結果說明此模型不適合Spray A條件下的高壓噴霧仿真計算。雖然RSM模型能夠較準確地預測噴霧體徑向發(fā)展,但軸向發(fā)展的誤差相對較大,而且RSM模型需要求解的方程較多,計算時間較長,對網格質量要求較高,所以也不建議采用RSM模型進行噴霧模擬。雖然SSTk-ω模型在噴霧體徑向發(fā)展的預測略優(yōu)于其他模型,但在對液相貫穿距的預測不如標準k-ε模型,對蒸氣的發(fā)展預測不如Realizablek-ε模型,也不推薦使用。標準k-ε模型和Realizablek-ε模型的預測準確性都比較高,如果想獲得更準確的液相貫穿距,優(yōu)先選用標準k-ε模型;如果想獲得更準確的氣相貫穿距,優(yōu)先選用Realizablek-ε模型。

2.2 液滴破碎模型對仿真結果的影響分析

通過標準k-ε模型結合不同破碎模型進行仿真,試驗和不同破碎模型對噴霧貫穿距的預測結果如圖13所示。

a)液相貫穿距 b)氣相貫穿距圖13 試驗和不同破碎模型對噴霧貫穿距的預測結果

由圖13可知:Wave模型與KH-RT模型計算得到的液相貫穿距與試驗結果相對接近,SSD模型的預測結果相對偏大;不同破碎模型的氣相貫穿距預測變化規(guī)律都與試驗結果一致,但KH-RT模型和SSD模型預測的氣相貫穿距相比Wave模型更接近試驗結果。由于SSD模型中的液滴破碎是離散隨機過程,液滴破碎的概率與液滴大小無關,導致部分液滴粒徑過大,液滴獲得的動量過大,軸向運動距離偏大,即液相貫穿距偏大。

不同破碎模型對氣相噴霧形態(tài)演變過程的預測結果如圖14～17所示,圖中白色點劃線表示噴霧軸向發(fā)展的最遠位置。不同破碎模型對氣相噴霧錐角的預測結果如表3所示。由圖14～17可知:3種破碎模型計算得到的噴霧輪廓形態(tài)與試驗結果基本相似,噴霧的徑向發(fā)展情況與試驗結果吻合良好,軸向發(fā)展與試驗結果存在一定差異。結合圖13中氣相貫穿距以及表3蒸氣錐角的仿真結果,KH-RT模型計算結果最接近試驗結果,其次是Wave模型,SSD模型計算得到的蒸氣形態(tài)在噴霧近場不連續(xù),這是由于SSD模型模擬的液滴粒徑較大,液滴蒸發(fā)速度慢,因此噴霧蒸氣在根部出現間斷。

表3 不同破碎模型對氣相噴霧錐角的預測結果

a)經過64 μs b)經過106 μs c)經過298 μs d)經過510 μs e)經過808 μs圖14 氣相噴霧形態(tài)演變過程試驗結果

由于SSD模型對噴霧蒸氣形態(tài)的預測存在根部間斷現象,所以SSD模型不適用于高壓噴霧過程的模擬研究;KH-RT模型在液相貫穿距的預測上與Wave模型的預測結果相差不多,但在氣相貫穿距和氣相噴霧錐角的預測都優(yōu)于Wave模型。所以,在對高壓噴霧霧化過程中的數值計算中破碎模型應優(yōu)先選擇KH-RT模型。

3 結論

基于歐拉-拉格朗日方法,在不同湍流模型以及液滴破碎模型下對Spray A條件下的噴霧霧化過程進行仿真研究,通過對比得到高壓噴霧霧化過程數值模擬研究中相對合適的湍流模型以及液滴破碎模型。

1)5種湍流模型中,RNGk-ε、SSTk-ω以及RSM模型預測的液相貫穿距和氣相貫穿距與試驗結果都相差比較大,這3種湍流模型不適合應用于高壓噴霧霧化過程模擬研究;標準k-ε模型與Realizablek-ε模型的模擬結果與試驗比較吻合,在噴霧液相貫穿距的預測中,標準k-ε模型的仿真結果更好,在氣相貫穿距的預測中,Realizablek-ε模型的仿真結果更好,這兩種模型都比較適合應用于高壓噴霧霧化過程數值模擬研究。

2)SSD模型預測的液相貫穿對比試驗結果明顯偏大,而且對噴霧蒸氣形態(tài)的預測存在根部間斷現象,不適用于高壓噴霧霧化過程的模擬研究;KH-RT模型與Wave模型的模擬結果與試驗結果都比較接近,但是KH-RT模型的模擬結果準確性更高。