納米流體液滴的耗散粒子動力學模擬

沈世元,周哲瑋

(上海大學上海市應用數學和力學研究所,上海　200072)

納米流體液滴的耗散粒子動力學模擬

沈世元,周哲瑋

(上海大學上海市應用數學和力學研究所,上海200072)

采用在介觀尺度下適用的耗散粒子動力學(dissipative particle dynamics,DPD)方法,對考慮靜電力作用的納米流體系統進行建模,研究了在納米顆粒帶電量和濃度影響下納米流體液滴接觸角的變化.通過對實驗過程進行的數值模擬,得到了定性相同的結果.

納米流體；介觀尺度；耗散粒子動力學；極化

2003年,You等[1]基于納米流體的沸騰熱傳遞(boiling heat transfer,BHT)研究,分析了臨界熱通量(critical heat flux,CHF)的變化,得到了臨界熱通量顯著升高的結論,但未給出臨界熱通量升高的原因.這種未能解釋的、有趣的現象吸引了大量的科研工作者對納米流體臨界熱通量的顯著變化進行研究.2005年,Vafaei等[2]制備了納米顆粒被巰基乙酸功能化的碲化鉍納米流體,并發現納米流體液滴的接觸角與納米顆粒的濃度及尺寸有關.2006年,Kim等[3]分別用3種不同的氧化物顆粒(氧化鋯、氧化鋁和氧化硅)制備了納米流體,并對其性質進行了研究,發現在干凈的固壁表面上,納米流體液滴的接觸角隨濃度變化不大,但是對于壁面上存在沉降的納米顆粒的情況來說,接觸角由于潤濕性的增大有明顯的減小.2010年,Ahn等[4]通過方腔流實驗發現,與純水相比,納米流體流過的固壁表面上的臨界熱通量顯著升高.通過對接觸角的測量和固壁表面性質的研究,確定臨界熱通量的升高主要是由納米顆粒沉降在固壁表面導致的潤濕性增大造成的.

對這一現象本課題組進行了研究,多體耗散粒子動力學(many-body dissipative particle dynamics,MDPD)方法的模擬顯示,對于懸浮在水中的SiO2納米顆粒,其濃度對液滴接觸角無影響；而沉降在固壁表面的SiO2納米顆粒數則影響液滴接觸角,該結果與Ahn等[4]的實驗結果定性符合.Vafaei等[2]指出,在給定的實驗條件下,納米流體液滴的接觸角與納米顆粒的濃度有關,原因是修飾納米顆粒表面的官能團在水中產生了極化,納米顆粒與流體的相互作用以及納米顆粒之間長程的靜電作用對納米流體的毛細管性能產生了影響.本工作主要參考膠體模型,通過帶有電荷的耗散粒子動力學(dissipative particle dynamics,DPD)方法模擬納米流體液滴來研究Vafaei等[2]的實驗.

1　DPD方法

1.1標準的DPD方法

1992年,DPD方法由Hoogerbrugge等[5]提出.1998年,Warren[6]認為耗散粒子可以表示分子、原子、高分子蛋白質聚團、帶電離子團等,指出了DPD方法的廣泛適用性.之后,學者們對于介觀尺度下的納米流體做了一些基于DPD方法的數值模擬工作.2008年,He等[7]利用能量守恒DPD方法模擬了納米流體的熱傳導,指出熱導率與流體中包含的納米顆粒數成正比關系,并認為產生這種現象的原因是納米顆粒的布朗運動.2012年,Yamada等[8]利用能量守恒DPD方法模擬了一維微通道納米流體的熱傳導,與實驗結果符合較好,并指出納米流體的流動速度、pH值、納米顆粒形狀都會影響納米流體的熱傳導率.目前,利用分子動力學(molecular dynamics,MD)方法研究納米流體的報道[9-13]較多,但本工作利用DPD方法對納米流體進行研究.

