劉學(xué)文，李金京，全曉軍，熊偉

（上海交通大學(xué)機(jī)械與動力工程學(xué)院，工程熱物理研究所，上海200240）

引 言

深入理解納米或微米顆粒和泡沫之間的相互作用在很多領(lǐng)域有重要意義。在粉塵抑制領(lǐng)域，更好地理解粉塵顆粒與泡沫相互作用的基本原理是生產(chǎn)具有更高粉塵抑制性能的泡沫的關(guān)鍵[1]，對顆粒氣泡或顆粒液膜相互作用的研究是膠體化學(xué)領(lǐng)域的研究熱點(diǎn)[2-3]。此外，隨著納米流體在沸騰換熱中的應(yīng)用，納米顆粒在沸騰現(xiàn)象中越來越重要[4-5]，核態(tài)沸騰中納米顆粒會影響氣泡穩(wěn)定性進(jìn)而影響沸騰換熱，Quan等[6-7]和Binks等[8]關(guān)于納米流體的沸騰實(shí)驗(yàn)中發(fā)現(xiàn)中等親水性納米顆粒可以吸附在氣泡的氣液界面上，因此可以減少氣泡聚結(jié)從而提高沸騰換熱性能。

當(dāng)氣泡接觸另一個(gè)氣泡或液體表面時(shí)，會形成液體薄膜，泡沫或氣泡的穩(wěn)定性與這些液膜的穩(wěn)定性密切相關(guān)[9]。因此對單個(gè)納米顆粒與薄液膜之間相互作用的研究是深入理解顆粒和泡沫之間的相互作用的基礎(chǔ)，目前已經(jīng)有很多文獻(xiàn)從理論和實(shí)驗(yàn)方面對該現(xiàn)象進(jìn)行了研究。

眾多的研究表明[1,10-20]，對于兩側(cè)皆為氣體的單層液膜，球形親水性顆粒會增強(qiáng)液膜的穩(wěn)定性，而球形疏水性顆粒會促進(jìn)液膜的破裂。然而，目前大多數(shù)研究注重親水顆粒對液膜穩(wěn)定性增強(qiáng)的研究，針對疏水顆粒導(dǎo)致液膜破裂的研究極少。Frye 等[21]提出了一種針對單個(gè)顆粒與單層薄膜相互作用的理論模型，發(fā)現(xiàn)球形顆粒的接觸角需要大于90°才能促進(jìn)膜破裂，但是該模型中疏水顆粒初始時(shí)刻便位于液膜中間，這顯然與實(shí)際不符。此外也有關(guān)于該現(xiàn)象的實(shí)驗(yàn)研究，Dippenaar[9]使用疏水玻璃球形微米顆粒對單個(gè)顆粒與單層薄液膜之間的相互作用進(jìn)行了實(shí)驗(yàn)研究，結(jié)果表明接觸角為102°的微米顆粒會促進(jìn)液膜破裂，在一定程度上證實(shí)了理論預(yù)測的液膜破裂機(jī)制。Zhu 等[1]則在納米尺度液膜上研究了單個(gè)球形玻璃顆粒與單層表面活性劑薄膜的相互作用，發(fā)現(xiàn)薄膜的壽命在加入親水性顆粒后延長，在添加疏水性顆粒時(shí)減小。這些研究中，液膜厚度都遠(yuǎn)小于顆粒半徑。對于理論模型，難以研究顆粒接觸液膜并浸入液膜的過程，而實(shí)驗(yàn)研究中難以對顆粒接觸角進(jìn)行精確控制。

對于單個(gè)納米或微米顆粒促進(jìn)液膜破裂的現(xiàn)象，其時(shí)間空間尺度小，實(shí)驗(yàn)難以準(zhǔn)確捕捉，單純的理論與實(shí)驗(yàn)研究都不能直觀展示顆粒與液膜相互作用的詳細(xì)過程，在這方面模擬研究有其獨(dú)特的優(yōu)勢。作為介觀方法，格子Boltzmann方法（LBM）已被證明是研究潤濕性對流體中固體顆粒運(yùn)動影響的有效數(shù)值方法。Ladd[22-23]提出用于模擬固液懸浮液的2D 格子Boltzmann 模型，Joshi 等[24]改進(jìn)Ladd 的模型[22-23]并與單組分/多相模型相結(jié)合，模擬懸浮在液體蒸氣表面上的2D顆粒的運(yùn)動。

