液液針?biāo)ǘ鄧娮卧獓婌F場數(shù)值模擬

王 凱,唐 亮,雷凡培,楊岸龍,周立新

(1.西安航天動力研究所 液體火箭發(fā)動機(jī)技術(shù)重點(diǎn)實(shí)驗(yàn)室,陜西 西安 710100;2.中國船舶工業(yè)集團(tuán)有限公司,北京 100044)

0 引言

自1970年首次被提出以來,針?biāo)ㄊ絿娮⒓夹g(shù)就成為變推力技術(shù)研究的熱點(diǎn),并被成功應(yīng)用于月球軟著陸探測和火星軟著陸探測。SpaceX的Tom Mueller首次將針?biāo)ㄊ絿娮⑵骷夹g(shù)應(yīng)用于獵鷹9號火箭的Merlin系列發(fā)動機(jī)上。不到十年間發(fā)展了3個改進(jìn)版本,使得獵鷹9成為商業(yè)航天市場上最具競爭力的火箭[1]。正是由于梅林發(fā)動機(jī)噴注器結(jié)構(gòu)簡單,零組件可快速更換實(shí)現(xiàn)改進(jìn),固有燃燒穩(wěn)定性好,變推力實(shí)現(xiàn)容易,為獵鷹火箭的可重復(fù)使用提供了堅(jiān)定的動力技術(shù)支持。

以往的噴注器推力變比最大為2.5∶1,變比再增大時由于低工況噴注壓降降低,會出現(xiàn)低頻燃燒不穩(wěn)定和燃燒效率顯著降低的情況。與以往的直流撞擊式及離心式噴注器不同,針?biāo)ㄊ絿娮⑵魍ㄟ^調(diào)節(jié)裝置改變氧化劑和燃料的噴注面積,保持各個工況的噴注壓降基本不變,進(jìn)而保持高的燃燒效率和良好的燃燒穩(wěn)定性,因而該噴注器技術(shù)成為深度變推力發(fā)動機(jī)的關(guān)鍵技術(shù)。針?biāo)ㄊ絿娮⑵骶哂歇?dú)特的幾何特性和噴注特性,它通過伸入推力室內(nèi)部的針?biāo)▏娮⑵鹘Y(jié)構(gòu),使軸向的液膜與徑向的液束(或液膜)垂直撞擊形成噴霧扇霧化混合,其工作原理結(jié)構(gòu)如圖1所示[2]。

圖1 針?biāo)ㄊ絿娮⑵髟韴DFig.1 Schematic diagram of pintle injector

目前針?biāo)ㄊ絿娮⑵麟m已成功應(yīng)用于工程,然而相關(guān)的研究多見于工程研制,基礎(chǔ)理論和數(shù)值研究很少。近年來,文獻(xiàn)[3-6]針對軸對稱針?biāo)▏娮⑵鏖_展了很低速度或很低動量比下的噴霧水試試驗(yàn)研究,然而僅能獲得描述噴霧場宏觀結(jié)構(gòu)的霧化角信息,難以獲得噴霧場內(nèi)部的詳細(xì)噴霧場結(jié)構(gòu),霧化角不能夠反映噴霧場內(nèi)部的詳細(xì)信息。這是因?yàn)獒標(biāo)ㄊ絿娮⑵鞯膰婌F場非常獨(dú)特,液霧極為稠密,且內(nèi)部呈現(xiàn)復(fù)雜的分區(qū)結(jié)構(gòu)特征。這給用現(xiàn)有光學(xué)測量設(shè)備進(jìn)行試驗(yàn)觀測帶來了極大挑戰(zhàn),能獲得的霧場特性參數(shù)非常有限。Sakaki等率先提出了平面針?biāo)▏娮⑵鲉卧臉?gòu)想,并通過試驗(yàn)清楚觀測到了噴霧場和燃燒火焰結(jié)構(gòu)[7]。后續(xù)Sakaki等又開展了同參數(shù)的軸對稱噴注器對比研究,結(jié)果表明兩者的特性一致,這為使用平面針?biāo)ㄊ絿娮⑵鞔孑S對稱針?biāo)ㄊ絿娮⑵鏖_展更深層次的研究提供了強(qiáng)有力的依據(jù),表明平面針?biāo)ǖ难芯糠椒ǚ浅Ｓ行8]。Kim等借助平面針?biāo)▏娮卧乃悸?針對徑向圓孔的平面針?biāo)ǘ鄧娮卧_展了自燃推進(jìn)劑的點(diǎn)火特性研究,結(jié)果表明阻塞率0.65的結(jié)構(gòu)在中等動量比下的點(diǎn)火延遲時間最短[9]。成鵬借助同樣的思路主要針對徑向圓孔的單噴注單元開展了較多的霧化特性試驗(yàn)研究及氣氧酒精的平面針?biāo)ㄈ紵^程光學(xué)觀測研究[10]。王凱等也借助平面針?biāo)▏娮卧乃悸烽_展了單噴注單元霧化過程高精度的數(shù)值模擬,首次分析了單噴注單元膜束撞擊變形特征、流場渦結(jié)構(gòu)特征及霧化混合分區(qū)結(jié)構(gòu)特征等,為揭示針?biāo)ㄊ絿娮⑵鲊婌F燃燒精細(xì)過程奠定了基礎(chǔ)[11]。

