微通道疏水表面滑移的耗散粒子動力學研究*

許少鋒 樓應侯 吳堯鋒 王向垟 何平

(浙江大學寧波理工學院,寧波 315000)

1 引 言

近年來,疏水表面在流動減阻、自清潔、防污等方面顯示了廣泛的應用前景[1,2].普遍認為,流體流經疏水表面時會發生滑移[3-5],尤其是在微觀尺度下,滑移更加明顯.經典流體力學中的無滑移假設適用于多數的宏觀黏性流體流動問題,但對于疏水表面,特別是在微納尺度下研究流體流動問題時,微小的滑移也會給流體的輸運產生重要的影響,需要利用滑移邊界才能準確獲得流體的流動行為.

隨著微納機電系統及微流控芯片的快速發展,疏水表面的滑移研究受到越來越多學者的關注.疏水表面的滑移已得到大量實驗研究的支持,如Choi等[6]以及Lee和Kim[7]利用流變儀研究了流體在疏水表面的滑移情況,測得在剪切速率為105sˉ1時的滑移長度接近30 nm；Tretheway 和Meinhart[8]采用粒子圖像測速技術(particle image velocimetry,PIV)測得含疏水涂層微管壁面的滑移長度為0.92 μm；Bhushan等[9]用原子力顯微鏡測得疏水和超疏水表面的滑移長度分別為43 和236 nm；Pit等[10]利用他們自己設計的實驗裝置直接觀察到了疏水表面的滑移.盡管各種新的實驗技術被用于疏水表面的滑移研究,但由于實驗材料各異,實驗條件也難以控制,而且微觀尺度下的實驗對精度要求苛刻,因此實驗研究結果還難以揭示疏水表面的滑移規律及其機理.而計算機數值模擬可以直接研究系統的微觀細節,能考慮到實驗研究和理論分析所忽略的細節,這使得數值模擬在研究微尺度流動方面具有明顯的優勢.很多學者采用分子動力學方法、格子Boltzmann方法與耗散粒子動力學方法等數值模擬技術開展疏水表面的滑移研究,試圖揭示其滑移規律(如滑移長度與接觸角之間的關系)及其機理.曹炳陽等[11]采用分子動力學方法模擬了液態氬在鉑納米通道內的流動,通過改變壁面與流體之間的勢能作用實現壁面的不同疏水性能,研究發現在疏水表面附近形成了一個低密度層,低密度層使流體在壁面發生了滑移.一般認為壁面接觸角越大,滑移長度越大,但Voronov等[12,13]的分子動力學模擬得到了不同的結論,他們模擬了類石墨烯平板間LJ (Lennard-Jones)流體的Couette流,發現接觸角在140°左右時滑移長度最大,并不是接觸角越大滑移長度越大.黃橋高等[14]利用格子Boltzmann方法研究了疏水表面微通道內的流體流動,通過改變吸附參數(反映壁面對流體粒子的作用強度)得到疏水表面,結果表明疏水作用在近壁區誘導了一個低密度層,表觀滑移則發生在低密度層,并發現滑移長度與吸附參數間存在近似的指數函數關系.Zhang等[15]也采用格子Boltzmann方法得到了與文獻[14]類似的結果,滑移長度與接觸角之間滿足指數函數關系.Cupelli等[16]利用耗散粒子動力學研究了毛細管的潤濕行為,改變固體壁面與流體粒子的相互作用參數得到了不同潤濕性能的壁面,模擬結果顯示接觸角與固液相互作用參數呈近似線性關系.關于疏水表面滑移機理,Tretheway和Meinhart[17]提出的理論認為在疏水表面存在不連續的氣體層或納米氣泡,這些不連續的氣體區域是表面產生滑移的原因.

從以上研究可以看到,疏水表面會產生滑移已獲得一致認同,但對其滑移機理以及滑移規律還沒有得到一致的結論,有待進一步研究.本文采用耗散粒子動力學(DPD)方法模擬微通道Couette流動,通過改變壁面粒子與流體粒子之間排斥力系數獲得不同疏水性能的壁面,探討疏水表面的滑移問題,以期揭示疏水表面的滑移規律及其機理.

