二維噴射液滴在重力場(chǎng)中的運(yùn)動(dòng)模型研究

劉展位,薄涵亮

(清華大學(xué) 核能與新能源技術(shù)研究院,北京 100084)

歐拉與拉格朗日是流體力學(xué)描述的兩種方法,拉格朗日方法著眼于描述流體質(zhì)點(diǎn)的位置隨時(shí)間變化的規(guī)律,跡線即為流體質(zhì)點(diǎn)在空間運(yùn)動(dòng)所描繪的曲線。歐拉方法著眼于空間點(diǎn),在空間每一點(diǎn)上面描述流體運(yùn)動(dòng)隨時(shí)間變化的情況,采用流線來(lái)描述流場(chǎng)。歐拉與拉格朗日方法的解耦問(wèn)題一直是流體力學(xué)研究的熱點(diǎn)和難點(diǎn)問(wèn)題。根據(jù)Yang[1]對(duì)多相流模擬的總結(jié),多相流模型包括歐拉-歐拉(EE)模型[2]、MUSIG模型、非齊次MUSIG模型[3]。其中EE模型計(jì)算代價(jià)較小但模擬精度相對(duì)比較有限。非齊次MUSIG模型擴(kuò)展了齊次的MUSIG模型,能夠使不同尺寸的離散相具有不同的速度場(chǎng),主要通過(guò)求解多組氣相的Navier-Stokes方程,以及改進(jìn)的連續(xù)性方程,從而獲得不同的離散相速度場(chǎng),其缺點(diǎn)是增加了計(jì)算成本,且一般速度組為2個(gè),因此對(duì)離散相速度組的描述受到一定限制。

基于粒子的方法包括歐拉-拉格朗日方法以及其他粒子法。歐拉-拉格朗日方法包括離散相模型(DPM)[4]、稠密離散相模型(DDPM)[5]、離散元法(DEM)[6]、光滑粒子法(SPH)[7]等。其他粒子法包括格子玻爾茲曼方法(LBM)[8]、耗散粒子動(dòng)力學(xué)(DPD)[9]。本課題組采用歐拉-拉格朗日方法的整體框架,對(duì)汽水分離器內(nèi)液滴行為進(jìn)行建模,從探究汽水分離器中液滴的微觀機(jī)理到獲取流場(chǎng)宏觀參數(shù),從而優(yōu)化設(shè)計(jì)汽水分離器結(jié)構(gòu)的研究路線,并通過(guò)自編程序與商用CFD軟件結(jié)合開發(fā)出一整套計(jì)算方法。這些行為包括:液滴的產(chǎn)生、運(yùn)動(dòng)、碰撞、聚合、破碎、蒸發(fā)、冷凝和消亡等[10]。歐拉-拉格朗日方法雖然較歐拉-歐拉法更能精確地模擬粒子運(yùn)動(dòng),但是計(jì)算量隨著粒子數(shù)目增加而明顯增加。因此,開發(fā)了歐拉網(wǎng)格逼近方法[11],歐拉網(wǎng)格逼近算法已被應(yīng)用于描述小液滴在重力空間中的運(yùn)動(dòng)特征,同時(shí)該方法也被應(yīng)用于多液滴的情況并驗(yàn)證了計(jì)算的可行性[11]。由于實(shí)際蒸汽發(fā)生器的汽水分離以及其他的汽水分離裝置空間內(nèi)具有部分大液滴,因此有必要將該方法拓展到大液滴的工況。

本文提出一種根據(jù)歐拉網(wǎng)格逼近方法來(lái)計(jì)算汽水分離器中大液滴運(yùn)動(dòng)的方法。介紹大液滴的計(jì)算公式和密度定義,驗(yàn)證歐拉網(wǎng)格逼近方法應(yīng)用于大液滴粒子計(jì)算的有效性,最后將歐拉網(wǎng)格逼近方法應(yīng)用于實(shí)際AP1000重力分離空間中大液滴的計(jì)算,分析3種尺寸的噴射液滴在重力分離空間中的行為,為汽水分離等裝置中的大液滴分析提供有效的計(jì)算分析工具。

