超聲速氣流中液體噴霧流動數值模擬研究

冼錦宇，陳錢，蔡樹峰

中山大學，廣東 深圳 518107

超燃沖壓發動機作為高超聲速飛行器的關鍵技術之一[1-5]，具有很高的工程應用價值。液態燃料是超燃沖壓發動機目前最主要的燃料，與氫氣相比，液態燃料具有能量密度高、環境友好且儲存簡便等特點。液態燃料在燃燒室內的霧化過程是燃料燃燒的第一步，而液態燃料霧化質量的好壞直接影響燃燒室內的燃燒過程。為了了解液態燃料的霧化特性，較多研究者對液體橫向射流霧化特性進行了研究。

Lin 等[6]對高速橫向氣流液體噴注的霧化特性和流動特性進行了系統的試驗研究，利用陰影成像技術和多普勒粒子分析儀（PDPA）對不同介質（水與酒精）垂直注入亞聲速氣體橫流的霧化特性與穿透高度進行研究，分析不同介質對穿透高度的影響并建立相應的經驗公式。之后，Lin等[7]進一步探索了液滴噴霧在超聲速流場中的流動特性，開展超聲速空氣氣流中的水射流試驗；通過PDPA 技術測量不同自由流位置的液滴和噴霧羽流中心線與半平面的特性，研究表明，射流噴霧過程受到噴嘴尺寸、氣液動量比和噴射角影響，并給出相當準確的液體射流穿透高度的關聯函數；PDPA測量結果顯示，不同橫截面的中心線上液滴平均直徑以及無量綱速度分別呈現出S形分布和鏡像S形分布。Perurena等[8]利用高速攝影和紋影技術對馬赫數為6的超聲速氣流中的液體橫向射流進行了試驗研究，分析對比了不同噴口形狀、氣液動量比對噴霧穿透高度以及橫向擴散的影響，并提出相應關聯式。

Ⅰm等[9]基于修正后的KH/RT混合破碎模型對馬赫數為1.92 的超聲速橫向射流水噴注試驗進行了數值模擬，并與Lin 等的試驗進行對比，計算通過與試驗中的噴霧特性（如噴霧穿透高度、液滴尺寸和流向速度）的定量對比，分析了新的混合破碎模型與試驗結果的差異，結果顯示，即使在流向速度上仍與實際結果有一定差異，但相較于TAB破碎模型，修正后的KH/RT混合破碎模型的噴霧高度更符合實際結果。楊東超等[10-11]基于Eulerian-Lagrangian方法研究了4種不同二次破碎模型對超聲速氣流中液體橫向射流的霧化效果，指出KH/RT的混合模型計算出的穿透高度與試驗更符合，之后他們還改進了KH/RT 混合模型，擴大了改進后的KH/RT混合模型的適用范圍，同時指出液滴實際追蹤數量對計算結果有一定影響。

總體而言，國內外研究者對超聲速氣流中液體橫向射流的許多方面都進行了較為細致的研究。然而，超聲速氣流中液體橫向射流的數值模擬研究中的初始液滴粒徑大多數采用固定值[12-14]，實際上，初始液滴粒徑并非均勻[15]，初始液滴粒徑分布可能對液滴霧化效果具有影響，需要開展較為細致的研究。因此，本文對超聲速氣流中不同初始液滴粒徑分布的橫向水射流進行了數值模擬研究，并與試驗進行了比較。

1 數值模擬方法

1.1 氣相控制方程

將基于歐拉框架的氣相大渦模擬方法與基于拉格朗日框架的液滴追蹤方法相耦合，對超聲速氣流中液體噴霧流動過程進行數值模擬。氣相由Navier-Stokes方程控制，通過Favre平均方法過濾，其流動的三維可壓縮氣相控制方程定義如下。

連續方程

式中，源項S?m，S?ui，S?h，S?Ym分別表示質量、動量、能量、組分輸運在氣液相間的相互作用。

1.2 液相控制方程

流場中液滴的運動根據拉格朗日液滴追蹤方法進行模擬，液滴位置矢量、速度矢量、溫度和質量的控制方程如下。

式中，ρf為氣相密度；Ap為液滴在氣流中的有效截面；CD為液滴的阻力系數；US為氣液兩相相對速度。

液滴的阻力系數公式定義為

1.3 幾何模型與物理模型

