孔徑比和動量比對雙股射流撞擊式霧化影響分析①

王宇奇，勾文進，陳明慧，張 帥，鄭 耀

(1.浙江大學 航空航天學院，杭州 310027；2.上海飛機設計研究院，上海 201210)

0 引言

雙股射流撞擊霧化噴嘴具有結構簡單、制造成本低、霧化性能較好、混合效率高等優點，已廣泛應用于液體火箭發動機中。噴嘴的霧化過程會影響后續的蒸發、燃燒過程，繼而影響發動機的燃燒效率與燃燒穩定性。因此，針對雙股射流撞擊式的霧化過程與霧化機理的研究對于發動機的設計具有重要意義。

已有雙股射流撞擊式噴嘴霧化特性研究以實驗手段為主。Heidmann等[1]對射流撞擊進行了大量實驗研究，揭示了不同射流速度、噴嘴直徑以及流體性質對噴霧結構的影響。Dombrowski等[2]對液膜破碎的原因進行了研究，研究表明撞擊產生的不穩定波導致了液膜的破碎，且只有在射流超過臨界韋伯數的時候才會出現，與雷諾數無關。Lai等[3]試驗研究了流體的物理性質對雙股撞擊霧化特性的影響，通過改變流體粘度和表面張力得到了多達10種的噴霧模式。張蒙正等[4-5]用激光全息及圖像處理技術研究了雙股撞擊式噴霧的撞擊夾角、孔徑比和動量比變化對霧化性能的影響，并通過對實驗數據的分析整理獲得了索特爾平均直徑的經驗公式。鄧寒玉等[6]實驗分析了射流自由長度對凝膠推進劑撞擊霧化的影響，結果表明,根據射流速度合理選取自由長度可以獲得更好的霧化效果。劉曉偉等[7]研究了魯泊數和孔徑比對直撞擊式噴注器性能的影響，試驗表明兼顧兩者才能使氧化劑和燃料達到最佳混合效果，提高燃燒效率。

近年來，霧化的數值模擬技術發展很快。與實驗研究相比，數值模擬能夠獲取霧化過程的更多細節，可以與實驗結果相互補充互相驗證，有助于更加深入地研究霧化機理。Inoue Chihiro等[8]運用CIP-LSM方法模擬了雙股射流撞擊霧化，發現液膜的動態特性對下游的液絲和液滴分布有重要影響。Ma Dongjun等[9]利用VOF方法結合AMR技術，對射流撞擊霧化形式和液膜破碎等現象進行了數值模擬。鄭剛等[10]基于CLSVOF方法對雙股互擊式噴嘴霧化進行了數值模擬，并詳細考察了動量比對液膜特性的影響。李佳楠等[11]基于開源程序Gerris完成了撞擊式噴注單元霧化過程的仿真模擬，與試驗數據進行對比，驗證了計算的有效性，給出了數值求解精度。但是，霧化過程的數值模擬往往涉及多相、多尺度流動問題，需要處理大變形、自有界面問題，因此計算量很大。

無網格粒子法是一種全Lagrange方法，在處理具有大變形、自由面等問題的時候，具有網格法無法比擬的巨大優勢。強洪夫等[12-13]將SPH方法探索性地應用于射流撞擊形成液膜和凝膠推進劑一次霧化仿真研究，解決了傳統網格法難以解決的自由面、大變形等問題，模擬結果與實驗結論基本一致。韓亞偉等[14]運用SPH方法對雙股液體射流撞擊霧化問題進行了三維數值模擬。勾文進等[15]將MPS方法應用于直流撞擊式噴嘴霧化模擬，成功模擬了霧化的三個模態，但受限于MPS方法較大的計算量，缺乏定量分析。國內還未見將MPS方法應用于雙股射流撞擊霧化孔徑比和動量比的分析相關文獻。

