液體燃料橫向射流霧化數值模擬研究

劉晶晶，葉桃紅

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

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

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

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

1 數學物理模型

1.1 歐拉框架下的控制方程

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

1)質量守恒方程

本文采取低馬赫數下不可壓縮假設，質量守恒方程為

(1)

2)動量守恒方程

(2)

(3)

本文采用Smagorinsky渦黏模型對亞格子應力進行模化，如式(4)所示：

(4)

(5)

式中：Δ為過濾尺度；Cs為Smagorinsky常數。Su,i為連續相和離散相之間的動量交換源項。

(6)

式中：nd為一個液滴包裹(parcel)的數密度；FD,i為液滴上的拖曳力；Sσ為表面張力在動量方程中產生的源項Sσ=σkα，采用連續表面張力模型[14]模化。

混合連續流體的密度和動力黏度根據相體積分數加權計算，下標g代表氣相，下標l代表液相：

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

(7)

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

(8)

體積分數輸運方程如下：

(9)

式中：Sα是與Eulerian/Lagrangian框架轉化有關的源項，計算方法為一個網格內的液相體積分數與時間步長Δt的比值：

(10)

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

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

(11)

(12)

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

(13)

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

液滴的二次破碎由Reitz Diwakar[16]二次模型模化，液滴的碰撞和融合由Nordin的碰撞算法[17]計算。

1.3 歐拉和拉格朗日框架之間的轉化

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

計算識別到的所有歐拉場中液團的體積、位置、速度和直徑：

(14)

根據計算得到的液團參數，設置以下準則：

(1)直徑小于某一臨界值dmax，一般為網格尺寸的4～6倍，因為直徑較大時，液團還會發生形變，直徑小于網格尺寸的4～6倍時，采用VOF方法已經不能準確描述液團的形狀和動力學特性。

如果滿足以上準則，歐拉場的液團就會被轉化為拉格朗日離散液滴注入計算域，同時將液團從歐拉場中刪除(即α=0)。

1.4 幾何模型與算例設置

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

圖1 計算域幾何結構示意圖

采用openfoam平臺網格工具snappyHexMesh對計算域進行網格劃分，計算時對氣液相界面采用動態自適應加密，最小網格尺度為0.03 mm。

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

表1 不同算例的詳細參數信息

表2 Jet A航空煤油的物性參數

1.5 數值計算方法與邊界條件

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

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

2 計算結果與分析

2.1 射流穿透深度

本文的模擬計算采用高溫高壓條件，選用的參考實驗未給出射流穿透深度，故選擇Li Lin等人[19]和Yogish Gopala等[20]提出的高溫高壓條件下的射流穿透深度經驗關系式，與模擬結果進行對比。Li Lin提出的射流穿透深度經驗關系式為

(14)

Yogish Gopala等人提出的射流穿透深度經驗關系式為

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

(15)

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

圖2 1#～4#射流穿透深度與經驗關系式的對比

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

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

圖3 1#在x=30 mm截面處統計結果與實驗結果的對比

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

圖4 1#～4#在x=30 mm截面處統計結果相互對比

2.3 氣相平均流場結構

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

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

3 結 論

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

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

(2)通過對液滴平均直徑和速度分布的分析，氣體韋伯數為33的模擬結果與實驗結果較為吻合，同時發現氣體韋伯數較小時(Weg為33)，射流初次破碎機制多模式破碎；氣體韋伯數較大時(Weg分別為470和800)，射流初次破碎機制為剪切破碎。

(3)在射流柱下游區域，氣相流場中存在回流區結構。