界面追蹤和拉格朗日粒子追蹤耦合計(jì)算方法

劉 靜，殷 義，時(shí)明偉

（北京航化節(jié)能環(huán)保技術(shù)有限公司，北京100176）

根據(jù)研究目的和需要可以將氣液兩相流動(dòng)的數(shù)值模擬方法分為兩大類(lèi)。 （1）拉格朗日粒子追蹤方法：這種方法發(fā)展歷史長(zhǎng)，應(yīng)用也比較廣泛。 拉格朗日粒子追蹤方法也稱(chēng)為離散顆粒模型， 顆粒被看作離散相，而氣相被視為連續(xù)相。這種方法對(duì)每個(gè)顆粒與氣體以及顆粒與顆粒間的作用都詳細(xì)考慮。 對(duì)于離散顆粒的破碎、 蒸發(fā)等物理過(guò)程都有相應(yīng)的子模型進(jìn)行描述。 （2）界面追蹤方法：這是一種較新的研究?jī)上嗔鲃?dòng)的方法。 該方法可以通過(guò)計(jì)算兩相交界面的運(yùn)動(dòng)，再現(xiàn)液相在氣相作用下的詳細(xì)霧化過(guò)程，從而達(dá)到研究氣液霧化機(jī)理的目的。 目前常用的有體積分?jǐn)?shù)方法和Level Set 方法等。

本文首先分別采用以上介紹的兩大類(lèi)方法對(duì)射流和液滴進(jìn)行了計(jì)算，針對(duì)兩種方法的優(yōu)缺點(diǎn)，提出了新的改進(jìn)的計(jì)算方法。 結(jié)合實(shí)驗(yàn)測(cè)量結(jié)果,對(duì)采用改進(jìn)方法的計(jì)算結(jié)果與兩種常用計(jì)算方法的結(jié)果進(jìn)行了對(duì)比分析。

1 計(jì)算方法

1.1 拉格朗日粒子追蹤方法

拉格朗日方法中，氣相控制方程為三維N-S 方程。液相采取拉格朗日顆粒跟蹤方法。液滴對(duì)周?chē)鷼庀嗟挠绊懕硎緸闅庀喾匠讨械脑聪郤。 而氣相對(duì)液相的影響則通過(guò)在計(jì)算液滴狀態(tài)的過(guò)程中， 當(dāng)用到氣相參數(shù)時(shí)直接采用液相所屬當(dāng)?shù)鼐W(wǎng)格的氣相參數(shù)來(lái)實(shí)現(xiàn)。

三維氣相N-S 方程為：

式中: U—守恒變矢量

F—x 方向的無(wú)粘矢通量

G—y 方向的無(wú)粘矢通量

H—z 方向的無(wú)粘矢通量

Fv—x 方向的粘性矢通量

Gv—y 方向的粘性矢通量

Hv—z 方向的粘性矢通量

t—時(shí)間變量，s

x—x 坐標(biāo)方向，m

y—y 坐標(biāo)方向，m

z—z 坐標(biāo)方向，m

S—液相影響源相矢量

液滴在運(yùn)動(dòng)過(guò)程中受到了空氣阻力和重力，液滴單位質(zhì)量受力用下式計(jì)算：

式中：下標(biāo)g—?dú)庀鄥?shù)

下標(biāo)d—液相參數(shù)

u—速度矢量，m/s

r—液滴半徑，m

g—重力加速度，m/s2

阻力系數(shù)CD計(jì)算公式為：

氣相雷諾數(shù)Re 為：

式中，μ 為氣體黏性系數(shù)。

單個(gè)液滴的運(yùn)動(dòng)速度計(jì)算公式為：

得到液滴的運(yùn)動(dòng)軌跡計(jì)算公式為：

液滴在運(yùn)動(dòng)過(guò)程中的霧化和蒸發(fā)由相應(yīng)的子模型來(lái)描述。本文中只考慮霧化過(guò)程，液滴的霧化采用K-H 和R-T 混合霧化模型。

計(jì)算流程如圖1 所示。

