基于歐拉-拉格朗日模型的內(nèi)燃機燃油噴霧數(shù)值模擬

李 亮，解茂昭，賈 明

(大連理工大學(xué)能源與動力工程學(xué)院，大連116024)

眾所周知，噴霧的初次霧化發(fā)生在液核周邊.在這一區(qū)域，流體以連續(xù)相的狀態(tài)射入，氣液兩相之間存在強烈的相互作用，而并非離散液滴法(DDM)中所描述的以眾多液體包的形式射入.而且，實驗中在該區(qū)域很少發(fā)現(xiàn)球形的液滴，而只能觀測到1個氣液界面為高度褶皺的連續(xù)相流體.因此，嚴格地說，DDM 框架下的多種亞模型[1-3]與真實噴霧存在較大差別，且在數(shù)值模型建立過程中其對入射粒子量、噴射錐角等參數(shù)的設(shè)定具有很強的人為經(jīng)驗性.對噴霧近場研究比較可行的方法是連續(xù)歐拉法.此方法中，把氣液界面的面積密度定義為一個函數(shù)，并導(dǎo)出精確的輸運方程.通過輸運方程的求解來描述氣液界面和噴霧下游離散液滴的行為.Vallet等[4]基于這種思想，首次提出了液體射流的 ELSA(Eulerian-Lagrangian spray and atomization)霧化模型；此模型主要針對的是高雷諾數(shù)和高韋伯?dāng)?shù)的流體計算，并逐步演化出多種應(yīng)用模型.

Blokkeel等[5]采用Vallet提出的方法對整個噴霧過程進行了數(shù)值研究，模擬了由近場連續(xù)相到遠場離散相的轉(zhuǎn)變，即假定兩個液滴間的平均自由間距達到平均液滴直徑的2倍時發(fā)生模型的切換，最終得出了非常滿意的結(jié)果，而且當(dāng)燃燒室氣體密度增加時模型同樣可以很精確地預(yù)測噴霧錐角的變化.Beau等[6]同樣采用了均相模型，假設(shè)在液相體積分數(shù)小于 0.1時發(fā)生歐拉到拉格朗日的轉(zhuǎn)換.Lebas等[7]對比ELSA模型和DNS模型的結(jié)果來確定模型的常數(shù)和參數(shù)值，并加入了蒸發(fā)模型，考慮了切應(yīng)力、破碎、碰撞、汽化等因素對界面密度輸運方程的影響，同時研究了不同環(huán)境氣體密度、溫度、射流壓力等條件下ELSA模型適用性.Ning等[8]在Vallet基礎(chǔ)上進一步增加了平衡蒸發(fā)模型，并同時修正了標準 κ-ε湍流模型，以考慮可壓縮性和渦管拉伸的影響.

Belhadef等[9]通過修改界面密度輸運方程中的生成和耗散項，利用歐拉均相模型對軸對稱旋流噴射進行模擬，結(jié)果顯示索特平均直徑與實驗值能很好地匹配.García-Oliver等[10]研究和驗證了蒸發(fā)和無蒸發(fā)條件下可壓縮Σ-Y模型，通過修改標準κ-ε湍流輸運方程中的耗散系數(shù)，可以精確地預(yù)測不同入射壓力和環(huán)境氣體密度下液體和噴霧的貫穿距、軸向速度和混合質(zhì)量分數(shù)分布.而且，隨著環(huán)境氣體密度和入射壓力的增加，預(yù)測精度也相應(yīng)增加.

從先前的研究可見，歐拉-拉格朗日模型在噴霧中的應(yīng)用相對較少.本文在先前基于 OpenFOAM 所構(gòu)建的歐拉模型[11]的基礎(chǔ)上，參照Vallet理論和前人的研究，構(gòu)建歐拉-拉格朗日模型.在Vallet的最初模型中，壓力是通過狀態(tài)方程直接求解[4]；在當(dāng)前構(gòu)建的模型中，則采用 PISO算法來迭代求解壓力泊松方程.由于該模型假設(shè)液相向氣體環(huán)境的彌散主要是由湍流擴散控制，故湍流模型的選擇將對計算產(chǎn)生影響.因此，本文對動量方程分別進行了RANS和LES封閉，以考察湍流模型的影響，其中RANS封閉模型見文獻[11].最后，應(yīng)用此模型對柴油噴射特性進行數(shù)值研究，模擬對象為 ECN(engine combustion network)噴霧 A[12]2106753#試驗，此噴霧已廣泛應(yīng)用于發(fā)動機各種噴霧模型的驗證.