為解決MPS方法計算量大的問題，可以將GPU并行加速技術應用于MPS方法。GPU是一種新的可應用于大規模并行計算的處理器和計算機集群的計算架構，最初主要作為游戲行業中的圖形計算處理器，隨著統一計算設備架構(Compute Unified Device Architecture，CUDA)的發布，GPU在并行計算中有著巨大優勢，被越來越多地應用于大規模科學計算。基于GPU的高性能計算主要采用異構架構，即CPU+GPU模式。計算過程中使用CPU執行串行工作，控制主程序的復雜流程，將需要批量處理的向量數據傳輸給GPU存儲器，由GPU的眾多處理器執行快速的并行計算。目前在國內公開發表資料中，采用GPU異構并行技術優化加速MPS方法的研究總體較少。

本文基于CUDA開發異構并行加速MPS霧化程序，實現了雙股射流撞擊霧化模擬。通過對典型工況下的霧化現象進行模擬，分析了孔徑比和動量比對霧化特性的影響規律。

1 方法

1.1 MPS方法控制方程與離散

MPS方法是一種完全拉格朗日粒子方法，界面始終清晰，不需要傳統網格方法在模擬過程中存在的界面追蹤或捕捉，在自由面模擬和大變形流動的模擬方面有自己的優勢。

MPS方法中，不可壓縮流體的控制方程為

▽·u=0

(1)

(2)

式中ρ為流體的密度；u為速度矢量；t為時間；F為體積力；p為壓力；ν為運動粘度。

MPS方法基于核函數對控制方程進行離散，其方程為

(3)

式中re為核函數影響域；r為兩個粒子之間的距離。

粒子數值密度為

(4)

式中rj和ri分別為粒子j和粒子i的坐標。

MPS方法采用梯度算子模型與拉普拉斯算子模型離散控制方程。梯度算子模型為

(5)

式中d為空間維數；n0為初始粒子數密度；φi為粒子i的物理量；φj為鄰居粒子j的物理量。

拉普拉斯算子模型為

(6)

其中

(7)

MPS方法采用了SMAC算法對控制方程進行求解。其壓力泊松方程PPE(Pressure Poisson Equation)為

(8)

式中n*為顯式階段粒子移動后的粒子數密度；Δt為時間步長；pn+1為下一時間步壓力值。

壓力求解后，速度修正量為

(9)

最后更新速度和位置：

(10)

(11)

本文在進行顯式計算時，加入表面張力的計算。本文采用一種基于自由能的表面張力模擬[16]，該模型非常易于應用到MPS中。

1.2 GPU并行加速

利用所述的梯度模型和拉普拉斯模型將流體的控制方程離散并求解，求解流程如圖1所示。

圖1 GPU加速MPS方法流程圖

粒子法雖然能夠較容易地模擬自由面和大變形流動等問題，但是計算量太大。一方面，不同于基于歐拉觀點的網格方法，MPS計算過程中粒子自由移動，每個時間步需要對領域粒子進行檢索；另一方面，不同于顯式的SPH方法，MPS方法需要求解大型稀疏矩陣。這兩部分在MPS方法的求解過程中占據時間較長，同時也是程序并行的難點。本文采用Simon[17]提出的Link-list算法實現領域粒子的搜索，該方法不僅可以高效搜索領域粒子，而且優化了后續步驟在GPU內存上的使用效率。為最大限度的發揮GPU的優勢，在程序開發時遵循了包括最小化GPU-CPU數據傳輸、最小化warp內線程分支、優化顯存訪問模式以及多使用共享顯存等原則。為減少設備端主機端數據傳輸，所有計算均在GPU上執行，數據只有需要保存時才會被傳輸到主機端。

GPU加速程序在典型算例中，最大加速比為16，取得了較好的加速效果[18]。

2 結果與討論

2.1 計算模型

根據直流撞擊式霧化試驗，本文數值模擬如圖2所示，夾角為2θ的兩股射流位于XOY平面，以初始速度V沿各自軸線方向運動，在撞擊點O處發生撞擊，進而在YOZ平面形成液膜、液絲和液滴。

圖2 雙股射流撞擊霧化計算模型

本文計算中涉及的物性參數主要有密度、粘性系數以及表面張力系數等，如表1所示。