圖1 計(jì)算流程示意圖

1.2 界面追蹤方法

假定在多介質(zhì)的界面上，壓力和速度為常數(shù)，而其它物理量在界面上可能存在著間斷。 以此假設(shè)為出發(fā)點(diǎn)， 可以得到多介質(zhì)可壓縮流動(dòng)的一維模型方程組，具體推導(dǎo)過(guò)程可以參考文獻(xiàn)[4]、[6]、[7]。

簡(jiǎn)便起見(jiàn)， 這里給出采用體積分?jǐn)?shù)形式描述的一維多介質(zhì)可壓縮流動(dòng)的歐拉方程組：

各介質(zhì)的狀態(tài)方程采用如下形式：

式中：ρ—密度，kg/m3

u—x 方向速度，m/s

p—壓強(qiáng)，Pa

e—單位質(zhì)量的內(nèi)能，kJ/kg

E=e+u2/2，單位質(zhì)量的總能量，kJ/kg

N 為介質(zhì)種類(lèi)總數(shù)，Y（i）為第 i 種介質(zhì)的體積分?jǐn)?shù)。 在一個(gè)網(wǎng)格中，N 種介質(zhì)的體積分?jǐn)?shù)總和為1。

若已知各組分的物理量 ρ（i），u（i），p（i），γ（i），p∞（i）。 各組分壓力與網(wǎng)格壓力相等 p=p（i），i=1，…，N。 網(wǎng)格中的平均物理量計(jì)算方法為各組分物理量的體積加權(quán)平均值。

狀態(tài)方程可表達(dá)為

求解方法采用的一維Lagrangian-Remapping 型PPM（Piecewise Parabolic Method）方法。 具體可分為以下四步：（1）各物理量的分段拋物差值；（2）近似黎曼問(wèn)題求解網(wǎng)格邊界上時(shí)間平均的壓力和速度；（3）拉格朗日方程組推進(jìn)求解；（4）將物理量轉(zhuǎn)換到靜止的歐拉網(wǎng)格上。 對(duì)于三維方程的求解采用Strang 分裂方法， 將多維的求解轉(zhuǎn)化為多個(gè)一維方向交替計(jì)算。 采用Strang 分裂計(jì)算簡(jiǎn)單而且有助于編制并行效率很高的程序，對(duì)開(kāi)展大規(guī)模并行計(jì)算很有幫助。

2 計(jì)算結(jié)果和分析

2.1 計(jì)算條件

算例一：矩形超聲速風(fēng)洞長(zhǎng)762 mm、寬152 mm、高 127 mm，來(lái)流 Ma 為 1.94，流場(chǎng)總溫 T* 為 533 K，總壓P* 為206 kPa，噴注器為直流形式，出口直徑d為0.5 mm，安裝在風(fēng)洞底部，動(dòng)壓比為10。

算例二：液滴在激波管中發(fā)生破碎，氣流速度為432 m/s，液滴直徑為2.6 mm。 該工況為文獻(xiàn)[9]中的實(shí)驗(yàn)工況，對(duì)于該工況文獻(xiàn)[9]還給出了實(shí)驗(yàn)結(jié)果，以便于對(duì)比分析。

2.2 拉格朗日粒子追蹤法計(jì)算結(jié)果

圖2 為三維噴霧流動(dòng)的液霧分布圖。 當(dāng)液體射流從噴嘴噴出后不久，在與高速氣流的相互作用下，一邊不斷地發(fā)生破碎生成新的液滴， 一邊隨著氣流向下游運(yùn)動(dòng)。

圖2 液霧分布圖

