納米通道粗糙內壁對流體流動行為的影響*

梅濤 陳占秀 楊歷 王坤 苗瑞燦

(河北工業大學能源與環境工程學院,天津 300401)

納米流動系統具有高效、經濟等優勢,在眾多領域具有廣泛的應用前景.因該類系統具有極高的表面積體積比,致使界面滑移效應對流動具有顯著影響.本文采用分子動力學方法以兩無限大平行非對稱壁面組成的Poiseuille流動為對象,分析了壁面粗糙度與潤濕性變化對通道內流體流動的影響.對于不同結構類型的壁面,需要通過水動力位置來確定固液界面位置,準確計算固液界面位置有助于更好地分析界面滑移效應.研究結果表明,上下壁面不對稱會引起通道內流場參數分布的不對稱,壁面粗糙度及潤濕性的變化會影響近壁面附近流體原子的流動特性,由于壁面凹槽的存在,粗糙壁面附近的數密度分布低于光滑壁面一側.壁面粗糙度及潤濕性的變化會影響固液界面位置,肋高變化及壁面潤濕性對通道中速度分布影響較大,界面滑移速度及滑移長度隨肋高和潤濕性的增大而減小;肋間距變化對通道內流體流動影響較小,界面滑移速度和滑移長度基本保持恒定.

1 引 言

對納米流體的深入研究是科技微型化導致的必然趨勢,納米流體具有優于常規介質的特性,包括運輸、傳熱、吸附等性能,對納米流體的深入研究將更有利于科技微型化的發展.納米量級流體在流動時存在許多微尺度效應,包括分子力效應、低雷諾數效應以及表面張力效應等.相比于宏觀尺度,微尺度條件下對流體流動規律的研究受實驗條件和測量精度的限制,使得數值模擬成為了微尺度流動的主要研究方法.而以微觀原子之間的作用模型為基礎的分子動力學,因不引入常規假設,能以原子級精度描述流體在近鄰壁面處的流動特性,是目前研究微尺度流動問題的主要數值模擬方法.

通常宏觀尺度下流體在通道內的流動特性都是基于無滑移邊界條件進行研究,因為壁面產生的表面效應對近壁處流體流動特性影響較小.然而,隨著流動通道尺寸的減小,比表面積急劇增加,微尺度下固體壁面對近壁處流體的影響直接關系到通道內整體的流動特性分布[1-3].Sun和Ebner[4]研究了固液界面間作用強度大小對流體流動特性的影響,發現當固液原子間作用力較強時,流體在壁面處幾乎無滑移,而當作用力較弱時,流體在壁面處會產生較大的邊界滑移.Voronov等[5]利用分子動力學對Couette流動進行模擬,結果顯示滑移長度與接觸角相關,對于疏水性表面,接觸角會增大,滑移長度也會隨之增大,但也會出現隨著接觸角增大滑移長度減小的情況.胡海豹等[6]利用分子動力學對超疏水納米通道內流體流動特性進行研究,結果表明當接觸角大于150°時,滑移速度和滑移長度出現隨接觸角增大而減小的反常現象,并進一步證明了改變固液原子間的勢能參數表征的潤濕性不能用來準確表示超疏水壁面對流體的影響.Nagayama和Cheng[7]對納米通道內的Poiseuille流動進行分析,發現改變壁面與流體間的勢能參數以及添加在流體部分的驅動力均會影響固液界面間流體原子的運動特性,固液間相互作用越強,邊界滑移速度越小,而添加的驅動力越大,滑移速度也越大.Barisik和Beskok[8]以及Shi等[9]利用智能壁分子動力學方法模擬研究了納米通道內氣體原子的流動特性,結果表明近壁區域內氣體的速度、密度和壓力的變化趨勢僅由壁面力場決定,不受通道高度以及密度的影響.張冉等[10]利用分子動力學分析了納米通道內近壁區域流體的流動特性,同樣發現近壁區域的氣體流動特性僅由壁面力場決定,壁面對氣體原子的勢能作用越強,氣體在近壁區域的密度越大,直至形成吸附層.Voronov等[5]利用分子動力學方法對Couette流動中流體的滑移現象進行了研究,發現固液間勢能強度以及平衡距離的減小均傾向于增大接觸角,而減小勢能強度會增加滑移長度,減小平衡距離會減少滑移長度,說明接觸角與滑移長度間的關系并不唯一.Cieplak[11]同樣利用分子動力學方法對Couette流動進行研究,主要探究了固液間的作用強度以及不同流體介質對流體滑移的影響,結果表明,滑移長度與固液間的作用有直接關系,而與流體介質無關.