本文的計算幾何模型與Lin 等[7]開展的超聲速橫向氣流中液體水射流噴注試驗一致。圖1給出了該模型的示意圖，模型是一個長方體，總長度為200mm，橫截面恒定為40mm×40mm，水射流通過距離空氣入口50mm中心底部直徑0.5mm 的圓形噴孔垂直噴入空氣中。整體計算域為200mm×40mm×40mm。

圖1 超聲速橫向氣流中水射流計算域 (單位:mm)Fig.1 Computation domain of water jet in supersonic transverse airflow

計算采用結構化網格，壁面第一層厚度為1×10-5m，對射流噴口以及底壁面附近進行網格加密，網格總數量約為105萬個。數值模擬的每個算例在96個CPU核心上運行4天。空氣來流和水射流的入口流動參數與Lin 等[7]的試驗一致，具體參數見表1。

表1 來流空氣和射流的入口流動參數Table 1 Inlet flow parameters for incoming air and jets

對于初始液滴粒徑均勻分布的算例，初始液滴以恒定100μm 注入超聲速氣流中；而對于采用Rosin-Rammler 分布的算例，其初始液滴最小直徑、最大直徑、平均直徑分別為15μm、150μm、100μm，擴散參數為3.9。

在超聲速氣流中的液體橫向射流中，液滴存在變形、振蕩、線性失穩過程，KH模型考慮氣液界面較強剪切作用導致的KH 不穩定波增長過程，能夠較好地模擬在噴口附近的小液滴從初始大液滴上剝離的過程；而當液滴的加速度較大時，RT 不穩定性則主導液滴破碎過程。基于上述分析，本文采用Patterson[16]建立的KH/RT 混合破碎模型進行霧化數值計算。

液滴阻力系數影響液滴破碎的計算，對霧化過程有重要影響。動態阻力模型能夠在液滴形狀發生變化時確定液滴阻力系數，并且能夠適應任何破碎模型，因此本文采用動態阻力模型計算液滴阻力系數。

2 計算結果與討論

2.1 噴霧流動基本特性

為了探索超聲速噴霧流動的基本特性，本節對由初始液滴呈均勻分布的算例計算得到的液霧場進行分析。該計算得到超聲速氣流中橫向射流液滴分布的三維液霧場如圖2所示，液滴由速度經過歸一化處理后著色。從圖2中可以看到，當液滴從噴孔噴出后受到垂直的超聲速氣流影響，液滴速度增加并發生方向變化。液滴在橫向和側向迅速擴散，外圍液滴由于與氣流相互作用更為強烈，其速度比液滴中心區域的更大，近壁面液滴由于受邊界層的影響，其速度比其他外圍液滴速度小。此外，采用大渦模擬方法能夠觀測到液霧場外邊界存在凸起結構（黑色箭頭標注），這與Lin等[17]的試驗結果一致。圖3給出噴孔附近區域液滴直徑分布，同時給出了氣相流場中心平面經過歸一化的壓力、橫截面經過歸一化的流向速度、水平面流向渦量。從圖3 中可以看到，噴注的初始液滴受到氣流影響液滴直徑逐漸減小。液滴與氣流相互作用后，液霧場前段形成激波，氣流在液滴經過區域，其流向速度降低，氣流在液霧場外圍近壁區形成了低流向渦量區，這些特性與Li等[18]的計算結果一致。

圖2 液霧場斜視圖Fig.2 Liquid spray field oblique view

圖3 噴孔附近區域液滴直徑分布Fig.3 Droplet diameter distribution in the area near the nozzle hole

本文計算結果為非定常的，且通過統計平均獲得了工程應用中關心的平均結果。統計平均包括瞬時空間平均和時間平均兩個步驟[19]。首先對不同時刻網格控制體內所有液滴進行空間平均，得到瞬時空間平均結果

式中，＜Un＞為液滴瞬時參數空間平均值；Nn為第n時刻控制體內液滴總數；Dkn表示第n時刻控制體內第k個液滴直徑。

由于不同時刻的網格控制體得到的平均結果可能會發生變化，因此基于空間平均值在時間上進行加權平均得到時間域上的平均值