然而以上研究未關(guān)注實(shí)際針?biāo)ㄊ絿娮⑵髦邢噜弴娮卧g的相互影響。實(shí)際針?biāo)ㄊ絿娮⑵飨噜弴娮卧g的間距較小,存在較強(qiáng)的相互影響,與單噴注單元形成的噴霧場結(jié)構(gòu)不同,這使得單噴注單元得到的結(jié)果不能直接應(yīng)用于實(shí)際針?biāo)ㄊ絿娮⑵?。因?需要對多噴注單元與單噴注單元噴霧場之間的差異性進(jìn)行深入的研究。

綜上所述,本文基于AMR技術(shù)和分相識別的PLIC VOF新方法,以液液平面針?biāo)ǘ鄧娮卧獮檠芯繉ο?通過對徑向和軸向兩路液體分別追蹤,開展液液針?biāo)ǘ鄧娮卧F化過程的高精度數(shù)值模擬。首先通過數(shù)值仿真結(jié)果獲得了多噴注單元噴霧場特有的結(jié)構(gòu)特征;然后深入分析了相鄰噴注單元噴霧扇的撞擊變過程,揭示了單元間的相互影響機(jī)制;最后對比分析了多噴注單元與單噴注單元噴霧場結(jié)構(gòu)及特性的差異。這對于將單噴注單元研究得到的結(jié)果應(yīng)用于實(shí)際針?biāo)ㄊ絿娮⑵骶哂兄匾囊饬x。

1 數(shù)值仿真方法

1.1 控制方程

樹狀結(jié)構(gòu)的自然分級特性非常適合多尺度網(wǎng)格的應(yīng)用,而且四叉樹/八叉樹網(wǎng)格的靈活性與簡單性在處理界面流動時有著巨大的優(yōu)勢?；谶@種考慮,Popinet提供了一種處理不可壓、變密度、帶有表面張力的N-S方程的數(shù)值方法,結(jié)合四叉樹/八叉樹離散、投影法與多級泊松求解器[12]。通過使用面向?qū)ο蟮腃語言編程實(shí)現(xiàn)后即為Gerris軟件。Gerris所采用的數(shù)值方法可表示為[13]

(1)

本文主要為了研究撞擊變形及流動過程,采用分三相VOF方法對兩種液體推進(jìn)劑和環(huán)境氣體的相界面分別進(jìn)行追蹤。為了分別識別兩路液體的變形運(yùn)動過程,認(rèn)為兩種液體互不相溶,且兩者之間存在液相界面。于是分別定義液相1的體積分?jǐn)?shù)c1(x,t)和液相2的體積分?jǐn)?shù)c2(x,t),對應(yīng)得到的流體密度和動力黏性系數(shù)為

工作倦怠發(fā)生的原因有很多種，然而我認(rèn)為員工的分配公平感是現(xiàn)階段最值得關(guān)注的問題之一。分配公平感，指的是職工自身對于資源的分配是否合理的個人性判斷，是以個人的公平觀為基礎(chǔ)的。不公平感是比較出來的，所以其具有相對性，但并無絕對標(biāo)準(zhǔn)。公平感是具有主觀性的，完全因個人的價值取向，這樣就可能造成因人而異的后果。不公平感是具有擴(kuò)散性的，職工在有不公平感時，往往會牽扯到整個營銷團(tuán)隊(duì)，負(fù)面情緒波及，這樣就會對公司的業(yè)績造成嚴(yán)重的影響，從而影響公司的正常發(fā)展。

(2)

式中:ρl1,ρl2和μl1,μl2分別是液相1和液相2的密度和黏度;ρg和μg分別為氣相的密度和黏度[16-18]。

