離子雙電層動態(tài)輸運特性及電場對液體油膜壓力的影響

馬瑞珂 梅冰昂 左正興

(北京理工大學機械與車輛學院 北京 100081)

隨著潤滑技術進入納米級領域，薄膜潤滑概念被提出。薄膜潤滑的潤滑膜膜厚范圍在幾十納米之間，在這種極小間隙的摩擦中，雙電層效應對潤滑膜的影響不可忽略。在潤滑過程中，摩擦副相對速度是變化的，所以需要對雙電層潤滑進行動態(tài)仿真。

在雙電層機制的研究中，HELMHOLTZ[1]首先研究雙電層現(xiàn)象，并提出平板雙電層模型，該模型認為雙電層類似于平板電容器，其正負離子分別在固液界面兩側整齊排列。但是在雙電層中出現(xiàn)的動電效應該模型無法解釋。于是GOUY[2]和CHAPMAN[3]提出了擴散雙電層模型，對平板雙電層模型進行改進。該模型認為溶液中的反離子以Boltsmann分布規(guī)律分布在液相中，靠近固液界面處離子濃度高，遠離固液界面時濃度降低。擴散雙電層模型能夠解釋雙電層的動電效應，但是無法解釋雙電層中電勢變號的問題。于是STERN[4]對擴散雙電層模型進行改進，考慮離子尺寸有限，在固液界面處緊貼一層反離子與特性吸附離子，能較好地反映雙電層的真實結構。Stern模型的雙電層由緊密層和擴散層兩部分組成，在固液滑移處與溶液本體之間的電勢差被稱為Zeta電勢。GRAHAME[5]對Stern模型進行改進，將Stern層分為內Helmholtz面和外Helmholtz面，其中內Helmholtz面包含特性吸附離子，外Helmholtz面包含表面吸附反離子，為擴散層起始面。此后，許多學者對以Stern模型和Poisson-Boltzmann理論為基礎的雙電層模型進行改進[6-11]。雙電層理論發(fā)展至今，主要應用于電化學儲能領域[12]，在潤滑領域應用尚不成熟。

針對潤滑問題中的雙電層研究中，BIKE和PRIEVE[13]在潤滑分析中考慮了雙電層的影響，在他們的模型中，外加電場使?jié)櫥瑓^(qū)內形成雙電層，雙電層和潤滑液流動產生的動電效應在固體上產生電動力，提高摩擦副承載能力，其結果表明，在其他條件相同時，更高的Zeta電勢可以提供更高的承載能力。ZHANG和UMEHARA[14]在研究工程陶瓷摩擦中考慮雙電層的影響，推導了考慮雙電層效應影響的Reynolds方程，得到了表觀黏度計算公式，并通過試驗驗證了其理論模型與試驗數據基本吻合。黃平、王新杰等[15-16]對雙電層對流體動力潤滑和彈流潤滑狀態(tài)參數的影響進行了相關的理論分析。結果表明：在膜厚較薄(幾十nm)時，雙電層的存在對潤滑膜的厚度有較明顯的影響，可以使最小膜的厚度明顯增加，隨著膜厚的增加這種影響將迅速減弱；雙電層的存在對潤滑膜的壓力影響不大，即使在膜厚很薄時，考慮和不考慮雙電層得到的最大壓力變化依然很小；雙電層產生了電黏度現(xiàn)象，隨著電場強度增加，電黏度現(xiàn)象加強，當電場強度超過某一值時電黏度現(xiàn)象減弱。白少先等[17-19]對雙電層的電黏度效應進行研究，改進了ZHANG和UMEHARA[14]提出的表觀黏度公式，提高了表觀黏度計算值的可靠性與準確性。其研究結果表明，雙電層對薄膜潤滑的有效黏度和摩擦因數影響顯著，使?jié)櫥ゐざ群湍Σ烈驍得黠@增加。左啟陽[20]對電解質的濃度與溫度對雙電層的影響進行了理論分析與數值分析，得出溫度升高使雙電層影響降低，而電解質濃度增加使雙電層效應增強。

