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

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

* (哈爾濱工業(yè)大學(xué)能源科學(xué)與工程學(xué)院,空天熱物理工信部重點(diǎn)實(shí)驗(yàn)室,哈爾濱 150001)

?(季華實(shí)驗(yàn)室,廣東佛山 528000)

引言

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

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

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

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

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

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

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

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

1 數(shù)學(xué)物理模型

1.1 宏觀控制方程

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

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

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

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

遠(yuǎn)離離子交換表面的上邊界為靜止?fàn)顟B(tài)的電解質(zhì)溶液,流動(dòng)速度為零,陰陽離子濃度相等

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

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

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

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

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

式(7)和式(8)中,外力F對應(yīng)于式(1)中的外力項(xiàng).

1.3 Poisson 方程的LB 求解格式

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

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

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

權(quán)重因子wi=?i=0,(i=0);1/4,(i=1-4) .平衡態(tài)分布函數(shù)為

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

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

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

式中 φk表示該對流擴(kuò)散方程的介觀分布函數(shù),M為從速度空間到矩空間的轉(zhuǎn)換矩陣,S為弛豫矩陣,是一個(gè)微分算子.平衡態(tài)分布函數(shù)為

其中

分布函數(shù)到矩空間的轉(zhuǎn)換關(guān)系由下式給出

關(guān)于弛豫矩陣和轉(zhuǎn)換矩陣的選取,以及對該數(shù)值格式更詳細(xì)的介紹可參閱文獻(xiàn)[29].

1.5 LB 邊界條件

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

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

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

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

圖4 多塊網(wǎng)格加密示意圖Fig.4 A sketch of multi blocks refinement

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

2.1 流場N-S 方程的LB 網(wǎng)格信息交換格式

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

由于上式左邊外力以及導(dǎo)數(shù)項(xiàng)與時(shí)空離散格式無關(guān),且平衡態(tài)分布函數(shù)只由宏觀量決定,而宏觀量應(yīng)在粗細(xì)網(wǎng)格中保持一致,因此粗細(xì)網(wǎng)格的分布函數(shù)只需滿足關(guān)系

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

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

根據(jù)式(22),分布函數(shù)的非平衡態(tài)部分可以寫為

將式(24)代入式(23),可得細(xì)網(wǎng)格向粗網(wǎng)格的分布函數(shù)信息交換格式

同理,可得粗網(wǎng)格向細(xì)網(wǎng)格的信息交換格式

式(26)右側(cè)的分布函數(shù)應(yīng)經(jīng)過時(shí)間和空間插值之后,再通過上式從粗網(wǎng)格轉(zhuǎn)換到細(xì)網(wǎng)格.

2.2 電勢Poisson 方程的LB 網(wǎng)格信息交換格式

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

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

由于 ?iRDφ-具網(wǎng)格無關(guān)性,為保證分布函數(shù)在粗細(xì)網(wǎng)格中保持一致,只需滿足

又因?yàn)?/p>

可得細(xì)網(wǎng)格到粗網(wǎng)格的信息交換方程

同理,分布函數(shù)由粗到細(xì)網(wǎng)格的交換格式為

2.3 離子輸運(yùn)方程的LB 網(wǎng)格信息交換格式

對于多松弛的對流擴(kuò)散模型,網(wǎng)格加密算法的實(shí)施依然應(yīng)從Chapman-Enskog 展開導(dǎo)出,式(13)的詳細(xì)信息及多尺度展開實(shí)施方法可參閱文獻(xiàn)[29],這里僅給出LB 方程在矩空間展開到不同階數(shù)并化簡之后得到的結(jié)果

其中

參考黃榮宗等[34]的MRT 多塊對流擴(kuò)散LB格式推導(dǎo),結(jié)合對矩空間分布函數(shù)與非平衡態(tài)矩的關(guān)系,有下式成立

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

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

分析上式,為滿足不同網(wǎng)格之間信息交換的網(wǎng)格無關(guān)性,應(yīng)滿足

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

又因?yàn)榕鲎埠蟮姆植己瘮?shù)矩可以寫為非平衡態(tài)格式

式(47)中,在粗網(wǎng)格到細(xì)網(wǎng)格的邊界信息傳遞過程中,粗網(wǎng)格邊界的值要先經(jīng)過一次空間插值,獲得對應(yīng)細(xì)網(wǎng)格所有節(jié)點(diǎn)的值,然后還要對矩空間分布函數(shù)實(shí)施時(shí)間插值以獲得信息交換的.

2.4 網(wǎng)格界面插值格式選取

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

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

以二倍加密系數(shù)為例,圖5 給出了多塊網(wǎng)格算法實(shí)施過程.

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

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

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

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

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

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

圖6 多塊網(wǎng)格排布示意圖及計(jì)算流線分布Fig.6 Schematics of multi blocks grids arrangement

3.1 電流隨電壓的變化

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