可見,目前科研工作者大多集中于進行固液間相互作用等因素對通道內流體流動影響規律的研究,且固體表面多是光滑壁面.實際上,任何固體表面都不可能是絕對光滑的,當宏觀尺度轉為納微尺度,流動通道尺度急劇減少,比表面積也隨之急劇增加,固體表面的粗糙程度對流體流動的作用也相應凸顯.因此,粗糙壁面對流動的影響機理研究已逐漸成為當前納微尺度傳熱傳質研究中的重點[12].張程賓等[13]利用分子動力學方法對含粗糙壁面納米通道內的流體流動進行了研究,發現粗糙壁面會限制近壁區流體原子的運動,導致流體流動速度及滑移速度降低.Rahmatipour等[14]利用分子動力學對含粗糙壁面的納米通道內的肋高變化進行了研究,結果顯示,相比于光滑壁面,隨著肋高的增大,壁面附近流體的密度波動范圍逐漸增大,但波動的峰值均低于光滑壁面.Toghraie等[15]對含納米顆粒的粗糙通道內的流體流動特性進行了研究,通過設定肋高與肋間距的比值來獲取不同壁面粗糙度,結果表明通道內流體的溫度及速度分布并不隨著粗糙度的增大而增大,而是由凹槽內限制的流體原子數量決定,限制的流體原子數量數量越多,流體流動速度越小,反之流動速度越大.但由該文獻的物理模型可知,在研究粗糙度變化時,粗糙壁面上肋的尺寸及數量也發生了變化,故不能很好地詮釋壁面肋間距及肋高變化對近壁區流體的影響.

綜上所述,目前研究壁面的粗糙度變化對流動的影響還不夠詳細,而針對粗糙壁面潤濕性變化的研究更是較少.為深入揭示納微尺度下粗糙壁面對流體流動的影響機理,本文擬構建含粗糙壁面的納米通道內流體流動的分子動力學模型,并與光滑壁面進行對比,分析粗糙壁面肋高及肋間距對流體流動特性的影響.在此基礎上,討論壁面潤濕性對壁面附近流體原子的影響,并揭示粗糙壁面與光滑壁面潤濕性變化的異同.

2 物理模型及模擬設置

納米通道流體流動的分子動力學模型如圖1(a)所示.模擬通道尺寸為Lx×Ly×Lz= 18.20σ× 19.68σ× 7.87σ,納米通道上壁為光滑壁面,下壁為粗糙壁面,壁面由銅原子構成,流體單原子氬均勻分布于上下壁面之間,且均按照面心立方結構進行排列.圖1(b)為納米結構示意圖,壁面長度為Lx= 18.20σ,壁面厚度D= 2.95σ,粗糙壁面肋間距為a,肋高為h,肋長度為b,通道整體高度為H,光滑壁面與切面間為平直流道,高度為H′.計算統計時,對通道高度H進行分層,其中統計流場中數密度分布時分為70層,對速度分布進行統計時分為44層.整個模擬系統在x,z方向上設置周期性邊界條件,在y方向上設置固定邊界條件.

流體原子間以及固液原子間的勢能作用均采用Lennard-Jones (12-6)勢能模型,其表達式為

式中rij為原子間的距離,ε為能量參數,σ為尺寸參數,rc為截斷半徑.流體氬原子間的能量參數ε=1.67×10-21J,尺寸參數σ= 0.3405 nm.壁面固體原子勢能參數εs=40ε,尺寸參數σs=0.69σ,截斷半徑取rc=2.5σ.