1 理論分析

1.1 噴射液滴形成機(jī)理

氣泡在自由表面處破裂主要過(guò)程為:首先,當(dāng)氣泡從水室上浮到蒸汽表面時(shí),氣泡上半部(液帽)突出汽空間,水膜不斷減薄最終破裂形成一些細(xì)小的水滴;之后,由于氣液兩相壓力差及流體表面張力的作用,氣泡周圍的水向中心聚集以填補(bǔ)空虛位置形成1個(gè)向上射流。射流以一定速度向上升起且在峰頂斷裂生成幾個(gè)大水滴。因此,可以歸納氣泡破裂產(chǎn)生液滴的兩種主要類型(圖1):第一類氣泡液膜厚度變薄后破裂產(chǎn)生的細(xì)小液滴被稱為膜液滴;第二類由于氣泡破裂的一瞬間在自由表面留下空洞,同時(shí)由于瑞利不穩(wěn)定性的作用,周圍液體涌入填補(bǔ)空缺位置便產(chǎn)生向上的不穩(wěn)定瑞利射流,進(jìn)而在表面張力作用下產(chǎn)生射流斷裂形成的液滴叫做噴射液滴[12]。文獻(xiàn)[13]的實(shí)驗(yàn)已觀察到這兩類液滴,關(guān)于噴射液滴的產(chǎn)生文獻(xiàn)[14]已進(jìn)行了大量的實(shí)驗(yàn)研究與數(shù)值模擬。

圖1 氣泡破裂產(chǎn)生液滴的兩種類型[12]Fig.1 Two types of bubble rupture to produce droplet[12]

文獻(xiàn)[15]對(duì)汽水分離器入口膜液滴的尺寸分布進(jìn)行了實(shí)驗(yàn)研究,其尺寸分布如圖2所示,而文獻(xiàn)[16]等對(duì)瑞利不穩(wěn)定性產(chǎn)生的噴射液滴開展了研究,得到了平均噴射液滴直徑為1.5 mm。

圖2 蒸汽發(fā)生器內(nèi)膜液滴尺寸的分布[15]Fig.2 Distribution of steam generator inner film droplet size[15]

1.2 基本假設(shè)

為了簡(jiǎn)化計(jì)算,假設(shè)如下:1) 外流場(chǎng)處于穩(wěn)態(tài);2) 單個(gè)液滴為剛性圓球,其質(zhì)量和大小都是固定的;3) 不考慮粒子的碰撞、聚并、破碎等作用;4) 液滴和蒸汽之間的傳熱和傳質(zhì)可以忽略不計(jì);5) 單向耦合假設(shè),認(rèn)為粒子對(duì)于流場(chǎng)的作用力可以忽略,只考慮流場(chǎng)對(duì)粒子的作用力。

2 數(shù)學(xué)描述

2.1 數(shù)學(xué)模型

圖3為網(wǎng)格中大液滴示意圖,大液滴的直徑大于網(wǎng)格尺寸,一部分網(wǎng)格能被大液滴完全覆蓋,而另一部分處于液滴邊緣的網(wǎng)格只能被大液滴部分覆蓋。認(rèn)為液滴每步的運(yùn)動(dòng)都是質(zhì)點(diǎn)運(yùn)動(dòng),大液滴覆蓋到的所有網(wǎng)格節(jié)點(diǎn)的速度均一致,而對(duì)于這兩種不同的網(wǎng)格需要采用不同的密度定義。其中D為大液滴直徑,d為網(wǎng)格長(zhǎng)度。

圖3 網(wǎng)格中大液滴示意圖Fig.3 Schematic diagram of large droplet in grid