噴霧的穿透高度是被用于衡量射流與超聲速氣流混合程度的重要參數。Lin 等通過PDPA 技術獲得水射流穿透高度經驗關系式，這里利用該關系式與計算的結果做對比。液霧穿透高度經驗關系式為式中，h為液霧穿透深度；d0為噴孔直徑大小；x則表示距射流噴孔中心的流向距離。

圖4 給出數值模擬與試驗的穿透高度，其中紅線表示試驗結果，黑線表示數值模擬液霧邊界的平均結果。可以看到，在x＜40mm時，平均液霧邊界數值計算結果與試驗結果完全吻合；在40mm＜x＜120mm 時，平均液霧邊界數值模擬與試驗結果非常相近，偏差較小；而在x＞120mm時，平均液霧邊界的數值模擬與試驗結果較為相近，出現較小程度的偏離。這說明數值計算的穿透高度結果符合試驗結果。下面將在不同橫截面上分析液相計算結果。

圖4 均勻分布算例的中心截面液滴穿透高度Fig.4 Droplet penetration height at the central plane of the case with the uniform distribution

圖5 是不同橫截面上水液滴的位置分布，圖中顯示的液滴大小是以液滴直徑按一定比例放大的。可以看到，液滴群整體呈現葫蘆形狀。隨著射流向下游發展，液滴群范圍由中心位置向四周延伸。兩個橫截面上直徑較大的液滴主要集中在液滴群中心；越到液滴群外圍，液滴粒徑越小，這是由于在邊緣處液滴與來流氣體相互作用更強烈，導致邊緣液滴破碎或蒸發。

圖5 均勻分布算例的橫截面液滴位置分布Fig.5 Position distribution of droplets in the cross-section of the case with the uniform distribution

在兩相流數值模擬的研究中，當液相體積分數大于10-3時，氣相湍流運動會受到液滴之間的相互碰撞影響[16]。圖6是經過統計平均后得到的不同橫截面上水液滴的體積分數。可以看到，液霧邊界附近的體積分數均小于10-3；液霧場上部分和中心截面上的體積分數較大，靠近壁面的液霧場體積分數較小；不同截面上整體變化趨勢相近，而在下半中心線上的液滴體積分數x=50mm處比x=100mm處高。

圖6 均勻分布算例橫截面液滴平均體積分數分布Fig. 6 Averaged volume fraction distribution of droplets in the cross-section of the case with uniform distribution

圖7 是不同橫截面上水液滴的平均速度分布，其中藍線是液霧邊界。可以看到，由于受到外邊緣氣流相互作用的影響，液滴平均速度由中心向邊緣不斷增大。隨著水射流向下游發展，液滴群范圍增大，橫截面液滴群受到的氣流作用增強，液滴群整體速度也增大，液霧場邊緣的平均速度增長明顯。

圖7 均勻分布算例橫截面液滴平均速度分布Fig.7 Averaged velocity distribution of droplets in the crosssection of the case with uniform distribution

Sauter 平均直徑（SMD）又稱表面積平均直徑，是一種衡量噴霧質量的重要參數，它的意義是與實際的顆粒具有相同表面積的球體的直徑，其定義可表示為

圖8是不同橫截面上水液滴的SMD值分布，其中藍線是液霧邊界。可以看到，不同截面之間的SMD值變化規律較為相近；SMD值最高接近35μm，主要集中在液霧場上部分和中心截面上，其余區域SMD值在10～20μm之間；同樣受氣流相互作用的影響，SMD值與平均速度的變化趨勢相反，由中心向邊緣處不斷減小，這說明中間部分液滴的破碎不夠充足。

圖8 初始液滴直徑均勻分布的算例的橫截面液滴SMD值分布Fig.8 SMD value distribution of droplets in the crosssection for the case with the uniform distribution

2.2 初始液滴直徑影響

本節將探索初始液滴直徑呈不同分布對液霧場的影響。圖9 是初始液滴直徑呈Rosin-Rammler 分布的算例的中心截面液滴瞬時穿透高度。可以看到，初始液滴直徑呈Rosin-Rammler 分布的算例與均勻分布的算例的穿透高度較為相近，但在x=140mm 附近Rosin-Rammler 分布的算例更接近試驗值。

圖9 Rosin-Rammler分布算例的中心截面液滴穿透高度Fig.9 Droplet penetration height at the center plane of the case with the Rosin-Rammler distribution