固液原子間的作用力大小決定了壁面潤濕性,流體在壁面鋪展得越充分,固液界面間的接觸角越小,可認為潤濕性越好.而對于光滑壁面而言,接觸角的理論公式為

圖1 (a)模擬系統圖;(b)納米結構示意圖Fig.1.(a) Simulation system;(b) schematic of nanostructure.

式中θ為光滑壁面接觸角,εsl為固液原子間的能量參數.而對于粗糙壁面,(2)式將不再適用,針對本文所研究的物理模型,粗糙壁面潤濕性可表示為

式中r表示粗糙壁面的粗糙程度;a,b及h為圖1(b)中所示參數.令,c為固液原子間的勢能系數.本文模擬工況與對應粗糙度下的接觸角如表1所列.

對于壁面原子與流體原子間的尺寸參數σsl,由Lorentz-Berthelot混合法進行計算：

平板間的Poiseuille流動是通過在x方向上施加外部驅動力來驅動流體流動的,并利用拋物線求解Navier-Stokes方程.Poiseuille流動被定義為在通道高度H(y方向上由0到H)內的流動,Navier-Stokes方程可簡化為[16]

式中μ為剪切黏度,ux為x方向上的流動速度,ρ為流體密度.對(6)式連續進行兩次積分,將中心對稱條件和邊界滑移條件ux|y=0=us或ux|y=H=us依次代入,其中us為邊界滑移速度,Poiseuille流動的速度場可表示為

固液界面間的滑移規律可以由滑移長度ls表示,表達式為[17-20]

式中 (?ux/?y)|y|=H為固液界面處流體的速度梯度.對于粗糙壁面滑移長度的確定,需要根據水動力位置來研究[21].圖2(a)和圖2(b)所示分別為Couette流動和Poiseuille流動示意圖,圖2(a)中yC為外推線性的速度直線達到壁面速度時對應的位置,圖2(b)中yP為外推拋物線的速度曲線達到壁面速度的對應位置.通過文獻[21]可知,粗糙壁面滑移長度可表示為,而水動力位置可表示為yH=yC+ls.

表1 不同模擬工況下對應的粗糙度與接觸角Table 1.Corresponding roughness and contact angle under different simulation conditions.

圖2 模型結構示意圖 (a) Couette流動;(b) Poiseuille流動Fig.2.Schematic of nanostructure：(a) Couette flow;(b) Poiseuille flow.

本文模擬過程中,改變壁面粗糙度會影響通道內流體密度與壓力的變化,但變化幅度較小,通道中流體原子的密度基本保持在(1.23 ± 0.01) g/cm3范圍內,造成的計算誤差可以忽略.粗糙壁面肋的長度b= 2.7σ,肋的數量為3.采用Velocity-Verlet算法對流體原子運動方程進行求解,時間步長為1 fs.首先利用正則系綜對初始模型進行弛豫,利用Nose-Hoover恒溫算法將系統溫度恒定在T=0.827ε/kB,流體原子的速度遵循高斯分布,進行50萬步使系統達到穩定.弛豫后對模型入口處厚度為2.0σ的流體區域施加沿x正方向的水平驅動力F= 0.05ε/σ,采用正則系綜對系統進行溫度控制,共運行400萬步,前200萬步使系統達到穩定,后200萬步對流場中所需參數進行計算統計.分析模型水動力位置時需要計算Couette流動的速度場分布,設置上下壁面在水平位置上以相同的速度反向運動,速度大小為.模擬采用LAMMPS程序編寫[22].

3 模擬方法驗證

為了驗證模型及參數設置的正確性,本文基于兩平行光滑平板間的Poiseuille流動,對不同壁面潤濕性下流體在通道內的密度分布進行驗證.具體參數設置如下：流體采用單原子氬,其中ε=1.67×10-21J,σ= 0.3405 nm;固液原子間的距離參數σsl=0.2872 nm,無量綱化后,取固液原子間的勢能系數c= 2.0,1.0,0.6,0.2;上下壁面設置為固體剛性壁面,截斷半徑rc= 2.5σ.