2 DPD方法與壁面邊界模型

2.1 DPD方法

DPD是由Hoogerbrugge和Koelman[18]于1992年提出的一種介觀尺度的模擬方法,后來Espanol和Warren[19]以及Marsh[20]建立了DPD統計力學理論模型,奠定了其理論基礎,Groot和Warren[21]的研究工作進一步推廣了DPD的應用.DPD模型中,每個粒子是由一團分子組成的,是一種粗粒化的粒子,粒子之間采用軟勢能作用,因此DPD模型模擬的時間和空間尺度比分子動力學大得多.DPD粒子之間具有保守力、耗散力和隨機力三種形式的作用力,三種作用力都是成對出現的,可以確保體系動量守恒,十分適合模擬流體動力學問題.DPD方法已被廣泛應用于液滴、多相流體、高分子溶液、血液等復雜流體的動力學模擬[22-24].

DPD方法中,粒子的運動規律通過牛頓第二定律得到.具體地說,對i粒子,其運動方程為

式中ri,vi分別為粒子的位置矢量和速度矢量.注意(1)式中,每個粒子的質量相等,并且質量采用無量綱單位1.Fext為施加的外場力.Fij是粒子j對粒子i的作用力,由保守力,耗散力,隨機力三部分組成：

粒子i和粒子j之間只在截斷半徑rc之內具有相互作用力,兩粒子之間的距離大于rc時則沒有相互作用力,三種作用力均沿粒子之間連線方向,都在rc范圍內才有相互作用.

式中aij為粒子i和j之間的最大排斥力系數,rij=ri-rj,rij=|rij|,eij=rij/|rij| .耗散力是依賴于粒子位置和速度的力,亦稱耗散阻力,其表達式為

式中γ為耗散力系數,wD(rij)為耗散力權重函數,vij=vi-vj.隨機力的表達式為

式中σ為隨機力系數,wR(rij)為隨機力權重函數,ζij是均值為0方差為1的隨機數,θij為高斯分布的隨機數,具有對稱性θij=θji,且具有以下隨機特性,對i≠j,k≠l：

式中 〈···〉 表示求平均.

耗散力阻止粒子間相對運動,使粒子速度減慢,其宏觀效應是使系統溫度降低,而隨機力宏觀效應是使系統升溫.這樣,在耗散力和隨機力的協同作用下,系統能保持溫度恒定,即引入的耗散力和隨機力相當于分子動力學模擬中的恒溫熱浴.

DPD模擬中,各物理量均采用無量綱單位,通常選取rc為特征長度,密度ρ為特征密度,kBT為特征能量,其他物理量單位由這個三個參數推導得到.在DPD模擬中,通常取無量綱截斷半徑rc=1,無量綱溫度kBT=1,粒子無量綱質量m=1,無量綱密度ρ= 3 — 6.通過定義DPD粒子的質量、平均熱能和截斷半徑,可以得到對應的實際物理單位.

2.2 壁面邊界模型

由于耗散粒子動力學中粒子之間作用為軟勢能的形式,固體壁面粒子的軟的排斥力不足以阻止流體粒子穿透壁面,導致施加固-液界面邊界條件比較困難.此外,壁面邊界結構處理不合理還會造成壁面附近一些假象的波動,如壁面附近流體密度和溫度波動.DPD固體壁面邊界條件是DPD模擬中的研究重點和熱點問題,過去20年很多學者提出了不同的施加固體壁面邊界條件的方法[25-31],其中一個最常見的方法是采用固定住的粒子并結合適當的反射機制構建固體壁面邊界,如Revenga等[25]提出了三種反射機制,即鏡像反射、向后反彈發射與Maxwellian反射,將穿過壁面的流體粒子重新移回流體區域,劉謀斌和常建忠[31]通過凍結隨機分布并且達到平衡狀態的DPD粒子,形成復雜的固體區域,結合反彈反射,該方法能有效地模擬復雜固體區域.本文模擬中采用固定住的粒子構建固體壁面,并通過一定的反射機制將穿透壁面的流體粒子重新移回流體區域.圖1為固體壁面結構示意圖,圖中黑色圓圈代表壁面粒子,空白圓圈代表流體粒子,空白圓圈上的箭頭表示粒子的速度方向.如圖1所示,采用三層固定住的粒子構建固體壁面,由Duong-Hong等[28]的研究工作可知,合理地調整這幾層固體壁面粒子與固-液界面之間的距離,可以消除壁面附近的密度波動并能阻止大多數流體粒子穿過固-液界面.本文模型中,三層固定住的粒子離固-液界面的距離分別為0.15rc,0.35rc,0.75rc,第一層壁面粒子離固-液界面的距離設置得比較小,能適當地增大壁面的排斥力,有助于阻止大部分流體粒子穿過固-液界面,另外兩層壁面粒子主要是為了給壁面附近的流體粒子提供均勻的壁面排斥力,進而消除壁面附近流體密度波動.相鄰粒子間的距離在x和y方向分別為0.5rc和rc.壁面粒子與壁面附近的流體粒子之間也受到(2)式中的三種作用力.

