低液氣密度比橫向射流一次霧化模擬研究

楊 帆，葉桃紅

(中國科學技術大學 熱科學和能源工程系，安徽 合肥 230001)

橫向射流霧化系統是傳統燃氣輪機燃燒裝置的重要組成部分，在航空、航天、艦船、農業等方面具有廣泛應用。一次霧化的機理研究起源于圓管射流霧化研究，Reitz[1]基于小擾動假設，提出了一次霧化和界面不穩定相關，并且最大增長速率的表面波決定了液滴的生成。但是這一理論僅適用于低速射流。在高速射流方面，Wu和Faeth[2]通過研究直射式噴嘴入注靜止空氣的霧化現象，提出了湍流誘導一次霧化機理，即湍流渦旋提供的能量與氣動壓降的能量之和大于生成的液滴所需要的表面能時，液滴將生成。Wu和Faeth 提出當液氣密度比大于500時，空氣動力學效應將被抑制，湍流成為驅動液滴形成的唯一機制。當密度比小于500，空氣動力和湍流都是一次霧化的重要驅動力。Sallam[3]在橫向射流實驗中也發現了類似的現象，管內流動為湍流的橫向射流的霧化同時受到湍流和空氣動力學因素的影響。Arienti[4]在混合歐拉-拉格朗日框架下，采用VOF方法描述連續的液核與氣相，采取拉格朗日粒子追蹤一次霧化生成的液滴，結合湍流是液滴生成的重要驅動力，研究了常溫常壓下液氣密度比為612和773的燃料的一次霧化，研究獲得的液滴直徑和液滴速度與實驗結果基本吻合。在現代燃氣輪機中，為了在短時間內實現高效的混合與燃燒，通常空氣經過進氣口壓縮機壓縮為高壓空氣噴注到燃燒室內[5]，此時液體燃料和橫流氣體之間的密度比要比常溫常壓下的密度比低，此時一次霧化不僅受到湍流的影響還受到空氣動力學影響。顯然Arienti的研究沒有考慮到空氣動力學的影響，不能直接應用于實際燃氣輪機工作工況下的霧化。本文在Arienti工作的基礎上，對低液氣密度比下的橫向射流一次霧化進行了模型修正，并開展了數值仿真計算及分析。

1 數學物理模型

1.1 基本控制方程

設生成的液滴的體積分數為，那么連續相的相體積分數為。設連續相中液相的體積分數為α1，密度為ρ1，連續相使用VOF方法追蹤連續相。

連續相的連續性方程為

(1)

式中：U為混合物速度，m/s；Sα為由一次霧化引起的連續相的質量損失，m/s。

質量交換的源項為

(2)

式中：Vcell為網格體積,m3；Δt時間對應網格因為霧化損失的質量為mp，kg。

連續相的動量方程為

-(1-ε)p+(1-ε)ρg+Mout+Mex+MS

(3)

式中：Reff為有效應力，包括黏性應力與雷諾應力的貢獻，N/s；g為重力加速度,m/s2；Mout為一次霧化引起的連續相動量損失,N/s；Mex為連續相與液滴之間的動量交換，N/s；MS為連續相之間的表面張力,N/s；表面張力采用CSF模型[6]計算。

一次霧化僅引起液相的動量損失：

Mout=ρ1SαU

(4)

一次霧化生成的液滴采用拉格朗日粒子追蹤，通過計算網格內液滴體積確定，液滴受到的曳力采用Schiller- Naumann模型[7]，由于氣體湍流引起的湍流分散力采用Shirolkar JS隨機漫步模型[8]，液滴的碰撞和聚合采用O’Rourke模型[9]，由于實驗中沒有觀測到明顯的二次破碎現象，因此在后續計算中并沒有采用二次破碎模型，此外由于蒸發的時間尺度遠大于液滴在計算域中的停留時間尺度，因此液滴也沒有采用蒸發模型。

1.2 一次霧化模型

在Arienti的研究基礎上，將同時受到液體湍流和空氣動力學因素影響的一次霧化拆分為兩個過程：

(1)第一階段僅僅受到液體湍流的作用，不受空氣動力學的影響，完全由湍流誘導的一次霧化經歷雷諾破碎[10]，根據理論分析其液絲大小滿足：

(5)

式中：τR為雷諾破碎的時間尺度，采用湍流渦旋到達離噴嘴距離為y的位置的時間尺度來估算，s;Uy為當地流速在射流方向的速度分量，m;C為經驗參數，這里采用Sallam[3]實驗中擬合的值0.4。

(2)第二階段完全受空氣動力學的影響，液絲破裂生成更多的小液滴，液滴直徑和液絲直徑之間的關系為[11]