Gerris中各項(xiàng)采用相應(yīng)的格式進(jìn)行離散。時間采用經(jīng)典的分裂投影法,達(dá)到二階精度?？臻g采用四叉樹/八叉樹進(jìn)行離散,達(dá)到二階精度,這使得自適應(yīng)加密算法實(shí)現(xiàn)更簡易靈活,在不損失計(jì)算精度的情況下顯著降低了計(jì)算量,非常適合處理存在多尺度流動的瞬態(tài)流場計(jì)算問題[18]。PLIC VOF幾何重構(gòu)方法非常適合用于包含破碎、聚合現(xiàn)象的多相流計(jì)算過程,圖2為霧化過程仿真自適應(yīng)生成的網(wǎng)格。Francois在采用連續(xù)表面張力模型(continuum surface force,CSF)計(jì)算液相表面張力的基礎(chǔ)上,通過將表面張力轉(zhuǎn)化為相界面附近的連續(xù)體積力源項(xiàng),結(jié)合高度函數(shù)曲率估計(jì)實(shí)現(xiàn)了表面張力的精確求解,并對壓力進(jìn)行修正,確保表面張力和壓力梯度之間的平衡[19]。采用單調(diào)集成大渦模擬MILES(又稱隱式大渦模擬ILES)[20]近似模擬亞格子SGS的能量傳遞。

圖2 基于PLIC VOF方法的AMR技術(shù)Fig.2 AMR based on PLIC VOF method

圖3 PLIC VOF方法Fig.3 PLIC VOF method

1.2 氣液相界面捕捉及表面張力

Gerris使用PLIC VOF方法對界面進(jìn)行重構(gòu)。PLIC VOF方法定義函數(shù)c作為第一相體積分?jǐn)?shù)來描述相界面(包含相界面的網(wǎng)格為0

Gerris采用的PLIC VOF方法對氣液相界面捕捉分兩步進(jìn)行:①界面重構(gòu),即如圖3所示;②幾何通量計(jì)算和界面對流的處理[13]。界面法向n可以通過考慮相鄰單元的體積分?jǐn)?shù),在八叉樹空間離散產(chǎn)生的3×3×3模板中使用Aulisa等提出的Mixed-Youngs-Centred(MYC)格式進(jìn)行計(jì)算[24]。一旦實(shí)現(xiàn)界面重構(gòu),在規(guī)則笛卡爾網(wǎng)格中采用交替方向方法很容易計(jì)算得到幾何通量,體積分?jǐn)?shù)對流項(xiàng)采用文獻(xiàn)[25]實(shí)現(xiàn)的方向交替格式計(jì)算,進(jìn)而對式(2)中體積分?jǐn)?shù)輸運(yùn)方程進(jìn)行離散計(jì)算。這種對流格式可以捕捉到銳利的界面,并且已在實(shí)際應(yīng)用中證明達(dá)到二階精度。

1.3 計(jì)算模型

對于平面針?biāo)ǘ鄧娮卧?建立如圖4所示計(jì)算域模型,至少選取3個噴注單元作為計(jì)算對象,徑向孔間距根據(jù)阻塞率計(jì)算確定。其中h為軸向液膜厚度,u1為軸向速度,d為徑向液束直徑,u2為徑向速度。計(jì)算域由10～46個L×L×L基本結(jié)構(gòu)Box構(gòu)成,根據(jù)計(jì)算需求選用不同的Box數(shù)量構(gòu)成不同大小的計(jì)算域。計(jì)算域左部和底部分別為軸向與徑向的速度入口邊界。壁面均為無滑移壁面條件,兩側(cè)設(shè)置為對稱面邊界,其余為出口,采用Outflow邊界條件,背壓為大氣環(huán)境。設(shè)置第一相和第二相均為水,對應(yīng)兩路液體,第三相為空氣。計(jì)算域Box的長度為L=10 mm,最高網(wǎng)格等級為9,對應(yīng)的最小網(wǎng)格約為19.5 μm,最低等級為6,同時絕大部分區(qū)域網(wǎng)格不超過100 μm。將兩種液相的體積分?jǐn)?shù)及其梯度設(shè)置為網(wǎng)格自適應(yīng)函數(shù),即網(wǎng)格會自動隨著體積分?jǐn)?shù)及其梯度大小而加密/粗化。具體結(jié)構(gòu)參數(shù)如表1所示,工況參數(shù)為u1=20 m/s,u2=25 m/s。

