冼錦宇，陳錢，蔡樹峰

中山大學(xué)，廣東 深圳 518107

超燃沖壓發(fā)動(dòng)機(jī)作為高超聲速飛行器的關(guān)鍵技術(shù)之一[1-5]，具有很高的工程應(yīng)用價(jià)值。液態(tài)燃料是超燃沖壓發(fā)動(dòng)機(jī)目前最主要的燃料，與氫氣相比，液態(tài)燃料具有能量密度高、環(huán)境友好且儲(chǔ)存簡(jiǎn)便等特點(diǎn)。液態(tài)燃料在燃燒室內(nèi)的霧化過程是燃料燃燒的第一步，而液態(tài)燃料霧化質(zhì)量的好壞直接影響燃燒室內(nèi)的燃燒過程。為了了解液態(tài)燃料的霧化特性，較多研究者對(duì)液體橫向射流霧化特性進(jìn)行了研究。

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

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

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

1 數(shù)值模擬方法

1.1 氣相控制方程

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

連續(xù)方程

式中，源項(xiàng)S?m，S?ui，S?h，S?Ym分別表示質(zhì)量、動(dòng)量、能量、組分輸運(yùn)在氣液相間的相互作用。

1.2 液相控制方程

流場(chǎng)中液滴的運(yùn)動(dòng)根據(jù)拉格朗日液滴追蹤方法進(jìn)行模擬，液滴位置矢量、速度矢量、溫度和質(zhì)量的控制方程如下。

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

液滴的阻力系數(shù)公式定義為

1.3 幾何模型與物理模型

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

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

計(jì)算采用結(jié)構(gòu)化網(wǎng)格，壁面第一層厚度為1×10-5m，對(duì)射流噴口以及底壁面附近進(jìn)行網(wǎng)格加密，網(wǎng)格總數(shù)量約為105萬個(gè)。數(shù)值模擬的每個(gè)算例在96個(gè)CPU核心上運(yùn)行4天?？諝鈦砹骱退淞鞯娜肟诹鲃?dòng)參數(shù)與Lin 等[7]的試驗(yàn)一致，具體參數(shù)見表1。

表1 來流空氣和射流的入口流動(dòng)參數(shù)Table 1 Inlet flow parameters for incoming air and jets

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

在超聲速氣流中的液體橫向射流中，液滴存在變形、振蕩、線性失穩(wěn)過程，KH模型考慮氣液界面較強(qiáng)剪切作用導(dǎo)致的KH 不穩(wěn)定波增長過程，能夠較好地模擬在噴口附近的小液滴從初始大液滴上剝離的過程；而當(dāng)液滴的加速度較大時(shí)，RT 不穩(wěn)定性則主導(dǎo)液滴破碎過程?；谏鲜龇治?，本文采用Patterson[16]建立的KH/RT 混合破碎模型進(jìn)行霧化數(shù)值計(jì)算。

液滴阻力系數(shù)影響液滴破碎的計(jì)算，對(duì)霧化過程有重要影響。動(dòng)態(tài)阻力模型能夠在液滴形狀發(fā)生變化時(shí)確定液滴阻力系數(shù)，并且能夠適應(yīng)任何破碎模型，因此本文采用動(dòng)態(tài)阻力模型計(jì)算液滴阻力系數(shù)。

2 計(jì)算結(jié)果與討論

2.1 噴霧流動(dòng)基本特性

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

圖2 液霧場(chǎng)斜視圖Fig.2 Liquid spray field oblique view

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

本文計(jì)算結(jié)果為非定常的，且通過統(tǒng)計(jì)平均獲得了工程應(yīng)用中關(guān)心的平均結(jié)果。統(tǒng)計(jì)平均包括瞬時(shí)空間平均和時(shí)間平均兩個(gè)步驟[19]。首先對(duì)不同時(shí)刻網(wǎng)格控制體內(nèi)所有液滴進(jìn)行空間平均，得到瞬時(shí)空間平均結(jié)果

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

