土俊鵬，劉漢濤,李海橋

(中北大學(xué)能源動力工程學(xué)院，山西 太原 030051)

隨著對自然界浸潤現(xiàn)象的深入研究，各種結(jié)構(gòu)表面的浸潤性能受到了廣泛關(guān)注，目前不同浸潤性表面仿生材料的開發(fā)與應(yīng)用成為科學(xué)研究的熱點和難點之一[1-3]。在表面材料浸潤性能的研究過程中發(fā)現(xiàn)，材料自身屬性起著非常重要的作用，多相間表面張力的協(xié)同作用是決定表面材料性能的最重要因素之一[4-6]。多相間表面張力的相互作用是一個復(fù)雜的動態(tài)過程，這個過程涉及到流體結(jié)構(gòu)大變形，且本質(zhì)上是一個跨尺度的復(fù)雜物理現(xiàn)象。對動態(tài)浸潤問題進行實驗研究是非常困難的，很難復(fù)現(xiàn)動態(tài)浸潤過程中的真實物理條件和流體動力學(xué)行為。數(shù)值模擬技術(shù)的發(fā)展給深入研究動態(tài)浸潤問題提供了便利，但基于網(wǎng)格的流體模擬方法很難捕捉多相間復(fù)雜的界面力學(xué)行為，在處理動態(tài)浸潤過程中涉及的流體結(jié)構(gòu)大變形還會遇到網(wǎng)格畸變等固有問題的困擾。無網(wǎng)格方法在處理上述問題時具有獨特的優(yōu)勢，其中光滑粒子動力學(xué)(smoothed particle hydrodynamics, SPH)作為最早的無網(wǎng)格方法之一，在研究多相流動等復(fù)雜流體動力學(xué)行為時受到了廣泛關(guān)注[7-8]。然而，SPH是一種宏觀尺度的計算方法，不能準(zhǔn)確刻畫復(fù)雜流體動力學(xué)行為中的多尺度機理[9]，因此Espaol等[10-11]將SPH和與耗散粒子動力學(xué)(dissipative particle dynamics, DPD)無網(wǎng)格方法相結(jié)合，提出了光滑耗散粒子動力學(xué)(smoothed dissipative particle dynamics, SDPD)方法。與DPD方法一樣，SDPD中考慮了粒子的熱運動，引入了隨機運動，可以描述介觀系統(tǒng)的隨機振蕩，因此SDPD能夠從介觀尺度上揭示浸潤問題內(nèi)部復(fù)雜的流體動力學(xué)行為。但欲真正揭示浸潤問題的內(nèi)部機理，準(zhǔn)確描述氣體、液體以及固體間復(fù)雜的界面張力是需要解決的最關(guān)鍵的問題之一[12]。鑒于SDPD方法自身的優(yōu)勢和特點，近年來已經(jīng)應(yīng)用到研究復(fù)雜多相流問題中[13-14]，其中Hu等[15]構(gòu)建的多相流模型受到了廣泛關(guān)注，雖然該模型在研究固體表面上的浸潤問題已有較多的應(yīng)用，但很少有該模型在通用性方面的研究。本文基于Hu等提出的SDPD模型，不僅研究了普通的親疏水問題，還研究了超親疏水問題，通過系統(tǒng)地對比研究，結(jié)合研究成果和動態(tài)浸潤過程中流體動力學(xué)行為，對利用SDPD模型研究浸潤問題的通用性進行討論，同時對該模型存在的問題進行剖析，并給出改進意見。本文的工作為揭示浸潤問題的本質(zhì)提供了理論依據(jù)，同時為SDPD方法的發(fā)展和改進打下了理論基礎(chǔ)。

1 數(shù)值研究方法

1.1 SDPD方法

SDPD方法演化方程的通用表達式可以寫成：

dri=vidt

(1)

(2)

(3)

式中：ρ，m分別為粒子的密度與質(zhì)量；σ為粒子的數(shù)量密度；W為核函數(shù)；下標(biāo)i和j為粒子的編號。

保守力FC的表達式為：

(4)

(5)

式中：η為動力黏度;上標(biāo)k和l表示不同的相；vij為從粒子i到粒子j的速度矢量。

根據(jù)GENERIC形式[16]中動量耗散形式，隨機力的表達式可以寫為：

(6)

其中

(7)

式中：kB和T分別為玻爾茲曼常數(shù)和系統(tǒng)溫度。本文采用五次樣條核函數(shù)[17]，其表達式為：

(8)

式中：W(R,h)為光滑核函數(shù)，其中R為兩粒子間的相對距離，h為光滑核函數(shù)W影響區(qū)域的光滑長度；αd為處理不同維數(shù)問題時的正則化參數(shù)，本文計算的是二維問題，αd=7/478πh2。粒子由于密度變化引起的壓力變化采用人工狀態(tài)方程[18]求解，其表達式為：