圖4 計(jì)算域模型Fig.4 Model of simulation zone

表1 徑向圓形孔的結(jié)構(gòu)參數(shù)

1.4 數(shù)據(jù)處理方法

圖5 測點(diǎn)設(shè)置示意圖Fig.5 Schematic diagram of measuring point arranged

2 仿真驗(yàn)證與結(jié)果分析

2.1 多噴注單元噴霧場結(jié)構(gòu)特征分析

為了降低多噴注單元數(shù)值仿真的計(jì)算量,同時又能說明多噴注單元特有的噴霧場結(jié)構(gòu),本節(jié)一開始先選用兩個徑向圓孔的噴注單元開展分析,重點(diǎn)關(guān)注噴注單元間相互影響造成的特征差異,后續(xù)將針對特征差異形成的原因開展分析。

從圖6所示的多噴注單元的噴霧場數(shù)值仿真結(jié)果可以看到,與單噴注單元噴霧場仿真結(jié)果(如圖7所示)相比,由于相鄰噴注單元之間存在相互影響,多噴注單元的噴霧場主要存在以下幾個特有結(jié)構(gòu):相鄰兩霧扇相撞背部呈脊?fàn)罱Y(jié)構(gòu),該結(jié)構(gòu)使得霧化區(qū)域大于霧化角;兩個霧扇相撞會在中間截面匯聚成一道薄液膜,使得液滴在整個霧化角范圍內(nèi)均有分布;相鄰兩孔之間會有一定寬度的液膜下漏,下漏液膜寬度和下漏率是重要的參數(shù);液滴空間分布及混合比分布也發(fā)生變化。

圖6 多噴注單元噴霧場結(jié)構(gòu)特征Fig.6 Spray field structures of multi-injector elements

圖7 單噴注單元噴霧場結(jié)構(gòu)Fig.7 Spray field structures of injector element

圖8 受相鄰噴注單元影響的噴霧扇變形過程Fig.8 Deformation process of spray fan affected by adjacent injector elements

圖9 相鄰兩個噴注單元的噴霧扇結(jié)構(gòu)Fig.9 Spray fan structure for adjacent two injector elements

圖10 多噴注單元的噴霧扇結(jié)構(gòu)Fig.10 Spray fan structure for multi-pintle injector

圖11和圖12分別為表1中多噴注單元工況參數(shù)下的數(shù)值仿真和高速攝影試驗(yàn)拍攝的結(jié)果,可以看到三者外形非常相似。

筆者在之前的研究結(jié)果中還對試驗(yàn)和數(shù)值仿真測得的霧化角及單噴注單元的液滴平均粒徑進(jìn)行了定量對比驗(yàn)證[26],相同工況下霧化角的試驗(yàn)測量值和仿真值分別為72.88°和70.28°,相對誤差為3.57%;單噴注單元液滴平均粒徑的試驗(yàn)測量值和仿真值分別為175 μm和158 μm,相對誤差為9.71%。上述結(jié)果表明該仿真方法的準(zhǔn)確性較高,也表明上述分析的變形過程及揭示的相鄰噴注單元間相互影響的機(jī)制是準(zhǔn)確合理的。

圖11 多噴注單元的噴霧扇數(shù)值仿真結(jié)果Fig.11 Numerical result of spray fan for multi-pintle injector

圖12 多噴注單元的噴霧扇試驗(yàn)拍攝結(jié)果Fig.12 Experimental result of spray fan for multi-pintle injector

另外,從多噴注單元的結(jié)果可以看出,相鄰噴霧扇由于受到相互影響被擠壓,最終形成的噴霧場呈扁平的多凹腔狀,每個噴注單元將整個噴霧場內(nèi)部分隔成一個個相對獨(dú)立的凹腔空間,而外圍拼接成一個大霧錐,與單噴注單元的噴霧場不同;同時,液膜破碎后的液霧在空間也呈現(xiàn)相應(yīng)的不規(guī)則曲面分布,從而使得噴霧場的液滴空間分布及混合比分布也發(fā)生相應(yīng)變化。

2.2 單噴注單元與多噴注單元的對比

