超音速橫向流作用下射流的二次破碎模型研究

楊東超，朱衛(wèi)兵，陳宏，郭金鑫，劉建文

(1.哈爾濱工程大學(xué)航天與建筑工程學(xué)院，黑龍江 哈爾濱 150001;2.中國航天科工集團(tuán)公司三十一所 高超聲速沖壓發(fā)動機(jī)技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京 100074)

隨著吸氣式高超聲速飛行器的發(fā)展，超聲速燃燒沖壓發(fā)動機(jī)的研究得到了世界各國的高度重視。與低速燃燒流場相比，超燃流場更為復(fù)雜，其內(nèi)部是包含各種波系結(jié)構(gòu)的復(fù)雜三維流動，對于液態(tài)燃料還包括破碎、碰撞、蒸發(fā)等物理過程。目前，超音速流場中的破碎機(jī)理尚未完全掌握，沒有統(tǒng)一的模型對這一過程進(jìn)行模擬。當(dāng)前的主要處理方法是在一定實(shí)驗(yàn)基礎(chǔ)上，提出帶有經(jīng)驗(yàn)常數(shù)的破碎模型。已經(jīng)建立的破碎模型都是在低速條件下提出的，對于超音速流中的霧化仍未有可靠的模型，普遍的做法是直接將已有的低速模型用于超音速流中。這種處理方式仍未得到理論與實(shí)驗(yàn)的充分驗(yàn)證，許多工作需進(jìn)一步完善。

目前所采用的低速半經(jīng)驗(yàn)霧化模型，主要有TAB(Taylor analogy breakup)模型、Reitz波不穩(wěn)定性模型(也稱WAVE模型，或K-H模型)和混合霧化模型。岳連捷等［1］利用KIVA程序，采用TAB模型對超聲速氣流中橫向煤油射流霧化過程進(jìn)行了數(shù)值模擬。Im Kyoung-Su等［2］采用修正的 K-H/R-T混合模型對超音速來流下水垂直噴注進(jìn)行了計算，得到液滴及馬赫數(shù)的分布。楊順華等［3］同樣采用混合破碎模型對Lin K C［4］的水噴注實(shí)驗(yàn)進(jìn)行了模擬，但與文獻(xiàn)［2］中不同的是，模型中的部分參數(shù)由實(shí)驗(yàn)給出，其結(jié)果更加準(zhǔn)確，但涉及的經(jīng)驗(yàn)參數(shù)不易給出。

本文利用Fluent分別選用4種二次破碎模型，對超音速條件下垂直射流過程進(jìn)行了研究，分析了氣相與液相的相互作用，同時對不同初始滴徑對計算的影響進(jìn)行了初步分析。計算中忽略了低溫時液體的蒸發(fā)，并將模擬結(jié)果與實(shí)驗(yàn)進(jìn)行了對比。

1 物理模型及控制方程

1.1 連續(xù)相控制方程

氣相的控制方程為三維N-S方程，并考慮了與液滴的相互作用項(xiàng)，其矢量形式為

式中:U= ρ ρu ρv ρw ρE[]為流動變量，F(xiàn)、G、H分別為對應(yīng)方向上帶有粘性項(xiàng)的通量矢量，S為氣相與液滴項(xiàng)相互作用的源項(xiàng)。

1.2 Eulerian-Lagrangian方法

為簡化計算，對液滴項(xiàng)的追蹤采用Lagrangian方法，不考慮液滴所占的體積?？紤]到本文所引用實(shí)驗(yàn)的真實(shí)工況，忽略液滴的蒸發(fā)過程。在Lagrangian坐標(biāo)系下，每個單獨(dú)被追蹤的粒子均帶有不同標(biāo)記，并且代表一定數(shù)量，具有相同尺寸、位置和速度的液滴(即DPM模型，discrete phasemodel)。氣液兩相間的質(zhì)量、動量和能量的相互作用，以源相的形式分別加到氣液兩相的控制方程中。其液滴的控制方程為:

1)液滴在任意時刻的瞬時速度通過求解動量方程得到

式中:FD=18μCDRep/(24ρpd2p);up、ρp和 dp是液滴的速度、密度和直徑;Rep是液滴雷諾數(shù);CD是阻力系數(shù)［4］，F(xiàn)為其他作用力。

2)液滴的軌跡方程:

3)液滴的能量方程:

在液滴溫度小于其蒸發(fā)溫度時的溫升方程是:

式中:mp、cp和Ap分別為液滴質(zhì)量、定壓比熱和表面積;T∞是氣相溫度;h是對流傳熱系數(shù)。

1.3 二次破碎模型

下面分別介紹本文采用的4種二次破碎模型:TAB 破碎模型［5］，WAVE 破碎模型［6］，K-H/R-T 破碎模型［7］，SSD 破碎模型［8］。

1.3.1 TAB模型

TAB模型是計算液滴裂化的傳統(tǒng)方法，是Taylor根據(jù)液滴的變形和振蕩與彈性質(zhì)量系統(tǒng)的相似性提出的，液滴表面張力類比恢復(fù)力，液滴阻力類比外力，液滴粘性力類比阻尼力。液滴在空氣阻力的作用下發(fā)生變形，設(shè)y=x/(Cbr)，y為液滴最大直徑無量綱變形量，x是與相對速度方向垂直的液滴最大直徑變形量，Cb是液滴變形判斷常數(shù)，則球形液滴控制方程為

式中:ρ、u和u＇分別表示液滴所在位置氣相密度、平均速度和湍流脈動速度;ρd、r和ud分別表示液滴密度、半徑和空間速度;σd是液滴的表面張力;μd是液滴的粘性系數(shù)。

假定方程(5)的系數(shù)為常數(shù)，可得方程的解為

式中:We是液滴Weber數(shù)，td是液滴中粘性阻尼的時間尺度，ω是液滴振蕩頻率。認(rèn)為只有當(dāng)液滴變形量y＞1時液滴才會發(fā)生二次破碎，子液滴數(shù)量通過質(zhì)量恒求得，子液滴的法向速度等于父液滴在破碎時刻的法向振動速度。因此使用TAB模型，噴霧錐角可以計算得出，而不需要預(yù)先設(shè)定。

TAB模型適宜于低韋伯?dāng)?shù)射流，對于特別高的韋伯?dāng)?shù)，噴霧液滴散落分布，用彈簧質(zhì)量系統(tǒng)類比是不適合的。

1.3.2 WAVE模型

WAVE模型基于表面波不穩(wěn)定理論，認(rèn)為液滴破碎的原因是液滴表面的不穩(wěn)定增長。作用在液體表面初始擾動的增長與液體表面不穩(wěn)定波的波長、液相和氣相的物性和動力學(xué)參數(shù)有關(guān)。液滴破碎分為低速噴射和高速噴射2種模式。在低速噴射下發(fā)生的Rayleigh-Taylor不穩(wěn)定會產(chǎn)生比原始滴徑更大的新液滴;在高壓噴霧下Kelvin-Helmholtz不穩(wěn)定波的增長會產(chǎn)生半徑小得多的液滴。Reitz提出的WAVE模型［6］認(rèn)為液體射流破碎是由氣液兩相之間的相對速度造成的，射流表面K-H不穩(wěn)定波的增長引起了液滴從液體表面剪切下來。

液滴半徑的變化率和破碎形成的子液滴的尺寸與最不穩(wěn)定波的頻率Ω和波長Λ有關(guān):

父液滴半徑按照下式變化:

式中:τ為破碎時間;B1為模型常數(shù)，Reitz取 1.73［6］。

1.3.3 K-H/R-T模型