此外,當流體粒子穿過固-液界面時,采用如圖2所示的修正過的向前反彈機制將穿透壁面的流體粒子重新移回流體區域.如圖2所示,t時刻粒子i處于位置1處,在t+Δt時刻粒子i穿過了固-液界面運動到位置2處,此時將粒子i重新移回到流體區域的位置3處,位置2與位置3關于固-液界面對稱,并且移回到位置3處的流體粒子速度方向與位置2處的速度方向相反(圖中從粒子中心引出的箭頭線表示速度分量).傳統的向前反彈機制是移回到位置3處的流體粒子速度方向與位置1處的速度方向相反,也即將t+Δt時刻的速度與t時刻的速度反向,兩個時刻的速度大小相同,而本文改進的反彈機制是將t+Δt時刻位置2處速度反向,得到位置3處的速度,這更符合物理運動規律.

圖1 固體壁面結構示意圖Fig.1.Sketch of the solid wall structure.

圖2 修正的向前反彈機制Fig.2.Sketch of the modified bounce-forward reflection.

2.3 模擬系統與模擬參數

本文以平板間的Couette流為例,探討微通道疏水表面的滑移問題.圖3所示為微通道(兩平板間)Couette流模擬系統示意圖,x和y方向采用周期性邊界條件,z方向采用上述壁面邊界模型.模擬區域的實際尺寸在x,y,z方向為Lx×Ly×Lz=60rc×5rc×30rc.DPD單位中,截斷半徑rc的無量綱單位為1,則模擬區域無量綱尺寸為Lx×Ly×Lz=60×5×30.笛卡爾坐標系的原點設在模擬區域的中心,上下固-液界面分別位于z=15與z=-15 處.所有流體粒子和壁面粒子的質量均相等,每個流體粒子或壁面粒子的質量設為1(m=1 ).所有的壁面粒子按照上述壁面邊界模型布置在上下壁面區域,則一共含有3600個壁面粒子,上下壁面各1800個粒子.所有流體粒子的初始位置按照面心立方體晶格布置在模擬區域,流體的密度取ρ=4.0,則一個單位立方體內包含4個流體粒子,流體粒子的總數為36000.所有的模擬中,系統的溫度取kBT=1,流體粒子的初始速度根據此溫度隨機設定.取耗散力系數γ=4.5,隨機力系數為σ=3.0 .

圖3 模擬系統示意圖Fig.3.Sketch of simulation system of Couette flow.

上下壁面沿x方向分別施加等值反向的速度v和—v,用于驅動形成Couette流.采用Groot和Warren[21]提出的修正Velocity-Verlet算法數值求解運動方程(1)—(9)式,獲得流體粒子的運動軌跡.數值求解的時間步長取 Δt=0.02 .為減少計算時間,采用文獻[32]描述的元胞分割法和臨位列表法計算粒子之間的作用力.

2.4 壁面疏水表征

本文模擬系統包含兩種粒子,即壁面粒子和流體粒子,壁面粒子間排斥力系數為aww,流體粒子間的排斥力系數為aff,其中下角標‘w’表示壁面粒子,‘f’表示流體粒子,根據文獻[27]定義壁面粒子與流體粒子之間的排斥力系數：

