艦船柴油機微通道內燃油空化現象的模型對比

李子銘, 劉振明, 劉景斌, 陳萍

(海軍工程大學 動力工程學院, 湖北 武漢 430033)

0 引言

作為現代艦船柴油機的主要發展方向,高速高功率密度柴油機具有高噴射壓力、高轉速和噴油持續期短等特點[1-2]。隨著高壓共軌系統中的噴射壓力不斷提高,噴油器噴孔內部的空化問題愈發嚴重。當噴孔進出口壓差較大時,噴孔內流道近壁面處會產生低壓區,致使該區域內的高速燃油壓力迅速降至其飽和蒸氣壓力值,導致空化現象的出現。噴孔內的空化不僅影響后續的噴霧燃燒效果,而且長期的空化腐蝕會導致噴孔內壁脫落,進而影響噴油精度和柴油機運行可靠性。目前,使用計算流體力學(CFD)方法驗證相關的噴孔燃油空化流動試驗,進而研究噴油器噴孔內的空化問題,是艦船柴油機領域中的研究熱點之一。

當前廣泛用于驗證數值模型計算可靠性的燃油空化試驗主要有Sou等[3]提出的二維噴孔試驗和Winklhofer等[4-5]提出的微通道試驗。由于微通道試驗更接近實際的噴孔工作情況,近年來國內外學者們將其廣泛應用于各類新型的噴孔數值模型驗證中。Dai等[6]將Winklhofer試驗用于弧形結構噴嘴模型可靠性的驗證,并利用該模型研究了噴嘴弧度對噴孔內瞬態空化特性的影響,發現噴嘴弧度可顯著降低孔內空化程度,弧度越大,空化程度越小,噴孔出口平均速度越小,但質量流量變化不大。通過結合Winklhofer試驗數據,Sa等[7]在計算中使用了Realizablek-ε(k為湍動能,ε為耗散率)湍流與Schnerr-Sauer(SS)空化組合模型(簡稱Realizablek-ε+SS模型),驗證了所構建多相流模型的空化分布與試驗一致,進而定量研究了噴油器針閥上螺旋反槽結構對孔內燃油湍流流動及后續燃油噴射霧化過程的影響。Cristofaro等[8]使用大渦模擬(LES)模型計算發現所構建的多相流模型的空化分布與試驗一致,并研究了液體燃料黏度對微通道節流流動時質量流量、速度曲線和空化分布的影響。Zhao等[9]使用網格尺寸為30 μm的模型驗證了模型的空化分布、出口質量流量和截面速度分布與試驗相近,最終研究了燃油的可壓縮性對燃油在噴孔中排放系數和臨界空化數的影響。Guo等[10]使用齊次松弛模型(HRM)和重整化群(RNG)k-ε湍流模型計算驗證了所用模型與試驗中的出口質量流量和空化分布誤差極小,研究了針閥運動對噴嘴內空化發展的影響。

目前在研究Winklhofer試驗時雖然能夠獲取到與試驗現象較接近的數值計算結果,但因為所使用的計算模型和網格劃分策略不盡相同,計算模型選取對計算結果和試驗驗證的影響尚無報道,影響機理尚未探明,導致模型驗證策略尚無規律可循。

本文使用CFD計算軟件,基于剪切應力輸運(SST)k-ω(ω為比耗散率)與Zwart-Gerber-Belamri(ZGB)組合模型(簡稱SSTk-ω+ZGB模型)、SSTk-ω與SS組合模型(簡稱SSTk-ω+SS組合模型)、Realizablek-ε與ZGB組合模型(簡稱Realizablek-ε+ZGB模型)、Realizablek-ε+SS模型、重整化群(RNG)k-ε與ZGB組合模型(簡稱RNGk-ε+ZGB模型)、RNGk-ε與SS組合模型(簡稱RNGk-ε+SS模型)6種湍流和空化模型組合,依次對Winklhofer微通道試驗中多種壓差工況下的壓力梯度、空化分布、出口質量流量和截面流速分布進行數值模擬研究,對比討論各模型組合之間的數值結果差異性,總結與評價模型選取對試驗驗證結果的影響。

1 計算模型理論

1.1 湍流模型

1.1.1 SSTk-ω模型

SSTk-ω模型由基線(BSL)k-ω模型改進而來,繼承了k-ω類模型在近壁面區域中的魯棒性、準確性,以及k-ε類模型在遠場區域中的自由流獨立性,并在湍流黏度中定義了湍流剪切應力的傳輸過程,其湍動能k和比耗散率ω的輸運方程分別為

(1)

(2)

