離子選擇性表面電對流的多塊LB 模擬1)

張 煜 李天富 羅 康 吳 健 易紅亮,2)

* (哈爾濱工業大學能源科學與工程學院,空天熱物理工信部重點實驗室,哈爾濱 150001)

?(季華實驗室,廣東佛山 528000)

引言

近幾十年來,隨著微機電加工技術的迅速發展,基于電流體動力學進行流動控制逐漸成為微尺度下介質操控的一種重要方式.當介電液體與離子選擇性材料(例如納米通道、離子交換膜等)相接觸,在外加電壓的作用下,會產生離子濃差極化(ion concentration polarization,ICP)現象,且在較高電壓下還會產生復雜的對流現象,即所謂的第二類電滲流或離子交換表面的電對流.這一類問題在海水淡化、離子濃度富集、燃料電池等諸多領域具有廣泛的應用前景[1-4].

向帶有離子交換表面的電解質溶液施加電壓,通過交換表面的電流會呈現出復雜的非線性變化過程.圖1 給出了典型離子交換膜系統的電流-電壓特性示意圖.隨著垂直施加在流體與離子交換膜表面的電壓逐漸增大,通過電解質溶液的電流會經歷三個階段: 歐姆區(Ohmic regime),電流先隨著電壓的增大而明顯增大;飽和區(limiting regime),電流不隨電壓發生明顯變化;過飽和區(over-limiting)[5]電壓繼續增大,電流隨電壓再一次出現增大的趨勢.

這些特性與系統內離子輸運機制隨電壓的變化有關[6].在離子交換膜表面施加電壓,異號離子被表面吸收,同號離子因電場力排斥遠離表面,形成雙電層(electric double layer,EDL),如圖2 中曲線I1.隨著外加電壓增大,電場對同性離子的排斥增強,雙電層外的離子濃度降低并形成濃差極化區.極化區的離子濃度主要由遠離表面的離子通過擴散作用進行補充,因此雙電層外部的離子濃度呈線性分布,形成電中性的外部擴散層(diffusion layer,DL),如圖2 中的曲線I2.當外加電壓增大到一定值時,離子膜表面同性離子被耗盡,如圖2 中曲線I3,介電液體表現出很高的電阻,電流幾乎不再隨電壓增加,與圖1 中的飽和區相對應所示.Levich 等[7]以帶有單種陰陽離子的稀溶液為模型進行分析,獲得這一飽和電流理論解Ilim=2zeDC0/L,其中z表示離子價,e為元電荷密度,D為離子的擴散系數,C0為邊界處于穩定狀態溶液的濃度,L為離子交換表面到與主體溶液濃度相等界面的距離.

圖1 離子交換膜表面電流-電壓特性示意圖Fig.1 A schematic voltage-current characteristics near ions selective surface

圖2 離子交換表面的濃差極化特性曲線[5]Fig.2 Ion concentration polarization on ions selective surface[5]

由于離子交換表面的存在,壁面附近同性離子濃度隨電壓的增大最終趨于零,而雙電層外部出現擴展空間電荷(extended space charge,ESC) 層.ESC 層外離子輸運仍保持擴散層特性.已有結果表明,在ESC 層與DL 層之間還存在一個過渡層(transition layer,TL)[6],如圖2 中曲線I4.當外加電壓超過某一臨界值,ESC 層變得不穩定,觸發電對流,離子濃度分布達到新的平衡.

關于離子交換表面的電對流問題研究已久.近年來隨著計算技術的發展,數值模擬成為一種有力的研究手段.多位國內外學者采用諸如譜方法[8]、有限體積[9]、有限元[10]等算法對此類問題進行數值模擬.格子Boltzmann 方法(lattice Boltzmann method,LBM)作為一種介觀方法,自誕生以來,以其易于編程、適合復雜形狀與并行計算、適合多相、多物理場耦合的模擬等優勢,被迅速擴展應用至復雜流動、傳熱傳質、多相及多場耦合流動等領域[11-14].本文采用LBM 數值求解這一問題.

