劉晶晶，葉桃紅

(中國科學(xué)技術(shù)大學(xué) 熱科學(xué)與能源工程系，安徽 合肥 230027)

液體燃料橫向射流是亞燃和超燃沖壓發(fā)動(dòng)機(jī)中一種常見且高效的燃料進(jìn)入方式[1]。液體燃料霧化過程對后續(xù)的蒸發(fā)、混合和燃燒過程有著重要影響。在液體燃料橫向射流霧化過程中，氣體韋伯?dāng)?shù)和射流動(dòng)量通量比是影響初次霧化不穩(wěn)定性和射流穿透深度的重要參數(shù)。為了揭示橫向射流霧化破碎機(jī)理，許多學(xué)者對不同條件下的液體燃料橫向射流霧化進(jìn)行了實(shí)驗(yàn)研究[2-4]。盡管對噴霧的實(shí)驗(yàn)研究取得了一定的進(jìn)展，但是在航空航天應(yīng)用中，高溫高壓環(huán)境使得實(shí)驗(yàn)研究非常困難和昂貴，所以數(shù)值模擬得到了發(fā)展。

目前液體射流霧化的模擬方法主要可以分為三大類：歐拉框架，歐拉-拉格朗日框架和混合歐拉-拉格朗日框架。第一類方法通過解析相界面直接求解射流破碎過程，常用的界面捕捉方法為流體體積(VOF)方法[5-6]，水平集(level-set，簡稱LS)[7]方法以及兩者的耦合(CLSVOF)[8-10]方法。但是這種方法要求網(wǎng)格極為精細(xì)以捕捉到液體表面產(chǎn)生液滴的過程，需要的計(jì)算量非常大。第二類方法氣相被看作連續(xù)相在歐拉框架下求解，液體射流被假設(shè)為與噴嘴直徑相同的Blobs大液滴在拉格朗日框架下求解，此種方法常常將二次破碎模型直接擴(kuò)展到初次破碎的求解中，典型的二次破碎模型有WAVE模型[11]，KH-RT(Kelvin-Helmholtz Rayleigh-Taylor)模型[12]，TAB(Taylor Analogy Breakup)模型[13]等，這種方法大大減小了計(jì)算量，但是忽略了真實(shí)的射流是連續(xù)的流體。第三類方法指在歐拉框架下求解連續(xù)的氣液相，在拉格朗日框架下求解離散相液滴，初次破碎過程由模型進(jìn)行?；?，并將該過程生成的液滴從Eulerian場中轉(zhuǎn)化到Lagrangian場中。此方法對氣液相界面直接求解，既可以較為準(zhǔn)確地求解液體射流和氣流之間的相互作用，又可以減少計(jì)算量。

本文基于混合歐拉-拉格朗日框架，采用VOF耦合Lagrangian算法對Lubarsky等人[4]實(shí)驗(yàn)中航空煤油(Jet A)噴入橫向氣流中的霧化過程進(jìn)行了大渦模擬研究，重點(diǎn)分析了射流穿透深度、液滴平均直徑和速度分布以及氣相平均流場結(jié)構(gòu)。

1 數(shù)學(xué)物理模型

1.1 歐拉框架下的控制方程

大渦模擬的基本思路是利用空間濾波方法消除湍流中起耗散作用的小尺度脈動(dòng)，并采用亞格子模型將其封閉，同時(shí)直接求解包含所有湍動(dòng)能的大尺度脈動(dòng)。根據(jù)流體連續(xù)性假設(shè)，氣液兩相連續(xù)流體大尺度運(yùn)動(dòng)的控制方程如下：

1)質(zhì)量守恒方程

本文采取低馬赫數(shù)下不可壓縮假設(shè)，質(zhì)量守恒方程為

(1)

2)動(dòng)量守恒方程

(2)

(3)

本文采用Smagorinsky渦黏模型對亞格子應(yīng)力進(jìn)行?；缡?4)所示：

