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

魯舟洋，邢 巖

（1.中國電建集團西北勘測設計研究院有限公司，陜西 西安 710065；2.交通運輸部天津水運工程科學研究所工程泥沙交通運輸行業(yè)重點實驗室，天津 300456）

0 引言

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

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

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

1 大密度比偽勢LBM 方法

本文采用的多松弛碰撞算子LBM偽勢模型，對應粒子分布方程為：

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

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

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

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

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

2 模擬結果與分析

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

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

式中：σ 為表面張力。

式中：H2為液膜厚度。

2.1 初始條件

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

圖1 液滴撞擊液面計算域示意圖

表1 液滴撞擊液面計算工況

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

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

2.2 不同Re 對液面演化的影響

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

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

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

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

2.3 不同We 對液面演化的影響

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

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

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

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

2.4 液膜厚度對液滴飛濺影響

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

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

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

圖7 不同液膜厚度H*條件下液冠半徑r/2R0 與無量綱時間之間關系

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

3 結論與展望

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

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

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

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