由(3)式可知,排斥力系數代表相互作用粒子間的最大排斥力,即awf為壁面粒子與流體粒子間的最大排斥力.因此,awf的大小可以反映出壁面與流體的相互作用強度,通過調整awf的大小即可實現壁面不同的親/疏水性能,awf越大,壁面與流體的排斥作用越強,疏水性能越強.文獻[16]和文獻[33]就是通過調整壁面粒子與流體粒子之間的排斥力系數awf,從而實現壁面的不同潤濕性能,文獻[16]研究結果表明,排斥力系數awf與接觸角呈近似線性關系,壁面與流體的排斥作用越強,流體的接觸角越大,調整awf大小可以實現從親水到疏水的轉變.本文模擬中取流體粒子間的排斥力系數為aff=18.75,為使壁面具有不同的親/疏水性能,實現無滑移到滑移的變化,設計了不同壁面粒子間的排斥力系數aww=5,10,15,20,25,對應的壁面粒子與流體粒子之間的排斥力系數為awf=9.68,13.69,16.77,19.36,21.65.

3 微通道疏水表面滑移

3.1 壁面無滑移

在研究壁面滑移之前,先討論壁面無滑移的情況.一方面有助于與壁面滑移進行比較,另一方面可以驗證上述壁面邊界模型的有效性以及編制的DPD程序的正確性.本文研究的表面滑移均為部分滑移(partial slip),圖4為壁面無滑移和部分滑移示意圖(壁面是靜止不動的).部分滑移邊界條件可以用Navier滑移邊界模型描述如下：

式中us為滑移速度,b為滑移長度,ux為x方向的速度,滑移速度和滑移長度的定義如圖4所示.

圖4 壁面邊界的無滑移和部分滑移示意圖Fig.4.No-slip and partial slip status at a solid-fluid interface.

在壁面無滑移模擬中,取壁面粒子間的排斥力系數aww=5,對應的壁面粒子與流體粒子間的排斥力系數為awf=9.68 .模擬先運行一段時間,使系統達到熱力學平衡狀態,然后給上壁面施加沿x正方向的速度=3.0,同時給下壁面施加沿x負方向的速度=-3.0,驅動形成Couette流動.模擬再運行一段時間,使系統達到穩定的Couette流后開始統計相關物理量.為獲得速度、密度、溫度以及剪切應力在通道z方向上的分布,通道在z方向上被劃分為200層,每 2×104個時間步長對每層的數據進行統計平均得到相關物理量.

圖5 DPD模擬的速度分布與Navier-Stokes(NS)分析解對比Fig.5.The velocity profile of DPD simulation result compares to the analytical solution.

圖6 DPD模擬的密度、溫度分布與Navier-Stokes(NS)分析解對比Fig.6.The density and temperature profiles of DPD simulation results compare to the analytical solutions.

圖7 DPD模擬的剪切應力分布與分析解對比Fig.7.The shear stress profile of DPD simulation result compares to the analytical solution.

圖5為DPD模擬的Couette流的速度分布與NS方程得到的分析解對比,圖6為DPD模擬得到的密度、溫度曲線與NS分析解比較,圖7為模擬的剪切應力分布與分析解對比.從圖5—圖7可以看到,DPD模擬得到速度、密度、溫度、剪切應力分布與理論分析解符合得很好,驗證了本文編寫的DPD程序的正確性.更重要的是,從圖5可知壁面處流體的速度與壁面的速度相同,即壁面無滑移.從圖6可以看到,模擬得到的密度與溫度分布比較均勻,密度和溫度曲線與理論曲線幾乎重合,在壁面附近密度、溫度幾乎沒有波動,圖7顯示模擬得到的剪切應力分布與分析解(Szx=μ?vx/?z,μ為動力黏度)符合得也很好,驗證了上述邊界模型的有效性.

由此可見,壁面粒子與流體粒子間的排斥力系數為awf=9.68 時,壁面實現了無滑移邊界條件,壁面附近密度分布均勻無波動.由模擬結果也可以看到,模擬得到的速度、溫度、密度以及剪切應力等物理量的分布與Navier-Stokes理論分析解符合得很好,這也說明了本文固體壁面邊界模型以及DPD模型具有很好的收斂性.