R-T不穩(wěn)定和K-H不穩(wěn)定在液滴破碎過程中是同時發(fā)生的。Patterson［7］結(jié)合這2種破碎模型，建立了K-H/R-T混合模型來同時處理這2種不穩(wěn)定因素。K-H/R-T模型認(rèn)為:在液滴破碎過程中K-H表面波和R-T擾動一直處于競爭關(guān)系。R-T機(jī)理適宜于描述由于液滴的快速減速而導(dǎo)致表面波在液滴的背風(fēng)面快速增長，引起變形導(dǎo)致破碎成小液滴，而K-H機(jī)理適宜于高相對速度和高環(huán)境密度的液滴破碎。用WAVE模型公式模擬K-H破碎，R-T擾動通過具有最大增長率的表面波的頻率ΩRT描述。

利用線性化的水動力學(xué)方程，可以得到液體表面具有最大增長速率的擾動波長ΛRT和擾動頻率ΩRT，具體表達(dá)式為

式中:ad是液滴在氣相中運(yùn)動的加速度。液滴加速度和阻力系數(shù)的關(guān)系為

式中:ur為液滴與氣相的相對速度。

破碎時間尺度τRT利用具有最大增長速率的擾動頻率計算得到

該時間尺度定義了R-T破碎發(fā)生的間隔。

液滴只有在方程(11)給出的擾動波長小于液滴直徑時才會發(fā)生破碎。當(dāng)破碎發(fā)生后，生成的新的子液滴的半徑為

R-T破碎只在液滴直徑大于液滴表面具有最大擾動波長的情況下才考慮。當(dāng)R-T擾動發(fā)展的累積時間超過了R-T破碎時間尺度時，液滴會根據(jù)RT破碎模型發(fā)生二次破碎，生成半徑更小的子液滴。

1.3.4 SSD模型

當(dāng)液滴噴入計算區(qū)域時，SSD模型會給定初始滴徑，并預(yù)測破碎過程中的時間、二次破碎液滴的數(shù)量及物性參數(shù)。當(dāng)液滴尺寸大于由SSD模型確定的臨界半徑時液滴將發(fā)生破碎，臨界半徑計算如下:

式中:Wecr是臨界韋伯?dāng)?shù)，韋伯?dāng)?shù)是需要具體指定的常數(shù)，默認(rèn)的韋伯?dāng)?shù)是6。破碎時間定義為

式中:B是一個用戶定義的常數(shù)，默認(rèn)取B=1.73。滴徑大于臨界半徑時，破碎時間會增加，當(dāng)液滴的破碎時間大于臨界破碎時間時，破碎就會發(fā)生。

當(dāng)液滴發(fā)生破碎時，父液滴發(fā)生破碎并且生成新液滴。這些子液滴的直徑是通過采用對數(shù)分布函數(shù)隨機(jī)獲得的，直徑范圍需預(yù)先設(shè)定，液滴破碎的計算過程將一直持續(xù)到父液滴完全破碎為止。

2 計算方法

2.1 計算模型及網(wǎng)格劃分

美國Taitech公司的Lin KC及其合作者對超音速橫向流作用下的液體噴射過程進(jìn)行了大量卓有成效的研究［4，9-13］。本文引用 Lin K C于2004年對超聲速流動中水射流結(jié)構(gòu)進(jìn)行的實(shí)驗(yàn)研究［4］。實(shí)驗(yàn)中直徑為0.5 mm的直流式噴注器安裝于矩形超聲速風(fēng)洞底部距入口下游139 mm處。水通過噴注器射入超聲速空氣來流中。為了節(jié)約計算工作量，沿來流方向只計算600 mm。圖1給出了計算模型示意圖和幾何尺寸。

圖1 實(shí)驗(yàn)?zāi)Ｐ虵ig.1 Experimental model

采用結(jié)構(gòu)化網(wǎng)格。對噴口附近進(jìn)行局部加密，第一層網(wǎng)格節(jié)點(diǎn)距離壁面1×10-4m，網(wǎng)格總數(shù)約為73萬，網(wǎng)格如圖2所示。