本文通過將Joshi 等[24]的顆粒模型與Gong 等[25]的單組分/兩相流模型結(jié)合，采用等溫格子Boltzmann 方法兩相流模型對單個(gè)顆粒與薄液膜之間的相互作用進(jìn)行了研究，并對疏水顆粒促進(jìn)液膜破裂的過程進(jìn)行了詳細(xì)研究。

1 格子Boltzmann方法

1.1 格子Boltzmann方法兩相流模型

本文采用Gong 等[25]提出的Gong-Cheng 單組分/兩相流模型。

格子Boltzmann 方法主要分為碰撞和遷移兩個(gè)過程，采用BGK 碰撞算子的流體粒子分布函數(shù)的碰撞過程為：

結(jié)合顆粒運(yùn)動模型之后其遷移過程有所不同，具體遷移過程的方程可參考文獻(xiàn)[26]，關(guān)于該二維兩相流模型各項(xiàng)的具體設(shè)置可參考文獻(xiàn)[25]。

1.2 顆粒運(yùn)動模型

本文顆粒運(yùn)動模型基于三個(gè)假設(shè)：(1)顆粒為球形，(2)忽略顆粒轉(zhuǎn)動，(3)忽略布朗運(yùn)動。根據(jù)牛頓第二定律，每個(gè)顆粒的運(yùn)動可以通過以下方程式描述：

式中，Up是粒子的速度，Mp是顆粒的質(zhì)量，二維情況下由式(3)給出：

式中,Rp是粒子的半徑，ρp是粒子的密度，此處取為液體密度的1.5倍。

式（2）中的Fp是作用在粒子上的總作用力，由于本文只研究單個(gè)顆粒的行為且忽略重力的作用，因此總作用力包括來自周圍流體的撞擊力Fimp，附著力Fadh和動量交換力Fmom三部分，可以表示為：

其中Fimp可參考文獻(xiàn)[26]，F(xiàn)adh是流體節(jié)點(diǎn)與粒子節(jié)點(diǎn)的相互作用力，表達(dá)式為[24]：

式中，b表示顆粒邊界相鄰的所有流體格點(diǎn)；Gsp是流固作用力強(qiáng)度系數(shù)，可用于調(diào)控顆粒表面潤濕性；sp(x+eiδt)為指示函數(shù)，顆粒內(nèi)部固體格點(diǎn)取1，流體格點(diǎn)取0。關(guān)于LBM 方法中固液接觸角的更多計(jì)算可以參考文獻(xiàn)[27]。關(guān)于Fmom的計(jì)算可參考文獻(xiàn)[24]。

Gong-Cheng 模型相比普通單組分/多相模型具有更好的數(shù)值精度和數(shù)值穩(wěn)定性[25,28]，將Gong-Cheng 單組分/兩相流模型與顆粒運(yùn)動模型相結(jié)合，固體顆粒單個(gè)時(shí)間步Δt（Δt=1）的位移為UpΔt，通過顆粒內(nèi)部格點(diǎn)和流體格點(diǎn)的相互轉(zhuǎn)化來實(shí)現(xiàn)顆粒位移。顆粒移動后，在t時(shí)刻由流體格點(diǎn)轉(zhuǎn)化為顆粒內(nèi)部格點(diǎn)的格點(diǎn)密度可直接取為ρp，而由固體顆粒內(nèi)部格點(diǎn)轉(zhuǎn)化為新流體格點(diǎn)的格點(diǎn)密度由周圍流體格點(diǎn)的平均密度來表示[26]：

式中，xa表示新流體格點(diǎn)周圍的流體點(diǎn)，Na表示周圍流體點(diǎn)的數(shù)量，ρa(bǔ)表示xa位置處流體格點(diǎn)的密度。由式（6）即可得到新的流體格點(diǎn)的密度。