3.2 疏水表面滑移

由前面分析知道,壁面粒子與流體粒子間的排斥力系數為awf=9.68 時壁面為無滑移壁面,這也能反映出壁面為親水性壁面.現增大awf,即增強壁面與流體的排斥作用,使壁面實現從親水到疏水的轉變,分析疏水表面的滑移機理,討論滑移長度與awf的關系.

圖8為壁面粒子與流體粒子間的排斥力系數awf分別為9.68,13.69,16.77,19.36和21.65時的速度分布.由于速度分布是關于坐標原點對稱的,為了顯示更清晰些,圖8只顯示了速度分布的一半,即z=-15—0 之間的速度曲線.awf=9.68 時,壁面滿足無滑移條件,而當增大awf,也即增強壁面與流體之間的排斥作用,使壁面向疏水轉變,速度分布發生了變化,也即當awf＞9.68 時,壁面產生了滑移.根據Navier滑移邊界的定義,將主流區的速度曲線線性延伸至壁面,得到的速度若小于壁面的速度則為滑移邊界.由圖8可以看到,awf=13.69,16.77,19.36,21.65時,壁面均產生了滑移,并且隨著awf的增大,即疏水性增強,壁面滑移越明顯.

圖8 排斥力系數 awf分別為9.68,13.69,16.77,19.36,21.65時的速度分布Fig.8.Velocity distributions due to different values of awf=9.68,13.69,16.77,19.36,21.65.

為探索排斥力系數awf(疏水強度)對壁面滑移的影響,下面分析滑移速度和滑移長度與排斥力系數awf的關系.根據Navier滑移邊界的定義,滑移速度和滑移長度可以通過如下方法獲得：對主流區的速度線性擬合得到直線形的速度分布,擬合的速度直線延伸至壁面處時得到速度與壁面速度的差值即為滑移速度,擬合的速度直線向外延伸出實際流場區域直至滿足無滑移邊界條件時離壁面的距離即為滑移長度.表1為awf= 13.69,16.77,19.36,21.65時通過線性擬合得到的速度分布以及對應的滑移速度和滑移長度,表中還列出無滑移時的情況(滑移速度和滑移長度均為0).根據表1中的數據,得到了圖9所示滑移速度us與排斥力系數為awf的關系,以及圖10所示滑移長度b與排斥力系數awf的關系.由圖9和圖10可以看到,隨著壁面粒子與流體粒子間的排斥力系數awf的增大,即壁面與流體之間排斥作用增強(疏水性增強),滑移速度和滑移長度均增大.分別對滑移速度us和排斥力系數awf、滑移長度b和排斥力系數awf進行曲線擬合,如圖9和圖10所示,擬合結果顯示滑移速度us或滑移長度b與排斥力系數awf之間存在近似的二次函數關系,函數表達式分別為：

表1 不同排斥力系數 awf時擬合的速度分布及對應的滑移速度和滑移長度Table 1.The linear fit velocity profiles,slip velocity and slip length with respect to different awf.

圖9 滑移速度 us與排斥力系數 awf的關系Fig.9.Slip velocity versus awf.

圖10 滑移長度b與排斥力系數 awf的關系Fig.10.Slip length versus awf.

表面的親/疏水性一般是通過接觸角表征的,表面疏水性越強,接觸角越大,表面滑移也越明顯,但目前關于滑移長度和接觸角之間還沒有明確的關系式.根據文獻[16]的DPD模擬結果可以看到,壁面與流體之間作用參數awf與接觸角呈近似線性關系,結合本文研究結果(13)式,可以認為滑移長度與接觸角之間存在近似的二次函數關系.接觸角是衡量液體對固體材料表面潤濕性能的重要參數,反映了液體與固體表面親和作用的強弱,也即接觸角的大小取決于液體與固體表面的相互作用強度.而疏水表面的滑移本質上是由液體與固體表面相互作用引起的,滑移長度的大小同樣取決于液體與固體表面的相互作用強度.因此,滑移長度與接觸角之間必然存在著某種函數對應關系.本文模擬中,壁面粒子與流體粒子之間的排斥力系數awf反映了流體與固體壁面的相互作用強度,由(13)式可以看到,滑移長度與接觸角之間滿足近似二次函數關系.