圖10～圖13 分別是初始液滴直徑呈Rosin-Rammler分布的算例的橫截面液滴位置分布、平均體積分數分布、平均速度分布和SMD值分布。可以看到液滴整體位置和平均速度變化較小，這說明初始粒徑的不同分布對液霧場的整體位置和平均速度影響不明顯；體積分數在不同截面的變化趨勢一致，但在液霧場上部分和中心截面上Rosin-Rammler 分布算例的體積分數比均值分布算例的小；x=50mm 橫截面上的SMD 值變化不明顯，而在x=100mm 橫截面上液霧場上部分的SMD 值明顯降低，說明采用Rosin-Rammler分布可以使液霧場上部分區域粒子破碎更充分。

圖10 Rosin-Rammler分布算例的橫截面液滴位置分布Fig.10 Position distribution of droplets in the cross-section of the case with the Rosin-Rammler distribution

圖11 Rosin-Rammler分布算例的橫截面液滴體積分數分布Fig.11 Volume fraction of droplets in the cross-section of the case with the Rosin-Rammler distribution

圖12 均勻分布的算例的橫截面液滴平均速度分布Fig.12 Average velocity distribution of droplets in the cross-section of the case with the Rosin-Rammler distribution

圖13 Rosin-Rammler分布算例的橫截面液滴SMD值分布Fig.13 SMD distribution of droplets in the cross-section of the case with the Rosin-Rammler distribution

圖14 給出中心截面不同流向位置的液滴平均速度法向分布曲線，其中液滴平均速度通過來流氣體聲速進行歸一化處理，y坐標通過當地穿透高度進行歸一化處理。可以看到，數值模擬結果整體呈現鏡像S 分布，與試驗結果的變化趨勢相一致，y/h＜0.4時，數值模擬結果與試驗值較為吻合；而當y/h＞0.4 時，數值模擬的值隨y/h的增大而緩慢提高，而試驗結果變化較為劇烈。如上述分析，不同初始液滴直徑分布對液滴速度的影響較小，因而值呈現相似變化。圖15 給出中心截面不同流向位置的SMD 法向分布曲線，y坐標通過當地穿透高度進行歸一化處理。從圖15 中可以看到，在y/h＜0.4 的近壁區域，計算結果與試驗變化趨勢一致；而在y/h＞0.4的中部區域，計算結果與試驗結果差異較大，這可能歸因于中部區域二次破碎不夠充分；而在液霧場上部，數值模擬結果與試驗結果相近。雖然計算的SMD 分布曲線與試驗給出曲線仍有一定差距，但仍然能夠看到初始液滴直徑不同分布對SMD 值的影響，且Rosin-Rammler分布比均勻分布能獲得與試驗更為接近的曲線變化趨勢。

圖14 中心截面液相速度平均分布Fig.14 Average velocity distribution of liquid phase in the central plane

圖15 中心截面液相SMD分布Fig.15 SMD distribution of liquid phase in the central plane

3 結論

本文采用數值模擬方法對超聲速氣流中液體噴霧流動進行了數值模擬研究，獲得以下主要結論:

（1） 基于Eulerian-Lagrangian兩相流計算架構、考慮氣液雙向耦合、采用KH/RT液滴二次破碎模型計算液滴霧化過程、采用大渦模擬計算氣相流動的數值模擬方法能獲得較準確的計算結果，如復現試驗中發現的外邊界凸起結構，獲得與試驗接近的液滴穿透高度、液滴平均速度分布等液霧特性。

（2） 液滴群之中直徑較大的液滴集中在中心區域，越向邊緣位置液滴直徑越小；而液滴速度呈相反趨勢，并且隨著液滴向下游發展整體速度逐漸提高；體積分數在不同截面的變化趨勢一致，但在液霧場上部分和中心截面上Rosin-Rammler 分布算例的體積分數比均值分布算例的小；在x=100mm 橫截面上液霧場上部分的SMD 值明顯降低，說明采用Rosin-Rammler分布可以使液霧場上部分區域粒子破碎更充分。

（3）不同液滴初始直徑分布對液滴的平均速度影響較小，而對SMD值分布影響較大，Rosin-Rammler分布比均勻分布獲得的SMD曲線變化趨勢與實驗更為接近。

（4）不同液滴初始直徑分布對霧化特性的影響規律可為超燃沖壓發動機液體燃料霧化質量的評估和改進提供重要參考。