圖3 為液霧霧化過(guò)程中的滴徑分布圖。 可以看出，當(dāng)?shù)扔趪娍谥睆降拇笠旱螐膰娍趪姵龊?，在KH 不穩(wěn)定波的作用下， 逐漸從上面剝離出直徑較小的液滴，并跟隨氣流向下游運(yùn)動(dòng)，形成液柱后的一片液霧， 此時(shí)大液滴在運(yùn)動(dòng)過(guò)程中隨著小液滴的逐漸剝離直徑也慢慢減小，這就類(lèi)似于液柱在剛噴出后，由于氣流的作用會(huì)從其上面剝離走一部分小液滴。而當(dāng)大液滴運(yùn)動(dòng)到其累計(jì)時(shí)間大于破碎特征時(shí)間時(shí)，R-T 不穩(wěn)定模型此時(shí)也加入競(jìng)爭(zhēng)。 在此之后，大液滴進(jìn)一步被剝離或破碎成較小的液滴， 液霧外圍的液滴直徑大于靠近壁面的液滴。

圖3 液霧霧化分布圖

采用拉格朗日粒子追蹤法，優(yōu)點(diǎn)是計(jì)算量較小，占用的計(jì)算存儲(chǔ)資源少。但在計(jì)算霧化子模型中，需要引進(jìn)一些計(jì)算經(jīng)驗(yàn)常數(shù)， 有時(shí)在計(jì)算中經(jīng)驗(yàn)參數(shù)的選擇對(duì)計(jì)算結(jié)果影響較大。

2.3 界面追蹤方法計(jì)算結(jié)果

圖4 為射流從底部垂直于來(lái)流方向噴入燃燒室后的形態(tài)圖。 可以看出，當(dāng)射流從下底面噴出后，在來(lái)流的壓力作用下， 圓柱流由剛噴出時(shí)的圓形截面迅速展開(kāi)呈現(xiàn)一個(gè)拱形的包絡(luò)面， 液柱在隨著氣流彎向下游的同時(shí)，由圓形截面展成了片狀截面?？梢钥吹皆谝褐L(fēng)面距離下底面一定距離處有明顯波動(dòng)的起伏。 而在液柱的兩側(cè)由于氣流對(duì)液體的剪切作用，呈現(xiàn)出一些條狀的凸起，有一些大的液塊已經(jīng)脫離主液柱隨流運(yùn)動(dòng)。

圖4 射流形態(tài)對(duì)比圖

對(duì)算例二液滴霧化過(guò)程也進(jìn)行了計(jì)算。 在高速氣體的作用下，液滴迅速霧化完畢。圖5 為液體體積分?jǐn)?shù)為0.01 的液滴霧化過(guò)程的等值面圖。 當(dāng)45 μs時(shí)，圓形的液滴被氣流逐漸壓扁，邊緣被氣流撕扯出一些不規(guī)則邊緣；當(dāng)75 μs 時(shí)，邊緣上的液膜和凸起進(jìn)一步被氣流作用， 液滴被壓得更扁；130 μs 時(shí)，液滴的凸起伸長(zhǎng)，呈現(xiàn)規(guī)則的八角形；當(dāng)時(shí)間達(dá)到170 μs 時(shí)，液滴已經(jīng)伸展成一片狀結(jié)構(gòu)，且迎風(fēng)面明顯出現(xiàn)一圈規(guī)則的凹坑；到290 μs 時(shí)，凹坑處破裂，液滴完全分解成塊狀。

通過(guò)對(duì)比計(jì)算結(jié)果圖與實(shí)驗(yàn)結(jié)果圖（見(jiàn)圖6）發(fā)現(xiàn)， 計(jì)算結(jié)果對(duì)液滴主體的追蹤基本上與實(shí)驗(yàn)觀(guān)測(cè)結(jié)果相一致， 但對(duì)于更加細(xì)小的液霧的追蹤用目前的界面追蹤方法就無(wú)能為力了。 因?yàn)榧?xì)密的液霧已經(jīng)達(dá)到了微米量級(jí)， 理論上來(lái)說(shuō)要想捕捉這些細(xì)小的液霧， 需要比液霧直徑還要小的網(wǎng)格分辨率才能實(shí)現(xiàn)， 而這對(duì)于目前的計(jì)算條件來(lái)說(shuō)是不可能實(shí)現(xiàn)的。而且本文所采用的體積分?jǐn)?shù)模型在計(jì)算過(guò)程中，本身對(duì)界面的捕捉并不是清晰的， 而是跨越了幾個(gè)網(wǎng)格，所以這也是目前計(jì)算存在的局限性。