目前關于疏水表面滑移現象的物理機理還沒有定論,比較常見的兩種解釋為：一是不連續氣體層或納米氣泡理論,認為在表面存在不連續氣體區域,因此固體表面同時存在液-固界面和液-氣界面,液-氣界面流體流動阻力小,產生了滑移；二是低密度層理論,由于排斥作用,疏水表面附近存在低密度區域,阻止了動量的傳遞,壁面發生了滑移.為探討疏水表面的滑移機理,圖11給出了排斥力系數awf=13.69,16.77,19.36,21.65時的密度分布和無滑移時的密度分布,由于遠離壁面區域的密度幾乎一樣,并且密度分布關于z=0 對稱,為了顯示更清晰,圖11只顯示了下壁面z=-15 至z=-10 之間的密度曲線.無滑移時理論密度分布是均勻的,當awf=13.69,16.77,19.36,21.65 時,即壁面產生滑移時,在緊鄰壁面附近流體密度很低,存在低密度區域,這些低密度區域阻礙了流體的動量傳遞,致使壁面產生滑移.從圖11可以看到,在緊鄰壁面處密度低,即存在一個低密度區域,緊鄰低密度區域又存在一個密度大于4.0的高密度區域.這是由于壁面附近存在低密度區域,導致高密度區域處的流體在垂直于壁面方向上受到的排斥力不平衡,進而出現了高密度區域.根據文獻[34,35]的研究可以知道,在DPD方法中,流體粒子受到的排斥力在垂直于壁面方向上的分力控制著密度的波動.本文中,在緊鄰壁面區域存在低密度區域,密度遠小于4.0,即流體粒子數目較少,而在遠離壁面區域密度分布均勻,密度為4.0.因此,低密度區域流體對高密度區域流體的排斥力小于密度均勻區域流體對高密度區域流體的排斥力,也即高密度區域流體粒子受到朝著壁面方向上的凈排斥力,該凈排斥力使流體粒子向低密度區域方向運動,流體粒子往該方向運動,使其與低密度區域一側流體粒子距離減小,由(3)式可知,粒子之間的距離減小,相互之間的排斥力增大,因此凈排斥力使流體粒子往低密度區域一側逐漸靠近,低密度區域一側流體粒子對該流體粒子的排斥力逐漸增大,必然存在一個平衡位置,當流體粒子運動到該平衡位置時,增大的排斥力與凈排斥力平衡,流體粒子會在該平衡位置堆積,故在緊鄰壁面低密度區域又會存在高密度區域.由圖11還可以看到,排斥力系數awf越大,即疏水性越強,緊鄰壁面附近的密度越低,阻礙動量傳遞的能力越強,滑移越明顯.

圖11 不同排斥力系數時和無滑移時的密度分布Fig.11.The density profiles for different awfvalues and for no slip condition.

4 結 論

本文利用耗散粒子動力學模擬了平板間的Couette流,研究了微通道疏水表面的滑移現象.采用三層固定住的粒子并結合修正的向前反彈機制構建了固體壁面邊界模型,通過調整壁面粒子與流體粒子之間的排斥力系數awf大小,壁面可以實現從親水到疏水的轉變,實現從無滑移到滑移的變化.當awf=9.68 時,壁面滿足無滑移邊界條件,壁面附近密度分布均勻；當awf＞9.68 時,壁面產生了滑移,并隨著排斥力系數awf增大(疏水性增強),壁面滑移越明顯,滑移速度和滑移長度越大.模擬結果還表明滑移速度和滑移長度與排斥力系數awf存在近似二次函數關系,也即與接觸角之間存在近似二次函數關系.模擬得到的密度分布顯示,無滑移時壁面附近的密度分布均勻,而疏水表面附近存在低密度區域,低密度區域阻礙了流體的動量傳遞,導致壁面產生了滑移.