表1和表2分別給出了氣液兩相的具體條件，此條件對應(yīng)的是水射流和氣流動壓比為10的情形

圖2 計算網(wǎng)格Fig.2 Computational grid

表1 空氣來流條件Table 1 Air flow conditions

表2 水噴霧計算條件Tab le 2 Water spray conditions

2.2 計算條件

計算利用Fluent軟件，采用基于密度的耦合隱式求解器，通量采用AUSM格式，空間離散采用二階迎風(fēng)格式，湍流模型采用超燃流場計算常用的可壓縮SST k-ω湍流模型，并采用代數(shù)多重網(wǎng)格技術(shù)加速收斂，使用ILU方法進(jìn)行光順處理。

目前多數(shù)的CFD模擬中，都忽略一次破碎過程，根據(jù)經(jīng)驗(yàn)給定液滴的初始尺寸、溫度和速度。本文給定的液滴按Rosin-Rammler分布，平均初始直徑為0.1 mm，流量根據(jù)實(shí)驗(yàn)噴注條件折算給出，并利用隨機(jī)游走模型考察湍流脈動效應(yīng)對液滴的作用。

3 計算結(jié)果及分析

對上述實(shí)驗(yàn)工況分別采用不同二次破碎模型進(jìn)行了計算，并對計算結(jié)果進(jìn)行對比研究。

3.1 不同二次破碎模型對比

圖3分別給出了不考慮液滴破碎和采用4種二次破碎模型得到的三維液霧流場中的水滴分布，圖中液滴大小代表水滴的相對大小，并給出了中間截面上的壓力分布圖。從圖中可以看出，水噴出后，在超音速氣流的作用下向氣流方向偏折，并向下游運(yùn)動。滴徑大的液滴由于自身動量較大，隨體性差，貫穿能力強(qiáng)，主要分布在液霧外圍，而小液滴主要分布在近壁區(qū)，即液滴的分布位置一定程度上反映了液滴的直徑。對比不考慮破碎和考慮破碎模型的結(jié)果圖3(a)可以看出，考慮液滴破碎時，液滴的分布范圍更廣更均勻，且在湍流脈動的作用下，液滴的分布呈現(xiàn)一定的隨機(jī)性。從壓力分布圖中可以看出，由于射流的阻礙作用，在射流口上游出現(xiàn)一道較弱的激波。對比4種不同破碎模型計算得到的液滴分布結(jié)果，TAB模型和SSD模型得到的液滴分布范圍較小，而另2種模型的結(jié)果呈現(xiàn)一定相似性;與其他3種模型相比，TAB模型得到的液滴更貼近壁面區(qū)。

圖3 液滴分布模擬結(jié)果Fig.3 The droplet distribution of simulation

圖4分別給出了4種不同的破碎模型得到的液滴穿透深度與實(shí)驗(yàn)值的比較。選取了2種不同測試手段得到的實(shí)驗(yàn)結(jié)果，分別是PDPA技術(shù)得到的穿透深度擬合公式［4］

式中:h是水霧穿透深度，d0是噴注器直徑，q是液氣動壓比，x是從噴注器位置開始的流向距離。另一個是劉靜［14］利用紋影法得到的擬合公式:

公式所表示的曲線為噴霧區(qū)的外包絡(luò)線。從圖中可以看出:1)利用TAB模型得到的噴注高度遠(yuǎn)低于實(shí)驗(yàn)值，說明該模型計算的滴徑較真實(shí)工況小，不適用于高速來流作用下的噴注過程;2)對于WAVE模型的結(jié)果，距離噴口較遠(yuǎn)的下游區(qū)液滴高度低于實(shí)驗(yàn)值，說明大滴徑的液滴數(shù)目少，噴注的貫穿力較小;3)K-H/R-T模型和SSD模型的結(jié)果與實(shí)驗(yàn)值符合得較好，相比而言，K-H/R-T破碎模型計算的液滴分布范圍更廣，其滴徑范圍也較SSD模型的大。