以上研究適用于穩(wěn)態(tài)數值分析和試驗驗證，而未考慮瞬態(tài)變工況下的離子輸運與潤滑特性。針對該問題，本文作者提出雙電層潤滑動態(tài)多場耦合模型，考慮離子輸運、流體流動與電場變化特性分析電壓和摩擦副相對運動速度對負載變化的影響。

1 理論模型

1.1 仿真區(qū)域與模型假設

文中采用Poisson-Nernst-Plank模型(PNP模型)探究離子輸運特性，采用Navier-Stokes方程探究流體運動特性。

圖1所示為仿真區(qū)域，上下邊界為壁面，左右兩側為潤滑液入口與出口。上壁面固定，下壁面沿x方向移動，帶動潤滑液流動。摩擦副兩壁面施加不同的電位使?jié)櫥瑓^(qū)存在電場，兩壁面形成雙電層。為簡化計算，模型有如下假設：(1)潤滑液中含有2種離子，其離子價相反；(2)2種離子的擴散系數相同，離子有效直徑相等；(3)潤滑區(qū)內溫度恒定，不考慮產熱傳熱；(4)無滑移邊界；(5)潤滑液為不可壓縮流體； (6)雙電層只考慮擴散層，不考慮緊密層電勢降低，故外加電勢等于雙電層滑移面電勢(Zeta電勢ζ)。

圖1 雙電層潤滑模型Fig 1 Electric double layer lubrication model

1.2 控制方程

離子的連續(xù)性方程[21]為

(1)

式中:ci(i=1,2)為i組分量濃度；Ji為i組分離子質量傳遞摩爾通量，其公式[22]為

(2)

式中:u為流體微團速度矢量；Di為i組分擴散系數；zi為i組分離子的代數價；F為法拉第常數；R為理想氣體常數；T為熱力學溫度；ψ為電勢；NA為阿伏伽德羅常數；a為離子直徑。

將式(2)代入式(1)得離子傳質微分方程

(3)

電場中的泊松方程[23]為

(4)

式中：ε0為真空介電常數；εr為相對介電常數；ρe為空間電荷密度。

溶液中空間電荷密度為

(5)

將式(5)代入式(4)得潤滑液中的泊松方程[24]

(6)

流場的動量方程[21]為

(7)

式中:ρ為溶液密度;p為壓力;μ為動力黏度；f為體積力。

潤滑區(qū)出口處流動電流[18]為

(8)

式中：h0為潤滑區(qū)出口處膜厚；ux為流場速度x方向分量；r為潤滑區(qū)寬度。

潤滑區(qū)膜厚h為

h=ωx+h0+Bω

(9)

式中：ω為上壁面斜率。

潤滑區(qū)電導L[18]為

(10)

式中：λ為潤滑液電導率。

在流動達到穩(wěn)定狀態(tài)時，流動電流等于傳導電流，即Is=Ic。假設流動電勢引起的電場為勻強電場，且只考慮其x方向的分量，則根據歐姆定律，x方向流動電勢引起的場強Ex為

(11)

在該模型中忽略重力，體積力只有電場力，流體微團的體積力為

(12)

式中：k為x方向單位向量;等式右邊第一項為外加電勢引起的電場力，第二項為流動電勢引起的電場力。

將式(12)代入式(7)中得

(13)

流體流動的連續(xù)性方程為

(14)

1.3 初始條件和邊界條件

為求解控制方程組式(3)、式(6)、式(13)、式(14)，每個方程都需要設置初始條件與邊界條件。

其中，初始速度u=0，初始壓力p=0，初始離子量濃度c0,i=1 mol/m3，初始電勢ψ0=0。

入口區(qū)(x=0)邊界條件為

p=0

(15)

ce,i=1 mol/m3

(16)

(17)

出口區(qū)(x=105nm)邊界條件為

p=0

(18)

(19)

(20)