1 理論模型

Vallet等[14]在 1998年提出了一種新的歐拉噴霧計算模型，此模型中同樣包含對離散液滴相的求解；因此，也可被稱作歐拉-拉格朗日模型.此模型的基本出發(fā)點是：高雷諾數(shù)和高韋伯?dāng)?shù)的高速液體湍流射流的霧化過程，與具有大密度差的氣體射流的湍流混合過程非常相似.從計算的觀點看，Vallet模型把噴霧視為液體和氣體混合物組成的一種單一的“有效流體”.液相向氣體環(huán)境的彌散主要是由于湍流擴散控制；對于小尺度的離散項，則通過一個局部的液體質(zhì)量分數(shù)以及氣液界面面積密度值來描述；此二者的值則通過求解其各自的輸運方程得出.

事實上，噴霧實驗中發(fā)現(xiàn)在噴嘴的出口存在一個連續(xù)的界面高度褶皺的液核，燃料的初次霧化發(fā)生在液核周邊，2次及多次霧化則在噴射的相對下游.正是基于這一現(xiàn)實，Vallet模型在接近噴嘴的稠密區(qū)采用歐拉方法描述氣液兩相行為，而在噴嘴下游的稀薄區(qū)則采用拉格朗日方法來描述.

模型中，歐拉相的基本控制方程的求解是采用均相模型方式，詳細表達式見筆者前期研究[11].湍流求解本文采用 LES方法，其中亞格子模型選用單方程湍動能模式，模型具體的描述見文獻[13].模型中兩相混合物密度按式(1)計算[4]：

式中：Y是當(dāng)?shù)匾合噘|(zhì)量分數(shù)，Y=1為純液相，Y=0為純氣相；ρl，ρg分別為液、氣相密度.氣相密度采用理想氣體狀態(tài)方程求解，液相密度采用線性的可壓縮方程求解.

其中φl，p0分別為液相可壓縮常數(shù)和初始參考壓力.

模型中質(zhì)量分數(shù)的輸運方程為

方程右側(cè)為液相的湍流擴散，即湍流引起的液、氣相之間的質(zhì)量彌散；其可由 Demoulin[14]模型進行封閉：

其中 Sct，Dp分別為 0.9 和 1.8.

對于拉格朗日離散相特性，需要借助氣液界面面積密度輸運方程[4]進行求解.具體方程如下：

此界面面積密度(Σ)則代表了每一個網(wǎng)格單元中液相的分裂程度的亞網(wǎng)格信息.輸運方程推導(dǎo)過程則采用了與火焰面密度輸運方程類似的方法；其中，方程右側(cè)的源項，考慮了表面張力、湍流拉伸、液滴碰撞和聚結(jié)等因素的影響.方程中各相關(guān)變量的含義如下：

DΣ為一擴散系數(shù)，其值等于湍流黏性除以施密特數(shù).

A為由于平均速度梯度引起的氣液界面的拉伸生成項即

a為考慮液滴破碎和碰撞而引起的界面的生成項，即

Vs為液滴碰撞和聚結(jié)引起的界面的耗散項：

式中：req為湍流特征時間尺度(lt)內(nèi)面密度達到平衡值時對應(yīng)的液滴的平衡穩(wěn)定半徑，即

拉格朗日液滴的尺寸和液滴數(shù)通過當(dāng)?shù)氐囊合噘|(zhì)量分數(shù)和界面面積密度來確定.具體如下：

2 物理模型及定解條件