式中:ρ為流體密度;ui為流體速度在i方向上的分量,i為ijk空間坐標系的i方向;xi為流線分量x在i方向上的分量;Γk、Γω分別為湍動能k和比耗散率ω的有效擴散率;Gk為湍動能k的生成項;Yk、Yω分別為湍流引起的湍動能k和比耗散率ω的耗散項;Dω為交叉擴散項;Gb為浮力引起的湍動能;Gωb為ω輸運方程中的浮力項。上述變量的具體表達式參見文獻[11]。

1.1.2 Realizablek-ε模型

Realizablek-ε模型基于均方渦度波動傳輸方程導出了耗散率ε的修正傳輸方程,并且在多數流動中尤其是分離流和具有復雜特征的二次流中具有優異的計算性能,其湍動能k和耗散率ε的輸運方程分別為

(3)

(4)

式中:μ為黏度;μt為湍流黏度;uj為流體速度在j方向上的分量,j為ijk空間坐標系的j方向;xj為流線分量x在j方向上的分量;YM為可壓縮湍流流動中脈動膨脹對總耗散率的貢獻程度;C1為方程中的一個系數項,

(5)

(6)

(7)

1.1.3 RNGk-ε模型

RNGk-ε模型在標準k-ε模型基礎上使用重整化群的統計方法改進而來,該模型的特點在于提高了渦流計算的準確性,能有效計算近壁面區域的低雷諾數效應,因此RNGk-ε模型相較于標準k-ε模型的精度更高,適用的流動情況更廣泛,其湍動能k和耗散率ε的輸運方程分別為

(8)

(9)

式中:μeff為有效黏度;αk、αε分別為湍動能k和耗散率ε的有效普朗特數的倒數;C1ε、C2ε、C3ε為耗散率ε輸運方程中的常數,C1ε=1.42、C2ε=1.68;Rε為耗散率ε輸運方程的附加項。上述變量的具體表達式參見文獻[13]。

1.2 空化模型

1.2.1 空化模型基本控制方程

空化模型可與控制氣液混合物的多相流模型、傳統的湍流模型一同構成多相流空化建模方法,其控制方程表述為

(10)

式中:fv為蒸氣質量分數;vv為蒸氣速度;Γ為擴散系數;Re為蒸氣生成速率;Rc為蒸氣凝結速率。

1.2.2 ZGB模型

ZGB模型假設多相流體系中的氣泡都具有相同的大小,忽略不凝結氣體對空化流的影響,并考慮使用單氣泡質量變化率和氣泡數密度來計算單位體積內的相間傳質速率[14]:

當p≤pv時,

(11)

當p≥pv時,

(12)

式中:p為流場局部壓力;pv為飽和蒸氣壓;Fv為蒸發系數,Fv=50;Fc為凝結系數,Fc=0.01;αn為氣核處體積分數,αn=5×10-4;αv為蒸氣體積分數;ρv、ρl分別為氣相密度和液相密度;RB為氣泡半徑,RB=1×10-6m。

1.2.3 SS模型

SS模型使用氣泡半徑和氣泡數密度來計算單位體積內的相間傳質速率,忽略不凝結氣體對空化流的影響,并假設當多相體系中沒有氣泡產生和湮滅時氣泡數密度應保持恒定[15]:

當p≤pv時:

(13)

當p≥pv時:

(14)

式中:蒸發系數Fv=1;凝結系數Fc=0.2;α為氣體體積分數。

2 數值模型構建及驗證

2.1 物理模型構建

本文選取Winklhofer可視化燃油試驗[5]中的U形管作為研究對象,其幾何結構如圖1[16]所示。 圖1(a)展示了可視化試驗使用的微通道組件,圖1(b)是微通道的放大示意圖。微通道模型幾何參數如下:管道厚度W=300 μm,微通道入口高度Hi=301 μm,微通道出口高度Ho=284 μm,微通道長度L=1 000 μm,微通道入口拐角半徑R=20 μm。

圖1 微通道幾何結構[16]Fig.1 Microchannel geometry[16]

在文獻[8-9,17]中,使用額外預置緩沖流道的微通道模型與Winklhofer試驗結果進行了對比驗證,目的是為了讓仿真時的微通道進出口壓力與試驗中的壓力條件保持一致,以便仿真結果更接近實際試驗測量結果。如圖2所示,本文參考了 Guan等[17]的仿真模型,在微通道模型前后分別預置了2 000 μm×300 μm×2 312 μm和3 000 μm×300 μm×2 312 μm的緩沖流道。

圖2 仿真計算模型Fig.2 Simulation calculation model

2.2 數值計算設置