圖7 離子交換表面的電流-電壓(I-V)特性曲線,藍(lán)色線條為模擬結(jié)果,橘色虛線為根據(jù)Yariv 的漸近分析獲得的理論解,紅色虛線為人為抑制流動(dòng)所獲得的數(shù)值解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 中點(diǎn)劃線表示通過人工抑制流動(dòng)發(fā)生從而得到的沒有對流輸運(yùn)的電流曲線,虛線代表Yariv等[32]通過漸近展開推導(dǎo)得到的一維情況下的電流解析解.一方面,本文的數(shù)值結(jié)果在有流動(dòng)和無流動(dòng)的情況下都與文獻(xiàn)[32]的理論解吻合得很好,定量證明了當(dāng)前多塊LB 模型的正確性;另一方面,后期電流的再次增大也定性驗(yàn)證了流動(dòng)對過飽和電流的貢獻(xiàn).

3.2 流動(dòng)未發(fā)生時(shí)的濃差極化現(xiàn)象

圖8(a)給出了沒有流動(dòng)發(fā)生時(shí)系統(tǒng)的正負(fù)離子濃度分布特性曲線,其中 25φ0的模擬結(jié)果是在人為抑制流動(dòng)的狀態(tài)下計(jì)算得到的.在離子交換表面附近,離子濃度迅速衰減,表現(xiàn)為雙電層擴(kuò)散分布的特性.擴(kuò)散層的離子濃度呈線性分布,此時(shí)離子輸運(yùn)由擴(kuò)散機(jī)制主導(dǎo).

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

圖8 流動(dòng)未發(fā)生時(shí)的離子濃度分布.(a) 正負(fù)離子濃度分布,正負(fù)離子濃度只在壁面附近有所不同,僅用點(diǎn)來區(qū)分不同電壓下的結(jié)果;(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 不同電壓下的流動(dòng)發(fā)展

進(jìn)一步以 20φ0下的穩(wěn)態(tài)結(jié)果為模擬初始條件,增加電壓,獲得不同電壓下的流動(dòng)演化.在本文所研究的電壓范圍內(nèi),離子交換表面的電對流狀態(tài)最終都會(huì)達(dá)到相似的穩(wěn)態(tài),但是流動(dòng)隨時(shí)間的演化特性會(huì)不同,下面對不同外加電壓下的流動(dòng)演化特性進(jìn)行分析.

如圖9 所示,在 22φ0的電壓下,系統(tǒng)需要經(jīng)歷較長時(shí)間才能產(chǎn)生流動(dòng).模擬初始階段,最大流動(dòng)速度呈指數(shù)增長,由于此時(shí)的流動(dòng)速度較小,不足以引起電流密度的顯著變化,整個(gè)系統(tǒng)的離子輸運(yùn)仍處于擴(kuò)散機(jī)制主導(dǎo)階段,流動(dòng)形態(tài)在初始階段即為對稱的渦結(jié)構(gòu).隨著流動(dòng)增強(qiáng),流動(dòng)對離子輸運(yùn)的貢獻(xiàn)逐漸顯現(xiàn),電流也隨之增大.流動(dòng)強(qiáng)度達(dá)到峰值之后稍有回落,進(jìn)入電流減小的過渡階段.

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

圖10 給出了系統(tǒng)到達(dá)穩(wěn)態(tài)之后的流動(dòng)形態(tài)和離子濃度分布,流場朝內(nèi)側(cè)聚集,離子濃度呈駝峰狀分布.在離子濃度耗盡程度較大的區(qū)域,電解質(zhì)溶液向上流動(dòng),攜帶離子遠(yuǎn)離表面,表現(xiàn)為對離子輸運(yùn)的加強(qiáng),并在兩側(cè)匯集到離子選擇性表面,抵消了一部分電遷移作用,表現(xiàn)為局部電流密度的減弱.由于流動(dòng)強(qiáng)化離子輸運(yùn)作用更強(qiáng),離子輸運(yùn)效率隨之增大,最終導(dǎo)致電流增大.

圖10 V=22φ0 穩(wěn)態(tài)流動(dòng)強(qiáng)度和和陰離子濃度分布.(a) 流場分布,(b) 陰離子濃度Fig.10 Contour of (a) fluid velocity and (b) anions concentration on steady state with V=22φ0

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

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

圖12 V=26φ0 不同時(shí)刻的流動(dòng)和離子濃度分布,(a),(b)對應(yīng)于圖11 中0.02 s,(c),(d)對應(yīng)于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 不同時(shí)刻的流動(dòng)和離子濃度分布,(a),(b)對應(yīng)于圖11 中0.02 s,(c),(d)對應(yīng)于0.03 s (續(xù))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 中,電流密度在流動(dòng)發(fā)生之后也出現(xiàn)了兩次轉(zhuǎn)變,分別對應(yīng)于小渦流動(dòng)增強(qiáng)離子輸運(yùn)、小渦轉(zhuǎn)變?yōu)榇鬁u流動(dòng).兩種流態(tài)的最大速度并沒有太大差距,但是大渦狀態(tài)下的離子輸運(yùn)效率更強(qiáng).該電壓下的最終的穩(wěn)態(tài)流動(dòng)與圖10 類似,但流動(dòng)強(qiáng)度明顯提高.

4 結(jié)論

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

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