基于格子玻爾茲曼方法的液滴撞擊液膜數(shù)值模擬

魯舟洋，邢 巖

（1.中國(guó)電建集團(tuán)西北勘測(cè)設(shè)計(jì)研究院有限公司，陜西 西安 710065；2.交通運(yùn)輸部天津水運(yùn)工程科學(xué)研究所工程泥沙交通運(yùn)輸行業(yè)重點(diǎn)實(shí)驗(yàn)室，天津 300456）

0 引言

液滴撞擊液膜現(xiàn)象廣泛存在于自然現(xiàn)象和工業(yè)制造中，該現(xiàn)象發(fā)生往往伴隨著復(fù)雜的傳熱傳質(zhì)過(guò)程。一系列研究表明，液滴撞擊液膜過(guò)程中液冠演化過(guò)程受到眾多參數(shù)的影響，這些參數(shù)通常包括雷諾數(shù)（Re），韋伯?dāng)?shù)（We），氣液密度比，氣液粘滯系數(shù)比，液膜高度等。Rioboo[1]等通過(guò)實(shí)驗(yàn)研究了液滴半徑對(duì)液冠演化的影響，并將該過(guò)程分為六個(gè)階段。Rein[2]、Yarin[3]等系列研究表明，液體飛濺過(guò)程中，會(huì)發(fā)生液滴從液冠脫落，即液滴飛濺現(xiàn)象。隨后Wang 等[4]研究了不同Re 數(shù)、We 數(shù)及液面厚度對(duì)液冠演化的影響。

隨著計(jì)算機(jī)技術(shù)的發(fā)展，計(jì)算流體力學(xué)研究結(jié)果可以有效展示復(fù)雜流場(chǎng)的演化過(guò)程，逐漸成為研究流體運(yùn)動(dòng)的重要手段。格子玻爾茲曼方法（LBM）作為介觀數(shù)值模擬方法，能夠有效模擬氣液交界面變化復(fù)雜的流動(dòng)現(xiàn)象。經(jīng)過(guò)近30 a 的發(fā)展，現(xiàn)有LBM多相流模型可以歸納為偽勢(shì)模型、顏色模型、自由能模型和相場(chǎng)模型[5]。Lee[6]等人首先利用相場(chǎng)模型模擬、不同Re條件下，液滴撞擊液膜的過(guò)程。Wang[7]等則利用基于有限體積法的LBM自由能模型模擬Re 從20～1000 條件下液滴飛濺的過(guò)程。黃虎[8]等則利用LBM相場(chǎng)模型模擬We=500 和Re=24～480條件下的液滴撞擊液膜的過(guò)程。以上研究均表明：不同LBM可以有效模擬液滴撞擊液面的過(guò)程。

總結(jié)前人研究發(fā)現(xiàn)，基于LBM 的數(shù)值模擬研究中，尚未利用偽勢(shì)LBM 模型對(duì)氣液密度比超過(guò)500 的液滴撞擊液膜現(xiàn)象進(jìn)行模擬。因此本文利用大密度比LBM 偽勢(shì)模型，結(jié)合可調(diào)節(jié)表面張力的外力模型，模擬兩相密度比為ρL/ρg，不同Re 數(shù)和We 數(shù)下液冠演化過(guò)程，分析這些參數(shù)對(duì)液冠半徑和液冠高度的影響，并進(jìn)一步研究液膜厚度對(duì)液滴飛濺過(guò)程的影響。

1 大密度比偽勢(shì)LBM 方法

本文采用的多松弛碰撞算子LBM偽勢(shì)模型，對(duì)應(yīng)粒子分布方程為：

式中：ε 用于調(diào)節(jié)熱力學(xué)一致性。其中ρ 為宏觀密度，ux，uy為宏觀速度。

其中式（5）為對(duì)應(yīng)宏觀量求解方式。F 為作用于流體上作用力，包括粒子間作用力Fm等，其中粒子間作用力：

pEOS為通過(guò)流體狀態(tài)方程求得壓強(qiáng)。本文采用Carnahan-Starling(C-S)狀態(tài)方程：