應用LBM 求解電流體力學問題由來已久.最初有學者利用Poisson-Boltzmann (PB)模型求解了微通道內的電滲流問題[15-18].2010 年,Wang 等[19]提出了完全求解Poisson-Nernst-Planck (PNP)模型的電流體力學LBM 求解模型.隨后,不斷有學者提出并改進基于LBM 的PNP 電場求解模型,用于研究單極注入電荷表面電(熱)對流[20]、電場強化相變[21]等電流體問題.

本文提出一種基于多塊網格的LB 模型來模擬離子選擇性表面的電對流問題.采用不同的LB 輸運方程對多物理場進行耦合,并給出了局部加密的多塊LB 信息交換格式.由于本文所考慮的物理問題在邊界雙電層內有很大的濃度梯度,且需要耦合求解離子輸運、流動以及電勢等控制方程,該局部加密的LB 格式適合本問題的求解.隨后,利用發展的LB 方法數值求解Navier-Stokes (N-S) 方程、PNP 方程互相耦合的電流體力學控制方程組,模擬了離子選擇性表面的電對流問題.

本工作有兩個難點: 一是表征離子輸運的NP方程的強對流占優使多個物理場的耦合變得不穩定;二是離子交換表面電荷濃度在很小的范圍內迅速衰減導致的大濃度梯度.針對前者,采用多弛豫時間(multi-relaxtion time,MRT) 模型以保證模型的穩定性;對于后者,采用多塊網格局部加密來克服標準LBM 均勻網格的不足[22].本文的結構安排如下: 第一節對物理及數值模型問題進行介紹;第二節給出了本文提出的多塊LB 電流體力學求解模型;第三節利用本文提出的多塊LB 模型模擬了離子交換表面的電對流問題;最后對文章內容做了總結與展望.

1 數學物理模型

1.1 宏觀控制方程

本文研究的物理模型如圖3 所示.求解區域內充滿電解質溶液,下方覆有離子交換表面并施加電壓.隨著電壓的增大,交換表面極性相同的離子濃度被耗盡,在壁面處積累大量的凈電荷,產生庫侖力,并引起流動.模擬區域中,離子交換表面被認為是無滑移壁面;上邊界為靜止溶液界面,數值處理格式與無滑移壁面一致;左右兩側數值上設為周期性邊界.

圖3 離子選擇性表面的電對流物理模型示意圖Fig.3 Illustration of physical model for electroconvection near ions selective surface

對應的控制方程包括以下N-S 方程組、Poisson方程以及NP 方程組[9]

其中,符號t,u,ρ ,μ,E,Ci,Di,ε,zi,T,φ 分別代表時間、流動速度矢量、流體密度、動力黏度、電場強度、離子濃度、電荷擴散系數、電解質溶液的介電常數、離子價、溫度和電勢,e和kb分別為元電荷量和玻爾茲曼常數,角標i表示第i種離子.本文考慮對稱二元單價離子的電解質溶液,i代表電荷的極性,z+=1,z-=-1.離子選擇性表面速度為無滑移邊界,施加固定負電勢.假設表面為理想離子膜,不考慮陰離子對膜內離子濃度的影響,設置陽離子濃度為定值[23],陰離子通量為零

遠離離子交換表面的上邊界為靜止狀態的電解質溶液,流動速度為零,陰陽離子濃度相等

下面對LB 求解以上方程式(1)～式(3)進行介紹.LBM 的基本建模思想是采用Chapman-Enskog 多尺度分析展開Boltzmann 方程,并滿足一定的宏觀物理量到介觀分布函數的對應關系,使得實際求解的格子離散Boltzmann 方程能夠恢復成原始的宏觀方程.方程求解的基本思路包括碰撞(弛豫迭代)和遷移,采用顯式時間步進的方程來依次求解不同的控制方程并耦合.本文對2D 問題做模擬,下面分別給出不同宏觀方程的具體LB 格式.

1.2 Navier-Stokes 方程組的LB 求解格式

采用單松弛的BGK 方程[24-25]求解控制流動的N-S 方程,方程中的外力項采用Guo 等[26]提出的外力模型

