電化學(xué)梯度作用下的poisson-nerns-一株運(yùn)移方程

電化學(xué)梯度作用下的poisson-nerns-一株運(yùn)移方程

1多離子運(yùn)移的數(shù)學(xué)模型在電位梯度的驅(qū)動下，離子通過sent-ml方程在多孔介質(zhì)中轉(zhuǎn)移。式中,ji為第i離子流量(mol/(s·m2));Du為i離子的擴(kuò)散系數(shù)(m2/s);▽ci為i離子濃度(mol/m3);為i離子的濃度梯度(mol/m4);▽φ為電勢梯度(V/m);F為法拉第常數(shù)(96,500Cmol-1);R為氣體常數(shù)(8.31Jmol-1K-1),T為絕對溫度(K);zi為i離子的電荷數(shù)。公式右邊第一項i為離子擴(kuò)散對流量的貢獻(xiàn);第二項為電遷移對流量的貢獻(xiàn)。可以方程(1)是在Fick定律的基礎(chǔ)上考慮了電勢對離子運(yùn)移的影響。方程(1)考慮質(zhì)量平衡,可以得到每個離子濃度隨時間的變化:式中,為離子的通量變化;t為時間,s;da為時間比例系數(shù),da=1為非穩(wěn)態(tài)運(yùn)移,da=0為穩(wěn)態(tài)運(yùn)移。方程(1)和(2)的未知數(shù)比方程數(shù)多一個,必須有另一個方程偶合離子運(yùn)移過程,為獲得方程的解,考慮各個離子之間的平衡條件(1)電中性條件(2)電位為常數(shù)假定(3)零電流條件式中,N為系統(tǒng)中離子總個數(shù);w為有效顆粒表面電荷密度。針對不同的研究問題,方程(3)~(5)僅在特定的條件下成立。最為嚴(yán)格的條件是電勢變化和顆粒表面電荷密度之間的關(guān)系。一般稱為高斯靜電理論(也叫Poisson方程)式中,ε0為介質(zhì)的介電常數(shù)(C2/(J·m));εw為液體的相對介電常數(shù)(無量綱);ρs為顆粒表面電荷密度(C/m2)。方程(2)、(6)即為多離子運(yùn)移Poisson-Nernt-Planck方程。當(dāng)方程(2)、(6)中zi和w均為0時,方程(2)、(6)簡化為Fick定律。2方程的求解2.1非穩(wěn)態(tài)優(yōu)化方法在過去的幾十年,不少學(xué)者對Nernst-Planck方程的求解進(jìn)行了研究。Planck本人于1980年給出了第一個解答。他研究了兩個單價離子(+,-)穩(wěn)態(tài)無量綱的擴(kuò)散方程,假定分析域兩個邊界施加固定離子濃度和電位。通過假定電中性和零電流條件,而獲得解析解。盡管幾何及邊界條件非常簡化和理想化,只有這樣才能獲得的解析解。文獻(xiàn)[1,4],用不同的方法獲得了穩(wěn)態(tài)無量綱的擴(kuò)散方程的解析解。用Planck同樣的邊界條件和同樣的電中性、零電流假設(shè),推導(dǎo)了解析解,該解析解可以計算任意個數(shù)不同離子價的離子。上述方法都是基于穩(wěn)態(tài)的擴(kuò)散試驗。同樣基于電中性和零電流條件,文獻(xiàn)推導(dǎo)了非穩(wěn)態(tài)問題的解析解。獲得了不同膜電位解析解,但卻不能得到濃度和電位的分布曲線。此外他們的應(yīng)用限定在單價離子。文獻(xiàn)基于常量的電位場假設(shè)(例如,假定電位在多孔介質(zhì)膜中線性分布)的解析解。通過該簡化假設(shè)后,可以直接求解Nernst-Planck方程。由于直接獲得方程組的解析解非常困難,研究者開始嘗試采用數(shù)值解法。文獻(xiàn)提出一種算法,可以解決非穩(wěn)態(tài)問題。采用逐步逼近的方法,求解時需要首先使用普朗克的解析解,然后逐步地推獲得最后結(jié)果,從而導(dǎo)致該算法只能用于分析非常簡單的例子。最近文獻(xiàn)發(fā)展了有限差分法,可以求解任意個數(shù)、不同離子價的離子、非穩(wěn)態(tài)問題。并非對方程(2)進(jìn)行直接差分,而是通過對電中性和零電流條件進(jìn)行變換,引入兩個新的變量作為迭代系數(shù)。然后,轉(zhuǎn)換后的方程組變得非常復(fù)雜,應(yīng)用到二維和三維空間問題時極為困難。同時轉(zhuǎn)換還包含了復(fù)雜的邊界條件。文獻(xiàn)引入膜內(nèi)可變的電流密度定義。對質(zhì)量平衡方程(2)直接采用有限差分法求解,通過采用顯式歐拉格式對時間離散。例如對于N個離子的方程,先計算N-1個離子濃度和電位,通過電中性假定計算第N個離子的濃度。由于引入的電流必須等于膜內(nèi)電流,因此可以確定新的電位分布。這個方法的主要優(yōu)點(diǎn)是可以非常方便的拓展到二維和三維分析,但是它仍然是基于電中性假定。前述的數(shù)值方法基本都是基于NERNST-PLANCK方程,耦合電中性條件或者常電位假定,已經(jīng)成功的解釋了離子擴(kuò)散機(jī)理。盡管這些假定對于某些時候是個不錯的選擇,但是文獻(xiàn)卻指出,對于某些實(shí)例,需要使用更為嚴(yán)格Possion方程。文獻(xiàn)指出,電中性和常電位假定都只是Possion方程的一個特例。電中性假定僅在離子濃度比較高時適用,而常電位假設(shè)則是在離子濃度相當(dāng)?shù)蜁r有效。一些研究者嘗試耦合Nernst-Planck方程和Poisson方程。例如,文獻(xiàn)耦合了PNP方程和Stokes方程,描述液體中包含帶電粒子在圓柱形孔中的運(yùn)移問題,分析僅限于穩(wěn)態(tài)問題。采用有限元法和伽遼金加權(quán)余量法求解質(zhì)量平衡方程。采用Picard迭代法,直到獲得收斂解。文獻(xiàn)提出了無量綱的穩(wěn)態(tài)PNP方程組的數(shù)值解法。一旦穩(wěn)態(tài)達(dá)到,擴(kuò)散通量應(yīng)該是常數(shù),采用有限差分法求解Nernst-Planck方程。該方法與文獻(xiàn)方法相似,主要不同是該方法先知道電位分布,然后計算濃度分布,其計算簡圖如1所示。文獻(xiàn)采用有限元法和伽遼金加權(quán)余量法求解Nernst-Planck方程,迭代采用收斂速度快的牛頓-拉夫遜方法,對穩(wěn)態(tài)和非穩(wěn)態(tài)的Nernst-Planck程進(jìn)行了求解。但是沒有指出如何具體實(shí)現(xiàn)該算法。2.2材料生態(tài)特性的解析法式中,ea為質(zhì)量系數(shù);da為質(zhì)量系數(shù)(或者阻尼系數(shù));c為擴(kuò)散系數(shù);a為吸收系數(shù);f為源項;α為保守通量對流系數(shù);β為對流系數(shù);γ為保守通量源項;q為邊界上的吸收系數(shù);g為邊界上的源項;h,r為方程系數(shù)。其中α、β可以是矢量,c可以是矩陣,用于表示各項異性的材料特性。方程(7)是求解域上的控制方程,方程(8)是邊界條件,分別為第三類邊界條件和第一類邊界條件,即羅賓邊界和狄利克邊界條件。當(dāng),方程(7)即簡化為瞬態(tài)的Nernst-Planck方程(2),如果da=0,進(jìn)一步簡化為穩(wěn)態(tài)的Nernst-Planck方程。,方程(7)簡化為泊松方程(6)。3分析程序驗證在以Comsol有限元分析軟件模擬之前,為了解本方法分析之可行性,參考文獻(xiàn)[2,12],對程序進(jìn)行驗證。文獻(xiàn)采用無量綱化分析,這里也給出PNP方程無量綱化。3.1無限概率的方程方程式的無因次化參數(shù)定義如下式中,“0”,是指一個相對參考值;L為擴(kuò)散空間的長度。方程可以重新寫為其中3.2多離子問題的穩(wěn)定擴(kuò)散3.2.1見表1:comol-pla多孔介質(zhì)孔隙中鹽溶液由兩個單價離子組成(即z1=+1,z2=-1)。兩個離子擴(kuò)散系數(shù)均為1,。邊界條件:(1)Nernst-Planck方程(2),左端;右端。(2)Possion方程(6):左端;右端。該問題的求解可以選擇自定義多物理場方程,也可以選擇Comsol已有的模塊,這里采用Comsol中已有的電動流模塊和泊松方程模型聯(lián)合求解該問題。選擇(1)化學(xué)工程模塊>質(zhì)量傳輸>電動流>穩(wěn)態(tài)分析,該模塊用于分析Nernst-Plankck方程;因變量為c1,c2。(2)COMSOLMultiphysics>PDE模式>古典偏微分方程>泊松方程,該模塊用于分析泊松方程(15),因變量為V。Comsol中電動流模塊和泊松模塊PDES方程如下電動流方程泊松方程其中,電動流方程中電位V和為系數(shù),其值由泊松方程計算獲得。泊松方程中右端匯源項,其中離子濃度c1,c2由Nernst-Planck方程計算獲得。可以看出NernstPlanck方程和泊松方程的耦合影響。Comsol具體的建模和求解步驟參見文獻(xiàn)。為和文獻(xiàn)結(jié)果比較,需要將計算結(jié)果導(dǎo)出,這里介紹Comsol二次開發(fā)工具進(jìn)行處理。通過文件>輸出>FEM結(jié)構(gòu),就可以通過COMSOLScript二次開發(fā)工具,方便地對結(jié)果進(jìn)行處理,例如下面程序可以導(dǎo)出指定位置的濃度到數(shù)據(jù)文件。x=0:0.1:1;%需要取點(diǎn)的x坐標(biāo)cc1=postinterp(fem,'c1',x);%獲得x坐標(biāo)處陽離子濃度c1值save'd:\jg.txt'cc1–ascii;%保存數(shù)據(jù)到文件本文計算結(jié)果和Samson計算結(jié)果對比見圖2,由圖可見兩者結(jié)果相當(dāng)吻合。3.2.2顆粒表面電中性假設(shè)該例子考察四個離子(SO42-,Mg2+,Na+,K+)在多孔介質(zhì)中的一維穩(wěn)態(tài)擴(kuò)散問題。四個離子的無量綱擴(kuò)散系數(shù)分別為:,均為1,為1.67。邊界條件:(1)Nernst-Planck方程(2),左端,;右端,。(2)Possion方程(6):左端;右端。有效顆粒表面電荷密度c0=w=1。α2=1×10-10。文獻(xiàn)基于電中性假設(shè),對于該問題獲得了解析解,文獻(xiàn)采用數(shù)值解法求解了該問題。本文的計算結(jié)果見圖3,本文計算結(jié)果和上述求解結(jié)果非常接近,誤差不超過3%。3.3數(shù)據(jù)分析及結(jié)果該例子考察離子在多孔介質(zhì)中的非穩(wěn)態(tài)擴(kuò)散問題。無量綱化的時間。多孔介質(zhì)孔隙中鹽溶液由兩個單價離子組成(即z1=+1,z2=+2)。兩個離子的無量綱化擴(kuò)散系數(shù)分別為:。擴(kuò)散域中兩種離子的初始濃度均為0。時,在左端處注入兩個無量綱化離子濃度均為1.0,且擴(kuò)散過程中保持不變;右端,濃度為0。α2分別取0.01和0.0001,即對于給定的擴(kuò)散域,后者濃度是前者濃度的100倍,借此可以考察背景濃度對離子運(yùn)移的影響。另外分析還和Fick定律進(jìn)行對比,考察電場對離子運(yùn)移的影響。Comsol求解方法及步驟基本同例1相同,不同處為Nernst-Plankck方程取瞬態(tài)分析模塊,且陽離子和陰離子時間比例系數(shù)均取1。后處理同樣可以通過輸出>FEM結(jié)構(gòu),就可以通過COMSOLScript二次開發(fā)工具,方便地對結(jié)果進(jìn)行處理,例如下面程序可以導(dǎo)出某一時刻time指定位置的濃度到數(shù)據(jù)文件。時刻,陽離子濃度c1和陰離子濃度c2本文計算結(jié)果和文獻(xiàn)計算結(jié)果如圖4所示,由圖可見陰離子濃度兩者計算結(jié)果極為一致,陽離子濃度略有偏差,最大不超過2%。說明兩者計算結(jié)果還是比較一致的。本文Fick定律計算結(jié)果和文獻(xiàn)計算結(jié)果對比列于圖5,由圖可以看出兩種方法計算的結(jié)果也較為一致。初始和邊界條件相同,PNP方程和Fick定律的計算結(jié)果比較列于圖6。從圖6(a)可以看出,離子之間電場的耦合影響還是比較顯著的。一方面,運(yùn)移快的陽離子由于電場的作用,速度變慢。相反,運(yùn)移慢的陰離子擴(kuò)散速度被加快。圖6(b)中這種現(xiàn)象更為明顯,兩種離子的濃度已經(jīng)非常接近。分析結(jié)果和已有的結(jié)論(隨著的遞減,電中性條件將逐步達(dá)到)是吻合的。3.4角形單元劃分分析實(shí)例取Nacl鹽溶液在兩維飽和粘土RVE單元擴(kuò)散問題。粘土微單元長度和寬度均為200nm,孔系率n為0.87。邊界條件為左邊c1=10.6mol/m3,右邊c2=10.0mol/m3,左右兩邊電位φ1=φ2=0,上下兩邊取軸對稱邊界。采用三角形單元劃分網(wǎng)格,共9876個單元,5182個節(jié)點(diǎn),60741個自由度。ε0=8.85×10-12,εw=78,離子擴(kuò)散系數(shù)DNa=1.334×10-9m2/s,DCl=2.032×10-9m2/s,顆粒表面電荷密度ρs=0.01。4離子和固相表面的電場耦合過程(1)采用ComsolMultiphysics求解Poisson-NernstPlanck方程,可以耦合分析多離子、多價位的穩(wěn)態(tài)和非穩(wěn)態(tài)離子在多孔介質(zhì)中的擴(kuò)散問題,模型可以方便的拓展到二維和三維空間。(2)在電化學(xué)位驅(qū)動下,多離子在多孔介質(zhì)中的運(yùn)移模型,可以考慮離子與離子之間,離子與固相顆粒表面電荷之間的電場的耦合影響。(3)分析表明,離子之間、離子與多孔介質(zhì)固相表面的電場之間相互耦合影響顯著,傳統(tǒng)的Fick定律僅考慮濃度梯度作用,不考慮離子之間