(9)

式中：p(ρ)為粒子壓力；p0,ρ0為參考值；γ通常為7。

1.2 表面張力模型

本文用連續(xù)表面張力模型[19]計算相與相之間的界面張力Π(1)，其表達式為：

(10)

(11)

在相內(nèi)部，色函數(shù)的梯度是為0的，而在相與相的界面處，支持域內(nèi)粒子分別屬于不同相時，色函數(shù)的梯度不為0，色函數(shù)梯度的計算式為：

(12)

(13)

式中：d為計算域的維度。粒子i總的表面張力為：

(14)

最終可以得到角動量和線動量均守恒的表面張力表達式：

(15)

固壁邊界處理是實現(xiàn)浸潤問題準(zhǔn)確計算的關(guān)鍵之一，本文采用文獻[15]中的方法，即利用鏡像粒子邊界[19]實現(xiàn)無滑移邊界條件，使用CFL[20](courant-friedrichs-levy)條件對時間步長進行估算。

2 結(jié)果及分析

2.1 計算模型及液滴變化過程

在不考慮重力的前提下，根據(jù)楊氏方程，固體表面上三相間的接觸角θ滿足α1w=α12cosθ+α2w，其中上標(biāo)1,2和w分別表示氣體、液體和固體，即α1w，α2w,α12分別為氣體與固體、液體與固體、氣體與液體間的表面張力。調(diào)節(jié)三相間的表面張力系數(shù)，并運用θ=2arctan(2H/L)來測量液滴的接觸角，如圖1所示。本文選擇超親水(理論接觸角接近10°)和超疏水(理論接觸角為170°)進行研究，深入研究該方法對超親、疏水的適應(yīng)性。

圖1 接觸角計算模型圖

本文計算模型的物理尺度為介觀尺度。在計算親、疏水時，若粒子間距離相同，則采用相同的光滑長度。在計算超親水時，由于液滴鋪展長度較大，因此計算區(qū)域選定為0 μm

圖2 初始粒子分布

計算超親水動態(tài)浸潤時的液滴形態(tài)以及流場的變化狀況如圖3所示。液滴頂部有向內(nèi)匯聚的流動，而在固體表面上，液體沿著壁面動態(tài)鋪展，尤其在三相接觸點的位置，液體的運動更加明顯，由于液滴頂部和底部運動方向相反，在液滴的腰部產(chǎn)生了明顯的渦，使得液滴形貌在很短的時間內(nèi)發(fā)生了劇烈變化。隨著時間的推移，液滴周圍的區(qū)域還有明顯的運動，而在內(nèi)部靠近壁面的液體在無滑移邊界條件的作用下，速度越來越小。當(dāng)液滴外部的曲率趨于一致時，表面張力趨于平衡，液滴在黏性力作用下，慢慢趨于穩(wěn)定。液滴穩(wěn)定后，根據(jù)θ=2arctan(2H/L)計算可得接觸角為27.3°，與理論值有較大的偏差。

圖3 超親水模擬變化過程

計算超疏水時，由于沒有長距離的鋪展，計算區(qū)域選定為0 μm

圖4 超疏水模擬變化過程

2.2 結(jié)果分析

通過上述研究，發(fā)現(xiàn)Hu等提出的方法在計算超親水或超疏水現(xiàn)象時有很大的誤差，為了更系統(tǒng)地研究該模型產(chǎn)生誤差的原因，設(shè)定α1w=1，α12=1，通過改變α2w的值來計算接觸角。從圖5可以看出，在親水性模擬中計算得到的液滴接觸角比理論角度都大，在疏水性模擬計算中比理論角度小，在90°附近誤差最小，為1%，60°～120°時的誤差在5%以內(nèi)，基本滿足計算精度的要求，但在計算超親水和超疏水角度時誤差比較大。由于液滴鋪展長度的變化比液滴高度的變化明顯，為進一步探討該模型的計算精度，通過改變不同條件，研究液滴鋪展長度的變化規(guī)律。其中液滴鋪展理論長度的計算方法采用Thmas[21]提出的方法，該方法中鋪展長度L與接觸角的關(guān)系為L=2Rsinθ，其中曲率半徑R是液滴的截面積A和接觸角θ的函數(shù)：

(16)

圖5 計算角度與理論角度對比圖