(6)

式中：ρf為液相密度，kg/m3；ρg為氣相密度,kg/m3；μf為液相黏度，Pa·s；li為液相長度尺度，m。

上述關系式僅僅適用于靜止空氣中的破碎，為了在橫向射流中使用，引入氣相和液絲的滑移速度Ur=Ul-Ug。Uinf采用當地的液相平均速度，Ug采用來流氣相速度近似。

一次霧化產生的液滴的質量采用Arienti提出的計算方式：

mp=ηρlμpSΔt

(7)

式中：S為界面的面積，m2。采用isoadvector[12]的Iso-Surface算法估計界面面積，已經證明該方法為二階精度。

不同于M. Arienti采用固定的霧化速率，這里采用Le[13]的實驗表達式，把霧化速率和破碎點高度結合在一起，越靠近破碎點，霧化速率越高。

(8)

式中：yb為橫向射流霧化破碎點高度，m。破碎點高度采用如式(9)估計：

(9)

式中：dj為噴嘴的直徑，m;q為動量通量比；Cb為一個和噴嘴設計、管內流動狀態相關的經驗參數。對于長徑比為10的噴嘴，Cb取為3.3[14]。

1.3 計算工況

數值模擬的幾何體的結構如圖1所示。圖1中入射孔直徑D=457 μm，噴嘴長度L/D=10，位于空氣入口下游5 mm中心線上，空氣入口寬度為46 mm，沿著射流方向的高度為30 mm，通道長度為74.6 mm。數值模擬的設置和Y.Gopala.E[15]的實驗設置一致，研究液氣密度比為312，馬赫數為0.2和0.35的預熱高壓空氣來流下，動量通量比為180時，液體燃料Jet A的霧化。具體的算例的主要參數見表1，Jet A密度為790，黏度為0.001 9，與空氣之間的表面張力系數為0.022。

圖1 計算域示意圖

表1 LJIC算例設置

1.4 數值方法和邊界條件

1)邊界條件

根據實驗測量在噴嘴上游5mm處的空氣流場屬于充分發展的湍流，且該位置的氣體流動基本不受射流流動的影響，氣體湍流強度為4%，根據冪次定律，空氣入口的速度分布滿足[4]：

δT>dw→u=uinf

(10)

式中：δT為速度邊界層厚度，mm，依據實驗測量，入射孔板上的邊界層厚度約為3 mm；dw為入口平面內任意點到壁面的最短距離,mm；Uinf為主流區流速。液體入口給定平均速度。出口采用壓力出口，假設出口各個物理量零梯度。壁面采用了無滑移壁面邊界條件。

2)數值方法

計算采用的求解器基于OpenFoam-7開發，采用有限體積方法離散控制方程，相體積分數方程的對流項采用Gauss vanLeer格式，其余方程對流項采用bounded Gauss Gamma格式，擴散項采用中心差分格式，湍流模型采用k-omega SST VLES模型[16]。速度-壓力修正采用PIMPLE算法。計算中控制最大庫朗數為0.6，計算中所有線性方程組的歸一化誤差設為10-8。所有的拉格朗日粒子統計時采用如下策略：計算射流穩定后的3個時間步內的粒子在邊長為1.5 mm的小正方體內的屬性平均值。

1.5 網格無關性驗證

為了排除網格對于計算結果的影響，采用了稀疏、中等、稠密的三套網格進行計算，網格的總數量分別為75萬、143萬、300萬，網格在壁面和噴嘴附近進行了加密。三套網格在壁面和噴嘴附近的最小尺度分別為0.044D、0.021D和0.01D，D為噴嘴直徑。

圖2(a)為三種網格在平面z=0內相體積分alpha為0.5的界面等值線，從圖中可以看到中等規模的網格下，計算得到的相界面和密網格下的結果差距不大，主要的差異處于破碎點附近，而稀疏網格計算的相界面的高度明顯偏低。為了檢查拉格朗日粒子統計值與網格大小的關系，統計了噴嘴下游20 mm附近一個長和寬為1.5 mm，沿著射流入射方向，高為30 mm的長方體內的液滴速度，圖2(b)展示了不同網格下，在噴嘴下游20 mm處液滴速度的概率密度統計，可以看到稀疏網格下液滴速度分布概率密度最大值略低于中等網格和密網格計算所得的最大值，中等網格和密網格計算得到的液滴速度分布基本一致。結合上述討論，在后面的計算中將采取中等網格進行計算。

圖2 不同網格下的計算結果

2 計算結果與討論