3 改進(jìn)方法和計(jì)算結(jié)果

3.1 改進(jìn)方法

圖5 液滴二次霧化體積分?jǐn)?shù)等值面圖

圖6 液滴二次霧化實(shí)驗(yàn)結(jié)果圖

前面分別用界面追蹤和拉格朗日粒子追蹤方法對(duì)霧化過(guò)程進(jìn)行了計(jì)算， 在對(duì)結(jié)果的分析過(guò)程中發(fā)現(xiàn)，界面追蹤方法能夠真實(shí)地再現(xiàn)霧化過(guò)程，是研究射流或液滴的變形破碎、 不同的因素對(duì)霧化結(jié)果的影響程度以及進(jìn)一步的霧化機(jī)理分析的有效工具。界面追蹤隨著計(jì)算機(jī)技術(shù)的發(fā)展而不斷發(fā)展起來(lái)，但它同樣存在自身的不足之處， 比如為了精確的捕捉氣液界面，需要在界面附近布置極細(xì)的網(wǎng)格。如果需要捕捉剝離的液滴，那么就需要更加細(xì)密的網(wǎng)格。在高速氣流中，霧化后的小液滴通常是微米量級(jí)的，所以如果想對(duì)整個(gè)霧化過(guò)程（一次霧化和二次霧化）進(jìn)行計(jì)算， 那么最小的網(wǎng)格尺度必須小于霧化后的小液滴的尺寸，可以想象，如果想捕捉整個(gè)液霧的分布，那么對(duì)網(wǎng)格的細(xì)密程度要求是極高的，這將造成巨大的計(jì)算量， 甚至可以說(shuō)這樣的計(jì)算是不經(jīng)濟(jì)也是幾乎不可能實(shí)現(xiàn)的。 而拉格朗日粒子追蹤方法是工程計(jì)算中常用的一種方法，它的優(yōu)點(diǎn)是計(jì)算量小，以追蹤多組典型的液滴分布來(lái)代替整個(gè)液霧的計(jì)算。它雖然不能對(duì)具體的霧化過(guò)程進(jìn)行描述，但對(duì)于更關(guān)心液霧運(yùn)動(dòng)和分布等宏觀(guān)現(xiàn)象的研究， 這種方法更具有優(yōu)勢(shì)。在拉格朗日粒子追蹤方法中，對(duì)霧化過(guò)程的模擬用一些霧化模型來(lái)代替。 由于對(duì)霧化機(jī)理研究的缺乏， 有限的霧化模型并不適用于所有的霧化過(guò)程。特別是對(duì)于一次霧化過(guò)程，到目前為止也沒(méi)有一個(gè)專(zhuān)門(mén)針對(duì)一次霧化的霧化模型可供使用。一般的做法是省略一次霧化， 認(rèn)為在出口處一次霧化已經(jīng)結(jié)束， 直接在出口處給出實(shí)驗(yàn)測(cè)量的液霧尺寸分布；或認(rèn)為出口處液滴尺寸近似等于噴孔直徑，即“Blob 模型”，總之都是一種近似。

為了更精細(xì)的模擬霧化過(guò)程， 綜合兩種方法的優(yōu)缺點(diǎn)， 嘗試把界面追蹤和拉格朗日粒子追蹤方法結(jié)合起來(lái)。 對(duì)于無(wú)法用霧化模型描述的一次霧化過(guò)程采用界面追蹤的方法來(lái)直接模擬， 而對(duì)之后的霧化過(guò)程采用拉格朗日粒子追蹤的方法。 這樣做的好處是：（1）采用界面追蹤的方法很好的解決了一次霧化過(guò)程的模擬；（2）用拉格朗日粒子追蹤方法對(duì)霧化后細(xì)密粒子的追蹤極大節(jié)約了計(jì)算網(wǎng)格， 提高了計(jì)算效率。這樣做結(jié)合了兩種方法的優(yōu)點(diǎn)，而且避免了各自的不足之處。表1 對(duì)幾種計(jì)算方法進(jìn)行了對(duì)比。