其中fi表示分布函數; δt表示格子單位下的離散時間步長;ci為格子離散速度,采用一般的D2Q9 格式[12];弛豫時間和平衡態分布函數以及外力離散由下式給出

宏觀量密度和速度由下式給出

式(7)和式(8)中,外力F對應于式(1)中的外力項.

1.3 Poisson 方程的LB 求解格式

采用LBM 求解Poisson 方程有多種思路.例如給Poisson 方程添加時間項,將橢圓型方程化為拋物型方程,步進迭代得到的穩態解即為原Poisson 方程的解.本文使用由Chai 等[27]提出的更接近橢圓型方程本質的一種LBM 求解格式,采用D2Q5 離散格式,碰撞項采用單松弛格式

其中 φi表示電勢的介觀分布函數,ci為格子離散速度;人工擴散系數Dφ、弛豫時間 τφ、權重因子?i的表達式如下

其中,α=1/2 .電勢宏觀量計算過程為

權重因子wi=?i=0,(i=0);1/4,(i=1-4) .平衡態分布函數為

關于這一Poisson 方程的求解格式,Liu 等[28]在工作中予以改進,本文采用基礎格式可滿足需求.

1.4 Nernst-Planck 方程的LB 求解格式

本文求解離子濃度分布的Nernst-Planck 方程的LB 格式選用Chai 等[29]提出的對流擴散方程的MRT 格式,電場輸運項的處理與Luo 等[20]一致.采用D2Q9 離散,格子速度離散方向與1.2 節NS 方程的LB 格式一致,具體的多松弛演化方程如下

式中 φk表示該對流擴散方程的介觀分布函數,M為從速度空間到矩空間的轉換矩陣,S為弛豫矩陣,是一個微分算子.平衡態分布函數為

其中

分布函數到矩空間的轉換關系由下式給出

關于弛豫矩陣和轉換矩陣的選取,以及對該數值格式更詳細的介紹可參閱文獻[29].

1.5 LB 邊界條件

流場的無滑移邊界條件采用Zou-He 邊界(非平衡反彈格式)[30],將分布函數分為平衡態和非平衡態兩部分,僅對非平衡態執行反彈.離子濃度和電勢采用Guo 等[31]提出的非平衡外推格式.對離子流通量邊界條件,可通過求解邊界離子濃度微分方程(4),獲得邊界處的離子濃度[32],進而將Robin 邊界轉化為Dirichlet 邊界條件.

2 多塊LB 局部加密算法介紹

如圖2 所示,在離子交換表面附近,離子濃度變化很大,在擴散層中線性分布,這樣劇烈的物理量的變化對計算網格要求很高.針對離子交換邊界附近很大的物理量變化,本文引入多塊網格來加密邊界附近的計算域,以減少網格消耗.

下面以圖4 為例進行說明.AB為細網格邊界,MN為粗網格邊界,模擬宏觀物理量在粗細網格交界面處應該保持一致.對于介觀分布函數,由于LBM算法的實施由碰撞和遷移兩步組成,碰撞過程只涉及當前節點的格子信息,而遷移會用到相鄰格點的分布函數,因此,分布函數在網格交界處的交換應在碰撞之后、遷移之前完成,保證在每個邊界處遷移進入多塊計算域的分布函數信息正確和網格無關性.因此,需要推導適用于三個不同的LB 方程的信息交換格式.

圖4 多塊網格加密示意圖Fig.4 A sketch of multi blocks refinement

下面用下標c表示粗網格,f表示細網格(下同).若粗細網格空間離散步長之比為m=δxc/δxf,為保證分布函數的同步遷移和物理量輸運的一致性,表征分布函數遷移速度的格子速度c=δx/δt以及表征物理量輸運的黏度和擴散系數應保持一致.根據黏度、擴散系數與弛豫時間的關系,以及格子時間步長滿足關系 δtc/δtf=m,細網格執行LBM 碰撞過程的弛豫時間也應該滿足 τf=1/2+m(τc-1/2)[22].至此,每次碰撞的弛豫參數得以保證,下面給出利用多尺度分析得到的信息交換格式.