(4)

(5)

式中：Δ為過濾尺度；Cs為Smagorinsky常數(shù)。Su,i為連續(xù)相和離散相之間的動(dòng)量交換源項(xiàng)。

(6)

式中：nd為一個(gè)液滴包裹(parcel)的數(shù)密度；FD,i為液滴上的拖曳力；Sσ為表面張力在動(dòng)量方程中產(chǎn)生的源項(xiàng)Sσ=σkα，采用連續(xù)表面張力模型[14]?；?/p>

混合連續(xù)流體的密度和動(dòng)力黏度根據(jù)相體積分?jǐn)?shù)加權(quán)計(jì)算，下標(biāo)g代表氣相，下標(biāo)l代表液相：

ρ=(1-α)ρg+αρl

(7)

μ=(1-α)μg+αμl

(8)

體積分?jǐn)?shù)輸運(yùn)方程如下：

(9)

式中：Sα是與Eulerian/Lagrangian框架轉(zhuǎn)化有關(guān)的源項(xiàng)，計(jì)算方法為一個(gè)網(wǎng)格內(nèi)的液相體積分?jǐn)?shù)與時(shí)間步長Δt的比值：

(10)

1.2 拉格朗日框架下的控制方程

離散相液滴運(yùn)動(dòng)規(guī)律遵循牛頓第二定律；采用雙向耦合，即除了考慮氣流對液滴的影響，還將考慮液滴對氣流的影響，通過動(dòng)量方程來描述。液滴的位置和速度通過以下方程得到：

(11)

(12)

式中：下標(biāo)p代表液滴；m，u和x分別為其質(zhì)量(kg)，速度(m/s)和位置(m)；下標(biāo)i代表矢量的三個(gè)分量；FG,i為重力和浮力的合力；FD,i為拖曳力，由式(12)計(jì)算：

(13)

式中：Dp為液滴直徑，mm；ρg為氣相密度，kg/m3；ug,i為在歐拉場中位置的氣相速度；CD為阻力系數(shù)，采用Schiller-Naumann阻力模型[15]計(jì)算。

液滴的二次破碎由Reitz Diwakar[16]二次模型?；旱蔚呐鲎埠腿诤嫌蒒ordin的碰撞算法[17]計(jì)算。

1.3 歐拉和拉格朗日框架之間的轉(zhuǎn)化

本文采用Martin等人[18]給出的算法將歐拉場中射流初次破碎產(chǎn)生的液團(tuán)轉(zhuǎn)化到拉格朗日場中的離散液滴進(jìn)行追蹤。該算法首先采用網(wǎng)格標(biāo)記關(guān)聯(lián)法對歐拉場中的液團(tuán)進(jìn)行識(shí)別，規(guī)定一個(gè)液相體積分?jǐn)?shù)臨界值αt(10-2數(shù)量級)，大于該臨界值的網(wǎng)格將被標(biāo)記為液相網(wǎng)格，否則被標(biāo)記為氣相網(wǎng)格，接下來該算法循環(huán)遍歷所有網(wǎng)格，并將每個(gè)液相網(wǎng)格都標(biāo)記為唯一的網(wǎng)格ID。然后，對于標(biāo)記過的所有液相網(wǎng)格，找到它們的相鄰網(wǎng)格，將它們的網(wǎng)格ID改為相鄰網(wǎng)格中的最小ID，這樣在歐拉場中具有相同ID的液相網(wǎng)格就構(gòu)成了一個(gè)液團(tuán)。

計(jì)算識(shí)別到的所有歐拉場中液團(tuán)的體積、位置、速度和直徑：

(14)

根據(jù)計(jì)算得到的液團(tuán)參數(shù)，設(shè)置以下準(zhǔn)則：