表1 幾種計(jì)算方法對(duì)比

界面追蹤和拉格朗日粒子追蹤方法在前面已經(jīng)詳細(xì)介紹過(guò)， 這里主要介紹一下兩種模型在結(jié)合過(guò)程中的一些具體方法。 主要的問(wèn)題是兩種方法的轉(zhuǎn)換方式，需要確定的重要參數(shù)如下。

（1）轉(zhuǎn)換位置。 掃描全場(chǎng)，轉(zhuǎn)換位置發(fā)生在同時(shí)滿(mǎn)足以下條件的網(wǎng)格中： ①該網(wǎng)格內(nèi)的液相體積分?jǐn)?shù)大于0，具體量值可根據(jù)具體情況輸入；②該網(wǎng)格正后方的網(wǎng)格非純液體網(wǎng)格， 即該網(wǎng)格正后方是氣體或氣液混合區(qū)。若滿(mǎn)足條件后，每個(gè)網(wǎng)格生成一個(gè)液滴組，其中所有的液體或部分液體轉(zhuǎn)化為液滴。

（2）轉(zhuǎn)化后液滴的初始速度。若相鄰網(wǎng)格有純液網(wǎng)格，則速度等于純液網(wǎng)格的速度中最小的那個(gè)；若相鄰沒(méi)有純液網(wǎng)格， 則按體積分?jǐn)?shù)取當(dāng)?shù)貧庀嗑W(wǎng)格速度的一部分。

（3）液滴初始位置。 生成該液滴的當(dāng)?shù)鼐W(wǎng)格位置。

（4）粒子大小。 由二次霧化模型求得，若已知液滴直徑，則液滴數(shù)也就可以得到。

（5）氣液相對(duì)速度的計(jì)算。由前面得到的液滴速度，可以計(jì)算出轉(zhuǎn)化液滴后當(dāng)?shù)貧庀嗟乃俣?，從而得到相?duì)速度。

（6）當(dāng)?shù)鼐W(wǎng)格計(jì)算考慮剝離液滴帶走的動(dòng)能，并設(shè)體積分?jǐn)?shù)為0，并相應(yīng)地修改密度。

3.2 計(jì)算流程

圖7 為耦合計(jì)算方法流程圖。 由于界面追蹤的計(jì)算時(shí)間步長(zhǎng)遠(yuǎn)小于拉格朗日粒子跟蹤計(jì)算時(shí)間步長(zhǎng)， 所以在進(jìn)行了若干步界面追蹤計(jì)算后進(jìn)行一次拉格朗日粒子追蹤計(jì)算， 使二者的推進(jìn)時(shí)間保持一致。拉格朗日粒子追蹤計(jì)算完成后，將體積分?jǐn)?shù)和氣相場(chǎng)源相的變化傳遞給界面追蹤流場(chǎng)計(jì)算。 如此交替計(jì)算，直至滿(mǎn)足計(jì)算條件為止。

3.3 計(jì)算結(jié)果

對(duì)氣流中的橫向射流破碎和液滴霧化進(jìn)行界面追蹤和拉格朗日粒子追蹤方法的耦合計(jì)算。 計(jì)算工況與前面相同。射流從噴注器噴出后，在氣流的作用下迅速破碎成液塊和小液滴。 圖8 是兩種方法耦合后的計(jì)算結(jié)果。 圖8（a）是xy 向液霧分布，三角形和圓形的標(biāo)記點(diǎn)分別表示的是實(shí)驗(yàn)測(cè)量擬合結(jié)果，三角形為陰影法的結(jié)果，圓形為PDPA 的結(jié)果。二者之所以出現(xiàn)差別， 是由于PDPA 捕捉到一些分布稀疏的用陰影法無(wú)法觀(guān)察到的外圍液滴， 所以PDPA 的穿透深度高于陰影法。 本文的計(jì)算結(jié)果對(duì)于液霧濃密區(qū)其穿透深度與陰影法接近。 而外圍擴(kuò)散的稀疏的液滴分布在液霧前端與PDPA 結(jié)果吻合， 后面有個(gè)別的液滴分布偏高于PDPA 的計(jì)算結(jié)果。 圖8（b）顯示的是xz 向俯視圖。 液霧分布呈現(xiàn)中間濃密兩側(cè)稀疏的趨勢(shì)。從整體來(lái)說(shuō)，計(jì)算結(jié)果與實(shí)驗(yàn)測(cè)量結(jié)果基本吻合。