由于計算模型中沒有考慮蒸發(fā)作用，因此選擇的試驗?zāi)Ｐ屯瑯右彩堑蜏叵碌腅CN Spray A 21067[12]，試驗條件詳見表 1.物理模型為圓柱體，長 80mm，直徑 50mm，噴嘴的直徑為 89.4μm.網(wǎng)格劃分為六面體正交網(wǎng)格，采用兩套網(wǎng)格系統(tǒng)進行無關(guān)性檢驗.網(wǎng)格劃分方面，噴嘴處沿徑向均勻劃分為 16個節(jié)點；腔室沿軸向和徑向按比例因子為 1.01、1.06劃分為300和200個節(jié)點，網(wǎng)格總量為300萬.

表1 實驗參數(shù)Tab.1 Experimental parameters

壁面采用無滑移絕熱壁面；出口為第 2類邊界，無回流；入口為第1類邊界，入口速度隨時間變化.

插值格式方面，時間項采用 1階 Euler隱式格式；對流項選用 Gauss Gamma差值格式；Laplacian項選用 Gauss linear corrected差值；梯度項為Gauss linear差值；其中，面法向梯度需進行非正交修正.

3 結(jié)果分析

由于ECN[12]實驗所得的數(shù)據(jù)和云圖是1.4ms以內(nèi)平均值；因此，此處的定性分析云圖和下面的定量分析數(shù)據(jù)需要將噴射時刻 1.5ms以內(nèi)的大渦模擬數(shù)據(jù)進行取樣平均(見圖1).體積分數(shù)與質(zhì)量分數(shù)的轉(zhuǎn)換如下：

圖1 噴霧質(zhì)量分數(shù)Fig.1 Mass fraction of spray

首先對噴霧中的歐拉特性進行分析.圖2反映了穩(wěn)定噴射后、軸距6mm以內(nèi)的燃料體積分數(shù)與實驗的對比.其中，實驗圖片為 ECN所得，詳見文獻[15].對比實驗圖片發(fā)現(xiàn)，在 3mm 內(nèi)模擬的預(yù)測結(jié)果較好，只是在末端附近模擬結(jié)果高于實驗值.然而，在 3～6mm 內(nèi)兩者的偏差增大，尤其是在軸距3mm 附近，實驗云圖的峰值較大，模擬值過低地預(yù)測了軸心燃料分布；然而，在下游6mm附近，軸心處的模擬結(jié)果顯示了比實驗值更高的體積分數(shù).

圖2 軸向體積分數(shù)Fig.2 Axial volume fractions

對于離散相的分布，圖3和圖4分別給出了t=0.5ms時刻粒徑D32和粒子量N的分布云圖.從圖3可見，在軸線處粒子直徑沿噴射方向先減小，后又逐漸增大.由此證實在噴嘴出口附近，由于氣液剪切和湍流作用，使得射流發(fā)生分裂破碎，粒徑減小、粒子量增加.同時，由圖3可見，在噴霧的徑向邊緣和末端下游附近，粒子的直徑較大；其增大的原因可能是在上述區(qū)域存在更為劇烈的氣液剪切作用和湍流擾動，造成液滴間的碰撞聚合作用加劇.對應(yīng)于圖4中的變化，是粒子量在噴霧下游減小.因此，可以認為在噴霧的下游，由于流速和湍流擾動的下降使得粒子間的碰撞融合加劇，從而導(dǎo)致粒徑的增加和粒子量的減小.而且，由圖4可見粒子量分布則是沿軸向先增加后減小，這與D32的整個云圖的分布正好相反.

圖5直觀地給出了兩相界面面積密度值的分布云圖.由噴嘴附近的放大云圖可見，在噴嘴出口附近，界面面積密度值接近于零，可以間接證實噴嘴出口附近射流液核的存在.同時，由于界面面積密度值代表的是當(dāng)前計算單元的液相的分裂程度，即其值越大說明此單元網(wǎng)格內(nèi)的分裂霧化效果越有效，平均液滴直徑越小；反之，其值越小，平均粒徑越大.因此，從圖5分布也可說明，在靠近液核的下游，液相的分裂霧化更強，產(chǎn)生的粒子直徑更小；然而，在噴霧的邊緣和下游，其分裂下降，粒徑增加，使得界面面積密度值減小.此結(jié)果與圖3的粒徑分布對應(yīng)；而且，在噴嘴出口的液核周邊，初次破碎的氣液界面為一種連續(xù)褶皺形態(tài)，且沿軸向逐漸向軸線收縮.