(1)直徑小于某一臨界值dmax，一般為網(wǎng)格尺寸的4～6倍，因?yàn)橹睆捷^大時(shí)，液團(tuán)還會(huì)發(fā)生形變，直徑小于網(wǎng)格尺寸的4～6倍時(shí)，采用VOF方法已經(jīng)不能準(zhǔn)確描述液團(tuán)的形狀和動(dòng)力學(xué)特性。

如果滿足以上準(zhǔn)則，歐拉場的液團(tuán)就會(huì)被轉(zhuǎn)化為拉格朗日離散液滴注入計(jì)算域，同時(shí)將液團(tuán)從歐拉場中刪除(即α=0)。

1.4 幾何模型與算例設(shè)置

根據(jù)Lubarsky等人的實(shí)驗(yàn)裝置，本文數(shù)值模擬采用的計(jì)算域及幾何結(jié)構(gòu)示意圖如圖1所示。燃料入射孔直徑為0.457 mm，燃料管道長度L=10D，位于空氣入口下游5 mm中心線上，空氣入口截面為46 mm×30 mm，坐標(biāo)軸原點(diǎn)位于射流出口圓心位置處，計(jì)算域總長71.2 mm。

圖1 計(jì)算域幾何結(jié)構(gòu)示意圖

采用openfoam平臺(tái)網(wǎng)格工具snappyHexMesh對計(jì)算域進(jìn)行網(wǎng)格劃分，計(jì)算時(shí)對氣液相界面采用動(dòng)態(tài)自適應(yīng)加密，最小網(wǎng)格尺度為0.03 mm。

根據(jù)實(shí)驗(yàn)條件，表1給出了本文開展的四個(gè)仿真算例的詳細(xì)參數(shù)，其中橫向來流為溫度555 K的預(yù)熱空氣，腔內(nèi)壓強(qiáng)維持在4atm。Case1-3的射流動(dòng)量通量比均為40，氣體韋伯?dāng)?shù)分別為33、470和800。此外，為了分析不同射流動(dòng)量通量比下的霧化特性，又設(shè)置4#的射流動(dòng)量通量比為10，氣體韋伯?dāng)?shù)為470。表2為燃料Jet A的物性參數(shù)。

表1 不同算例的詳細(xì)參數(shù)信息

表2 Jet A航空煤油的物性參數(shù)

1.5 數(shù)值計(jì)算方法與邊界條件

采用有限體積方法離散控制方程，相體積分?jǐn)?shù)方程的對流項(xiàng)離散格式為Gauss vanLeer格式，動(dòng)量方程中的對流項(xiàng)為二階Gauss LUST grad(U)格式離散，其余方程對流項(xiàng)為高斯迎風(fēng)(Gauss upwind)格式，其他項(xiàng)(面法向梯度，梯度項(xiàng)，拉普拉斯項(xiàng)等)采用高斯線性(Gauss linear)格式。采用PIMPLE算法進(jìn)行壓力速度的耦合求解，采用動(dòng)態(tài)調(diào)整時(shí)間步，保持庫朗數(shù)Co<1。

對于初始場和邊界條件，給定來流空氣入口參數(shù)作為初始條件，空氣和燃料入口邊界條件均為均勻邊界條件加上4%的湍流脈動(dòng)，壁面邊界采用無滑移條件，計(jì)算域出口采用壓力出口邊界，假設(shè)出口處各個(gè)物理量零梯度。

2 計(jì)算結(jié)果與分析

2.1 射流穿透深度

本文的模擬計(jì)算采用高溫高壓條件，選用的參考實(shí)驗(yàn)未給出射流穿透深度，故選擇Li Lin等人[19]和Yogish Gopala等[20]提出的高溫高壓條件下的射流穿透深度經(jīng)驗(yàn)關(guān)系式，與模擬結(jié)果進(jìn)行對比。Li Lin提出的射流穿透深度經(jīng)驗(yàn)關(guān)系式為

(14)

Yogish Gopala等人提出的射流穿透深度經(jīng)驗(yàn)關(guān)系式為