假設(shè)液滴能完全覆蓋的網(wǎng)格節(jié)點(diǎn)為Nc,不完全覆蓋的網(wǎng)格節(jié)點(diǎn)為Ne,可以得到完全覆蓋的網(wǎng)格節(jié)點(diǎn)所占的面積Ac為:

Ac=Ncd2

(1)

大液滴所占面積Ao為:

(2)

每個(gè)網(wǎng)格面積Am為:

Am=d2

(3)

不完全覆蓋的網(wǎng)格節(jié)點(diǎn)所占面積Ae為:

(4)

不完全覆蓋的網(wǎng)格節(jié)點(diǎn)的系數(shù)k為:

(5)

其中,k的范圍為(0,1)。

2.2 密度公式

歐拉網(wǎng)格逼近方法的基本方程參考文獻(xiàn)[11],認(rèn)為在二維情況下完全覆蓋的網(wǎng)格節(jié)點(diǎn)的密度是一致的,而不完全覆蓋的網(wǎng)格節(jié)點(diǎn)的密度需要乘以1個(gè)修正系數(shù)。因此,對(duì)于大于網(wǎng)格直徑的液滴噴射密度定義如下。

1) 在完全覆蓋的網(wǎng)格中

當(dāng)運(yùn)動(dòng)速度|v(i)|≤Df時(shí),則:

ρl(i)=ρf

(6)

當(dāng)運(yùn)動(dòng)速度|v(i)|>Df時(shí),則:

(7)

2) 在非完全覆蓋的網(wǎng)格中

當(dāng)運(yùn)動(dòng)速度|v(i)|

ρl(i)=kρf

(8)

當(dāng)運(yùn)動(dòng)速度|v(i)|>Df時(shí),則:

(9)

其中:ρl為網(wǎng)格節(jié)點(diǎn)密度;ρf為水密度;md為液滴質(zhì)量;N為液滴數(shù)目;ΔV為網(wǎng)格體積;v為液滴運(yùn)動(dòng)速度;f為液滴入射頻率。

2.3 方程求解

大液滴的方程求解與小液滴的方程求解類似,主要求解步驟如下。

1) 計(jì)算x和y方向路程

對(duì)于x方向速度vx和y方向速度vy均為正方向,且xacc≥0、yacc≥0時(shí),則:

Δx(i)=dx-xacc(i)

Δy(i)=dy-yacc(i)

(10)

對(duì)于vx為正方向、vy為負(fù)方向,且xacc≥0、yacc≤0時(shí),則:

Δx(i)=dx-xacc(i)

Δy(i)=-dy+yacc(i)

(11)

其中,xacc、yacc分別為x方向和y方向的路程。

2) 計(jì)算運(yùn)動(dòng)時(shí)間

(1) 求解x方向方程和y方向方程得到tx和ty

對(duì)于vx和vy均為正方向,則:

(12)

其中,ax和ay分別為x和y方向的加速度。

(2) 判斷第i步計(jì)算的tx(i)與ty(i)是否為虛根,修正tx(i)和ty(i)

當(dāng)粒子運(yùn)動(dòng)到y(tǒng)方向最高點(diǎn)時(shí),可能會(huì)出現(xiàn)ty等于虛根的情況,原因在于粒子不能到達(dá)y正方向下一個(gè)網(wǎng)格節(jié)點(diǎn),而是經(jīng)歷了節(jié)點(diǎn)間的上升和下降過(guò)程。因此需要采用新的y方向求解方程,以保證ty為實(shí)根。對(duì)于vx為正方向、vy為正方向,且xacc(i)≥0、yacc(i)≥0時(shí),y方向產(chǎn)生了虛根,此時(shí)修正方程為式(13),得到的tx(i)和ty(i)均為實(shí)根。

Δx(i)=dx-xacc(i)

Δy(i)=-yacc(i)

(13)

(3) 比較tx和ty,并將其中較小值作為時(shí)間參數(shù)t

t(i)=min{tx(i),ty(i)}

(14)

3) 更新速度、累積量、密度和位移