圖3所示為不同固液原子間勢能系數c對應的微通道內流體的無量綱密度分布.由圖3可知,模擬數據與文獻[23]所給數據基本相符,在近壁面區域內,由于表面效應的存在使得流體密度分布均出現了有序振蕩現象,流體密度不均勻;而在通道中心的主流區,流體受壁面的影響較小而趨近于平緩.因此,可認為本文建立理論模型、選用算法及編寫的程序準確可靠.

4 結果與分析

4.1 壁面粗糙度變化對流動特性的影響

本節以固液原子間勢能系數c= 0.75,對應的光滑接觸角θ= 60°時的工況為例.通過改變粗糙壁面上肋高及肋間距研究壁面粗糙度對通道內流體數密度及速度分布的影響.

4.1.1 數密度分布

為研究肋高變化對近壁區流體數密度的影響,固定肋間距a= 3.6σ不變,分別取肋高h=0.5σ,1.0σ,1.5σ,2.0σ,研究微通道內流體流動數密度分布曲線規律.如圖4所示,橫坐標為沿y方向的高度,縱坐標為無量綱數密度ρ?=ρσ3.由圖4可知,由于壁面效應使得近壁區流體原子分布不均勻,近壁區域流體數密度分布出現明顯的振蕩現象,而通道主流區域流體受壁面影響較小,數密度分布基本保持恒定.由于通道壁面形狀的不對稱,導致流場中數密度分布的不對稱,由圖4(a)和圖4(b)所示,光滑壁面附近流體的數密度波動幅度要大于粗糙壁面,且呈現逐漸衰減趨勢.通過改變肋高h來改變壁面粗糙度,實現不同粗糙表面構造.圖4(a)結果顯示,當h較小時,近壁區流體數密度波動呈現逐漸衰減趨勢,但隨著肋高的增大,近壁區流體數密度分布出現一次回升現象,這是由于凹槽內的壁面與切面處的壁面對流體均有影響,導致數密度分布呈現衰弱、增強、再衰弱的趨勢.為研究肋間距a變化對通道內流體數密度分布的影響,固定肋高h=2.0σ不變,分別取長度a=2.7σ,3.6σ,4.5σ,5.4σ,結果如圖5所示.由圖5(a)和圖5(b)可知,納米通道內流體在不同肋間距下數密度的振蕩周期一致,肋間距a變化對近壁區數密度影響較小.因此,可認為肋間距a的變化基本不影響壁面粗糙度對流體數密度分布.

圖3 不同勢能系數c下流體沿z方向的密度分布 (a)c = 2.0;(b)c = 1.0;(c)c = 0.6;(d)c = 0.2Fig.3.Density profiles in thez-direction with different energy coefficientc：(a)c = 2.0;(b)c = 1.0;(c)c = 0.6;(d)c = 0.2.

圖4 肋高h對壁面附近流體數密度分布的影響 (a)粗糙壁面;(b)光滑壁面Fig.4.Effect of rib heighthon the distribution of fluid number density near wall surface：(a) Rough wall surface;(b) smooth wall surface.

4.1.2 速度分布

圖5 肋間距a對壁面附近流體數密度分布的影響 (a)粗糙壁面;(b)光滑壁面Fig.5.Effect of rib spacingaon the distribution of fluid number density near wall surface：(a) Rough wall surface;(b) smooth wall surface.

圖6 不同肋高h下流體沿y方向的速度分布 (a) Couette流動;(b) Poiseuille流動Fig.6.Velocity profiles in they-direction with different rib heighth：(a) Couette flow;(b) Poiseuille flow.