圖3 粒徑D32分布云圖Fig.3 Nephogram of D32

圖4 粒子量N分布云圖Fig.4 Nephogram of the number of droplets N

圖5 兩相界面面積密度∑分布云圖Fig.5 Nephogram of two-phase interfacial area density

圖6給出了不同噴射時刻噴霧的貫穿距與實驗值對比；其中，RANS模型為利用兩相密度修正后的湍流模型[11].對比RANS模型和實驗數(shù)據(jù)發(fā)現(xiàn)，在噴霧的初始時刻，LES模型具有更長的貫穿距；然而在0.2ms之后，即噴霧達到穩(wěn)定噴射后，LES與實驗結(jié)果更加吻合.

圖7顯示了軸心線的燃料體積分數(shù)分布.在噴嘴附近，LES模擬與實驗對比顯示了良好的發(fā)展趨勢；而且，展現(xiàn)了噴嘴出口處穩(wěn)定完整液核的存在；同時，實驗值的體積分數(shù)保持在 0.9上下波動，這可能與實驗中噴嘴內(nèi)部的氣穴現(xiàn)象或者周圍的湍流氣相卷吸有關(guān)；然而，在本文的模擬中并沒有考慮噴嘴內(nèi)的空化因素.

圖6 噴霧貫穿距Fig.6 Spray penetrations

圖7 噴霧軸線體積分數(shù)分布Fig.7 Distribution of volume fractions on the axis

在液核下游燃料的體積分數(shù)迅速衰減，LES和RANS模型均很好地捕捉了此衰減過程；但是，LES在軸距 6mm附近卻存在一個明顯的小峰值；這說明LES的衰減過程中在此處存在燃料的部分聚集，這與圖2的結(jié)果相對應(yīng).

為了進一步分析噴嘴出口附近RANS和LES結(jié)果的不同，圖8給出了軸向速度分布.由圖可見，LES在6mm內(nèi)的高速度值區(qū)間較RANS結(jié)果更長，尤其是在6mm處出現(xiàn)一個小的速度峰值；這對應(yīng)于LES模擬 6mm附近燃料分布的波動.此后，LES模擬的速度下降速率過快，速度值低于 RANS模擬結(jié)果，這使得 LES模型噴霧末端的噴霧貫穿距低于RANS模型的預(yù)測結(jié)果.

圖8 噴霧軸心線速度分布Fig.8 Distribution of axial velocity in the spraying process

圖9對不同橫截面處體積分數(shù)的徑向擴散分布進行對比驗證，由于圖2、圖7顯示預(yù)測值與實驗值的差別主要位于軸距 6mm 內(nèi).因此，分別選取軸距為 0.1mm、2mm和 6mm 共 3個不同軸向位置.在軸距 0.1 mm 處，由于位于噴嘴出口處存在液核，體積分數(shù)在液核邊緣存在斷崖式下降；通過對比實驗值可見，LES的結(jié)果與邊緣體積分數(shù)的下降趨勢吻合較好；然而，RANS結(jié)果則呈現(xiàn)出更加尖銳的斷崖式變化.同時，在中心的液核區(qū)，LES的預(yù)測結(jié)果同樣與實驗測量值更接近；然而，RANS結(jié)果略微偏高；同時，在軸心處的 RANS模型曲線的水平線性的跨度更寬.在軸距 2mm 處，同樣在體積分數(shù)下降的徑向邊緣，LES的結(jié)果與實驗值更接近；然而，在軸線處，LES的預(yù)測值略高；RANS的預(yù)測值略低.對于軸距6mm 處，雖然 LES在徑向邊緣處的預(yù)測較好，但是在軸心處的狹小區(qū)間內(nèi)卻存在過高的體積分數(shù)，其值大約為實驗值的 2.5倍；然而，RANS的預(yù)測值雖然在徑向邊緣略高，但在軸心處與實驗值更吻合.

圖9 燃料體積分數(shù)的徑向分布Fig.9 Radial distribution of volume fractions