式中：a=0.4963R2Tc2/pc，b=0.1873RTc/pc,其中Tc，pc分別為該氣體的臨界溫度與臨界壓強(qiáng)。本研究中狀態(tài)方程中相關(guān)參數(shù)為a=1，b=4，Rg=1。同時(shí)本文引入Li提出表面張力調(diào)節(jié)項(xiàng)，可以在偽勢(shì)模型中不通過(guò)調(diào)節(jié)密度比達(dá)到調(diào)節(jié)表面張力的目的：

式中：k 為調(diào)節(jié)表面張力系數(shù)，取值在0～1 之間。

2 模擬結(jié)果與分析

由于液滴撞擊液膜受到液滴撞擊速度，液相粘滯系數(shù)，液體表面張力等多種因素的影響，為更系統(tǒng)的研究概述這些因素的影響，本研究依舊沿襲前人的分析方法，引入多個(gè)無(wú)量綱參數(shù)分析液面的演化過(guò)程。引入無(wú)量綱參數(shù)包括雷諾數(shù)（Re）、韋伯?dāng)?shù)（We）、以及相對(duì)厚度和無(wú)量綱時(shí)間。

式中：vL為液體粘滯系數(shù)，U 為液滴撞擊速度。

式中：σ 為表面張力。

式中：H2為液膜厚度。

2.1 初始條件

為驗(yàn)證模型的準(zhǔn)確性，本文采用1001×351 矩形計(jì)算域，計(jì)算域底部固壁采用標(biāo)準(zhǔn)反彈邊界，左右均采用周期邊界，而計(jì)算域頂部采用非平衡外推邊界，計(jì)算區(qū)域見(jiàn)圖1。其中液滴初始半徑為R0，液滴初始垂向速度為U，液滴中心距液膜高度為H1，液膜高度為H2，計(jì)算中各工況選取參數(shù)見(jiàn)表1。

圖1 液滴撞擊液面計(jì)算域示意圖

表1 液滴撞擊液面計(jì)算工況

本研究中液滴撞擊速度設(shè)為U=0.125，通過(guò)調(diào)節(jié)松弛系數(shù)τv和k 來(lái)獲得不同的Re 數(shù)和We 數(shù)。為保證熱力學(xué)一致性，式（4）中ε 取值為0.11。同時(shí)為保證計(jì)算時(shí)數(shù)值計(jì)算穩(wěn)定性，減小初始界面附近密度梯度，液滴周圍密度場(chǎng)初始化見(jiàn)式（14）。

其中（x0，y0）為液滴中心，初始界面厚度w=5。而液膜附近密度初始化公式為：

2.2 不同Re 對(duì)液面演化的影響

本研究在保證撞擊速度不變的前提下，通過(guò)改變液體粘滯系數(shù)，獲得四組不同的Re 數(shù)，分別為Re=40、200、500、1000，用來(lái)觀察不同Re 數(shù)對(duì)液冠形態(tài)的影響。圖2 展示了不同Re 對(duì)液冠形態(tài)的變化的影響，圖中可以看到Re=40 時(shí)，并未形成有效液冠，但其他三個(gè)工況形成液冠均非常明顯。從圖中可以看到當(dāng)t*=1.5 時(shí)刻，不同Re 條件下數(shù)液冠均未發(fā)生液滴飛濺現(xiàn)象。但當(dāng)t*=2.5 時(shí)，而當(dāng)Re 數(shù)增大到1000 后，液體尖端受到Rayleigh-Plateau 不穩(wěn)定性的影響，從而發(fā)生液滴飛濺。

圖2 不同Re 條件下液冠形態(tài)演化過(guò)程

圖3 展示了不同Re 條件下，無(wú)量綱液冠半徑與無(wú)量綱時(shí)間之間的關(guān)系。從圖中可以看出，不同Re 數(shù)條件下r/2R0和均成正比關(guān)系，與理論分析一致。同時(shí)斜率k 隨著Re 數(shù)的增大而增大。同時(shí)在圖上我們可以看出，當(dāng)Re=40 時(shí)斜率與Re=200、500、1000 三個(gè)工況斜率相差較大，這是由于Re=40 時(shí)受到液體粘滯系數(shù)的影響，并未形成有效液冠。

圖3 不同Re 數(shù)液冠半徑r/2R0 與無(wú)量綱時(shí)間之間關(guān)系

2.3 不同We 對(duì)液面演化的影響