y/dj=1.393q0.482log(1+1.267x/dj)

(15)

圖2為1#-4#模擬得到的射流穿透深度與上述兩個(gè)經(jīng)驗(yàn)關(guān)系式的對比。由圖2可知，四個(gè)算例模擬得到的射流穿透深度均與兩個(gè)經(jīng)驗(yàn)關(guān)系式對比較好，且隨著射流動(dòng)量通量比增加，射流穿透深度明顯增加。對于相同射流動(dòng)量通量比的1#-3#，1#在X/D>40時(shí)射流穿透深度較大于2#和3#，這是由于1#中來流空氣和射流雷諾數(shù)較小，導(dǎo)致了較小的湍流作用和空氣動(dòng)量，射流破碎程度減弱，射流穿透深度增加。由此看來，盡管射流穿透深度主要受動(dòng)量通量比的影響，但是由于來流空氣和液體射流的雷諾數(shù)影響著液滴破碎、分散和混合，因而會(huì)影響射流穿透深度[21]。此外，在距離射流出口較遠(yuǎn)的下游，模擬的射流穿透深度略高于經(jīng)驗(yàn)關(guān)系式。

圖2 1#～4#射流穿透深度與經(jīng)驗(yàn)關(guān)系式的對比

2.2 液滴平均直徑和速度分布

圖3為1#在x=30 mm截面處統(tǒng)計(jì)結(jié)果與實(shí)驗(yàn)結(jié)果的對比，模擬結(jié)果與實(shí)驗(yàn)基本保持一致，但液滴穿透深度未達(dá)到實(shí)驗(yàn)值，由于模擬結(jié)果在高射流穿透深度處只存在非常少的大液滴，不具備統(tǒng)計(jì)意義，故這里將其忽略。從圖3(a)和圖3(b)可以看出，隨著噴霧穿透深度增加，液滴索特平均直徑(SMD)和算術(shù)平均直徑(AMD)均增加。在較高穿透深度處，液滴SMD和AMD均大于90 μm，表明射流頂部發(fā)生了柱破碎現(xiàn)象，該機(jī)制生成的液滴尺寸較大；在低穿透深度處，液滴的平均直徑在20～50 μm，表明射流柱初次破碎生成液滴的機(jī)制為剪切破碎，該機(jī)制生成的液滴尺寸較小；綜上可知Weg=33時(shí)射流初次破碎機(jī)制為多模式破碎機(jī)制，即包含柱破碎和剪切破碎。在圖3(d)中，壁面附近液滴y方向速度分量Vy為負(fù)值，說明靠近壁面處存在渦旋結(jié)構(gòu)[4]。

圖3 1#在x=30 mm截面處統(tǒng)計(jì)結(jié)果與實(shí)驗(yàn)結(jié)果的對比