與MD方法相同,DPD方法采用Lagrange坐標描述粒子,體系中的粒子運動滿足牛頓運動方程：

ωC,ωD,ωR為各力的權重函數；rij為i,j粒子間的距離；eij為和rij方向相同的單位矢量；aij為保守力系數,反映粒子i和j之間的最大排斥力；γ和σ分別為耗散力和隨機力幅值；vij=vi?vj,θij為高斯白噪聲項,滿足以下統計學特性：

Espanol等[13]將耗散漲落定理引入耗散粒子動力學方法中,為耗散力和隨機力提供了統計熱力學上的平衡條件,即

本工作采用Groot等[14]提出的修正的Velocity-Verlet算法：

式中,λ為可調參數.當λ=0.5時,該算法和標準的Velocity-Verlet算法[15]一致.目前,絕大部分DPD模擬都采用這種修正的Velocity-Verlet算法.

在DPD體系中,將單個DPD粒子的質量mi、截斷半徑Rc和能量kBT選取為系統的參考單位.例如,以水為例,將N0個質量為m的水分子粗粒化為一個DPD粒子,粒子的質量為mi=N0·m,繼而可以得到無量綱單位1對應的真實物理長度：

式中,ρD為DPD粒子的數密度,ρ0為水的密度.

通過長度、質量、能量共3個無量綱基本量,可以分別導出速度、時間和力的無量綱量表達式：

式中,kB是玻爾茲曼常數,一般來說,kBT=1.0；M?=N0×m,即單個DPD粒子的質量.如無特別說明,本模擬中采用的參數均以DPD無量綱單位給出.

1.2改進的DPD算法

在DPD模型中最重要的參數是各種粒子間的保守力參數.Groot等[14]通過維里方程得到了狀態方程(equation of state,EOS),并通過算例在數密度足夠高(ρD＞2)的前提下得到了保守力參數的近似形式：

式中,a為排斥力系數.標準的DPD方法在介觀模擬方面已經取得了很大的成功,而MDPD方法使狀態方程中出現密度的高次項,可以解決標準的DPD方法所無法解決的模擬氣液共存的問題[14,16].在本工作中,為了模擬氣液共存的狀態,采用MDPD方法來對液滴進行模擬. MDPD方法是在保留了耗散力、隨機力的基礎上,通過引入吸引相互作用,對保守力進行了修正.修正后的保守力的表達形式為

rc和rd分別為吸引和排斥截斷半徑.

Warren[16]改進了狀態方程的形式.改進后的狀態方程為

式中,α為一個平均場的系數,數值在0.101左右；A,B為系數矩陣；c為密度相關常數,d由10%×cρ2D得出.Warren[16]通過一系列的數值實驗,選取了兩組參數開展進一步的研究,最終確定了系數矩陣A和B的具體數值.

2　靜電相互作用

當DPD粒子帶有電荷時,靜電相互作用隨粒子間距增大而衰減的速度遠低于范德華作用,如果對靜電勢能函數采取直接截斷的方式,其截斷半徑要遠大于范德華作用的截斷半徑,否則將會產生較大的誤差.同時,由于DPD粒子的軟粒子性質,當一對異性電荷無限接近時,電場力將出現奇異性,導致人工離子對的產生.為了避免上述這些問題,當采用DPD模擬靜電相互作用時,粒子所帶電荷需要用電荷分布來表示.本工作基于標準Ewald求和法和分布電荷來求解靜電力.

在求解靜電相互作用時,相關參數的取值參考文獻[17].對于本工作中長方體的計算域的算例,采用Ewald求和法求解DPD模型中兩個帶電粒子之間的作用力,結果如圖1所示.可以看出,計算結果與理論結果吻合較好.同時也可以看出,DPD模型中的靜電力也是軟作用力,在兩個DPD粒子質心距離r＜1的部分有定義說明DPD粒子之間可以相互侵入.

圖1　DPD模型中兩個帶電粒子間的作用力Fig.1　Interaction force between two particles in DPD model