液滴的索特爾平均直徑(Sauter Mean Diameter，SMD)是燃燒流場中常用的性能評價指標[19]。對液滴粒徑的統計方法為在某一時刻，統計位于撞擊點下游10 mm和20 mm處高度為1 mm的長方體區域內的液滴數目Ni和尺寸Di，如圖3所示。根據索特爾平均直徑的計算公式：

(12)

表1 物性參數總結

得到該區域內的索特爾平均直徑，并近似認為是該平面該時刻的索特爾平均直徑。對不同時該平面的索特爾平均直徑進行算術平均即認為該平面的索特爾平均直徑。液膜破碎長度Lp定義為自撞擊點到液絲從液膜上完全剝離時的長度，見圖3。

圖3 液膜破碎長度和SMD統計方法

噴霧擴散角β定義為垂直于撞擊平面內，在撞擊點下游以撞擊點為頂點的具有一定夾角的霧化區域，如圖2所示。

計算模擬了雙股射流噴注單元在典型工況下的霧化過程，并與試驗數據進行了對比分析。噴嘴的初始直徑為1 mm，孔徑比為1，動量比為1，撞擊夾角為60°，單股射流質量流量為10 g/s，則可計算得到射流速度為12.75 m/s。

如圖4所示，計算得到了兩股射流撞擊形成液膜，液膜在撞擊波作用下發生破碎形成液絲，液絲進一步破碎形成液滴的全過程。圖4給出了霧場達到穩定狀態后的圖像與相同條件下的試驗結果[11]的對比圖，模擬結果與試驗結果吻合較好。

2.2 孔徑比的影響

孔徑比是雙股射流撞擊式霧化的重要工作參數。實際工程中，由于燃料與氧化劑的密度不同，當要求同時滿足推進劑最佳混合效果以及燃燒效率的時候，會將燃料孔與氧化劑孔設計成不同的大小，即雙股射流孔徑比不為1。本文模擬了一系列不同孔徑比的雙股射流撞擊霧化過程。撞擊夾角均為60°，射流的動量比為1，射流孔的尺寸、孔徑比以及射流的質量流量之比如表2所示。

圖4 霧場模擬和試驗對比圖

表2 射流孔尺寸、孔徑比和射流質量流量

孔徑比為1.6時的正面霧場圖片如圖5所示，與相同條件下的試驗拍攝的霧場[11]結構一致：都觀察到了霧場存在明顯的偏斜；射流孔徑比越大，霧場的偏斜程度就越大；液膜的形狀是凹形的，有將孔徑較小的射流包圍的趨勢。

圖5 孔徑比為1.6時模擬結果與試驗結果對比

由表2可知，孔徑較小的射流速度較大，撞擊時的穿透能力較強。雖然兩股射流的噴射動量相同，孔徑較大的射流動量沒有全部用來與孔徑較小的射流發生撞擊，實際發生撞擊的流體動量是不相同的，這導致了合成射流偏向撞擊動量較小的射流，即偏向于孔徑較大的射流。孔徑較大射流的邊緣區域由于未參與撞擊沿著射流方向繼續運動，因此形成了凹形液膜。

孔徑比的變化對噴霧擴散角的影響如圖6所示。可見，孔徑比的增大會引起噴霧擴散角的增大，但是增大幅度不明顯，驗證了文獻[4]的試驗結論。

圖6 不同孔徑比下噴霧擴散角的變化

由圖7可知液膜破碎長度隨孔徑比的增大而增大，而且在大孔徑比的時候增加較快，這是因為孔徑加大，大孔徑射流沒有參與撞擊的部分增加，參與撞擊的射流的相對動量減小，造成了不充分的霧化，因此液膜會在更長的距離內達到破碎臨界點。

圖7 不同孔徑比下液膜破碎長度的變化

對不同孔徑比的霧化液滴的索特爾平均直徑進行了統計，結果見圖8，隨著孔徑比的增加，霧化場的索特爾平均直徑也會增加。射流孔徑比對液滴尺寸的影響是正變的，這與文獻[20-21]結論一致。