液滴撞擊液膜過(guò)程中，液冠飛濺效應(yīng)會(huì)隨著We 數(shù)增大而變強(qiáng)。圖4 展示了Re=500，H*=2.5 條件下不同We 數(shù)條件下，液冠在t*=0.625 和t*=2.5 時(shí)流場(chǎng)形態(tài)。當(dāng)t*=0.625 時(shí)，不同We條件下，液滴兩側(cè)均形成液冠水柱，可以看到液冠水柱厚度隨著We 的增大而減小。而當(dāng)t*=2.5 時(shí)，We=87.8、139.4、1165.5 三種條件下均能觀測(cè)到液滴飛濺的過(guò)程，但飛濺液滴大小隨著We 的增大而顯著減小，液冠厚度也減小。原因在于表面張力減小，導(dǎo)致液柱受到Rayleigh-Plateau 不穩(wěn)定性更強(qiáng)，液冠難以維持穩(wěn)定形態(tài)，末端液滴更容易發(fā)生脫落。

圖4 不同We 條件下液冠形態(tài)變化

圖5 展示了不同We 條件下，無(wú)量綱液冠半徑與無(wú)量綱時(shí)間之間的關(guān)系。本研究表明，在We=69.0、87.8、139.4 三種工況下，液冠半徑r/2R0和均成正比關(guān)系，與理論分析一致，但當(dāng)We 增加到1165.5 時(shí)，液冠半徑與無(wú)量綱時(shí)間并不成線性關(guān)系，這是由于液冠末端液滴脫落過(guò)早，液滴不受液冠影響而導(dǎo)致的。圖中可以得到在We 由69.0 增加到139.4 過(guò)程中對(duì)斜率k 值影響并不大。

圖5 不同We 數(shù)液冠半徑r/2R0 與無(wú)量綱時(shí)間之間關(guān)系

2.4 液膜厚度對(duì)液滴飛濺影響

液膜厚度對(duì)液冠形狀的演化具有重要意義，Wang 實(shí)驗(yàn)研究表明[4]，隨著液膜厚度減小，液滴撞擊液膜過(guò)程中更容易發(fā)生散裂。圖6 展示了Re=500 時(shí)六種不同液膜高度條件下，t*=1.875 時(shí)液冠形狀。液冠演化過(guò)程中，液冠液柱與液膜之間的夾角隨著液膜厚度的增加而增大，這與Wang 實(shí)驗(yàn)一致。

圖6 不同H*條件下液冠形態(tài)變化

不同H*條件下，液冠半徑與無(wú)量綱時(shí)間之間的關(guān)系見(jiàn)圖7。本研究表明，在不同液膜厚度條件下，液冠半徑r/2R0和均成正比關(guān)系，與理論分析一致。不同液膜厚度H*也會(huì)影響到液冠高度。圖8 展示了不同H*條件下t*=2.5 時(shí)液冠高度的變化，當(dāng)t*=2.5 時(shí)，H*由0.1 增加到0.5 過(guò)程中液冠高度隨著H*的增加而增加，原因在于液滴撞擊過(guò)程中液柱與液膜之間的夾角減小了，導(dǎo)致水平方向液冠速度增加，而垂直方向上速度減小，液冠半徑增長(zhǎng)速度也隨之增大，但垂直增長(zhǎng)速度隨之減小。

圖7 不同液膜厚度H*條件下液冠半徑r/2R0 與無(wú)量綱時(shí)間之間關(guān)系

圖8 不同液膜厚度t*=2.5 時(shí)液冠高度的變化

3 結(jié)論與展望

本文采用大密度比LBM偽勢(shì)模型，結(jié)合可調(diào)節(jié)表面張力的外力項(xiàng)，模擬液滴撞擊液膜的過(guò)程，并分析不同Re、We 和液膜高度H*條件下，液冠的演化過(guò)程，得到結(jié)論如下：

1）在相同的液膜高度，高Re 數(shù)或高We 數(shù)條件下，液冠末端會(huì)出現(xiàn)液滴脫離現(xiàn)象，大密度比LBM偽勢(shì)模型可以有效模擬出此類現(xiàn)象。

2）隨著液相粘滯系數(shù)增大或液膜厚度H*增加，液冠末端變得越不容易散裂。

3）液滴撞擊形成液冠半徑與無(wú)量綱時(shí)間存在線性關(guān)系，其斜率隨著Re、We 的增大而增大，隨著液膜厚度H*的減小而增大。