關(guān)于顆粒運(yùn)動模型的更多細(xì)節(jié)可參考文獻(xiàn)[24]。

2 顆粒親疏水性

本節(jié)討論在等溫條件下單個(gè)顆粒漂浮在水面上的穩(wěn)態(tài)情況，格子Boltzmann 方法中顆粒在水蒸氣-水界面處的平衡位置取決于顆粒由于固液作用力引起的潤濕性[29]，LBM 中通過調(diào)節(jié)系數(shù)Gsp可得到不同顆粒接觸角?，F(xiàn)在通過格子Boltzmann 方法(LBM)模擬單個(gè)粒子在無重力的情況下漂浮在液體表面的情況。圖1(a)給出了格子單位為Nx×Ny =400 × 400 的計(jì)算域的2D 示意圖，通過網(wǎng)格無關(guān)性測試[圖1(b)]發(fā)現(xiàn)這一格子計(jì)算區(qū)域可以保證計(jì)算精度，高度為Ny/2的液態(tài)水充滿下方區(qū)域，固體顆粒（半徑Rp=20 個(gè)格子單位，實(shí)際單位約為30.6 nm）最初放置在氣-液界面的中心，液體和水蒸氣都處于飽和溫度T0=0.85Tcr，在此溫度下，液體和水蒸氣的飽和密度分別約為6.63 和0.34。在x方向上施加周期性邊界條件，并且在計(jì)算區(qū)域的頂部（y=400）和底部（y= 0）施加無滑移壁面邊界條件。關(guān)于格子單位與實(shí)際物理單位轉(zhuǎn)換可參考Wang 等[30]的方法，本文中1格子單位約對應(yīng)1.53 nm。

接觸角大于90°的疏水顆粒與接觸角小于90°的親水顆粒在液面不同位置如圖1 中所示，作用于氣-液界面處半徑為Rp的球形顆粒上的界面力（Fi）可表示為[11]：

圖1 液面上漂浮單個(gè)粒子的計(jì)算域示意圖(a);不同網(wǎng)格數(shù)下Gsp與顆粒接觸角關(guān)系圖(b)Fig.1 Schematic of the computation domain for a single particle floating on liquid surface (a);Relationship between Gsp and contact angle with different grid numbers(b)

式中，σ為液體的表面張力，J/m2；θ為固體顆粒的接觸角，(°)；x為顆粒浸入液體中的相對深度，等于其絕對浸入深度h與顆粒半徑之比：x=h/Rp。如果Fi> 0，則向量Fi指向液相。在沒有其他力的情況下，如果Fi= 0，則顆粒將在氣-液界面處于平衡位置。此時(shí)可得顆粒平衡位置浸入液體深度為：

如式（8）所示，接觸角唯一地確定了顆粒在液-氣界面處的平衡位置。

3 結(jié)果與討論

3.1 計(jì)算域設(shè)置

3.2 顆粒接觸角對液膜破裂的影響

圖2 單個(gè)固體顆粒與薄液膜相互作用的計(jì)算域示意圖Fig.2 Schematic diagram of the computational domain simulating the interaction between a single solid particle and thin liquid film

式中，R1和R2是描述顆粒附近彎曲的氣-液界面的兩個(gè)曲率半徑，σ是水蒸氣-水的表面張力。上下的氣液界面會在該毛細(xì)力驅(qū)動下會沿著顆粒表面移動，當(dāng)上下的氣液界面接觸時(shí)，顆粒與液膜接觸極不穩(wěn)定，液膜便會脫離顆粒表面而破裂[21]。

圖3 液面彎曲產(chǎn)生毛細(xì)力示意圖Fig.3 Schematic diagram of capillary force generated by liquid surface bending

模擬發(fā)現(xiàn)顆粒導(dǎo)致液膜破裂的整個(gè)過程如圖4所示，液膜破裂現(xiàn)象與理論分析[32]和實(shí)驗(yàn)現(xiàn)象[9]基本一致。模擬發(fā)現(xiàn)疏水納米顆粒導(dǎo)致液膜的破裂可以細(xì)分為兩個(gè)階段，階段1 是顆粒接觸液膜并在界面力（Fi）的作用下向液膜內(nèi)部移動直至顆粒恰好接觸下氣液界面的過程，階段2 是上下氣液界面被毛細(xì)力驅(qū)動在顆粒表面移動的過程。