Elsam[17]和Yogish Gopala[18]通過實驗測量了高溫高壓下的射流軌跡，圖3為Case 1和Case 2模擬獲得的射流軌跡和Elsa和Yogish Gopala獲得的經驗關系式對比圖。模擬獲得的Case 1和Case 2射流軌跡高度與上述經驗公式相比都偏低，Case2更加明顯。可能的原因是Elsam和Yogish Gopala的測量的動量通量比的最大值為40，模擬的Case的動量通量比為180，動量通量比代表了液體射流和氣相橫流的慣性力之比，對射流軌跡的影響很大。另外模擬的射流的邊界條件與實驗真實情況存在一定的偏差，也是造成射流軌跡偏低的重要原因。

圖3 射流軌跡對比圖

圖4 Case 1和Case 2在y=0截面上的液滴速度分布

圖4為Case1和Case 2 在y=0截面上的液滴速度分布圖，液滴大小為實際大小的24倍，并采用液滴的速度大小著色，圖中灰色部分為液柱，可以清晰地看到液滴從液柱背面剝離的表面霧化現象，且越遠離壁面的液滴其直徑越大，這可以通過液體湍流是驅動一次霧化的動力因素來解釋，噴嘴出發到達更遠的地方的湍流渦旋的時間尺度更大，相應的其長度尺度更大，而產生的液滴直徑又和液體湍流渦旋的長度尺度成正比，因此離噴嘴越遠的液滴的直徑也就越大。大液滴主要集中于液柱破碎點后，但是最大液滴出現破碎點下游，這說明這些更大的液滴不是由于霧化破碎產生而是液滴聚合產生。另外Case2的液滴直徑明顯小于Case1的液滴直徑，雖然Case1 和Case2有相同的動量通量比，但是Case2的液體入射雷諾數更大，液體湍流的長度尺度更小，在湍流驅動霧化假設下產生的液滴直徑相對就較小，Case2的氣體韋伯數相對于Case1的更大，這意味著Case1中空氣動力學作用更強，在空氣動力學作用下，將產生的更小的液滴。從液滴速度來看，靠近噴嘴所在的壁面的液滴的速度往往較小，而遠離噴嘴所在壁面的液滴的速度較大。這是由于表面破碎產生的液滴和液柱背面的The Wake[29]區域進行動量交換，速度迅速和當地氣流速度一致，在The Wake區域內保持較小的速度向下游運動，而遠離噴嘴所在壁面的大液滴直接受到來流氣流的影響，液滴和來流氣流進行動量交換，液滴速度逐漸和來流速度相同，因此這些大液滴的速度更大。

圖5 Case 1和Case 2在x=20 mm處的液滴D32、D10和液滴速度統計圖

圖5(a)和圖5(b)給出了Case1和Case2在噴嘴下游20 mm處的液滴的D32和D10的統計信息與實驗測量值的對比，液滴直徑分布和圖4中提及的現象一致，隨著液滴遠離壁面液滴直徑逐漸增大，且Case 2產生的液滴直徑比Case 1小。圖5(c)給出了Case 1和Case 2在噴嘴下游20 mm處液滴速度大小的統計信息與實驗測量值對比，為液滴速度，和實驗數據一樣僅考慮空氣流動方向和射流入射方向的速度，為來流空氣速度。可以看到液滴的速度大小和來流空氣速度大小有明顯的區別，其最大的特點為液滴的速度相對于氣體來流速度有較大的滯后，Case1 和Case 2在噴嘴下游20 mm處液滴的速度大小的分布十分相似，Case 2的液滴的速度滯后率相對較小，最大滯后位置相對于Case 1偏后。這說明動量通量比對液滴的速度分布起決定性的影響。在y/D大于40后，模型預測的D32明顯比實驗值低，可能是由于射流穿透深度預測偏低引起。Case 2中液滴的速度最大滯后位置預測相對于實驗值偏后，除了液滴直徑在y/D大于30后與實驗有一定差別外，液滴直徑和液滴速度分布與實驗基本一致。

3 結 論

本文在混合歐拉拉格朗日框架下，針對低液氣密度比同時受到液體湍流和空氣動力學影響的橫向射流，引入了一次霧化的液滴直徑的修正，發展了新的一次霧化模型。對液氣密度比為312、空氣入口馬赫數為0.2和0.35以及液體射流和空氣主流的動量通量比為180的直射式噴嘴的橫向射流霧化開展了模擬，結果表明模擬射流軌跡與實驗經驗關系式相比偏低，除了液滴直徑在y/D大于30后與實驗有一定差別外，液滴直徑和液滴速度分布與實驗基本一致。

致謝

感謝中國科學技術大學超級計算中心對本文數值模擬計算的支持。