數值計算時使用的燃油熱物性參數如表1[17]所示,設置燃油為可壓縮流體。此外,微通道進出口分別定義為壓力入口條件和壓力出口條件,保持100 bar的入口壓力不變,為計算模型設置不同的出口壓力條件(19～85 bar)。

表1 燃油熱物性參數[17]Table 1 Fuel thermophysical properties[17]

數值模型的離散基于有限體積法,對于氣-液兩相流使用Mixture模型描述,壓力-速度耦合方法使用壓力隱式分裂算子(PISO)算法,壓力求解格式和梯度求解格式分別設為PRESTO!格式、Least Squares Cell-Based格式,并將時間插值格式設置為2階隱式格式以提高求解精度,其余格式設為QUICK格式。

2.3 網格無關性驗證

對于微通道流域及微通道進出口附近的緩沖區流域,使用ICEM軟件進行局部網格加密,以便后續的模型計算得到精確的流動參數結果,網格劃分如圖3所示。

圖3 仿真模型網格分布Fig.3 Simulation mesh model

對于上述網格劃分策略,依次采用量級為5萬、20萬、40萬、60萬、80萬的網格進行網格無關性驗證。在進出口壓差為80 bar的工況下,Winklhofer流動試驗已能觀測到微通道中出現了超空化現象。因此,網格無關性驗證基于上述壓力工況,5種網格數下的微通道出口質量流量如圖4所示。從圖4中可見,當網格數量超過60萬時,微通道出口處質量流量趨于一穩定值?？紤]到后續數值計算的精度和計算資源分配,最終確定計算模型的網格數量為60萬,其中近壁面處最小網格尺寸為0.5 μm,中央流道處最大網格尺寸約為4 μm。

圖4 網格無關性驗證Fig.4 Grid independence verification

3 計算結果和討論

當前后緩沖區流域具有一定的壓差,燃油流經微通道時,流道截面突然變小,流速突增,壓力突降至燃油飽和蒸汽壓力以下,此時會發生流動分離、空化等現象。上述現象直接影響了微通道流域的壓力與速度分布,并反映在出口質量流量、空化區域分布、流速分布和壓力梯度分布等結果參數上。本文分別研究以上物理參數的分布特性,對比分析多種數值模型應用情景下的結果差異,最終對可壓縮燃油在多模型應用下的數值計算準確性進行討論。

3.1 壓力梯度結果對比

為了闡明通道內空化初生到空化發展的內在機理,首先結合Winklhofer等[4-5]試驗數據,分析壓差為70 bar情況下模型中心線處沿流向的壓力梯度變化情況,如圖5所示。其中,圖5(a)中的x軸刻度對應圖5(b)中的模型中心線標識。在該壓力工況下,微通道內的燃油流動剛好達到臨界空化條件。但由于此時的空化區域并未完全發展,燃油壓力值先是在通道入口至空化區域中部迅速降低,并在局部流域橫截面最小的空化區域中部降至極小值,然后流域橫截面有所增加,壓力值回升。在燃油液體流過空化區域末端并重新附著回壁面后,隨著流體在不斷接近通道出口的過程中,流域橫截面逐漸減小,流速逐漸增加,導致壓力也逐漸降低。

圖5 壓力梯度變化圖Fig.5 Change in pressure

圖5展示了多種湍流和空化模型應用下的壓力梯度曲線變化差異。在燃油自上游緩沖區流域流至微通道入口附近,不同模型組合所得壓力梯度曲線幾乎重合,并與試驗結果相近。隨著燃油先后流經空化區域和無空化區域,除含有RNGk-ε模型的2種組合外,其余4種模型組合都有壓力值先突降至一個極小值,再有所回升,最后逐漸減小的變化趨勢。由表2可知:當模型中心線位置點處于 -1.00～-0.25 mm和-0.25～0.22 mm兩個區段時,6種模型組合所得結果與試驗數據接近,誤差均在7%以內;處于0.22～0.50 mm和0.50～2.00 mm兩個區段時,只有RNGk-ε+ZGB模型、RNGk-ε+SS模型這兩種組合所得壓力梯度平均值與試驗數據相近,誤差在6%以內。對于其他模型組合所得結果,其計算誤差主要在10%～18%范圍內,仿真效果較差。

表2 壓力梯度平均值結果對比Table 2 Comparison of average pressure gradient results bar

3.2 空化分布結果對比

表3為試驗和多種模型組合仿真的空化分布結果對比,其中Δp為前后緩沖區流域壓差。由于前后緩沖區域存在壓差,同時燃油在微通道中流動增加,使得局部區域的燃油壓力降低至燃油蒸發壓力,區域中的燃油瞬間蒸發為氣相燃油,空化產生。