顆粒接觸角越大，液膜在顆粒表面的彎曲程度越大，所產(chǎn)生的毛細(xì)力越大，液面線在顆粒表面移動的速度越大，上下氣液界面會更快相遇從而導(dǎo)致液膜破裂，因此階段2 時(shí)間大大縮短。但是顆粒接觸角越大，由式（7）可知使顆粒向液體內(nèi)部運(yùn)動的界面力Fi越小，因此顆粒從接觸液膜到接觸液膜下界面的移動速度減小，而這卻會導(dǎo)致階段1 時(shí)間大大延長，階段1 所耗費(fèi)的時(shí)間遠(yuǎn)長于階段2。圖5 是不同接觸角顆粒導(dǎo)致定厚度液膜破裂兩個(gè)階段的時(shí)間分布模擬結(jié)果，其中破裂時(shí)間為顆粒接觸液膜直到液膜恰好脫離顆粒表面的時(shí)間間隔，圖5表明，在顆粒接觸角為106.7°附近時(shí)兩個(gè)階段綜合耗費(fèi)的時(shí)間最短，因此在此接觸角附近顆粒導(dǎo)致液膜破裂時(shí)間最短。

圖4 球形疏水顆粒導(dǎo)致液膜破裂的兩個(gè)階段（接觸角θ=106.7°，液膜厚度=0.4）Fig.4 Two stages of liquid film rupture caused by spherical hydrophobic particle(contact angle θ=106.7°,liquid film thickness=0.40)

圖5 不同接觸角顆粒導(dǎo)致定厚度液膜破裂的兩個(gè)階段的時(shí)間分布(液膜厚度 =0.4)Fig.5 Time distribution of two stages of liquid film rupture caused by particle with different contact angle at a given thickness(liquid film thickness=0.4)

3.3 液膜厚度的影響

當(dāng)顆粒接觸液膜時(shí)，顆粒會在界面力的作用下向液膜內(nèi)部移動，此時(shí)會有兩種情況，一種是液膜的厚度小于顆粒在氣液界面平衡時(shí)進(jìn)入液體的深度，另一種是大于。

當(dāng)?shù)谝环N情況發(fā)生時(shí)，顆粒會在界面力的作用下直接穿透液膜并接觸液膜下界面，然后開始圖4所示的液膜破裂的第二個(gè)階段。

第二種情況多發(fā)生于顆粒接觸角過大即過于疏水的情況，此時(shí)界面力的作用不足以使顆粒穿透液膜并接觸下界面，模擬結(jié)果如圖6所示，此時(shí)顆粒下部液體不斷向兩側(cè)排開，隨著液體不斷向兩側(cè)排開下液面則會不斷向上移動，直到顆粒下部液面線接觸顆粒，此后的現(xiàn)象與第一種情況相同。在這種情況下，前期顆粒排開液體向兩側(cè)運(yùn)動的過程極為緩慢，因此顆粒導(dǎo)致液膜破裂的時(shí)間大大延長。

圖6 展示的是當(dāng)顆粒半徑網(wǎng)格數(shù)為20、30 和40時(shí)導(dǎo)致液膜排水的示意圖，液膜厚度均為=0.45，可見隨顆粒半徑網(wǎng)格數(shù)增加依然會出現(xiàn)排水現(xiàn)象，因此可以排除網(wǎng)格數(shù)量的影響。

圖6 疏水顆粒導(dǎo)致液膜排水示意圖（接觸角θ=150.5°，液膜厚度=0.45）Fig.6 Schematic diagram of liquid film drainage caused by particle(contact angle θ=150.5°,liquid film thickness=0.45)

