汪鳳山，毛曉芳，虞育松，張 榛，王 平

(1.北京控制工程研究所，北京100190;2.北京交通大學，北京100044)

雙組元離心式噴注器霧化性能的大渦模擬數(shù)值研究

汪鳳山1，毛曉芳1，虞育松2，張 榛1，王 平1

(1.北京控制工程研究所，北京100190;2.北京交通大學，北京100044)

基于氣液兩相體積混合分數(shù)建立某型雙組元離心式噴注器內(nèi)部流場及霧化場的數(shù)學模型，并采用大渦模擬(LES)方法對其霧化性能進行研究.計算結(jié)果表明，噴注器外路噴射速度低于內(nèi)路噴射速度，內(nèi)外兩路同時噴射時，兩路相互影響，霧化錐角增加;隨著噴注器出口直段長度的增加，內(nèi)外路噴射速度和流量系數(shù)均呈單調(diào)下降趨勢;隨著噴注壓降的增加，噴注器流量系數(shù)和出口噴射速度隨之增加，而霧化粒度不斷降低，霧化錐角不斷增大.

同軸離心式噴注器;混合分數(shù)模型;大渦模擬;數(shù)值模擬;霧化粒度;霧化錐角

隨著航天科技的發(fā)展，新一代大型、長壽命衛(wèi)星對雙組元姿控發(fā)動機的穩(wěn)態(tài)及脈沖工作性能要求提出了更高的要求.液體火箭發(fā)動機工作過程涉及到流動、霧化、燃燒及傳熱等多個過程，其中推進劑霧化是直接影響推力器比沖、可靠性與壽命的關(guān)鍵過程.為提高發(fā)動機的工作性能，通常需要對推進劑霧化單元(即噴注器)進行良好設(shè)計，尤其要對推進劑的霧化性能進行詳細的研究.

液體霧化現(xiàn)象包括液膜形成、液膜破碎和霧滴形成，是一種復雜的物理過程，工質(zhì)物性、噴注器結(jié)構(gòu)設(shè)計、出口背壓條件及入口壓力等均對其霧化性能產(chǎn)生影響.一直以來，國內(nèi)外學者試圖通過采用流體力學原理建立對霧化過程的數(shù)學描述，并建立多種數(shù)學模型[1].目前，針對液體在入口壓力作用下的霧化現(xiàn)象，通常采用歐拉近似方法與歐拉—拉格朗日迭代方法(即ALE方法)兩種方數(shù)值方法進行仿真分析.第一種方法完全忽略噴霧粒子的初速和在流場中的受力，液滴分布完全由氣流決定，通過平均計算，以擬流體的方法計算噴霧過程中的氣流和噴霧粒子的運動規(guī)律.該方法的優(yōu)點是計算工作量較小，速度快，能夠在較短的時間內(nèi)獲得噴霧粒子的運動規(guī)律，特別適用于計算粒徑較小、初速較低的噴霧輸運現(xiàn)象，如飽和蒸氣運動、煤粉輸運過程等.第二種方法中噴霧粒子的運動不僅與初速有關(guān)，而且還受氣動力的影響，運動規(guī)律通過拉格朗日方法進行描述，更符合噴霧過程的實際情況，能夠捕捉到霧化過程中的一些細節(jié)現(xiàn)象，精度也較高，因此在實際中也應用較廣.

針對雙組元液體火箭發(fā)動機，噴注器出口推進劑霧滴速度通常較高，且在氣流中的受力情況復雜，因此，不宜采用歐拉近似方法對其霧化過程進行仿真計算，通常采用歐拉-拉格朗日迭代方法.但由于噴霧粒子尺寸較小，受氣流場中湍流運動關(guān)系密切，而通常采用的雷諾應力湍流模型空間分辨率較低，很難得到噴霧場的具體細節(jié).如20世紀90年代，Jeng等[2-4]運用歐拉-拉格朗日迭代方法追蹤壓力旋流噴嘴的氣液分界面，同時運用數(shù)值模擬和實驗相結(jié)合的方法研究了結(jié)構(gòu)參數(shù)對液膜厚度和噴霧錐角的影響.周立新等[5]采用基于單流體模型及氣液兩相混合分數(shù)概念對離心式噴嘴內(nèi)部流場進行數(shù)值計算，并與Jeng的ALE方法進行對比.但由于采用了雷諾應力湍流模型，未能給出噴霧場的具體分布規(guī)律.

