非理想體系多離子非穩(wěn)態(tài)擴散模型及數(shù)值模擬

張志紅， 陳 楊， 鄭九州， 田改壘， 楊 凡

(北京工業(yè)大學 城市與工程安全減災教育部重點實驗室， 北京 100124)

1 研究背景

擴散是海水入侵、地熱開發(fā)、油氣開發(fā)、核廢料泄露、污染物運移等眾多工程領域中普遍存在的現(xiàn)象[1-3]，而控制擴散過程是環(huán)境修復和廢水處理基本技術中最有效且廣泛應用的手段之一。實際工程中不可避免地遇到多種離子共存的現(xiàn)象，流體或溶液呈現(xiàn)出非理想特征，為了合理準確地預測實際擴散過程，亟需構建非理想體系中多離子擴散數(shù)學模型，并對其影響機理開展深入探討。

目前多離子擴散模型主要有Fick、Maxwell-Stefan(MS)、Onsager-Fuos(OF)及Nernst-Planck(NP)擴散模型。Fick定律基于機械動力學理論，擴散通量與濃度梯度呈線性關系，形式簡單，但模型建立未體現(xiàn)溶液非理想性特征[4]。多數(shù)研究均是基于其線性的擴散通量方程，將溶液非理想性對擴散的影響體現(xiàn)在有效擴散系數(shù)上，對有效擴散系數(shù)進行了預測或者試驗數(shù)據(jù)的反演[5-9]。MS方程采用流體動力學理論，基于離子驅動力等于摩擦力，建立了各離子交叉耦合的多離子擴散模型，由化學勢表征的驅動力能夠考慮溶液非理想性效應，但擴散通量與離子濃度間的非線性關系使得模型求解困難，大多數(shù)MS模型在實際應用中將溶液環(huán)境視為理想狀態(tài)，并未獲得溶液非理想性對離子擴散影響的相關定量或定性分析[10-13]。NP模型考慮了離子運動過程中由于擴散速度不同導致的離子間靜電力作用[14-19]，忽略了非理想溶液中離子強度梯度對擴散的影響。1923年Debye和Huckel提出了離子間的相互作用會造成電解質(zhì)溶液偏離理想溶液的觀點[20-21]，王仁遠等[22]基于電解質(zhì)溶液微觀結構熵效應的概念建立電解質(zhì)溶液理論，闡明了離子間靜電作用是電解質(zhì)溶液偏離理想溶液的主要根源。因此為了獲得離子真實擴散運動特征，在擴散建模時應充分考慮溶液非理想性特征的影響。

Li Hang等[23]基于化學勢為驅動力的非線性形式，推導得到了廣義線性擴散通量方程，其形式簡單，方便運用和求解。本文在化學勢表征驅動力的基礎上，運用廣義線性擴散通量方程，引入Guntelberg修正的活度系數(shù)表示理想溶液與實際溶液中離子濃度的偏差，結合質(zhì)量平衡方程建立了適用于非理想體系中多離子非穩(wěn)態(tài)擴散模型。所建模型是經(jīng)典Fick擴散模型的拓展，實現(xiàn)了對多離子體系非理想性特征的描述。采用有限元軟件Comsol Multiphysics進行數(shù)值求解，分析了擴散機理的主要影響因素，對多離子并存條件下的擴散時空分布規(guī)律開展了定量分析和探討，從而揭示多離子體系在多孔介質(zhì)中的真實擴散規(guī)律。

2 非理想體系多離子非穩(wěn)態(tài)擴散模型

2.1 多離子擴散通量方程

對多孔介質(zhì)孔隙液相而言，當溶液呈理想狀態(tài)時，可采用經(jīng)典的Fick定律描述離子的擴散，其擴散通量表達式為：

Ji=-nτDici

(1)

式中：Ji為離子i的擴散通量，mol/(m2·s);n為多孔介質(zhì)孔隙率;τ為彎曲因子;Di為擴散系數(shù),m2/s;為梯度算子；ci為離子i的濃度，mol/L。

已有研究結果表明多離子共存時離子間靜電作用將導致溶液偏離理想狀態(tài)，需引入活度系數(shù)表示真實溶液和理想溶液之間的偏差。本文采用考慮溶液非理想性的廣義線性擴散通量方程描述實際溶液中離子擴散通量[23]：