由以上分析可知，孔徑比的增加將會導致兩股射流混合的不均勻，同時也導致了霧化變差。因此，考慮雙股射流撞擊式霧化的時候，孔徑比取值不宜過高。

圖8 不同孔徑比下SMD的變化

2.3 動量比的影響

動量比也是雙股射流撞擊式霧化的重要工作參數。當孔徑比為1時，在實際工程中，為滿足需要的混合比時，由于燃料與氧化劑的密度以及射流速度的差異，兩股射流的動量比可能不為1。因此，本文模擬了一系列不同動量比的霧化過程，研究不同動量比對霧化結果的影響。計算中，撞擊夾角為60°，射流孔直徑為1 mm且孔徑比為1，射流動量比以及射流質量流量之比如表3所示。

表3 射流動量比和射流質量流量統計

動量比為1.5、1.75和2.0時霧場的正面模擬結果圖像分別如圖9(a)～(c)所示。圖9的霧場模擬結果與相同條件下文獻[11]試驗拍攝霧場變化趨勢一致，都觀察到了霧場向動量較小的一側出現了明顯的偏斜。且隨著動量比的增加，偏斜程度有所增加。霧場的偏斜角α定義為以液膜和液絲為代表的的霧場主體所在的平面與經過撞擊點且與雙股射流平面垂直的平面之間的夾角。當動量比為2.0時，模擬得到的偏斜角為21°，根據文獻[11]，試驗所得的偏轉角為24°，模擬得到的偏斜角誤差為12.5%。

(a)Momentum (b)Momentum (c)Momentum

對模擬得到的霧場偏斜角度進行測量統計，結果如圖10所示。可以看到隨著動量比的增加，霧場的偏斜程度增大，偏斜角度相應的增大。當動量比由1.5增加到2.0時，霧場的偏斜角度由9°增加到了21°。

動量比的變化對噴霧擴散角的影響如圖11所示。隨著動量比的增加，噴霧擴散角也會增加，但增加幅度不明顯，因此動量比對噴霧擴散角的影響較小，這與文獻[4]試驗結論一致。

圖10 不同動量比下霧場偏斜角的變化

圖11 不同動量比下噴霧擴散角的變化

由第2.2節分析可知，孔徑比變化的影響實質上是雙股射流參與撞擊的動量比發生了變化。因此，當孔徑比保持不變的時候，隨著動量比的增加，未參與撞擊部分的動量一直在增加，導致了霧化撞擊不充分，霧化結果變差。因此，動量比的增加會導致液膜破碎長度增加，如圖12所示。圖13統計了不同動量比下的霧化液滴索特爾平均直徑。由圖可知，霧化液滴的索特爾平均直徑隨著射流動量比的增加而增加。

圖12 不同動量比下液膜破碎長度的變化

圖13 不同動量比下SMD的變化

3 結論

本文基于GPU加速的移動粒子半隱式法，對雙股射流霧化過程進行了直接數值模擬，成功模擬了液膜變成液絲，液絲破碎成液滴的一次霧化過程。分析模擬結果，得到結論如下：

(1)孔徑比通過改變參與撞擊的射流動量的變化來影響霧化過程。孔徑比不為1的時候，霧場會存在明顯的偏斜，同時會形成凹形液膜。

(2)動量比主要影響霧場的偏斜程度，且動量比越大，霧場偏斜程度越大。本文計算中，撞擊夾角為60°、孔徑比為1時，動量比由1.5提高到2.0，霧場偏斜角由9°提高到21°。

(3)孔徑比和動量比的增大都會造成噴霧擴散角的增大，但是變化幅度較小。

(4)孔徑比和動量比的增大都會造成霧化撞擊不充分，因此液膜會在更長距離內到達臨界破碎點，所以使得液膜破碎長度和霧化液滴的索特爾平均直徑增大。

后續將進一步改進模擬的精度，對多股射流霧化開展研究，分析各噴注單元之間的干涉關系，并對多股射流霧化進行數值的研究。