大渦模擬(LES)方法是通過將耗散尺度的湍流脈動過濾掉，只求解大尺度湍流脈動的一種新型納維-斯托克斯(N-S)方程數(shù)值求解方法.與傳統(tǒng)的雷諾時均模型(RANS)方法相比，具有分辨率高、普適性高、更能真實反映流動微觀結(jié)構(gòu)等優(yōu)點，但計算量要比雷諾時均模型方法的要大.該方法能夠得到較為精細和準確的流場細節(jié)，為通過采用數(shù)值計算的方法獲得噴霧場的細節(jié)提供了可能.隨著近年來計算機運行速度和計算能力的迅速進步，采用大渦模擬方法已經(jīng)可以處理部分工程問題的湍流流動，在國內(nèi)外噴霧學界已經(jīng)得到了越來越廣泛的應用.

本文通過采用大渦模擬數(shù)值方法對某型雙組元離心式噴注器的噴霧流場進行了數(shù)值仿真分析，得到了不同條件下噴注器的霧化性能參數(shù)和液滴分布特性，探索噴注器霧化特性數(shù)值求解的可行性與方法，所得結(jié)果可為未來雙組元噴注器的設(shè)計提供參考.

1 數(shù)值仿真模型與算法

1.1 氣液兩相流數(shù)學模型

針對噴注器出口噴霧區(qū)域內(nèi)同時存在氣液兩種不相容工質(zhì)，氣液兩相之間存在界面的情況.為對氣液兩相共存情況進行描述，通過引入體積分數(shù)的概念，可實現(xiàn)采用一組納維-斯托克斯(N-S)方程描述氣液兩相及相界面，這種方法被稱之為體積分數(shù)法(VOF)，單元體積分數(shù)εk的計算公式為

通過引入體積分數(shù)，對于液體噴霧特性可用如下基于壓力、速度和體積分數(shù)的N-S方程描述:△

對于牛頓流體，在斯托克斯假設(shè)下:

由公式可得:

式中，V為流體速度矢量，ρ為流體密度，p為流體各向同性壓力，F(xiàn)為體積力(表面張力)，g為重力，μ為流體動力黏度，τ為剪應力，I=δij為克羅內(nèi)克符號，S為速度梯度張量.對于牛頓流體，粘性系數(shù)為常數(shù)，則:

考慮到:

將方程(6)代入方程(5)，則有:△

相界面的跟蹤是不相溶多相流數(shù)值模擬的關(guān)鍵問題之一，通過基于異相流體混合分數(shù)的假設(shè)，VOF方法可以建立每一相流體體積分數(shù)的輸運方程，對相界面進行跟蹤:

式中，Sεk是描述相間質(zhì)量轉(zhuǎn)移的物理量，由于霧化過程模擬不考慮相變，因此相間不存在質(zhì)量轉(zhuǎn)移，即Sεk=0.對于一定體積分數(shù)下的密度和黏度滿足:

式中，下標l和g分別代表液相和氣相.

氣、液相界面結(jié)構(gòu)在流場作用下變得非常復雜，精確求解三維空間的界面曲率和表面張力相當困難.Brackbill等[6]利用較為簡潔的連續(xù)表面張力模型(CSF)解決了表面曲率和表面張力求解的問題，具體的表達式為:

1.2 方程的數(shù)值解法

N-S方程的基本特點是非線性和耦合性，在一般情況下不能用解析方法得到封閉形式的解，因此需要采用一定的數(shù)值解法對N-S方法進行求解.在本文中，通過采用LES方法來完成對控制方程的數(shù)值求解.

LES方法求解N-S方程的思想是對控制方程進行濾波處理，將流場中的大小湍流渦團進行分離.大尺度湍流渦由控制方程直接求解;小尺度湍流渦通過引入亞網(wǎng)格尺度(SGS)湍流模型進行模擬求解.由于小尺度湍流渦結(jié)構(gòu)比大尺度湍流渦具有更好的各向同性，因此LES對SGS模型的要求不像雷諾應力模型(RANS)那樣高.當大渦模擬的網(wǎng)格間距逐漸減小時，SGS黏度趨于零，預測結(jié)果逐漸逼近直接數(shù)值模擬(DNS)結(jié)果.

在采用LES數(shù)值方法對N-S方程求解過程中，為了使方程封閉，需要引入模型對該亞網(wǎng)格尺度應力項進行計算.目前，使用最廣泛且最簡單的亞網(wǎng)格尺度應力模型是Smagorinsky渦黏性模型[7]，該模型通過建立未知亞網(wǎng)格尺度應力與已知的變形速率張量之間的簡單函數(shù)關(guān)系來實現(xiàn)對亞網(wǎng)格尺度應力項的模擬.假設(shè)亞網(wǎng)格應力張量正比于平均應力張量，則:

式中，δij為克羅內(nèi)克符號，是已知的應變張量，υt為湍流渦黏度(m2/s)，為濾波尺度(定義為(△χ△y△z)1/3).LES的湍流渦黏度υt∝uλ，其中，u