從表3中可以看出,不同模型組合所得的云圖與試驗結果之間存在明顯的差異:由于實際燃油流動中存在著不凝結氣體,導致燃油抗拉強度下降,燃油液體在壓力下降過程中更容易產生氣核,空化更容易發生[18]。同時,ZGB模型和SS模型皆未考慮不凝結氣體對空化域的影響。表3中所示壓差為60 bar工況下,6種模型組合都能模擬得到通道入口處的附著空化,與試驗結果一致;壓差為70 bar和80 bar工況下,6種模型組合所得空化域長度都比試驗結果中的小。與試驗結果對比可知,70 bar的壓差工況下,6種模型組合模擬到的空化域長度遠小于試驗結果;80 bar工況下,SSTk-ω+ZGB模型組合所得空化域長度稍小于試驗值,Realizablek-ε+ZGB模型組合所得結果與70 bar壓差工況下的試驗結果相當,其他4種模型組合所得空化域長度皆稍大于或相近于60 bar壓差工況下的試驗結果。

表3 氣相體積分數云圖

從試驗結果可以看出,壓差從60 bar增加至 80 bar,氣相燃油不再局限于微通道入口處的回流區域,而是逐漸沿著壁面延伸至出口處。然而,試驗中在微通道觀察到的空化域比模擬所獲取到的空化域大,不同空化模型所得到的空化分布情況也不同。Yu等[19]的研究結果表明,對同一微通道結構使用不同的空化模型,觀測到的空化分布不同。 Giannadakis等[20]指出,空化域發展情況不依賴于空化泡初始直徑。Altimira等[16]也發現試驗中空化域比模擬中的大,并指出燃油壓縮性會影響蒸氣冷凝速率,這可能是造成模擬中的空化域偏小的原因。同時他進一步研究發現,使用不同密度的網格分布以及調節空化模型中影響空化泡發展的控制參數大小,對近壁面處的空化分布沒有顯著影響。

綜上所述,隨著進出口壓差的增大,依附于微通道壁面的空化域不斷發展,沿壁面逐漸延伸至通道出口處;在相同的進出口壓差下,造成不同模型組合之間的空化分布差異性的原因是不同的湍流和空化模型在數學表述上存在差異。

3.3 出口質量流量結果對比

圖6揭示了微通道出口質量流量與微通道進出口壓差的變化特性,同時展示了6種模型組合所得模擬結果的差異。本文研究了13個壓差工況下(Δp1=19 bar、Δp2=45 bar、Δp3=58 bar、Δp4=60 bar、Δp5=63 bar、Δp6=65 bar、Δp7=67 bar、Δp8=69 bar、Δp9=70 bar、Δp10=71 bar、Δp11=75 bar、Δp12=80 bar、Δp13=85 bar)的出口質量流量,與試驗所采用的工況條件完全一致。

圖6 出口質量流量圖Fig.6 Outlet mass flow rate diagram

由試驗結果可知,在壓差從19 bar逐漸增加至85 bar的過程中,出口質量流量首先隨著進出口壓差的增大而增大,隨后質量流量曲線在70 bar壓差下出現拐點,在這之后的出口質量流量保持一定值,不再隨壓差的增大而變化。如圖6所示,19～70 bar壓差范圍內,6種模型組合都能夠獲得與試驗數據吻合較好的數值結果。表4為壓力梯度平均值結果對比,由表4可知,將模型結果按試驗數據劃分成出口質量流量變化趨勢不同的兩個區段(19～70 bar和70～85 bar),6種模型所得結果與試驗的平均值誤差在10%以內。其中,Realizablek-ε+ZGB模型、RNGk-ε+ZGB模型兩種組合的仿真誤差較小,在兩個區段內的誤差均在4%以內。在 70 bar壓差之后,Realizablek-ε+ZGB模型、RNGk-ε+ZGB模型2種組合計算到的結果雖然稍大于試驗數據,但可以復現出口質量流量不再增加的趨勢,其余4種模型組合模擬效果較差,無法反映出這一趨勢。

3.4 流速結果對比

圖7所示為55 bar和67 bar壓差下微通道內2個橫截面(x1=0.053 mm;x2=0.17 mm)沿y軸的速度分布情況。由于當壓差為67 bar時,氣化燃油完全占據了空化區域,試驗過程中難以測量近壁面處的速度變化[5]。依照Winklhofer試驗中截面流速分布的變化趨勢,可將截面x1、x2在兩種壓差下的仿真結果分區段計算流速平均值,列表匯總,如表5～表8所示。