3　建模與計算

3.1建模

在帶電膠體體系中存在雙電層結構,該結構對膠體顆粒的穩定性起到決定性作用.20世紀40年代,前蘇聯學者Derjaguin等[18]與荷蘭學者Verwey等[19]分別獨立地提出了各種形狀的膠體顆粒間相互吸引能與雙電層排斥能的計算方法,并對膠體體系的穩定性進行了定量描述,即描述膠體體系穩定性的DLVO(Derjaguin-Landau-Verwey-Overbeek)理論.DLVO理論認為,膠體顆粒間同時存在著范德華吸引作用和雙電層引起的靜電排斥作用,并且體系的穩定性取決于膠體顆粒之間吸引作用和排斥作用的相對大小.

MDPD程序采用的參數如表1所示,其中參數均以無量綱單位給出.

表1　MDPD程序的參數取值Table 1　Parameters selection of MDPD program

3.2計算

在實驗中,被巰基乙酸功能化的Bi2Te3納米流體,可通過改變Bi2Te3顆粒的濃度來改變接觸角的大小.在數值模擬中,不僅要考慮改變納米顆粒的濃度大小,由于體系極化程度與納米顆粒帶電量的關系還不清楚,因此還需要考慮改變納米顆粒帶電量大小.

對于帶電量Q的設置,只需要在程序中賦值即可.對于濃度的改變,有兩種方案：①調整納米顆粒數,從50開始,增量為50,增加到350為止,以使濃度逐漸增加；②調整初始液滴半徑R,R越大,納米流體的濃度就越小,但是缺點是計算成本會相應增加.

黃汝佳等[22]對于純水與固壁組成的系統,給出了MDPD方法中液體與固壁相互吸引力作用系數A12與接觸角θ的函數關系：θ=5.5046A12+228.98(見圖2),并介紹了計算接觸角的圓擬合算法.本工作采用圓擬合算法來計算接觸角.由于吸引力作用系數A12對接觸角的影響比其他Aij要大得多,在文獻[22]給出的“自上而下”方法中,可以由實際測量得到的液滴的靜態接觸角來決定A12,以處理非單純物質組成的較復雜的實際流體和實際固壁.

圖2　純水與固壁的吸引力作用系數與接觸角的關系Fig.2　Relationships between attractive force coefficients and their contact angles

首先選取兩個參數：初始接觸角θ0、單個納米顆粒帶電量Q,計算納米顆粒濃度由低到高變化時納米流體液滴接觸角的變化,并與文獻[2]中的實驗數據進行對比,結果如表2和圖3所示.由圖3可以看出,隨著納米顆粒濃度V的增大,納米液滴接觸角θ呈先增大后減小的趨勢.

表2　數值模擬結果與文獻[2]中實驗數據的對比(R=14,Q=8,25)Table 2　Comparisons of numerical simulation with experimental results in Ref.[2](R=14,Q=8,25)

接下來,分析各種參數情況下納米顆粒濃度和帶電量的改變對接觸角的影響.由于圓擬合算法對于接近鈍角的計算相對更加準確,所以采用85°的初始接觸角.

(1)初始液滴半徑R=5.

對于較小的液滴半徑進行計算,可在納米顆粒濃度較高、結果改變明顯時分析定性的趨勢.當納米顆粒帶電量為1,3,5時,接觸角隨納米顆粒濃度的變化如表3和圖4所示.

圖3　納米流體液滴的接觸角隨納米顆粒濃度的變化Fig.3　Relationships between nanofluid droplet contact angles and nanoparticle concentrations

表3　納米顆粒帶電量Q=1,3,5時,不同納米顆粒濃度下液滴的接觸角(R=5)Table 3　Droplet contact angles for different nanoparticle concentrations when Q=1,3,5(R=5)

當納米顆粒濃度為1.54%時,納米流體液滴接觸角隨納米顆粒帶電量Q的變化如表4和圖5所示.