時(shí)間參數(shù)t(i)后的vx和vy為:

vx(i+1)=vx(i)+ax(i)t(i)

vy(i+1)=vy(i)+ay(i)t(i)

(15)

對(duì)于時(shí)間參數(shù)t=ty(i),則:

xacc(i+1)=xacc(i)+vx(i)t(i)+

0.5ax(i)t(i)2

yacc(i+1)=0

(16)

對(duì)于時(shí)間參數(shù)t=tx(i),則:

yacc(i+1)=yacc(i)+vy(i)t(i)+

0.5ay(i)t(i)2

xacc(i+1)=0

(17)

4) 更新液滴覆蓋的節(jié)點(diǎn)密度和液滴覆蓋的節(jié)點(diǎn)速度

大液滴運(yùn)動(dòng)時(shí)可能會(huì)出現(xiàn)網(wǎng)格節(jié)點(diǎn)重疊的問(wèn)題,需要比較此步節(jié)點(diǎn)密度和原有節(jié)點(diǎn)密度,并取其中的較大值作為節(jié)點(diǎn)密度。其中,ρn為此步大液滴的節(jié)點(diǎn)密度,ρo為節(jié)點(diǎn)原密度。同時(shí),如果節(jié)點(diǎn)密度發(fā)生改變,節(jié)點(diǎn)速度也改變?yōu)閷?duì)應(yīng)的值。

ρl(i+1)=max(ρn(i+1),ρo)

(18)

3 模型應(yīng)用

3.1 模型驗(yàn)證與比較

為了展現(xiàn)模型在工程應(yīng)用上的意義,本文以AP1000蒸汽發(fā)生器的汽水分離系統(tǒng)為研究對(duì)象。AP1000汽水分離器的壓力p=5.76 MPa,外流場(chǎng)為純重力場(chǎng)且蒸汽流速v=1.45 m/s。噴射液滴的速度、半徑、高度信息由文獻(xiàn)[16]得到。

首先根據(jù)壓力確定臨界特征液滴半徑,液滴受到的力包括重力、浮力、曳力、附加質(zhì)量力、Magnus力、Saffman力[17]。圖4為不同尺寸噴射液滴的網(wǎng)格節(jié)點(diǎn)示意圖,其中黑色節(jié)點(diǎn)代表完全被液滴覆蓋的節(jié)點(diǎn),紅色節(jié)點(diǎn)代表不完全被液滴覆蓋的節(jié)點(diǎn)。由圖4可知,液滴尺寸越大,其覆蓋的黑色節(jié)點(diǎn)越多,且總的覆蓋節(jié)點(diǎn)也更多。液滴的受力達(dá)到平衡時(shí),即得到懸浮液滴的半徑r=1 356 μm,此時(shí)液滴為臨界尺寸液滴,臨界尺寸液滴減速為0且懸浮在某一位置。

液滴尺寸:a——臨界尺寸(r=1 356 μm);b——小于臨界尺寸(r=1 250 μm);c——大于臨界尺寸 (r=1 500 μm)圖4 網(wǎng)格節(jié)點(diǎn)示意圖Fig.4 Schematic diagram of grid node

分別對(duì)3種類型尺寸的液滴進(jìn)行了對(duì)比驗(yàn)證,兩種方法的速度相對(duì)誤差達(dá)到0.1%(圖5),密度相對(duì)誤差達(dá)到0.6%(圖6)。

液滴尺寸:a——臨界尺寸(r=1 356 μm);b——小于臨界尺寸(r=1 250 μm);c——大于臨界尺寸(r=1 500 μm)圖5 歐拉網(wǎng)格逼近方法和歐拉-拉格朗日方法的速度對(duì)比Fig.5 Comparison of velocity between Euler grid approximation and Euler-Lagrange methods