為了直觀地顯示和分析 LES模擬值在此處出現(xiàn)較高體積分數(shù)值的原因，圖10給出了燃料1.5ms時刻的瞬時質(zhì)量分數(shù)、渦量云圖，以及整個噴射時間平均化處理后的質(zhì)量分數(shù)、渦量云圖.由圖可見，LES模擬的噴霧在 6mm附近存在高的質(zhì)量分數(shù)和強的渦量分布.而且，此處的不穩(wěn)定波動增強，促使噴霧產(chǎn)生初次分離和破碎，形成高密度液體包；由公式(12)知，質(zhì)量分數(shù)與體積分數(shù)相對應(yīng)，同時，考慮到高密度梯度對射流的流體動力學(xué)穩(wěn)定作用，即所謂固壁效應(yīng)，使得射流的進一步不穩(wěn)定性破碎產(chǎn)生延遲.然而，一旦失穩(wěn)發(fā)生，大的稠密的流體包裹將再次分離破碎，產(chǎn)生快速移動的塊狀結(jié)構(gòu)或韌帶，增強了射流末端的局部湍流和不穩(wěn)定波動，從而對氣液混合產(chǎn)生顯著影響，進而影響體積分數(shù)的分布.

圖10 LES模型給出的燃料質(zhì)量分數(shù)Y和渦量Ω 分布云圖Fig.10 Nephograms of mass fraction and vortices given by the LES model

圖11給出了噴霧D32隨噴射時間的發(fā)展變化和分布.隨著時間的發(fā)展，噴霧的下游和邊緣的液滴直徑逐漸增大，液滴的大直徑范圍也逐漸增加.圖12給出了平均化處理后的不同軸距處 D32的徑向分布，同時與García-Oliver的實驗進行對比.圖12(a)顯示在噴霧的下游，液滴直徑的徑向分布沿軸向逐漸增加；同一徑向曲線的最大液滴直徑均位于噴霧的徑向邊緣.圖12(b)中García-Oliver實驗條件為噴射壓力80MPa，環(huán)境密度ρ為 10～40kg/m3；本文中的模擬算例的噴射壓力為 150MPa，環(huán)境密度為 22 kg/m3；由實驗數(shù)據(jù)可見，隨著環(huán)境密度的增加，液滴徑向的直徑逐漸增加，直徑徑向的分布也更加均一.對比實驗值與模擬曲線可見，本文和 García-Oliver的模擬值均小于測量值.在保證環(huán)境密度值近似的條件下，相同軸距處，本文模擬的液滴直徑較 García-Oliver的模擬值更小，徑向分布更寬；這是完全正確的，因為高的噴射壓力必然導(dǎo)致高的流速和湍流強度，液滴的分裂破碎加劇，產(chǎn)生較小的液滴直徑.

圖11 軸向D32外部輪廓Fig.11 Exterior profile of axial D32

圖12 D32的徑向分布Fig.12 Radial distribution of D32

4 結(jié) 論

(1) 對于噴霧中的歐拉特性，除了在軸距 6mm附近軸心線的體積分數(shù)的實驗值存在一定差別外，LES模擬的預(yù)測值與實驗吻合較好，尤其是燃料的徑向分布；同時，也可以捕捉到噴嘴出口附近恒密度液核的存在.對于軸距 6mm 處的偏差，這與大渦模擬下此處高的質(zhì)量分數(shù)和渦量分布有關(guān)，使得液核末端不穩(wěn)定波動產(chǎn)生高密度的液體包，這些高密度液體包的進一步破碎受到氣液密度差的抑制而產(chǎn)生延遲，進而影響液核末端的體積分數(shù)分布.

(2) 對于離散相而言，在噴嘴出口處粒子直徑基本等同于噴嘴直徑；而且，界面面積密度值接近零，這進一步證實了噴嘴出口連續(xù)液核的存在.同時，由于界面面積密度值可表征液體的分裂程度，因此，由其分布可知：在液核末端的下游，液相的分裂破碎增強，產(chǎn)生較小的粒子直徑；而在噴霧的下游，液滴粒子間的碰撞融合加劇，使得液滴直徑略顯增加.對比液滴直徑實驗值與模擬值可見，相同軸距處，本文模擬的液滴直徑較小，徑向分布較寬，這與高的噴射壓力直接相關(guān).