下壁面(y=0)處電勢為ζ1，切向速度為v0，法向離子通量為0

n·Ji=0

(21)

上壁面(y=ωx+h0-Bω)處電勢為ζ2,無切向速度，法向離子通量為0

n·Ji=0

(22)

1.4 參數取值

為數值求解上述控制方程組及其初始條件與邊界條件，一共需要15個參數，包括離子擴散系數、溫度、相對介電常數、溶液黏度、離子直徑、潤滑液電導率、陰離子與陽離子化學價、下壁面速度、上下壁面電勢、入口處離子濃度、出口處膜厚、仿真區(qū)域長度、上壁面斜率。其取值或取值范圍如表1所示。

表1 參數取值

1.5 仿真計算

文中使用COMSOL商用軟件進行仿真，使用稀物質傳遞模塊、層流模塊和靜電模塊3個物理場接口耦合，得出仿真區(qū)域速度場、壓力場、離子濃度場、電勢場的分布。仿真區(qū)域主體采用自由三角形網格，網格最大允許尺寸為20 nm，仿真結果具有網格無關性。在靠近固-液界面處設置邊界層網格，其初始尺寸為1.5 nm，網格增長率為1.02，可以精確反映固-液界面處物理量變化。

2 結果與討論

2.1 PNP模型驗證

圖2所示為潤滑區(qū)中段(x=5×104nm)電勢沿y方向的分布曲線，其中上下壁面雙電層Zeta電勢差為0.1 V?？芍?，電勢變化發(fā)生在2個雙電層中，2個雙電層之間的潤滑液處于靜電平衡狀態(tài)。當下壁面速度接近0時(v0=0.1 m/s)，PNP模型與經典的PB模型結果一致。在v0=0.1 m/s時，雙電層電勢分布沿y方向對稱分布，這是由于速度趨近于0時，2種離子出入口處摩爾流量大致相等，使?jié)櫥瑓^(qū)內陰陽離子濃度相等，上下壁面雙電層電荷密度分布一致，電勢分布同Possion-Boltsmann模型的電勢分布一致，由此可以驗證模型的可靠性。而當壁面速度不為0時，上、下壁面電勢分布不對稱，并且當v0增大時，不對稱程度增加。

圖2 不同模型的電勢沿y方向分布(x=5×104 nm)Fig 2 Potential distribution along the y direction ofdifferent models(x=5×104 nm)

為闡明電勢不對稱性的原因，圖3展示了入口處與出口處陰陽離子的流量隨壁面速度變化的規(guī)律?？梢?，當上壁面固定、下壁面向右移動時，靠近下壁面的溶液流動速度遠大于靠近上壁面的溶液流動速度。由于存在外加電勢，上壁面雙電層中大量聚集陽離子，下壁面雙電層大量聚集陰離子，離子隨溶液流動，垂直于流動方向的截面上陰離子流量始終大于陽離子流量。隨著壁面速度增大，入口處2種離子流量大致相等，而在出口處陽離子摩爾流量小于陰離子摩爾流量，使陽離子在潤滑區(qū)內積累，其潤滑區(qū)內平均濃度大于陰離子平均濃度。

圖3 離子出入口摩爾流量對比Fig 3 Comparison of ion molar flow at inlet and outlet(a)flow at inlet;(b) flow at outlet