液滴尺寸:a——臨界尺寸(r=1 356 μm);b——小于臨界尺寸(r=1 250 μm);c——大于臨界尺寸(r=1 500 μm)圖6 歐拉網(wǎng)格逼近方法和歐拉-拉格朗日方法的密度對(duì)比Fig.6 Comparison of density between Euler grid approximation and Euler-Lagrange methods

在計(jì)算效率方面,采用傳統(tǒng)歐拉-拉格朗日方法在時(shí)間步長(zhǎng)為0.000 1 s時(shí)對(duì)于單液滴粒子的計(jì)算需要2 000 s,而采用歐拉網(wǎng)格逼近算法計(jì)算時(shí)間為10 s,能將計(jì)算效率提升1～2個(gè)數(shù)量級(jí)。

3.2 模型應(yīng)用

根據(jù)文獻(xiàn)[12]計(jì)算得到單位面積、單位時(shí)間產(chǎn)生的氣泡數(shù)目。下部重力空間的直徑為5 m,本文采用的網(wǎng)格尺寸為0.5 mm。根據(jù)二維截面得到氣泡產(chǎn)生總頻率f=250 s-1,即每隔20 mm的節(jié)點(diǎn)產(chǎn)生氣泡頻率f=1 s-1。假設(shè)每個(gè)氣泡能產(chǎn)生8個(gè)噴射液滴,即可得到每個(gè)噴射液滴產(chǎn)生的頻率。

圖7～9為不同尺寸噴射液滴局部圖,可看出,小于臨界尺寸的液滴會(huì)隨蒸汽流場(chǎng)向上流動(dòng),而大于臨界尺寸的液滴會(huì)先向上做減速運(yùn)動(dòng)到最高點(diǎn),然后下落回自由表面。

圖7 半徑1 250 μm的噴射液滴局部圖Fig.7 Local view of jet droplet with radius of 1 250 μm

圖8 半徑1 356 μm的噴射液滴局部圖Fig.8 Local view of jet droplet with radius of 1 356 μm

圖9 半徑1 500 μm的噴射液滴局部圖Fig.9 Local view of jet droplet with radius of 1 500 μm

計(jì)算得到的重力空間中不同尺寸噴射液滴速度場(chǎng)和密度場(chǎng)如圖10、11所示。

液滴半徑:a——1 250 μm;b——1 356 μm;c——1 500 μm 圖10 重力空間中的不同尺寸噴射液滴速度場(chǎng)Fig.10 Velocity field of jet droplet with different sizes in gravity space

計(jì)算得到的重力分離空間中的不同尺寸噴射液滴合并密度場(chǎng)如圖12、13所示,密度較大的位置是可能的碰撞點(diǎn)位置,實(shí)際重力分離空間的液滴總密度圖需要再疊加膜液滴的密度分布。

圖12 不同尺寸噴射液滴合并密度局部圖Fig.12 Local plot of combined density for jet droplet with different sizes

圖13 不同尺寸噴射液滴合并密度圖Fig.13 Combined density diagram for jet droplet with different sizes

4 結(jié)論

本文基于歐拉網(wǎng)格逼近理論基礎(chǔ)提出了一套用于大液滴計(jì)算的模型,采用歐拉-拉格朗日方法進(jìn)行了模型驗(yàn)證,兩種方法的速度相對(duì)誤差達(dá)到0.1%,密度相對(duì)誤差達(dá)到0.6%。同時(shí)歐拉網(wǎng)格逼近方法的計(jì)算效率相較于歐拉-拉格朗日方法能提升1～2個(gè)數(shù)量級(jí)。本文將該模型應(yīng)用于AP1000汽水分離器的下部重力分離空間中噴射液滴的模擬,分別得到達(dá)到懸浮的臨界尺寸噴射液滴、小于臨界尺寸和大于臨界尺寸噴射液滴的運(yùn)動(dòng)圖,且給出了合成密度圖,合成密度圖中密度較大的位置處是可能的碰撞點(diǎn)區(qū)域。本工作可為汽水分離等裝置中的大液滴分析提供有效的計(jì)算分析工具。