2.1 流場N-S 方程的LB 網格信息交換格式

對求解N-S 方程的LBM 分布函數做Chapman-Enskog 展開

由于上式左邊外力以及導數項與時空離散格式無關,且平衡態分布函數只由宏觀量決定,而宏觀量應在粗細網格中保持一致,因此粗細網格的分布函數只需滿足關系

將非平衡態引入式(6),格子碰撞過程可以寫為

其中,上標pc代表碰撞后的量,結合式(21) 與式(22),可得

根據式(22),分布函數的非平衡態部分可以寫為

將式(24)代入式(23),可得細網格向粗網格的分布函數信息交換格式

同理,可得粗網格向細網格的信息交換格式

式(26)右側的分布函數應經過時間和空間插值之后,再通過上式從粗網格轉換到細網格.

2.2 電勢Poisson 方程的LB 網格信息交換格式

Poisson 方程的LB 格式已由式(9) 給出,對不同網格尺度的LB 演化方程,應保證格子聲速和人工擴散系數一致,相應的弛豫時間的表達式已由前文給出.對Poisson 方程的分布函數做Chapman-Enskog 展開

再對式(9)左側第一項做泰勒展開并結合式(27),其中 ξ 作為一個小量,可令 ξ=δt[33],有

由于 ?iRDφ-具網格無關性,為保證分布函數在粗細網格中保持一致,只需滿足

又因為

可得細網格到粗網格的信息交換方程

同理,分布函數由粗到細網格的交換格式為

2.3 離子輸運方程的LB 網格信息交換格式

對于多松弛的對流擴散模型,網格加密算法的實施依然應從Chapman-Enskog 展開導出,式(13)的詳細信息及多尺度展開實施方法可參閱文獻[29],這里僅給出LB 方程在矩空間展開到不同階數并化簡之后得到的結果

其中

參考黃榮宗等[34]的MRT 多塊對流擴散LB格式推導,結合對矩空間分布函數與非平衡態矩的關系,有下式成立

式中m(0)即為平衡態分布矩,滿足

將式(39)代入式(36),可得

分析上式,為滿足不同網格之間信息交換的網格無關性,應滿足

結合式(35)和式(38),可得

又因為碰撞后的分布函數矩可以寫為非平衡態格式

式(47)中,在粗網格到細網格的邊界信息傳遞過程中,粗網格邊界的值要先經過一次空間插值,獲得對應細網格所有節點的值,然后還要對矩空間分布函數實施時間插值以獲得信息交換的.

2.4 網格界面插值格式選取

由于LBM 具有時間和空間二階精度,因此網格界面需要選取更高精度的插值格式.本文中時間插值選用三點拉格朗日插值,空間插值采用自然邊界條件的三次樣條插值[35]

其中,式(48)代表三點拉格朗日插值,式(49)為三次樣條插值,系數選取參見式(50),σi代表二階導數,具體的實施可查閱相關參考文獻[35]

以二倍加密系數為例,圖5 給出了多塊網格算法實施過程.

圖5 多塊網格 LBM 算法求解步驟Fig.5 Solving steps of multi blocks LBM

3 離子交換表面的電對流模擬

本節采用前文提出的多塊LB 方法對離子交換表面電對流進行數值模擬.考慮對稱二元電解質氯化鉀(KCl)溶液為工質,并假設正負離子擴散系數相同,其他所選用的物性參數如表1 所示.

表1 模擬選用物性參數Table 1 Physical properties used in simulations

根據物性表1 計算得到,飽和限制電流密度為Ilim=7.15m2/s,德拜長度 λD=13.48 nm,規定特征電壓(熱電壓) φ0=kBT/(ze)= 25.25 mV.