由于不同時(shí)刻的網(wǎng)格控制體得到的平均結(jié)果可能會(huì)發(fā)生變化，因此基于空間平均值在時(shí)間上進(jìn)行加權(quán)平均得到時(shí)間域上的平均值

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

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

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

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

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

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

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

圖7 是不同橫截面上水液滴的平均速度分布，其中藍(lán)線是液霧邊界?？梢钥吹?，由于受到外邊緣氣流相互作用的影響，液滴平均速度由中心向邊緣不斷增大。隨著水射流向下游發(fā)展，液滴群范圍增大，橫截面液滴群受到的氣流作用增強(qiáng)，液滴群整體速度也增大，液霧場(chǎng)邊緣的平均速度增長明顯。

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

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

圖8是不同橫截面上水液滴的SMD值分布，其中藍(lán)線是液霧邊界?？梢钥吹剑煌孛嬷g的SMD值變化規(guī)律較為相近；SMD值最高接近35μm，主要集中在液霧場(chǎng)上部分和中心截面上，其余區(qū)域SMD值在10～20μm之間；同樣受氣流相互作用的影響，SMD值與平均速度的變化趨勢(shì)相反，由中心向邊緣處不斷減小，這說明中間部分液滴的破碎不夠充足。

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

2.2 初始液滴直徑影響

本節(jié)將探索初始液滴直徑呈不同分布對(duì)液霧場(chǎng)的影響。圖9 是初始液滴直徑呈Rosin-Rammler 分布的算例的中心截面液滴瞬時(shí)穿透高度?？梢钥吹?，初始液滴直徑呈Rosin-Rammler 分布的算例與均勻分布的算例的穿透高度較為相近，但在x=140mm 附近Rosin-Rammler 分布的算例更接近試驗(yàn)值。

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

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

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

圖11 Rosin-Rammler分布算例的橫截面液滴體積分?jǐn)?shù)分布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 給出中心截面不同流向位置的液滴平均速度法向分布曲線，其中液滴平均速度通過來流氣體聲速進(jìn)行歸一化處理，y坐標(biāo)通過當(dāng)?shù)卮┩父叨冗M(jìn)行歸一化處理。可以看到，數(shù)值模擬結(jié)果整體呈現(xiàn)鏡像S 分布，與試驗(yàn)結(jié)果的變化趨勢(shì)相一致，y/h＜0.4時(shí)，數(shù)值模擬結(jié)果與試驗(yàn)值較為吻合；而當(dāng)y/h＞0.4 時(shí)，數(shù)值模擬的值隨y/h的增大而緩慢提高，而試驗(yàn)結(jié)果變化較為劇烈。如上述分析，不同初始液滴直徑分布對(duì)液滴速度的影響較小，因而值呈現(xiàn)相似變化。圖15 給出中心截面不同流向位置的SMD 法向分布曲線，y坐標(biāo)通過當(dāng)?shù)卮┩父叨冗M(jìn)行歸一化處理。從圖15 中可以看到，在y/h＜0.4 的近壁區(qū)域，計(jì)算結(jié)果與試驗(yàn)變化趨勢(shì)一致；而在y/h＞0.4的中部區(qū)域，計(jì)算結(jié)果與試驗(yàn)結(jié)果差異較大，這可能歸因于中部區(qū)域二次破碎不夠充分；而在液霧場(chǎng)上部，數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果相近。雖然計(jì)算的SMD 分布曲線與試驗(yàn)給出曲線仍有一定差距，但仍然能夠看到初始液滴直徑不同分布對(duì)SMD 值的影響，且Rosin-Rammler分布比均勻分布能獲得與試驗(yàn)更為接近的曲線變化趨勢(shì)。

圖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 結(jié)論

本文采用數(shù)值模擬方法對(duì)超聲速氣流中液體噴霧流動(dòng)進(jìn)行了數(shù)值模擬研究，獲得以下主要結(jié)論:

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

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

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

（4）不同液滴初始直徑分布對(duì)霧化特性的影響規(guī)律可為超燃沖壓發(fā)動(dòng)機(jī)液體燃料霧化質(zhì)量的評(píng)估和改進(jìn)提供重要參考。