壁面粗糙度變化引起的通道內流體數密度的變化,也會導致通道內流體流速分布的變化.為確定模型水動力位置,分別計算Couette流動及Poiseuille流動時通道內的速度分布,更深一步地分析邊界滑移機理.圖6(a)和圖6(b)為不同肋高對兩種流動狀態下速度分布的影響,橫坐標為沿y方向上的通道高度分布,縱坐標為無量綱速度u?=u/(ε/m)1/2.圖6結果顯示,微通道內流體速度分布呈現不對稱性,壁面粗糙度的存在使近壁區流體剪切產生了額外黏性耗散,其次粗糙壁面凹槽內流體原子的運動受到限制,難以從凹槽內逃脫進布差異越明顯.通過兩種流動狀態下的速度分布可推算出固液邊界位置及滑移長度,本文主要分析Poiseuille流動時的滑移機理,通過施加不同外部驅動力F= 0.03ε/σ,F= 0.05ε/σ,F= 0.07ε/σ來獲取滑移長度計算的標準差η,并取不同驅動力下滑移長度的平均值作為分析結果.圖7所示為不同肋高h下滑移長度標準差分布.圖8所示為不同肋高h對Poiseuille流動中通道內滑移速度及滑移長度的影響.分析結果顯示,光滑壁面一側的滑移速度及滑移長度要高于粗糙壁面一側,隨著肋高h的增大,粗糙壁面一側滑移速度及滑移長度逐漸減小,與Schmatko等[24]研究粗糙高度對界面處滑移速度的影響結論一致.圖9為肋間距a變化對兩種流動狀態下速度分布的影響.由圖9可知,兩種流動狀態下肋間距a的變化對流場中速度分布影響較小,這是由于肋間距a變化不會產生額外的入主流區域,導致Couette流動中壁面更容易帶動流體運動,速度大小分布高于光滑壁面一側,而在Poiseuille流動中則剛好相反,流體的運動受到壁面的阻礙,速度大小分布低于光滑壁面一側.隨著肋高h的增大,凹槽內限制的流體原子數量增多,產生的黏性耗散越大,不同結構壁面附近的速度分黏性耗散,且凹槽限制流體原子的數量不會產生太大差異.圖10為不同肋間距a下滑移長度標準差分布.圖11(a)和圖11(b)為肋間距a變化對Poiseuille流動時通道內滑移效應的影響,可以看出肋間距a變化對滑移速度及滑移長度的影響較小.通過分析壁面粗糙度對速度分布的影響,發現肋高h的變化會影響水動力位置的變化,對邊界滑移影響較大;而肋間距a變化對通道內速度分布影響較小,水動力位置變化不大,故對邊界滑移影響較小.

圖7 不同肋高h下滑移長度標準差分布Fig.7.Standard deviation distribution of slip length with different rib heighth.

圖8 (a)肋高h對滑移長度的影響;(b)肋高h對滑移速度的影響Fig.8.(a) Effect of rib heighthon the slip length;(b) effect of rib heighthon the slip velocity.

圖9 不同肋間距a下流體沿y方向的速度分布 (a) Couette流動;(b) Poiseuille流動Fig.9.Velocity profiles in they-direction with different rib spacinga：(a) Couette flow;(b) Poiseuille flow.

圖10 不同肋間距a下滑移長度標準差分布Fig.10.Standard deviation distribution of slip length with different rib spacinga.

4.2 壁面潤濕性變化對流動特性的影響

4.2.1 壁面潤濕性變化對數密度的影響

壁面潤濕性決定了固體壁面與流體間的相互作用,不僅會影響固液界面處動量的傳遞,也會改變近壁區流體原子的分布狀態[25].為研究壁面潤濕性對微通道內流動的影響,固定凹槽高度h=2.0σ,長度a= 3.6σ,對勢能系數c= 1.0,0.75,0.5,0.25下流體在通道內的流動特性分別進行研究.圖12為不同潤濕性下對應的通道內流體數密度分布,由圖12可知,無論是粗糙壁面還是光滑壁面,壁面與流體間的作用力越強,壁面潤濕性越好,壁面附近吸附的流體原子數量越多,近壁區流體的數密度均隨著勢能系數的增大而增大.對于粗糙壁面,當固液界面間的勢能系數c≤ 0.5時,隨著固液間勢能系數c的減小,凹槽內流體原子的數密度有所下降,但下降幅度小于光滑壁面,且凹槽內流體數密度要大于切面附近.而當固液間的勢能系數c= 0.25時,對應的接觸角最大,凹槽內流體原子的數密度明顯減少,如圖12(a)所示,凹槽內的數密度波動峰值要低于切面附近,說明此時凹槽對流體原子幾乎處于排斥狀態.