圖6 給出了不同網格加密級數下的網格配置(1:2:2 與1:2:4)示意圖以及部分相應的模擬結果,其中最粗的網格密度為200.從流線可以看出模擬結果在不同網格的加密處有很好的連續性,說明該局部加密算法的可靠性.此外,分別對格子密度為160,200,250 (均代表最粗網格的格子密度,細網格部分按比例增加,下同),模擬結果均有很好的一致性.下文的計算結果均采用200 個格子(以粗網格計)計算得到,計算域設為1:1.

圖6 多塊網格排布示意圖及計算流線分布Fig.6 Schematics of multi blocks grids arrangement

3.1 電流隨電壓的變化

圖7 給出了隨著電壓增加,溶液中電流密度隨電壓的演化特性曲線.可以看出電流密度隨電壓的變化趨勢表現為先增長后趨于飽和,在電壓到達一個閾值后繼續增加.電壓低于 7φ0時,隨著電壓的增加,通過系統的電流密度迅速增大;電壓大于 7φ0時,電壓繼續增大,電流密度的增長幾乎停滯,趨于飽和;當電壓達到 21φ0時,電流密度繼續增大,表現為非線性增長趨勢.

圖7 離子交換表面的電流-電壓(I-V)特性曲線,藍色線條為模擬結果,橘色虛線為根據Yariv 的漸近分析獲得的理論解,紅色虛線為人為抑制流動所獲得的數值解Fig.7 The relationship of voltage-current near the ions selective surface.The blue line represent the numreical result,the orange dotted line is the analytical solution form Yariv's asymptotic analysis,the red dashed line is the numerical solution obtained by suppressing the flow artificially

圖7 中點劃線表示通過人工抑制流動發生從而得到的沒有對流輸運的電流曲線,虛線代表Yariv等[32]通過漸近展開推導得到的一維情況下的電流解析解.一方面,本文的數值結果在有流動和無流動的情況下都與文獻[32]的理論解吻合得很好,定量證明了當前多塊LB 模型的正確性;另一方面,后期電流的再次增大也定性驗證了流動對過飽和電流的貢獻.

3.2 流動未發生時的濃差極化現象

圖8(a)給出了沒有流動發生時系統的正負離子濃度分布特性曲線,其中 25φ0的模擬結果是在人為抑制流動的狀態下計算得到的.在離子交換表面附近,離子濃度迅速衰減,表現為雙電層擴散分布的特性.擴散層的離子濃度呈線性分布,此時離子輸運由擴散機制主導.

圖8(b)給出了離子濃度C=C+-C-的分布特性.電壓達到飽和限制區域之后,離子交換表面附近的負離子被耗盡,形成了ESC 層.ESC 層積累的凈電荷量遠大于EDL 內部,而外部的擴散層仍保持電中性.從圖8(b)中還可以看出,在擴散層與擴展空間電荷層之間還存在一個凈離子濃度增大的過渡層TL,這一模擬結果與相應的理論分析一致.

圖8 流動未發生時的離子濃度分布.(a) 正負離子濃度分布,正負離子濃度只在壁面附近有所不同,僅用點來區分不同電壓下的結果;(b) 凈離子濃度分布Fig.8 The ions concentration distribution near the ions selective surface without convection.(a) Cations and anions concentration.The cations and anions are different near the wall while are the same in DL,so we only use markers to distinguish the results under different voltages.(b) Net charge density

3.3 不同電壓下的流動發展

進一步以 20φ0下的穩態結果為模擬初始條件,增加電壓,獲得不同電壓下的流動演化.在本文所研究的電壓范圍內,離子交換表面的電對流狀態最終都會達到相似的穩態,但是流動隨時間的演化特性會不同,下面對不同外加電壓下的流動演化特性進行分析.

如圖9 所示,在 22φ0的電壓下,系統需要經歷較長時間才能產生流動.模擬初始階段,最大流動速度呈指數增長,由于此時的流動速度較小,不足以引起電流密度的顯著變化,整個系統的離子輸運仍處于擴散機制主導階段,流動形態在初始階段即為對稱的渦結構.隨著流動增強,流動對離子輸運的貢獻逐漸顯現,電流也隨之增大.流動強度達到峰值之后稍有回落,進入電流減小的過渡階段.