式中，k是表面曲率，σ是表面張力系數(shù)，n是垂直表面的法向單位向量，εl為液相體積分數(shù)，其中:是流動速度，λ為最大的耗散湍流渦泰勒尺度(Taylor尺度).

通過對上述微分方程組離散化，將連續(xù)變化形式描述的微分方程轉(zhuǎn)化為離散形式的代數(shù)方程，采用有限體積法在計算機上進行數(shù)值求解.

2 仿真分析

2.1 噴注器的幾何結(jié)構(gòu)

本文研究的噴注器采用雙旋渦離心噴嘴結(jié)構(gòu)，該結(jié)構(gòu)具有霧化效果好，噴霧場分布均勻，以及便于在燃燒室壁面組織液膜冷卻等優(yōu)點，如圖1所示.為了盡量減小數(shù)值計算的網(wǎng)格數(shù)量，取網(wǎng)格平均間距為100μm，總網(wǎng)格單元數(shù)約5×106，噴嘴外區(qū)域的張角根據(jù)試算預先確定.

圖1 雙旋渦離心噴注器工作原理圖Fig.1 Principle of coaxial centrifugal injector

2.2 邊界條件

在噴注器入口處根據(jù)給定壓力入口邊界條件，湍動能和耗散率根據(jù)經(jīng)驗公式給出，出口設(shè)置壓力出口邊界條件，壁面設(shè)置無滑移固體邊界條件，近壁函數(shù)采用標準壁面函數(shù)法求解.

3 仿真結(jié)果與分析

圖2為噴注器在不同時刻的噴霧形態(tài)數(shù)值仿真結(jié)果，從圖中可以看出，液體離開噴嘴后，液膜迅速分裂為數(shù)根細長的液絲，隨后連續(xù)液絲會分裂為較長的離散液滴，從而完成整個霧化過程.對于外路噴嘴，由于出口液滴噴射速度相對較低(約4m/s)，其霧化過程相對較慢;而對于內(nèi)路噴嘴，由于出口液滴噴射速度較高(約9m/s)，其霧化過程進度相對較快，但噴霧形態(tài)具有明顯的周向不對稱性，噴霧邊緣區(qū)域存在較長的連續(xù)液絲;對于兩路同時噴射的情況，內(nèi)外噴嘴出口液滴相互之間發(fā)生強烈的相互作用，霧化速度介于內(nèi)外路獨立噴射之間，且基本消除了內(nèi)路單獨噴射時的周向不對稱現(xiàn)象.

圖2 噴注器不同時刻的噴霧形態(tài)Fig.2 Shape of the propellant spray varied with time

圖3 為噴注器在不同時刻噴霧錐角變化示意圖.從圖中可以看出，對于噴注器，大概經(jīng)過4ms左右，其噴霧錐已基本成形.由于外路速度較低，其霧化錐角與實際數(shù)值(110°)相差較大，而內(nèi)路霧錐角度卻略小于實際數(shù)值(90°);內(nèi)外兩路同時噴射時，由于同向旋流疊加的加強作用，其霧化角度比內(nèi)外路單獨噴射的時候均略有增加.

圖3 噴注器霧化錐角隨時間的變化Fig.3 Spray angle varied with time

圖4 為在0.5MPa壓差下，不同出口直段長度下噴注器噴射速度與流量系數(shù)的變化示意圖.從圖中可以看出，隨著出口直段長度的增加，內(nèi)外路噴射速度和流量系數(shù)均呈單調(diào)下降趨勢.可見隨著出口直段長度的增加，在摩擦力的作用下，噴注器總流阻不斷增加，從而導致噴嘴旋流速度不斷降低，進而影響到噴注器噴霧特性.

圖4 內(nèi)外路噴射速度和流量系數(shù)隨出口直段長度的變化Fig.4 Spray velocity ejected from the flow passage varied with the injector exit length

圖5 為不同噴注壓差下內(nèi)外路噴射速度和流量系數(shù)的變化.從圖中可以看出，隨著噴注壓差的增大，噴注器噴射速度和流量系數(shù)也隨之增大.

圖6為外路噴嘴霧化粒度和霧化錐角隨噴注壓降的變化.從圖中可以看出，隨著噴注壓降的增加，外路霧化粒度不斷降低，霧化錐角不斷增大.可見，隨著噴注壓降的提高，液體出口速度不斷升高，液膜與氣流之間的剪切作用更為強烈，液膜的不穩(wěn)定性增強，霧化粒度不斷降低;另外，由于噴注壓降的升高，增強了出口液膜的旋流速度，從而導致形成的霧錐角增大.