計算區(qū)域不變，當(dāng)粒子間距為0.083 μm和0.042 μm時，對應(yīng)液滴內(nèi)粒子數(shù)目分別為216和648，兩種分辨率下的計算結(jié)果如圖6所示。當(dāng)粒子間距減小，粒子數(shù)目增加時，在計算親水角度時，液滴鋪展距離更長，在計算疏水角度時，液滴收縮程度更高，說明隨著粒子數(shù)目的增加，計算精度會提高。

圖6 不同分辨率計算結(jié)果

選定與理論角度出現(xiàn)較大偏差的超親水(10°)，計算液滴內(nèi)粒子數(shù)目從216到1 080時液滴長度與理論值的相對誤差，結(jié)果如圖7所示，由圖可知，隨著液滴內(nèi)粒子數(shù)目的增加，其接觸角與理論值的偏差逐漸減小，說明在接觸角附近，通過提高分辨率可以更細致地刻畫液滴的鋪展過程。

圖7 相對誤差曲線

SDPD方法中，用隨機力來模擬粒子熱波動，由式(7)可以看出，計算尺寸對隨機力的計算有顯著影響。設(shè)系統(tǒng)溫度為300 K，計算模型不變，粒子數(shù)目一定(液滴內(nèi)粒子數(shù)為216)，通過改變粒子單元的尺寸來改變計算區(qū)域的大小，粒子單元尺寸分別為0.10 μm×0.10 μm，0.25 μm×0.25 μm和1.00 μm×1.00 μm。在這3個計算尺寸下計算多個接觸角對應(yīng)的液滴鋪展長度，并對應(yīng)到同一坐標(biāo)系中。從圖8的對比結(jié)果可以看出，不同尺度的液滴總體變化規(guī)律相近，但對應(yīng)的鋪展長度會因為粒子熱波動呈現(xiàn)輕微的差異。

圖8 不同計算尺寸的影響

圖9所示是3個不同尺寸液滴達到穩(wěn)定狀態(tài)時的粒子分布(60°)，由圖可知，粒子單元尺寸越小，液滴外輪廓越不規(guī)則，內(nèi)部粒子分布無序度越高；粒子單元尺寸越大，外輪廓越光滑，內(nèi)部粒子分布越具有一定的規(guī)律性。由此可證明，由隨機力引起的熱波動隨著粒子單元尺寸的減小而顯著增大，表明了SDPD方法能描述介觀尺度不同物理尺寸的熱波動特性。

圖9 不同單元尺寸液滴平衡狀態(tài)時粒子分布圖

就目前的研究結(jié)果可以看出，鏡像邊界處理方法中壁面粒子的位置會隨著流體粒子位置的改變而發(fā)生改變，當(dāng)邊界處壓力梯度較大時，產(chǎn)生較大的排斥力，會導(dǎo)致邊界粒子與流體粒子的距離增大，致使局部區(qū)域粒子變得稀疏，影響計算結(jié)果，采用固定的邊界粒子，可以消除由此引起的誤差。在計算超親水和超疏水時，通過兩種接觸角的求解方法，均沒有達到預(yù)期的結(jié)果，這表明：1)通過穩(wěn)定時液滴外輪廓的曲線方程，求出三相點處的曲線斜率的方式求出接觸角，可以更準(zhǔn)確地計算出接觸角的大??；2)除了表面張力的作用外，還有固體微結(jié)構(gòu)與液體間的相互作用，可以通過引入特殊的液固粒子間作用力模型來處理相關(guān)的問題，減小超親水和超疏水接觸角計算誤差；3)本文在計算過程中，并未考慮能量方程，而在超親水或超疏水現(xiàn)象中，能量的變化也有一定的影響，考慮能量方程可以更加準(zhǔn)確地刻畫超親水或超疏水現(xiàn)象。通過減小粒子間的距離可以在一定程度上提高計算精度，證明在研究超親水或超疏水問題時，通過提高分辨率可以刻畫浸潤過程中的更多細節(jié)，引入粒子自適應(yīng)劃分技術(shù)，細化三相接觸線及其附近的區(qū)域，可以更好地研究各類浸潤問題。

3 結(jié)論

本文基于SDPD方法研究超親水和超疏水兩種特殊的浸潤問題，結(jié)果表明：

1)SDPD方法在計算超親水或超疏水等涉及流體結(jié)構(gòu)大變形問題時，方法自身的粒子自適應(yīng)分布能夠捕捉多相間復(fù)雜的界面力學(xué)行為，為研究各類復(fù)雜的介觀尺度浸潤問題提供了可靠的數(shù)值方法；

2)Hu等提出的多相流模型和鏡像邊界條件在處理浸潤問題時，當(dāng)接觸角為30°～135°時基本能滿足計算精度的要求，但在研究超親水和超疏水現(xiàn)象時，計算誤差較大，需要引入其他方法來緩解計算中誤差較大的問題。