表4 壓力梯度平均值結果對比Table 4 Comparison of average pressure gradient results g/s

在55 bar壓差工況下,x1處的模擬結果可以反映出試驗數據中的速度變化趨勢,如圖7(b)所示;圖7(d)則表明在x2處所有模型組合均無法在遠離微通道上下壁面約50～100 μm的范圍內再現復雜的流速變化。由表5和表7可知,距離上下壁面約0～50 μm的近壁面流域內,仿真結果與試驗數據仍有一定差距。在y軸位置點處于約50～250 μm范圍內時(即截面x1、x2的中心流域),RNGk-ε+SS模型所得結果與試驗數據最接近,即兩處截面中的計算誤差均在10%以內。SSTk-ω+ZGB模型的計算誤差最大,在兩處截面中的計算誤差都超過了12%。

結合圖7(c)和圖7(e),在67 bar壓差工況下,6種模型組合所得結果均與截面x1處試驗數據的變化趨勢吻合,但這仍難以在截面x2處復現遠離壁面約50～100 μm范圍內的流速分布。由表5和表7可知,相較于55 bar壓差工況,67 bar壓差工況下的近壁面流速分布的計算誤差略有改善,但仿真結果與試驗數據尚有一定差距。對于截面x1、x2的中心流域,RNGk-ε+SS模型的計算誤差依舊最小,即兩個截面處的計算誤差均在6%以內。而SSTk-ω+ZGB的計算誤差最大,在兩處截面中的計算誤差都超過了9%。

4 結論

本文分別使用SSTk-ω、Realizablek-ε和RNGk-ε湍流模型以及ZGB和SS空化模型,對不同進出口壓差工況下的Winklhofer微通道試驗模型中的空化現象進行了數值計算研究。得到主要結論如下:

1)在壓差為70 bar時,雖然RNGk-ε+ZGB模型、RNGk-ε+SS模型相較于其余4種模型組合,反映出微通道中心線處沿流向的壓力先突降至一個極小值,再有所回升,最后逐漸減小的變化趨勢的效果較差,但RNGk-ε+ZGB模型、RNGk-ε+SS模型這兩種組合所得壓力梯度平均值與試驗數據較相近?？傮w上看,上述兩種模型組合的計算誤差都在7%以內。在相同的壓差工況下,模擬所得空化域長度小于試驗觀測所得的空化域長度。不同模型組合之間的空化分布存在差異性,原因是不同的湍流和空化模型在數學表述上存在著差異。

2)出口質量流量的數值計算方面,6種模型組合均可在19～70 bar的壓差范圍內計算得到與試驗數據吻合較好的數值結果,計算誤差控制在10%以內。在壓差超過70 bar后,只有Realizablek-ε+ZGB模型、RNGk-ε+ZGB模型2種組合的計算結果仍然可以反映試驗數據的變化趨勢,且計算誤差在4%以下,而其余4種模型組合均無法體現出實際測量下的出口質量流量變化情況。

表5 55 bar壓差下截面x1流速平均值結果對比Table 5 Comparison of average flow rate results of sectionx1 at 55 bar differential pressure m/s

表6 67 bar壓差下截面x1流速平均值結果對比Table 6 Comparison of average flow rate results of sectionx1 at 67 bar differential pressure m/s

表7 55 bar壓差下截面x2流速平均值結果對比Table 7 Comparison of average flow rate results of sectionx2 at 55 bar differential pressure m/s

3)沿y軸的流速分布數值計算方面,在距離上下壁面約0～50 μm的近壁面流域內,6種模型組合仿真所得結果與試驗數據仍有一定差距。在截面x1、x2的中心流域內,RNGk-ε+SS模型的計算結果變化趨勢與試驗數據最吻合,并且55 bar和 67 bar 兩種壓差工況下的計算誤差分別在10%以下和6%以下。SSTk-ω+ZGB模型組合的計算誤差最大,兩種壓差工況下的計算誤差分別超過12%和9%。>

表8 67 bar壓差下截面x2流速平均值結果對比Table 8 Comparison of average flow rate results of sectionx2 at 67 bar differential pressure m/s

4)與Winklhofer試驗數據對比,RNGk-ε+ZGB模型、RNGk-ε+SS模型適用于計算微通道中心線處的壓力梯度變化;6種模型組合未能較好地反映處試驗時的空化分布情況;Realizablek-ε+ZGB模型、RNGk-ε+ZGB模型可用于計算噴油器微通道的出口質量流量情況;RNGk-ε+SS模型則可以用于計算微通道在發生空化時的某一截面的流速分布情況。