圖7是不同接觸角顆粒導(dǎo)致不同厚度的液膜的破裂時(shí)間的模擬結(jié)果。對于接觸角為90°和106°的顆粒而言，其平衡時(shí)進(jìn)入液體的深度極大，并不會產(chǎn)生圖6 的排水現(xiàn)象，因此液膜破裂時(shí)間隨液膜厚度近似線性增長。對于接觸角為125.7°的顆粒，當(dāng)液膜厚度大于等于0.45 時(shí)，顆粒導(dǎo)致液膜破裂的時(shí)間大大延長，這是因?yàn)榻佑|角為125.7°的顆粒平衡時(shí)能夠浸入液體的深度為=1+cosθ=0.42，當(dāng)液膜厚度大于0.42 時(shí)，顆粒僅依靠界面力無法穿透液膜，達(dá)到穩(wěn)定的深度heq后(圖6)，顆粒下部的液體不斷向兩側(cè)排開，下方液相線逐漸靠近顆粒，直到接觸顆粒表面，這個(gè)過程遠(yuǎn)遠(yuǎn)慢于顆粒依靠界面力運(yùn)動的過程，因此耗費(fèi)的時(shí)間呈現(xiàn)突然的增長。而對于接觸角為150.5°和161.2°的顆粒而言，其能深入液膜的深度極小，因此在液膜厚度為0.13 的時(shí)候便已經(jīng)是圖6破裂機(jī)理，因此隨液膜厚度增加，其導(dǎo)致液膜破裂的時(shí)間呈近似線性增長，并不會因?yàn)槠屏涯Ｊ降霓D(zhuǎn)變而導(dǎo)致破裂時(shí)間的突變。當(dāng)液膜厚度為0.5 時(shí)，接觸角為125.7°、150.5°以及161.2°的顆粒無法依靠界面力浸入液膜，因此會對液膜產(chǎn)生圖6的排水現(xiàn)象，此時(shí)接觸角越大排水越快，因此接觸角161.2°的顆粒反而比125.7°和150.5°的顆粒更快導(dǎo)致液膜破裂。從圖7 可以看出，當(dāng)液膜厚度小于顆粒半徑時(shí)，接觸角為106.7°的疏水顆粒導(dǎo)致液膜破裂的時(shí)間總是最短。

圖7 液膜破裂時(shí)間與液膜厚度及顆粒接觸角關(guān)系Fig.7 Time distribution of liquid film rupture of different thickness caused by particle with different contact angles

4 結(jié) 論

采用等溫格子Boltzmann 方法兩相流模型結(jié)合顆粒運(yùn)動模型對單個(gè)顆粒導(dǎo)致液膜破裂的現(xiàn)象進(jìn)行了數(shù)值模擬研究，得到如下結(jié)論。

（1）疏水性球形顆粒會導(dǎo)致液膜破裂，顆粒促進(jìn)液膜破裂可以細(xì)分為兩個(gè)階段，一是顆粒接觸液膜并在界面力的作用下向液膜內(nèi)部移動直至顆粒恰好接觸液膜下界面的階段，二是上下氣液界面在毛細(xì)力驅(qū)動下沿顆粒表面移動并接觸的階段。綜合兩個(gè)階段的時(shí)間，發(fā)現(xiàn)顆粒接觸角為106.7°附近時(shí)導(dǎo)致液膜破裂時(shí)間最短。

（2）液膜厚度與顆粒接觸角的相對關(guān)系對液膜破裂時(shí)間有很大影響，當(dāng)液膜的厚度大于顆粒在氣-液界面平衡時(shí)浸入液體的深度時(shí)，此時(shí)界面力的作用不足以使顆粒穿透液膜并接觸液膜下界面，液膜將會發(fā)生排水現(xiàn)象直至顆粒接觸液膜下氣液界面，在這種情況下，疏水顆粒導(dǎo)致液膜破裂的時(shí)間會延長。當(dāng)不同接觸角下均處于排水狀態(tài)時(shí)，接觸角越大則液膜排水越快。

由于缺少通用的理論分析與詳細(xì)的實(shí)驗(yàn)數(shù)據(jù)，僅能對該現(xiàn)象進(jìn)行模擬分析，模擬現(xiàn)象與理論和實(shí)驗(yàn)基本一致，但缺少定量對比，這將是下一步的研究重點(diǎn)。