表4　納米顆粒濃度為1.54%時,液滴接觸角隨納米顆粒帶電量Q的變化(R=5)Table 4　Droplet contact angles for different nanoparticle charge quantities when V=1.54%(R=5)

圖4　納米顆粒帶電量Q=1,3,5時,液滴接觸角隨顆粒濃度的變化(R=5)Fig.4　Droplet contact angles for different nanoparticle concentrations when Q=1,3,5(R=5)

圖5　納米顆粒濃度V=1.54%時,液滴接觸角隨納米顆粒帶電量Q的變化(R=5)Fig.5　Droplet contact angles for different nanoparticle charge quantities when V=1.54%(R=5)

(2)初始液滴半徑R=7.4.

當納米顆粒帶電量為1,3,5時,液滴接觸角隨納米顆粒濃度的變化如表5和圖6所示.

表5　納米顆粒帶電量Q=1,3,5時,不同納米顆粒濃度下液滴的接觸角(R=7.4)Table 5　Droplet contact angles for different nanoparticle concentrations when Q=1,3,5(R=7.4)

圖6　納米顆粒帶電量Q=1,3,5時,液滴接觸角隨納米顆粒濃度的變化(R=7.4)Fig.6　Droplet contact angles for different nanoparticle concentrations when Q=1,3,5(R=7.4)

當納米顆粒濃度為0.62%時,納米流體液滴接觸角隨納米顆粒帶電量Q的變化如表6和圖7所示.

表6　納米顆粒濃度V=0.62%時,液滴接觸角隨納米顆粒帶電量Q的變化(R=7.4)Table 6　Droplet contact angles for different nanoparticle charge quantities when V=0.62%(R=7.4)

圖7　納米顆粒濃度為0.62%時,液滴接觸角隨納米顆粒帶電量Q的變化(R=7.4)Fig.7　Droplet contact angles for different nanoparticle charge quantities when V=0.62%(R=7.4)

(3)初始液滴半徑R=14.

當納米顆粒帶電量為1,20時,液滴接觸角隨納米顆粒濃度的變化如表7和圖8所示.當納米顆粒濃度為0.07%時,液滴接觸角隨納米顆粒帶電量Q的變化如表8和圖9所示.

表7　納米顆粒帶電量Q=1,20時,不同納米顆粒濃度下的液滴接觸角(R=14)Table 7　Droplet contact angles for different nanoparticle concentrations when Q=1,20(R=14)

圖8　納米顆粒帶電量Q=1,20時,液滴接觸角隨納米顆粒濃度的變化(R=14)Fig.8　Droplet contact angles for different nanoparticle concentrations when Q=1,20(R=14)

表8　納米顆粒濃度V=0.07%時,液滴接觸角隨納米顆粒帶電量Q的變化(R=14)Table 8　Droplet contact angles for different nanoparticle charge quantities when V=0.07%(R=14)

圖9　納米顆粒濃度為0.07%時,液滴接觸角隨納米顆粒帶電量Q的變化(R=14)Fig.9　Droplet contact angles for different nanoparticle charge quantities when V=0.07%(R=14)

綜上所述,可以得出以下結論：

(1)納米顆粒帶電量較大時(Q＞20),液滴接觸角最大值在顆粒濃度為0.5%附近.當納米顆粒濃度小于0.5%時,隨著顆粒濃度的增加,納米流體液滴的接觸角緩慢增大,趨勢較不明顯；當納米顆粒濃度大于0.5%時,隨著顆粒濃度的增加,納米流體的接觸角逐漸減小.

(2)納米顆粒帶電量較小(Q=1)時,無論顆粒濃度如何變化,液滴接觸角的變化趨勢都不明顯.

(3)納米顆粒濃度較高(V＞0.5%)時,液滴接觸角最大值在Q=4附近.隨著顆粒帶電量的增大,液滴接觸角呈現先增大后減小的趨勢.