圖4 液滴穿透深度模擬結(jié)果與實(shí)驗(yàn)對比Fig.4 Comparison of drop let penetration depth between simulation and experiment

圖5分別給出了流場出口視角的4種破碎模型及無破碎計算得到的液霧分布。從圖中可以看出，無破碎時，液滴噴注最高，且分布范圍集中。TAB模型和SSD模型得到的液霧在展向的分布范圍較小，而WAVE模型和K-H/R-T模型結(jié)果的展向分布范圍相似。K-H/R-T模型和WAVE模型的結(jié)果呈“圓形”，表明破碎產(chǎn)生多種滴徑的液滴。而SSD模型的結(jié)果呈“樹形”，即液霧沿展向的運(yùn)動主要集中在距離壁面較遠(yuǎn)的區(qū)域內(nèi)，且分布相對集中。

圖5 液滴分布模擬結(jié)果(流場出口視角)Fig.5 Drop let distribution of simulation(outlet perspective)

圖6給出了4種破碎模型計算得到沿流向的液滴SMD(Sautermean diameter)分布。從圖中可以看出，液滴在離開噴口后便發(fā)生了破碎，從初始的100μm破碎成10～30μm的液滴。TAB模型得到的破碎后小液滴的平均尺度最小，沿流向幾乎不變，而其他3種模型得到的破碎后液滴的平均尺度接近，且隨著向下游流動而呈下降趨勢。實(shí)驗(yàn)給出的全流場內(nèi)統(tǒng)計的SMD為15.2μm，4種模型的結(jié)果均有一定偏差。

圖6 液霧沿流向的SMD分布Fig.6 SMD distribution of droplets along the direction of flow

3.2 初始滴徑的影響

液滴的初始直徑為一次破碎后的滴徑，由于計算中忽略了一次破碎過程，初始滴徑根據(jù)經(jīng)驗(yàn)給出，本文初步取為0.1 mm，下面給出了平均初始滴徑分別為0.2 mm和0.5 mm，基于K-H/R-T模型的部分計算結(jié)果。

圖7和圖8分別給出了軸側(cè)視角和出口視角的滴徑分布圖。結(jié)合平均初始滴徑為0.1 mm的結(jié)果(圖3(d)和圖5(d))，可以看出，初始滴徑對液滴的分布有較大影響:隨著初始滴徑的增大，大液滴比例和直徑均增加，且液滴分布范圍增大，其噴注高度較0.1 mm初始液滴的結(jié)果也相應(yīng)增加。

圖7 不同初始滴徑的液滴分布(軸側(cè))Fig.7 Droplets distribution for different initial droplet diameter

圖8 不同初始滴徑的液滴分布(出口視角)Fig.8 Droplet distribution for different initial droplet diameter

圖9給出了3種不同初始滴徑條件下，液滴沿流向SMD分布。從圖中可以看出，SMD分布的3個結(jié)果相近，說明該參數(shù)對初始滴徑的改變并不敏感，初始滴徑的變化只改變了破碎后的滴徑分布規(guī)律，而不是統(tǒng)計平均值。

圖9 不同初始滴徑的液滴沿流向SMD分布Fig.9 SMD distribution of droplets along the direction of flow

4 結(jié)論

本文采用4種不同二次破碎模型，對超聲速氣流中水射流進(jìn)行了模擬。模擬得到了不同模型下的水射流破碎后噴霧結(jié)構(gòu)，穿透深度及液滴平均尺寸沿流向分布，并將計算結(jié)果與實(shí)驗(yàn)擬合的穿透深度進(jìn)行了對比。結(jié)果表明:

1)TAB模型得到的滴徑最小，其穿透深度也最小，不適用于高Weber數(shù)下的超音速流場條件;