圖11 (a)肋間距a對滑移長度的影響;(b)肋間距a對滑移速度的影響Fig.11.(a) Effect of rib spacingaon the slip length;(b) effect of rib spacingaon the slip velocity.

圖12 勢能系數c對壁面附近流體數密度分布的影響 (a)粗糙壁面;(b)光滑壁面Fig.12.Effect of energy coefficientcon the distribution of fluid number density near wall surface：(a) Rough wall surface;(b) smooth wall surface.

4.2.2 壁面潤濕性變化對速度分布的影響

圖13 不同勢能系數c下流體沿y方向的速度分布 (a) Couette流動;(b) Poiseuille流動Fig.13.Velocity profiles in they-direction with different energy coefficientc：(a) Couette flow;(b) Poiseuille flow.

圖13為勢能系數c變化對兩種流動狀態下通道內速度分布的影響.結果顯示,隨著勢能系數c的增大,固液原子間的作用力逐漸增強,兩種流動狀態下速度分布呈現相反的變化趨勢.如圖13(a)所示,Couette流動中通道內的速度分布隨著勢能系數c的增大而增大,而圖13(b)顯示Poiseuille流動中通道內的速度分布隨勢能系數c的增大而減小.另外值得一提的是,對于Couette流動來說,粗糙壁面附近流體速度變化幅度要大于光滑壁面一側,而Poiseuille流動則剛好相反.圖14為不同勢能系數c下滑移長度標準差分布.圖15(a)和圖15(b)分別為Poiseuille流動中通道內的滑移速度及滑移長度分布,可以發現,隨著勢能系數c的增大,無論是光滑壁面還是粗糙壁面,滑移速度和滑移長度分布均逐漸降低.通過分析可知,改變勢能系數c會影響通道內的速度分布,使水動力位置發生變化,并對邊界滑移影響較大.

圖14 不同勢能系數c下滑移長度標準差分布Fig.14.Standard deviation distribution of slip length with different energy coefficientc.

圖15 (a)勢能系數c對滑移長度的影響;(b)勢能系數c對滑移速度的影響Fig.15.(a) Effect of energy coefficientcon the slip length;(b) effect of energy coefficientcon the slip velocity.

5 結 論

本文采用分子動力學方法研究了含粗糙壁面納米通道內流體的流動特性,探討了不同壁面粗糙度以及壁面潤濕性對通道內流體的數密度和速度場分布的影響,所得結論如下.

1)相比于光滑壁面,粗糙壁面附近流體的數密度分布較低,隨著肋高的增大,密度波動呈現一次回升現象;改變肋間距對近壁區流體影響較小,數密度波動趨勢基本一致;無論是光滑壁面還是粗糙壁面,增大壁面潤濕性均會使近壁區數密度波動增大.

2)通過分析Couette流動和Poiseuille流動時通道內的速度場分布確定模型的固液邊界位置.分析結果表示,改變肋高及壁面潤濕性會明顯影響通道內的速度分布,使固液邊界位置發生偏離,而肋間距變化對固液邊界位置影響較小.

3)以Poiseuille流動為主,分析了壁面粗糙度及潤濕性對界面滑移的影響.結果表明,粗糙壁面一側滑移速度及滑移長度均小于光滑壁面一側,且隨著肋高及壁面潤濕性的增大,滑移速度和滑移長度逐漸減小;肋間距變化對界面滑移影響較小,滑移速度及滑移長度分布基本恒定.