3.3進一步探討

進一步分析液滴半徑R=5的算例.為了觀察納米流體液滴(考慮靜電相互作用)內部的粒子分布情況,不顯示阻礙觀察的水的DPD粒子.系統中納米顆粒和正離子DPD粒子的分布情況如圖10所示,其中黑色的是固壁,也就是實驗中的基板初始設置固壁帶負電.圖10(a)中,淺色小球表示水的DPD粒子；圖10(b)中,深色小球和淺色小球分別代表液滴內部的納米顆粒以及正離子的DPD粒子,禁止顯示水粒子.由圖10可以看出：由于同性相斥,帶負電的納米顆粒受到了來自固壁相斥的作用力而遠離固壁；由于異性相吸,帶正電的正離子被固壁的庫侖力所吸引,沉降到了固壁.

因此,納米液滴的接觸角增大很有可能是因為初始濃度比較低時,納米流體液滴的頂端被懸浮的納米顆粒撐高,固壁附近的納米顆粒被推離壁面,攜帶附近的流體一起,使得液滴表面朝著接觸角變大的方向運動.

圖10所示的基板上布置了一層帶負電的羥基,可以將其看作雙電層中帶負電的固壁.根據雙電層的機理可知,這層帶負電的固壁會吸附帶正電的粒子并于附近形成一個Stern層.如果假設固壁帶電量是一定的,那么吸附正離子的數量僅取決于納米流體內部本身的極化程度.納米顆粒濃度越大,液滴內部游離的正離子也就越多,容易被吸附到固壁上的正離子也就越多.當一定數量的正離子被吸附到固壁上之后,固壁表面的物理性質就會發生變化,主要表現在三相點附近與液滴的相互作用發生了變化.三相點附近物理性質的變化也會影響接觸角,使得接觸角變小.

圖10　納米流體液滴中納米顆粒和正離子的分布Fig.10　Distribution for nanoparticles and positive ions in nanofluid droplet

4　結束語

Vafaei等[2]的實驗表明,納米流體液滴的接觸角與納米顆粒的濃度有關,并猜想造成這個現象的原因是修飾納米顆粒表面的官能團在水中產生了極化,納米顆粒與流體的相互作用以及納米顆粒之間長程的靜電作用對納米流體的毛細管性能產生了影響.本工作參考膠體模型,用帶有電荷的DPD系統來模擬納米流體液滴.研究結果表明,該模型可以定性地描述納米流體液滴的靜態接觸角隨納米顆粒濃度的增加而呈現先增大后減小的變化趨勢.由于顆粒濃度與實驗結果還有數量級的差別,說明僅用納米顆粒攜帶負電荷和水中布置相應正離子的模型來模擬納米顆粒被巰基乙酸功能化的碲化鉍納米流體還過于粗糙,需要進一步細化.最后通過分析納米顆粒在液滴中的分布,探討了納米顆粒影響接觸角的物理機制.

[1]YOU S M,KIM J H,KIM K H.Effect of nanoparticles on critical heat flux of water in pool boiling heat transfer[J].Applied Physics Letters,2003,83(16)：3374-3376.

[2]VAFAEI S,BORcA-TAScIUc T,PRODOwSKI M Z,et al.Effect of nanoparticles on sessile droplet contact angle[J].Nanotechnology,2006,17(10)：2523-2527.

[3]KIM S J,BANG I C,BUONGIORNO J,et al.Effects of nanoparticle deposition on surface wettability influencing boiling heat transfer in nanofluids[J].Applied Physics Letters,2006,89(15)：153107.

[4]AHN H S,KIM H,JO H J,et al.Experimental study of critical heat flux enhancement during forced convective flow boiling of nanofluid on a short heated surface[J].International Journal of Multiphase Flow,2010,36(5)：375-384.

[5]HOOGERBRUGGE P J,KOELMAN J M V A.Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics[J].Europhysics Letters,1992,19(3)：155-160.