圖7 耦合計(jì)算方法流程圖

圖8 液柱破碎液霧分布圖

圖 9（a）是圖 8（a）中方框區(qū)域放大圖。和圖 9（b）純界面追蹤計(jì)算結(jié)果相比， 液柱的部分轉(zhuǎn)化為液滴并隨流運(yùn)動(dòng)。這也是耦合方法的優(yōu)勢(shì)所在，純界面追蹤在液柱部分模擬得很好， 但對(duì)于射流后部濃密液霧的計(jì)算卻無(wú)能為力。采用耦合方法，能自然地根據(jù)射流特性計(jì)算出其穿透深度， 而不再需要像使用混合霧化模型時(shí)， 須給出液柱破碎時(shí)間判斷標(biāo)準(zhǔn)或破碎點(diǎn)的判斷標(biāo)準(zhǔn)。

圖9 耦合法和單純界面追蹤法對(duì)比圖

通過(guò)耦合計(jì)算， 對(duì)于液滴霧化的整個(gè)流場(chǎng)進(jìn)行了數(shù)值模擬。前段母液滴采用界面追蹤法進(jìn)行計(jì)算，之后的液霧區(qū)采用拉格朗日粒子追蹤方法進(jìn)行計(jì)算。 圖10 為計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的對(duì)比。45 μs 時(shí)，開(kāi)始有液霧從液滴上剝離生成；75 μs 時(shí)， 液霧繼續(xù)發(fā)展，且由實(shí)驗(yàn)結(jié)果看到，液滴上下兩側(cè)的液霧濃度較中間區(qū)域偏高；隨著時(shí)間的推移，母液滴生成的液霧越來(lái)越多，濃度也越來(lái)越大。由圖11 可見(jiàn)，從始至終液霧始終產(chǎn)生于整個(gè)母液滴的邊緣區(qū)， 也就是變形最劇烈的地方，液霧分布呈現(xiàn)環(huán)狀。

圖10 二次霧化計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比圖

圖11 yz 向界面分布圖

4 結(jié)論

為對(duì)霧化過(guò)程進(jìn)行全面而準(zhǔn)確的模擬， 結(jié)合界面追蹤和拉格朗日粒子追蹤的優(yōu)點(diǎn)， 采用兩種方法耦合的思路對(duì)射流破碎和液滴霧化進(jìn)行了數(shù)值模擬。 射流破碎的計(jì)算過(guò)程在初始段采用界面追蹤方法，比采用簡(jiǎn)單的一次霧化簡(jiǎn)化模型（Blob 模型）的計(jì)算結(jié)果更符合實(shí)際。與之前采用的混合模型相比，優(yōu)點(diǎn)在于減少了經(jīng)驗(yàn)參數(shù)的使用， 不再需要在使用混合霧化模型時(shí)給出經(jīng)驗(yàn)破碎時(shí)間或破碎點(diǎn)的判斷標(biāo)準(zhǔn)。液滴霧化的結(jié)果更顯示了耦合方法的優(yōu)越性，對(duì)于單個(gè)液滴霧化過(guò)程的模擬， 如果單獨(dú)采用界面追蹤或拉格朗日粒子追蹤方法都不能夠得到真實(shí)的霧化現(xiàn)象， 用耦合方法得到的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果定性符合較好。 耦合方法結(jié)合了兩種計(jì)算方法的優(yōu)勢(shì)，摒棄了不足，有助于研究一些對(duì)于兩相界面形狀和液霧分布都很重要的問(wèn)題。