2)SSD模型計算的滴徑尺寸較為均勻，該模型涉及臨界破碎直徑這一經(jīng)驗(yàn)參數(shù)，該值的選取對超音速條件下霧化的影響還需進(jìn)一步驗(yàn)證;

3)WAVE模型與K-H/R-T模型的結(jié)果相近，而K-H/R-T模型得到的穿透深度與實(shí)驗(yàn)更為相符;

4)對于K-H/R-T模型，隨著初始滴徑的增加，較大液滴所占比例和直徑均增加，而流向SMD分布幾乎不變。

此項(xiàng)工作為今后深入研究超燃沖壓發(fā)動機(jī)中的燃料霧化特性奠定了基礎(chǔ)。下一步的工作是對超燃沖壓發(fā)動機(jī)燃燒室的燃料霧化、燃燒全過程進(jìn)行數(shù)值模擬。

［1］岳連捷，俞剛.超聲速氣流中橫向煤油射流的數(shù)值模擬［J］.推進(jìn)技術(shù)，2004，25(1):11-14.YUE Lianjie，YU Gang.Numerical simulation of kerosene spray in supersonic cross flow［J］.Journal of Propulsion Technology，2004，25(1):11-14.

［2］IM K S，LIN K C，LAIM C.Spray atomization of liquid jet in supersonic cross flows［J］.AIAA Paper，2005:732.

［3］楊順華，樂嘉陵.超聲速氣流中液體燃料霧化數(shù)值模擬［J］. 推進(jìn)技術(shù)，2008，29(5):519-522.YANG Shunhua，LE Jialing.Numerical simulation of liquid fuel atomization in supersonic cross flow［J］.Journal of Propulsion Technology，2008，29(5):519-522.

［4］LIN K C，KENNEDY P J，JACKSON TA.Structures of water jets in a Mach 1.94 supersonic crossflow［J］.AIAA Paper，2004:971.

［5］O’ROURKE P J，AMEDEN A A.The TAB method for numerical calculation of spray droplet breakup.No.LA-UR-87-2105-Rev［R］.(s.l.):Los Alamos National Lab，1987.

［6］LIU A B，MATHER D，REITZR D.Modeling the effect of drop drag and breakup on fuel sprays.No.TP-930072［R］.(s.l.):Madison Engine Research Center，1993.

［7］PATTERSON M A.Modeling the effects of fuel injection characteristics on diesel combustion and emissions［D］.Madison:University of Wisconsin-Madison，1997.

［8］REITZ R D.Modeling atomization processes in high pressure vaporizing sprays［J］.Atomization and Spray Technology，1987，3:309-337.

［9］LIN K C，KENNEDY P J，JACKSON T A.Spray penetration heights of angled-injected aerated-liquid jets in supersonic cross flows［J］.AIAA Paper，2000:0194.

［10］LIN K C，KENNEDY P J，JACKSON T A.Spray structures of aerated-liquid jets in subsonic cross flows［J］.AIAA Paper，2001:0330.

［11］LIN K C，RAJNICEK C，MCCAL J，et al.Structures of aerated-liquid jets injected from various nozzle contours［J］.AIAA paper，2011:232.

［12］OLINGER D S，LEE J，LIN K C，et al.Digital holographic analysis of near-field aerated liquid jets in cross flow partI:algorithm development［J］.AIAA Paper，2013:0164.

［13］OLINGER D S，LEE J，LIN K C，et al.Digital holographic analysis of near-field aerated，liquid jets in cross flow.part II:measurements［J］.AIAA Paper，2013:0167.

［14］劉靜，王遼，張佳，等.超聲速氣流中橫向射流霧化實(shí)驗(yàn)和數(shù)值模擬［J］.航空動力學(xué)報，2008，23(4):724-729.LIU Jing，WANG Liao，ZHANG Jia，et al.Experimental and numerical simulation of atomization of liquid jet in supersonic cross flow［J］.Journal of Aerospace Power，2008，23(4):724-729.