針對表1相同結(jié)構(gòu)參數(shù)(d=0.8 mm,h=0.45 mm)和工作參數(shù)(u1=20 m/s,u2=25 m/s)的單噴注單元和多噴注單元,通過數(shù)值仿真獲得了噴霧場液滴粒徑信息。獲得的單噴注單元與多噴注單元的噴霧場計(jì)算結(jié)果對比如圖13所示,對應(yīng)的全場液滴粒徑概率密度分布如圖14所示,液滴粒徑沿徑向的空間分布如圖15所示,相對流強(qiáng)(即流強(qiáng)概率密度分布)及混合比沿徑向的空間分布分別如圖16和圖17所示。通過對比可以看到,正如2.1節(jié)所述,單噴注單元與多噴注單元的噴霧場存在顯著不同。

圖13 噴霧場計(jì)算結(jié)果Fig.13 Simulation results of spray field

圖14 全場液滴粒徑概率密度分布Fig.14 Droplet diameter probability density distribution in the whole spray field

從圖14可以看出,與單噴注單元相比,多噴注單元液膜路和液束路的液滴分布峰值均偏向更大粒徑方向,大液滴數(shù)目占比明顯增加,液滴平均粒徑S均明顯增大,分別從119.38 μm和127.44 μm增大至170.29 μm和171.73 μm,分別增大了42.6%和34.8%。

圖15 液滴粒徑沿徑向的空間分布Fig.15 Radial distribution of droplet diameter

如圖16所示,實(shí)際的多噴注單元噴霧扇僅有較小半徑處的部分發(fā)生重疊,液膜厚度加倍。在較小半徑位置采用虛線液膜厚度加倍的理論結(jié)果,在大半徑位置采用實(shí)線單噴注單元原來液膜厚度的理論結(jié)果,中間半徑位置介于兩種結(jié)果之間。由于中間過渡區(qū)霧場結(jié)構(gòu)十分復(fù)雜,目前還沒更好的理論模型進(jìn)行解析描述,還不能形成一個統(tǒng)一的解析表達(dá)式。目前只能用實(shí)線和虛線組合起來描述多噴注單元。這樣的模型雖存在不足,但相對簡單,可以滿足一定精度要求下的理論預(yù)估使用需求。

圖16 相鄰霧扇的相互影響Fig.16 Interaction between adjacent spray fans

圖17 相對流強(qiáng)沿徑向的空間分布Fig.17 Radial distribution of relative mass flow rate flux

圖18 混合比沿徑向的空間分布Fig.18 Radial distribution of mixture ratio

3 結(jié)論

為了研究相鄰噴注單元間相互影響對針?biāo)ㄊ絿娮⑵鲊婌F場的影響,基于AMR技術(shù)和分相識別的PLIC VOF新方法,實(shí)現(xiàn)了針?biāo)ǘ鄧娮卧F化過程的高保真數(shù)值模擬,并對比了其與噴注單元噴霧場結(jié)構(gòu)及特性的差異,揭示了相鄰噴注單元間的相互影響機(jī)制,得出以下結(jié)論。

1)數(shù)值仿真獲得多噴注單元噴霧扇結(jié)構(gòu)與高速攝影試驗(yàn)拍攝的結(jié)果對比,兩者外形非常相似,試驗(yàn)與數(shù)值仿真的霧化角相對誤差為3.57%,液滴平均粒徑相對誤差為9.71%,表明新的仿真方法在精細(xì)研究針?biāo)ㄊ絿娮⑵鲊婌F場方面具有較好的準(zhǔn)確性,也表明受相鄰噴注單元影響的噴霧扇變形過程分析是準(zhǔn)確合理的。

2)與單噴注單元相比,多噴注單元噴霧場主要存在以下特殊結(jié)構(gòu):相鄰兩霧扇相撞背部呈脊?fàn)罱Y(jié)構(gòu),使得霧化區(qū)域大于霧化角;兩霧扇相撞在中間對稱面匯聚形成薄液膜,使整個霧化角范圍內(nèi)均有液滴分布;相鄰兩孔之間形成一定下漏率和下漏液膜寬度;液滴空間分布均發(fā)生顯著變化。

3)相鄰噴注單元間的相互影響機(jī)制為:相鄰噴霧扇相撞后原先各自向外展開的霧扇被擠回中心對稱面,變形成薄液膜,其厚度是原霧扇的兩倍,其他未發(fā)生撞擊位置的液膜厚度保持不變,最終形成的噴霧扇結(jié)構(gòu)呈扁平的多凹腔狀。

4)單噴注單元與多噴注單元的噴霧場特性存在顯著不同。受相鄰噴注單元間相互作用影響,多噴注單元液膜路和液束路的液滴粒徑均顯著增大了約35%,流強(qiáng)和混合比沿徑向分布更趨于均勻。