Ji=-nτDi(rici)

(2)

式中：ri為離子i的活度系數(shù)，用以體現(xiàn)多孔介質(zhì)孔隙液相的活性特征。

為了獲得i組分的擴散通量Ji，需確定活度系數(shù)的具體表達式，本文選用Debye-Huckel活度系數(shù)模型[24]，其表達式為：

(3)

(4)

式中：常數(shù)A和B為取決于水的介電常數(shù)、密度和溫度；常數(shù)a與離子有效直徑(水化半徑)有關；zi為離子的電荷數(shù)；I為離子的離子強度。

2.2 多離子非穩(wěn)態(tài)擴散模型

考慮擴散過程中離子濃度和活度系數(shù)隨時空的變化，根據(jù)質(zhì)量守恒定律，建立多孔介質(zhì)孔隙液相中離子i的非穩(wěn)態(tài)擴散方程：

(5)

將公式(2)代入公式(5)進一步擴展為：

(6)

Guntelberg對Debye-Huckel活度系數(shù)模型參數(shù)進行了具體取值，公式(3)可寫為：

(7)

公式(6)及(7)共同構成了適用于非理想體系的多離子非穩(wěn)態(tài)擴散模型，模型能夠體現(xiàn)離子存在環(huán)境的非理想性特征且溶液活性隨時空發(fā)生改變，符合真實的多離子擴散現(xiàn)象規(guī)律。

3 模型驗證

為驗證所建模型的正確性和有效性，采用Xi Yonghui等[25]獲得的試驗數(shù)據(jù)進行對比分析。試驗采用自制的三區(qū)域容器裝置分別對ZnSO4和CuSO4溶液中鋅(Zn2+)和銅(Cu2+)離子的濃度進行了擴散研究。

圖1和2分別為在相同的邊界條件和參數(shù)取值條件下Zn2+和Cu2+濃度擴散的試驗數(shù)據(jù)、本文模型及Fick模型預測的對比結果。結果表明：經(jīng)典的Fick定律預測結果與試驗數(shù)據(jù)存在偏差，所建模型預測結果與實測試驗數(shù)據(jù)吻合度較高。

4 數(shù)值計算

4.1 定解條件

初始條件和邊界條件如下：

初始條件：ci(z,0)=0

上邊界條件：ci(0,t)=ci0(t≥0)

所選下邊界條件為Cauchy邊界條件，介于透水與不透水邊界之間，更適用于實際工況[26]。

4.2 模型參數(shù)

選用鈉(Na+)、鋅(Zn2+)和 鉻(Cr3+)3種離子作為數(shù)值模擬特征離子，輸入濃度設為恒定源，孔隙液中化學成分為可溶無機物且不與黏土顆粒發(fā)生化學反應。模型參數(shù)取值如表2所示。

4.3 數(shù)值模擬及分析

4.3.1 離子價位對擴散的影響 圖3、4、5分別為鈉(Na+)、鋅(Zn2+)和鉻(Cr3+)3種離子不同模擬年限的濃度隨擴散距離的變化曲線。

圖1Zn2+濃度隨擴散距離的變化曲線 圖2Cu2+濃度隨擴散距離的變化曲線

表2 模型參數(shù)

圖3Na+濃度隨擴散距離的變化曲線 圖4Zn2+濃度隨擴散距離的變化曲線

圖5 Cr3+濃度隨擴散距離的變化曲線

由圖3、4、5可知Na+、Zn2+、Cr3+離子濃度均隨擴散距離的增大呈現(xiàn)逐漸衰減的趨勢，在同一土層深度處離子的積累濃度隨擴散時間的延長而不斷增大且所建模型預測結果明顯滯后于Fick模型預測結果，究其原因主要是多離子共存的環(huán)境下,離子間的靜電作用導致溶液偏離理想狀態(tài)，使得各離子濃度得不到完全有效的發(fā)揮，對比結果表明由于多離子共存導致孔隙液相呈現(xiàn)的非理想性特征，對離子的擴散范圍及擴散穿透時間具有顯著的延滯效果。另外值得注意的是：擴散時間相同時，本文模型與Fick模型模擬結果的差異隨著離子價位升高而逐漸增大。