進而，圖4展示了不同壁面速度下陰陽離子濃度沿y方向的分布曲線。可知，當下壁面速度接近0(v0=0.1 m/s) 時，上下壁面雙電層反離子量濃度大致相等，當下壁面速度較大(v0=1.1 m/s)時,上壁面陽離子量濃度(8 mol/m3)明顯高于下壁面陰離子量濃度(3 mol/m3)。根據泊松方程(見式(4))，電勢的梯度僅受電荷密度影響，而空間電荷密度由溶液中離子濃度和離子電荷量決定，文中模型中陰、陽離子電荷數相同，故而電勢的梯度僅受離子濃度的影響。下壁面雙電層中電荷密度為負值，電勢的二階導數為正值，電勢沿y方向分布的函數為凹函數，負離子濃度增大時，該雙電層電勢的二階導數增大，下凹程度增大，電勢變化更快；上壁面雙電層中電荷密度為正值，電勢的二階導數為負值，電勢沿y方向分布的函數為凸函數，正離子濃度增大時，該雙電層電勢的二階導數降低，上凸程度增大，電勢變化更快。故得出雙電層中反離子濃度越大，電勢變化幅度越大。上述分析中，上壁面雙電層反離子濃度大于下壁面雙電層反離子濃度，所以上壁面雙電層電勢變化幅度大于下壁面電勢變化幅度，液相平衡電勢偏向下壁面Zeta電勢，形成圖2中所示y方向電勢分布不對稱現(xiàn)象。

圖4 不同速度下陰陽離子沿y方向分布Fig 4 The distribution of anions and cations along the y directionat different speeds (a) x=5×104 nm,v0=0.1 m/s;(b) x=5×104 nm,v0=1.1 m/s

2.2 壓力與承載能力

圖5 不同Zeta電勢下壓力沿x方向分布Fig 5 Pressure distribution along the x directionunder different Zeta potentials

圖6所示為不同壁面運動速度下最大承載能力隨電勢差的變化趨勢?？芍?，承載能力隨Zeta電勢的增加而先增大后減小。這是由于雙電層的動電效應引起表觀黏度變化，在潤滑區(qū)流體流動時，雙電層反離子隨流體運動，形成流動電流；反離子沿同一方向運動，沿流動方向積累，產生流動電勢；流動電勢阻礙反離子沿流體流動方向的運動，表現(xiàn)為阻礙了潤滑區(qū)內流體的流動，降低了流速，表現(xiàn)為黏度提高，在摩擦副相對速度不變的情況下，承載能力得到提高，在假設潤滑膜膜厚不變的情況下，最大壓力也相對于無雙電層時增加。而最大壓力和承載能力降低的原因是雙電層疊加引起電黏度降低[20]。

圖6 不同壁面運動速度下雙電層Zeta電勢對承載能力的影響Fig 6 The influence of Zeta potential of electric double layeron bearing capacity at different wall speeds

為闡述雙電層疊加現(xiàn)象，圖7截取x=8.4×104nm處電勢沿y方向分布?？芍妱莶?ζ1-ζ2)增大使雙電層厚度增加，電勢差增大到0.3 V時電勢變化趨于線性變化，雙電層出現(xiàn)疊加現(xiàn)象;當電勢差繼續(xù)增大時，雙電層疊加程度增大。雙電層的存在可以使?jié)櫥け碛^黏度提高，但當雙電層產生疊加時，由左啟陽[20]的結論得雙電層疊加使雙電層電黏度效應降低。

圖7 不同Zeta電勢下電勢沿y方向分布Fig 7 The distribution of potential along y directionunder different Zeta potentials

3 結論

(1)為解決雙電層動態(tài)仿真問題，綜合考慮離子輸運特性與流體流動及電場效應耦合作用機制，對雙電層潤滑進行動態(tài)仿真，通過將文中模型與經典Poisson-Boltsmann模型對比，分析電勢分布與電荷密度分布，驗證了模型結果的可靠性，為雙電層潤滑變工況的分析研究提供理論基礎。

(2)探究摩擦副相對運動速度對潤滑區(qū)域電勢分布的影響，結果表明，摩擦副相對運動速度對雙電層電勢分布有顯著影響，滑動壁的雙電層電勢變化比固定壁的雙電層電勢變化更大，隨著相對速度增加，兩個雙電層電勢變化程度的差距增大。

(3)探究雙電層Zeta電勢對雙電層壓力分布和承載能力大小的影響，結果表明，雙電層影響潤滑膜壓力和承載能力，提高上下壁面電勢差，壓力和承載能力先增加后減少，這是由于雙電層厚度增加使電黏度增加，但是雙電層過厚產生疊加時電黏度又降低。