圖5 內(nèi)外路噴射速度和流量系數(shù)隨噴注壓差變化計算結(jié)果Fig.5 Spray velocity from the flow passage and flow coefficient varied with the injector pressure drop

圖6 外路霧化錐角和霧化粒度隨噴注壓差的變化Fig.6 Spray velocity and SMD ejected from the outer flow passage varied with the injector pressure drop

4 結(jié) 論

本文通過引入氣液兩相流體積分數(shù)近似方法描述氣液兩相及相界面特性，并采用LES方法對某型雙組元離心噴注器的霧化性能進行了數(shù)值研究，結(jié)果如下:

1)噴注器外路噴射速度低于內(nèi)路噴射速度，從而導致外路達到霧化穩(wěn)定的時間要比內(nèi)路更長;

2)內(nèi)外兩路同時噴射時，兩路相互影響，由于旋轉(zhuǎn)方向相同的緣故，噴注器霧化錐角較兩路均有所增加;

3)隨著噴注器出口直段長度的增加，內(nèi)外路噴射速度和流量系數(shù)均呈單調(diào)下降趨勢;

4)隨著噴注壓降的增加，噴注器流量系數(shù)和出口噴射速度隨之增加，而霧化粒度不斷降低，霧化錐角不斷增大.

對計算結(jié)果進行分析，可見通過引入體積分數(shù)(VOF)對氣、液兩相霧化流動納維-斯托克斯(N-S)方程進行簡化，并采用大渦模擬(LES)的方法對其進行數(shù)值求解，能夠較為準確地捕捉到霧化過程中的細節(jié)，計算精度也較高，是一種值得推廣應用噴注器霧化過程數(shù)值仿真方法.如能通過與實驗測試數(shù)據(jù)進行充分比較和分析，近一步優(yōu)化計算模型，將能為液體火箭發(fā)動機噴注器的設(shè)計、優(yōu)化提供更準確的參考依據(jù).

[1] Lefebvre A H.Atomization and sprays[D].Hemisphere Publishing.New York.1989

[2] Jeng SM，Joe M A，et al.Computational and experimental study of liquid sheet emanating from simple fuel nozzle[R].AIAA-97-0796

[3] Holtzclaw D，Sakman T，Jeng SM，et al.Investigation of flow in a simplex fuel nozzle[R].AIAA-97-2970

[4] Rizk N K，Chin J S.Comprehensive fuel nozzle model [R].AIAA-94-3278

[5] 周立新，雷凡培.離心式噴嘴內(nèi)流場特性的數(shù)值模擬[J].推進技術(shù)，2002，23(6):480-484 Zhou L X，Lei F P.Numerical simulation of internal flow field of swirl nozzle[J].Journal of Propulsion Technology，2002，23(6):480-484

[6] Brackbill J U，Kothe D B，Zemach C.A continuum method for modeling surface tension[J].Journal of Computational Physics，1992，100:335-354

[7] Smagorinsky J.General circulation experiments with the primitive equations.I.the basic experiment[J]. Monthly Weather Review，1963，91(3):99-164

Num erical LES Study of Sp ray Perform ance in a Bi-Propellant Coaxial Centrifugal In jector

WANG Fengshan1，MAO Xiaofang1，YU Yusong2，ZHANG Zhen1，WANG Ping1
(1.Beijing Institute of Control Engineering，Beijing 100190，China; (2.Beijing Jiaotong Uniυersity，Beijing 100044，China)

A numerical study of spray performance in a bi-propellant coaxial centrifugal injector is proposed by using the large eddy simulation(LES)method，based on themultiphase volume fraction model of the gas and liquid.The results show that the spray velocity ejected from the inner flow passage is larger than the outer flow passage，and the spray angle increases slightly when both of the two flow passages simultaneously inject.The spray velocity and flow coefficient increase with the decrease of the injector exit length and increase of the injector pressure drop.The spray droplet size(SMD)decreases and the spray angle increases with the increases of the injector pressure drop.

coaxial centrifugal injector;multiphase fraction model;large eddy simulation(LES);numerical study;spray droplet size;spray angle

V43

A

1674-1579(2012)06-0013-05

汪鳳山(1981—)，男，高級工程師，研究方向為雙組元推力器技術(shù);毛曉芳(1971—)，女，研究員，研究方向為航天器推進技術(shù);虞育松(1979—)，男，副教授，研究方向為燃燒科學與技術(shù);張 榛(1983—)，男，工程師，研究方向為雙組元推力器技術(shù);王 平(1981—)，男，高級工程師，研究方向為雙組元推力器技術(shù).

2011-08-10

DO I:10.3969/j.issn.1674-1579.2012.06.003