圖4為1#～4#在x=30 mm截面處統(tǒng)計(jì)結(jié)果之間的相互對比，對于射流動(dòng)量通量比相同(q=40)的1#～3#，三個(gè)算例的液滴平均直徑分布和速度分布趨勢較為一致，4#的穿透深度明顯低于1#～3#，與2.1節(jié)的分析一致。從圖4(a)和圖4(b)可以看出，2#～4#的SMD和AMD值均在10～60 μm，表明了在較高氣體韋伯?dāng)?shù)時(shí)，空氣動(dòng)力作用較強(qiáng)，破碎生成的液滴直徑較小，液滴的形成機(jī)制為剪切破碎。在靠近射流壁面處，液滴平均直徑較小，這些小液滴在空氣動(dòng)力的作用下從液柱表面剝離[3]，隨著距壁面距離的增加，液滴平均直徑逐漸增大。由于4#的噴霧穿透深度較小，在同一深度，4#的平均直徑較大于1#～3#，對于1#～3#，隨著氣體韋伯?dāng)?shù)增大，同一深度的液滴平均直徑減小。圖4(c)中，y/D在15～20附近時(shí)，存在一個(gè)液滴平均速度較小的區(qū)域，這個(gè)區(qū)域?yàn)槲擦鲄^(qū)，主要是由于液柱表面波與空氣的動(dòng)量交換形成的，尾流區(qū)內(nèi)的液滴會(huì)失去部分的橫向動(dòng)量[3]。此外，在同一深度處，2#～3#的Vx/Vair值大于4#，這主要是由于液滴直徑的不同導(dǎo)致的，當(dāng)剪切破碎形成的液滴較小時(shí)，易被周圍氣體加速，故沿x方向的速度較大。但是這個(gè)規(guī)律卻不適用于1#，1#的平均直徑較2#、3#大，但1#的Vx/Vair值反而大于2#～3#，這主要是由于液滴形成機(jī)制的不同導(dǎo)致的[4]，1#為多模式破碎機(jī)制，2#～4#為剪切破碎機(jī)制。圖4(d)中，在射流壁面附近均出現(xiàn)負(fù)的y速度，說明1#～4#壁面附近均存在渦旋結(jié)構(gòu)，靠近壁面處的小液滴跟隨渦旋結(jié)構(gòu)運(yùn)動(dòng)，一部分向著壁面方向即-y方向運(yùn)動(dòng)。

圖4 1#～4#在x=30 mm截面處統(tǒng)計(jì)結(jié)果相互對比

2.3 氣相平均流場結(jié)構(gòu)

圖5示出了液柱與周圍空氣相互作用的流線圖，流線由氣體渦量的y方向分量進(jìn)行著色，觀察到兩個(gè)算例在射流柱下游均產(chǎn)生了對稱的回流區(qū)結(jié)構(gòu)，2#中回流區(qū)尺度大于4#，這是由于2#的射流動(dòng)量通量比較4#大，射流柱較高，對來流空氣的阻礙作用更強(qiáng)；除此之外，在靠近射流柱底部的下游處，來流空氣繞過液柱，2#中在靠近壁面處形成了尺度相對較小的回流區(qū)，而在4#中并沒有觀察到，這是因?yàn)?#中較大尺度的回流區(qū)產(chǎn)生的渦結(jié)構(gòu)位置較靠近壁面，從而導(dǎo)致渦旋與壁面相互作用，阻礙了壁面附近小尺度回流區(qū)的產(chǎn)生。

圖5 液體射流周圍的空氣平均流線圖

3 結(jié) 論

采用流體體積耦合拉格朗日方法對航空煤油(Jet A)噴入橫向氣流中的霧化特性和流場結(jié)構(gòu)進(jìn)行了大渦模擬(LES)研究，得到結(jié)論如下：

(1)不同氣體韋伯?dāng)?shù)Weg(分別為33、470和800)和不同射流動(dòng)量通量比q(分別為10和40)下的射流穿透深度均與經(jīng)驗(yàn)關(guān)系式對比良好，且隨射流動(dòng)量通量比增加，射流穿透深度增加。盡管射流穿透深度主要受動(dòng)量通量比的影響，但是若來流空氣和射流射流雷諾數(shù)較大，則會(huì)增強(qiáng)液滴的破碎，分散和混合過程，導(dǎo)致液滴穿透深度降低。

(2)通過對液滴平均直徑和速度分布的分析，氣體韋伯?dāng)?shù)為33的模擬結(jié)果與實(shí)驗(yàn)結(jié)果較為吻合，同時(shí)發(fā)現(xiàn)氣體韋伯?dāng)?shù)較小時(shí)(Weg為33)，射流初次破碎機(jī)制多模式破碎；氣體韋伯?dāng)?shù)較大時(shí)(Weg分別為470和800)，射流初次破碎機(jī)制為剪切破碎。

(3)在射流柱下游區(qū)域，氣相流場中存在回流區(qū)結(jié)構(gòu)。