由公式(7)可知，當離子價位升高時，表征孔隙流體活性特征的活度系數(shù)逐漸減小，則由孔隙溶液非理想性造成的對離子擴散的延滯效果逐漸凸顯并不斷增強。

4.3.2 離子初始濃度對擴散的影響 針對不同初始濃度對離子擴散的影響展開了數(shù)值模擬。圖6所示為初始濃度0.001、0.01、0.1、1 mol/L的情況下，鋅(Zn2+)離子擴散的時空分布規(guī)律。

圖6 Zn2+濃度隨擴散距離的變化曲線

由圖6(a)～6(d)可以看出，各模擬工況下Zn2+濃度均隨擴散距離的增大逐漸衰減，在同一土層深度處，所建模型預測結果明顯小于Fick模型預測結果。尤為重要的是，隨著初始源輸入濃度的增大，所建模型與Fick模型所得結果的差異也隨之不斷增大。主要是由于離子初始濃度值越大，導致多孔介質(zhì)孔隙流體活性越弱，活度系數(shù)越小，宏觀表現(xiàn)為對擴散離子的阻滯效果越強，造成離子運動減慢，離子濃度隨時空的分布規(guī)律受到明顯的影響。

4.3.3 分析與討論 由離子價位對擴散影響的模擬結果(圖3～5)和離子初始濃度對擴散的影響結果(圖6)可知，活度系數(shù)作為多組分離子共存的重要表征參數(shù)，對擴散趨勢的發(fā)展具有極其重要的延緩作用。

活度系數(shù)可取值為常數(shù)和非常數(shù)兩種。若r=1，則表明離子所處的溶液環(huán)境為理想狀態(tài)；若r≠1的常數(shù)，則多孔介質(zhì)孔隙液相中離子的非穩(wěn)態(tài)擴散公式(5)為：

(8)

由于ri為常數(shù)，公式(8)可簡化為經(jīng)典Fick定律形式，即：

(9)

基于公式(9)得到的數(shù)值模擬結果如圖3～5和6中Fick模型預測結果所示，離子濃度隨時空的擴散分布規(guī)律不受活度系數(shù)的影響。

但實際工程中，污染源一旦滲入多孔介質(zhì)中，離子濃度將隨時空發(fā)生改變，則與離子特性(化合價及濃度)相關的活度系數(shù)也必將隨之發(fā)生變化，即真實的活度系數(shù)是隨污染物離子擴散發(fā)生變化的變量，由此導致多孔介質(zhì)孔隙液相中離子的非穩(wěn)態(tài)擴散方程形式無法退化為經(jīng)典的Fick定律描述方程，可由公式(6)方程形式描述。通過圖3～5和6采用本文所建模型得到的預測結果可知，擴散機制中將活度系數(shù)視為變量所得結果與Fick擴散模型所得結果間的差異，隨離子化合價的升高和初始濃度的增大而逐漸增強，考慮溶液的非理想性特征對離子擴散進程具有顯著的遲滯效果。由此表明考慮離子擴散的溶液環(huán)境非理想性特征對預測結果的影響不可忽視，建議各類工程屏障設計及服役性能評估時應充分重視多孔介質(zhì)孔隙液相活性的影響。

5 結 論

在經(jīng)典Fick定律的基礎上考慮溶液非理想性特征對離子擴散的影響，引入Guntelberg修正后的活度系數(shù)表示理想溶液與實際溶液中離子濃度的偏差，并結合廣義線性非穩(wěn)態(tài)擴散通量描述方程，建立了適用于非理想體系的多離子非穩(wěn)態(tài)擴散模型。主要結論如下：

(1)隨著擴散距離的增大，模擬的特征離子濃度逐漸衰減，且所建模型模擬結果呈現(xiàn)出滯后于Fick模型模擬結果的整體趨勢。

(2)溶液非理想性特征對離子擴散具有顯著的阻滯作用。當離子價位一定時，初始濃度越高對離子擴散阻滯效果越明顯；當離子初始濃度相同時，對離子擴散約束效果隨著離子價位增高而增大。

(3)活度系數(shù)隨擴散離子濃度和價位的改變而變化，與實際工程更相符，預測結果更有效和準確。