圖9 V=22φ0 最大橫向速度與電流密度的演化曲線Fig.9 The maximum lateral velocity and current density development curve when V=22φ0

圖10 給出了系統到達穩態之后的流動形態和離子濃度分布,流場朝內側聚集,離子濃度呈駝峰狀分布.在離子濃度耗盡程度較大的區域,電解質溶液向上流動,攜帶離子遠離表面,表現為對離子輸運的加強,并在兩側匯集到離子選擇性表面,抵消了一部分電遷移作用,表現為局部電流密度的減弱.由于流動強化離子輸運作用更強,離子輸運效率隨之增大,最終導致電流增大.

圖10 V=22φ0 穩態流動強度和和陰離子濃度分布.(a) 流場分布,(b) 陰離子濃度Fig.10 Contour of (a) fluid velocity and (b) anions concentration on steady state with V=22φ0

當電壓增大為 26φ0時,流動隨時間的演化表現出不同的特性.如圖11 所示,相比于V=22φ0,這一電壓下流動發生了兩次明顯的變化,最終達到穩態.較強的電場力在流動開始階段就導致了明顯的擾動,流動在初始階段表現為兩對小渦,強度隨著擾動的積累而增大,最終影響離子濃度分布.如圖12(a)和圖12(b)所示,此時離子耗盡區的濃度分布與流動形態一致.這種小渦流動并不穩定,會隨著時間的演化逐漸過渡到大渦.如圖12(c)和圖12(d)所示,右側的兩個渦逐漸增大,合并了左側的小渦,相應的離子濃度分布也隨流動而改變,對流渦對離子輸運的影響繼續增大.

圖11 V=26φ0 最大橫向速度與電流密度的演化曲線Fig.11 The maximum lateral velocity and current density development curve when V=26φ0

圖12 V=26φ0 不同時刻的流動和離子濃度分布,(a),(b)對應于圖11 中0.02 s,(c),(d)對應于0.03 sFig.12 Contour of velocity and anions concentration at different evolution time when V=26φ0,the subfigures (a),(b) corresponds to 0.02 s and (c),(d) corresponds to 0.03 s in Fig.11

圖12 V=26φ0 不同時刻的流動和離子濃度分布,(a),(b)對應于圖11 中0.02 s,(c),(d)對應于0.03 s (續)Fig.12 Contour of velocity and anions concentration at different evolution time when V=26φ0,the subfigures (a),(b) corresponds to 0.02 s and (c),(d) corresponds to 0.03 s in Fig.11 (continued)

圖11 中,電流密度在流動發生之后也出現了兩次轉變,分別對應于小渦流動增強離子輸運、小渦轉變為大渦流動.兩種流態的最大速度并沒有太大差距,但是大渦狀態下的離子輸運效率更強.該電壓下的最終的穩態流動與圖10 類似,但流動強度明顯提高.

4 結論

針對離子交換表面的電對流問題,本文首先提出了一種多塊LB 電流體求解算法.該方法采用多松弛格式計算電場分布,提高了程序的穩定性;用多塊網格局部加密格式,較好地解決了計算資源在模擬區域不同范圍內的分配問題,可有效處理離子交換表面邊界離子分布的強烈變化.隨后,利用新發展的模型對離子交換表面的濃差極化和流動現象進行了模擬,分析了不同電壓下的流動隨時間的演化特性.模擬結果與漸近理論解吻合良好,捕捉到了不同電壓下的電流特性曲線.在相對低的電壓下,流動傾向于直接形成大渦,而更大的電壓則會在初期就激發小渦流動,并逐漸合并為大渦流動結構.相比于小渦,大渦流動具有更高的離子輸運效率.

本文提出的多塊LB 電流體求解格式給出了三個對應于不同宏觀方程的介觀分布函數在不同精細度網格界面的信息交換格式,也可用于數值求解電荷單極注入電對流、電滲流等其他電場與流動的耦合問題.該多塊局部LB 算法具有良好的可拓展性,后續可在此基礎上發展自適應網格的求解算法[34],以適應更加靈活的電流體動力學問題的求解.