[6]WARREN P B.Dissipative particle dynamics[J].Current Opinion in Colloid&Interface Science, 1998,3(6)：620-624.

[7]HE P,QIAO R.Self-consistent fluctuating hydrodynamics simulations of thermal transport in nanoparticle suspensions[J].Journal of Applied Physics,2008,103(9)：094305.

[8]YAMADA T,ASAKO Y,GREGORy O J,et al.Simulation of thermal conductivity of nanofluids using dissipative particle dynamics[J].Numerical Heat Transfer,2012,61(5)：323-337.

[9]EvANS W,FISH J,KEBLINSKI P.Role of Brownian motion hydrodynamics on nanofluid thermal conductivity[J].Applied Physics Letters,2006,88(9)：093116.

[10]VLADKOv M,BARRAT J L.Modeling transient absorption and thermal conductivity in a simple nanofluid[J].Nano Letters,2006,6(6)：1224-1228.

[11]EApEN J,LI J,YIp S.Mechanism of thermal transport in dilute nanocolloids[J].Physical Review Letters,2007,98(2)：028302.

[12]EApEN J,WILLIAMS W C,BUONGIORNO J,et al.Mean-field versus microconvection effects in nanofluid thermal conduction[J].Physical Review Letters,2007,99(9)：095901.

[13]ESpANOL P,WARREN P.Statistical mechanics of dissipative particle dynamics[J].Europhysics Letters,1995,30(4)：191-196.

[14]GROOT R D,WARREN P B.Dissipative particle dynamics：bridging the gap between atomistic and mesoscopic simulation[J].The Journal of Chemical Physics,1997,107(11)：4423-4435.

[15]ALLEN M P,TILDESLEy D J.Computer simulation of liquids[M].Oxford：Clarendon Press, 1987.

[16]WARREN P B.Vapor-liquid coexistence in many-body dissipative particle dynamics[J].Physical Review E,2003,68(6)：066702.

[17]王永雷.靜電相互作用的計算在耗散粒子動力學模擬中的實現及應用[D].長春：吉林大學,2012.

[18]DERjAGUIN B,LANDAU L.Theory of the stability of strongly charged lyophobic sols and the adhesion of strongly charged particles in solutions of electrolytes[J].Acta Physicochim,1941, 14(1)：30-59.

[19]VERwEy E J W,OvERBEEK J T G.Theory of the stability of lyophobic colloids[J].Journal of Colloid Science,1955,10(2)：224-225.

[20]MIAO J,DONG M,REN M,et al.Effect of nanoparticle polarization on relative permittivity of transformer oil-based nanofluids[J].Journal of Applied Physics,2013,113(20)：204103.

[21]HwANG J G,ZAHN M,PETTERSSON L A.Bipolar charging and discharging of a perfectly conducting sphere in a lossy medium stressed by a uniform electric field[J].Journal of Applied Physics,2011,109(8)：084331.

本文彩色版可登陸本刊網站查詢：http://www.journal.shu.edu.cn

Dissipative particle dynamics simulation of nanoparticles droplet

SHEN Shiyuan,ZHOU Zhewei
(Shanghai Institute of Applied Mathematics and Mechanics,Shanghai University, Shanghai 200072,China)

Mesoscale dissipative particle dynamics was used to model a nanofluid system considering electrostatic interaction.The effects on the droplet static contact angle of the concentration and charge quantity of the nanoparticles were investigated.Computational results were qualitatively agreed with experiments.

nanofluid；mesoscale；dissipative particle dynamics；polarization

O 35

A

1007-2861(2016)05-0560-13

10.3969/j.issn.1007-2861.2015.02.015

2015-06-05

國家自然科學基金資助項目(11372175)

周哲瑋(1950—),男,教授,博士生導師,博士,研究方向為流體穩定性理論及其應用. E-mail:zhwzhou@shu.edu.cn