基于SDPD的動(dòng)態(tài)浸潤(rùn)研究

土俊鵬，劉漢濤,李海橋

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

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

1 數(shù)值研究方法

1.1 SDPD方法

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

dri=vidt

(1)

(2)

(3)

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

保守力FC的表達(dá)式為：

(4)

(5)

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

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

(6)

其中

(7)

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

(8)

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

(9)

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

1.2 表面張力模型

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

(10)

(11)

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

(12)

(13)

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

(14)

最終可以得到角動(dòng)量和線動(dòng)量均守恒的表面張力表達(dá)式：

(15)

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

2 結(jié)果及分析

2.1 計(jì)算模型及液滴變化過程

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

圖1 接觸角計(jì)算模型圖

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

圖2 初始粒子分布

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

圖3 超親水模擬變化過程

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

圖4 超疏水模擬變化過程

2.2 結(jié)果分析

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

(16)

圖5 計(jì)算角度與理論角度對(duì)比圖

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

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

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

圖7 相對(duì)誤差曲線

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

圖8 不同計(jì)算尺寸的影響

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

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

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

3 結(jié)論

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

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

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