王 哲，王玉平，李江波，張振龍，馮喜楊，羅 瑩，張家千，曾秋平，劉 艷，易發(fā)成，張 灝，趙 雨，王振雨，吳亞東

1.西南科技大學(xué) 核廢物與環(huán)境安全國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，四川 綿陽 621010；2.中國(guó)科學(xué)技術(shù)大學(xué) 地球與空間科學(xué)學(xué)院，安徽 合肥 230026；3.宜賓學(xué)院 國(guó)際應(yīng)用技術(shù)學(xué)部，四川 宜賓 644000；4.中國(guó)工程物理研究院 材料研究所，四川 江油 621700；5.四川輕化工大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，四川 宜賓 644000

隨著核技術(shù)的廣泛應(yīng)用，不可避免地產(chǎn)生大量的放射性廢物，特別是那些具有放射性強(qiáng)、毒性大、半衰期長(zhǎng)、釋熱率高和處理難度大的高放廢物，需要對(duì)這些放射性廢物進(jìn)行有效地安全處置，才能夠避免其進(jìn)入生物圈而對(duì)環(huán)境和人類健康構(gòu)成嚴(yán)重威脅。目前，對(duì)于高放廢物的安全處置，國(guó)際上普遍接受的、唯一可行的處置方式是深地質(zhì)處置[1-5](圖1[5])。 它是一項(xiàng)以核素的固化、阻滯為核心內(nèi)容，以多重屏障為主要手段，以1萬年～10萬年時(shí)間尺度內(nèi)人類健康和環(huán)境不受其破壞為目標(biāo)的處置庫(kù)系統(tǒng)工程[6]。高放廢物深地質(zhì)處置庫(kù)的長(zhǎng)期穩(wěn)定性取決于多重屏障體系中各子系統(tǒng)的長(zhǎng)期安全性，即天然屏障的完整性及工程屏障的長(zhǎng)期有效性，而前者與場(chǎng)址的地層巖性、地質(zhì)構(gòu)造與巖體的力學(xué)穩(wěn)定性及地下水賦存情況有關(guān)，可通過選取有利工程地質(zhì)與水文地質(zhì)條件及完整性好的圍巖來實(shí)現(xiàn)，后者可通過優(yōu)化處置單元設(shè)計(jì)和選取性能優(yōu)良的工程屏障材料(選取有利的固化體、廢物罐與緩沖材料)來實(shí)現(xiàn)[7-8]。

在高放廢物深地質(zhì)處置的多重屏障系統(tǒng)中，緩沖材料可作為最后一道人工屏障，對(duì)放射性核素的吸附滯留及阻滯特性將直接影響到放射性廢物隨地下水逐漸滲出地質(zhì)屏障而返回到生物圈的能力。因此，從這個(gè)意義上說，研究放射性核素在緩沖材料中的吸附、阻滯特性對(duì)高放廢物深地質(zhì)處置庫(kù)的長(zhǎng)期安全性具有重要的意義。放射性核素遷移并返回生物圈的主要?jiǎng)恿κ亲鳛楹怂貍鬟f載體的地下水的流動(dòng)和核素濃度差引起的擴(kuò)散，而阻滯核素遷移的因素主要是天然及人工工程屏障對(duì)放射性核素的吸附作用。因此，核素阻滯效應(yīng)研究的內(nèi)容就是放射性核素在天然地質(zhì)體和各種人工屏障材料中的吸附、滯留、擴(kuò)散以及在地下水中的遷移行為特征[9]。

(a)——豎向處置單元， (b)——水平處置單元圖1 高放廢物深地質(zhì)處置庫(kù)示意圖[5]Fig.1 Schematic diagram of high-level radioactive waste deep geological disposal repository[5]

目前國(guó)內(nèi)外關(guān)于放射性核素遷移的研究主要是通過室內(nèi)模擬實(shí)驗(yàn)來測(cè)定各種放射性核素在多種屏障材料中的吸附分配系數(shù)Kd、阻滯因子Rd、擴(kuò)散系數(shù)D和容量因子α，來表征放射性核素在屏障材料中的吸附、滯留、擴(kuò)散性質(zhì)以及在地下水中的遷移行為、遷移規(guī)律和遷移模型[10-12]。

20世紀(jì)70 年代，美國(guó)、瑞典、英國(guó)等國(guó)在這方面作了大量的工作，其中具有代表的是英國(guó)的Bradbury 等[13]的研究，他們研究了86Sr、99Tc、137Cs等核素在花崗巖、砂巖等介質(zhì)中的吸附及一維擴(kuò)散行為，證明了核素在巖石中的擴(kuò)散也符合Fick 擴(kuò)散定理，并由實(shí)驗(yàn)曲線獲得了核素?cái)U(kuò)散系數(shù)、巖石容量因子和吸附系數(shù)等參數(shù)。進(jìn)入20世紀(jì) 80 年代后，我國(guó)也開展了核素遷移方面的相關(guān)研究工作，對(duì)Sr、Cs、I、Co、Ru、Ce等元素在花崗巖、灰?guī)r、泥土與地下水等介質(zhì)中的遷移進(jìn)行了較多的研究，測(cè)得一系列滲透、擴(kuò)散和滯留參數(shù)，并利用這些參數(shù)作出了放射性核素的濃度隨遷移距離和遷移時(shí)間的變化關(guān)系圖[14]。

當(dāng)多重屏障體系的高放廢物深地質(zhì)處置庫(kù)在關(guān)閉并投入運(yùn)營(yíng)后，在一定的時(shí)期內(nèi)，地下水通過浸潤(rùn)將使緩沖層逐漸飽和并與包裝容器直接接觸。隨著時(shí)間推移，廢物包裝容器在輻照、地下水腐蝕及應(yīng)力等作用下，最終包裝容器銹蝕和破裂，核素將隨著地下水的滲流在整個(gè)多重屏障系統(tǒng)內(nèi)進(jìn)行遷移，其遷移過程主要有三個(gè)階段：第一階段為核素釋放-擴(kuò)散遷移，即核素在工程屏障中的遷移；第二階段為核素在天然地質(zhì)屏障中的對(duì)流-彌散遷移；第三階段為核素在生物圈中的分區(qū)傳遞[15-16]。而緩沖材料對(duì)核素的阻滯效應(yīng)，將直接影響到核素在第二、第三階段中的遷移進(jìn)程及處置庫(kù)系統(tǒng)的長(zhǎng)期穩(wěn)定性與安全性。

由于放射性衰變，核廢物隨著時(shí)間不斷放熱(幾千年～一萬年)，并在人工屏障中傳輸，改變了緩沖材料的溫度，進(jìn)而引發(fā)緩沖材料體積膨脹和熱應(yīng)力，改變流體的性狀和運(yùn)動(dòng)，破壞水-緩沖材料間的化學(xué)平衡并產(chǎn)生新的地球化學(xué)作用；而地球化學(xué)場(chǎng)具有特別重要的作用，表現(xiàn)在因溫度、濕度變化導(dǎo)致處置庫(kù)近場(chǎng)的化學(xué)作用引發(fā)的化學(xué)成分變化對(duì)固化體的溶解與核素的析出具有重要作用，進(jìn)而對(duì)核素在緩沖材料中的對(duì)流、擴(kuò)散及吸附等物理化學(xué)過程產(chǎn)生影響[17]?？梢姡谔幹脦?kù)的近場(chǎng)環(huán)境條件下，核素在緩沖材料中遷移擴(kuò)散受控于溫度場(chǎng)(T)、滲流場(chǎng)(H)、膨脹應(yīng)力場(chǎng)(M)和化學(xué)吸附場(chǎng)(C)的耦合作用，并且這種復(fù)雜的THMC耦合作用會(huì)直接影響到以核素為核心的溶質(zhì)遷移。為了揭示處置庫(kù)環(huán)境中多場(chǎng)耦合作用下緩沖材料行為變化的規(guī)律，瑞典、加拿大、比利時(shí)、日本、捷克、西班牙、法國(guó)、韓國(guó)等國(guó)研究人員開始對(duì)本國(guó)高放廢物處置庫(kù)預(yù)選膨潤(rùn)土進(jìn)行了多場(chǎng)耦合條件下的長(zhǎng)期性能實(shí)驗(yàn)[18-27]。此外，為了對(duì)高放廢物處置庫(kù)長(zhǎng)期穩(wěn)定性進(jìn)行評(píng)價(jià)，國(guó)際上已經(jīng)開展了數(shù)項(xiàng)重大合作項(xiàng)目，先后開展了DECOVALEX項(xiàng)目、FEBEX項(xiàng)目和China Mock-up項(xiàng)目等，為推動(dòng)高放廢物地質(zhì)處置安全評(píng)價(jià)作出了重要貢獻(xiàn)[28-33]。要充分理解這些耦合過程，數(shù)值建模是一個(gè)不可或缺的工具，因?yàn)樘幹脦?kù)的安全性必須保證10萬年或更長(zhǎng)時(shí)間。然而，在緩沖材料中模擬THMC耦合過程是具有挑戰(zhàn)性的，一方面因?yàn)轳詈线^程的本構(gòu)關(guān)系復(fù)雜且難以獲得，另一方面本構(gòu)關(guān)系的選擇和參數(shù)對(duì)模擬結(jié)果有很大影響。THMC多場(chǎng)耦合研究的數(shù)值軟件主要有TOUGHREACT、FLAC3D和COMSOL。Nardi等[34]采用Java開發(fā)了一款能夠有效協(xié)同多物理場(chǎng)耦合分析軟件COMSOL Multiphysics與地球化學(xué)模擬軟件PHREEQC的接口(iCP)，可以有效地處理復(fù)雜的地球化學(xué)反應(yīng)過程與其他物理過程的耦合，為有效模擬大量與地球化學(xué)反應(yīng)過程有關(guān)的多物理耦合問題提供了一個(gè)數(shù)值模擬平臺(tái)。Nasir等[35]基于COMSOL Multiphysics有限元和PHREEQC，并利用MATLAB編寫一個(gè)連接兩個(gè)軟件的代碼，開發(fā)了一種新的模擬工具，實(shí)現(xiàn)了兩種軟件的耦合，完成了模擬核廢料深部地質(zhì)庫(kù)(DGRs)近場(chǎng)發(fā)生的熱-水-力-地球化學(xué)(THM-GeoC)耦合過程及其對(duì)巖石孔隙度和滲透率演化的影響。張志紅等[36]基于流體的守恒方程，建立了某一個(gè)溶質(zhì)在飽和的土體中遷移的THMC全耦合模型，并使用COMSOL軟件研究了滲濾液在不同環(huán)境溫度下對(duì)溶質(zhì)遷移的影響。郭志光等[37]建立了非均衡吸附問題理論模型并利用COMSOL數(shù)值分析方法討論污染物遷移規(guī)律。張玉軍[38-40]將滲透遷移行為引入THM耦合控制方程中，建立了THMC耦合模型，將其運(yùn)用到假想的核素泄漏后的地下處置庫(kù)中，并對(duì)產(chǎn)生泄露的處置庫(kù)進(jìn)行多場(chǎng)耦合數(shù)值計(jì)算，研究其溫度、核素濃度、正應(yīng)力隨時(shí)間的變化和分布情況。林文勝等[41]分析了世界各國(guó)對(duì)高放廢物處置設(shè)施THMC耦合效應(yīng)的研究進(jìn)展。劉澤佳等[42]在溫度影響的前提下，將化學(xué)函數(shù)添加到孔隙水中，用來模擬有機(jī)污染物對(duì)多孔介質(zhì)力學(xué)性質(zhì)產(chǎn)生的影響，建立了一個(gè)關(guān)于多孔介質(zhì)的THMC數(shù)學(xué)模型。Yasuhara等[43]建立了一種THMC耦合數(shù)值模型，通過模擬高放廢物深埋地下情景，應(yīng)用模擬圍壓和溫度條件，預(yù)測(cè)巖石滲透性的長(zhǎng)期演化，模型預(yù)測(cè)清楚地顯示出壓力溶解對(duì)滲透率隨時(shí)間變化的顯著影響。Changsoo等[44]為了提高韓國(guó)高放廢物的處置密度，提出了一種多層儲(chǔ)存庫(kù)的概念來替代單層儲(chǔ)存庫(kù)的概念，采用數(shù)值模擬的方法比較了現(xiàn)有韓國(guó)高放廢物處置庫(kù)的THMC耦合行為。Cinzia等[45]利用地下水模擬程序PMWIN，以非均勻地下水流場(chǎng)為例，模擬了放射性核素在地質(zhì)處置庫(kù)的遠(yuǎn)場(chǎng)遷移過程并討論了分布系數(shù)及水力梯度兩個(gè)重要參數(shù)對(duì)地下水流場(chǎng)的影響。上述關(guān)于高放廢物地質(zhì)處置庫(kù)的長(zhǎng)期穩(wěn)定性評(píng)價(jià)工作主要側(cè)重于天然地質(zhì)體和開挖擾動(dòng)區(qū)內(nèi)裂隙巖體的多場(chǎng)耦合研究，較多且較全面，其研究主要集中在熱-水-力等多場(chǎng)耦合作用下對(duì)巖體的物理力學(xué)性能響應(yīng)特征及數(shù)值模擬。而關(guān)于緩沖材料在多場(chǎng)耦合作用下的性能研究，主要關(guān)注于壓實(shí)膨潤(rùn)土基緩沖材料在多場(chǎng)耦合條件下的物理力學(xué)特征的變化規(guī)律，并未對(duì)多場(chǎng)耦合作用下核素在緩沖材料中的遷移擴(kuò)散行為進(jìn)行深入研究。因此，本工作擬通過建立飽和緩沖材料溫度場(chǎng)、滲流場(chǎng)、膨脹應(yīng)力場(chǎng)和化學(xué)吸附場(chǎng)的耦合數(shù)學(xué)模型，基于模型實(shí)驗(yàn)獲得的參數(shù)，采用多物理場(chǎng)仿真軟件COMSOL中的PDE接口和內(nèi)置接口，將推導(dǎo)出的控制方程和本構(gòu)關(guān)系進(jìn)行導(dǎo)入，建立多因素耦合作用下緩沖材料對(duì)核素的長(zhǎng)期阻滯效應(yīng)的數(shù)值模擬方法，為工程實(shí)踐和放射性廢物的安全處置提供有價(jià)值的依據(jù)，同時(shí)也能夠?yàn)槲覈?guó)高效廢物深地質(zhì)處置庫(kù)地下實(shí)驗(yàn)室開展1∶1工程尺度的工程屏障設(shè)計(jì)與安全性能評(píng)價(jià)提供具有參考價(jià)值的科學(xué)數(shù)據(jù)。

1 飽和緩沖材料中核素遷移擴(kuò)散的熱-水-力-化耦合模型建立

1.1 模型假設(shè)

為構(gòu)建建立在彈性力學(xué)基礎(chǔ)上的多孔介質(zhì)熱-水-力-化耦合數(shù)學(xué)模型，需作出如下假設(shè)[9]：

(1) 非飽和的介質(zhì)被視為多相系統(tǒng)(固體、液體和氣體)；固體骨架的孔隙部分被液態(tài)水填充，部分被氣體填充；氣相假定為由干燥空氣和水蒸氣組成的理想氣體混合物；液相假定為水和溶解的空氣組成；

(2) 多相介質(zhì)被認(rèn)為是一種混合物；每一相都是連續(xù)的，混合物中的每一個(gè)空間點(diǎn)都被認(rèn)為同時(shí)被每一相的一個(gè)物質(zhì)點(diǎn)占據(jù)；混合物作為一個(gè)整體的平衡定律與單相物質(zhì)平衡定律的方程形式相同；

(3) 氣相密度是溫度和壓力的函數(shù)，且滿足理想氣體狀態(tài)方程；

(4) 流體滲流運(yùn)動(dòng)滿足廣義的Darcy 定律；氣相氣體中包括水蒸氣和干空氣，兩者壓力分布符合 Dalton 分壓定律，飽和蒸氣壓力滿足 Kelvin 適度定律；水蒸氣的擴(kuò)散符合廣義 Fick 定律；

(5) 多孔介質(zhì)的變形滿足小變形假設(shè)，且假定無論飽和、還是非飽和狀態(tài)下，修正的Terzaghi有效應(yīng)力原理均能適用；

(6) 熱傳導(dǎo)過程滿足廣義Fourier定律，并且假定固相、液相、氣相，三相保持局部熱平衡狀態(tài)，即在空間內(nèi)的任意一點(diǎn)，固相、液相、氣相，三相均具有相同的溫度，并且兩相之間的熱交換只存在于相變過程中；

(7) 緩沖材料對(duì)放射性核素的吸附過程是線性的、可逆的等溫吸附；

(8) 放射性核素遷移到緩沖材料之前就已達(dá)到水飽和狀態(tài)，且固體骨架不發(fā)生變化；

(9) 放射性核素在緩沖材料遷移的過程中不考慮由于化學(xué)反應(yīng)而引起核素質(zhì)量的變化，遵循質(zhì)量守恒；

(10) 除了從高放廢物固化體中浸出的放射性核素外，不存在其他源項(xiàng)，核素在固相和液相中濃度達(dá)到平衡；

(11) 不考慮放射性核素衰變對(duì)源項(xiàng)濃度的影響。

1.2 飽和緩沖材料中核素遷移擴(kuò)散的熱-水-力-化耦合控制方程

耦合過程的分析不是單獨(dú)的，物理場(chǎng)變量與其它物理場(chǎng)變量之間相互影響[46]。多孔介質(zhì)的多場(chǎng)耦合研究了溫度場(chǎng)、應(yīng)力場(chǎng)、滲流場(chǎng)和化學(xué)場(chǎng)耦合作用下單相流或兩相流體在孔隙中的傳輸，以及固相骨架和流體溫度的分布和固相骨架形變的規(guī)律。因此，高放廢物處置庫(kù)中的緩沖材料作為一種多孔介質(zhì)，滿足多場(chǎng)耦合作用的理論體系，其中就包括：連續(xù)介質(zhì)理論、混合物理論、滲流力學(xué)、傳熱學(xué)理論、普遍的守恒原理及污染物擴(kuò)散Fick定律等。

為了建立核素在緩沖材料這種多孔介質(zhì)中遷移擴(kuò)散的多場(chǎng)耦合數(shù)學(xué)模型，需要結(jié)合混合物理論和連續(xù)介質(zhì)理論以及三大守恒定理(質(zhì)量守恒、動(dòng)量守恒、能量守恒)與溶質(zhì)擴(kuò)散的Fick定律來推導(dǎo)飽和緩沖材料中核素遷移擴(kuò)散的熱-水-力-化耦合控制方程。

(1) 飽和緩沖材料動(dòng)量守恒方程

當(dāng)沒有慣性力的情況下，非飽和緩沖材料動(dòng)量守恒方程為式(1)。

(1)

式中：α代表不同的相，即s表示固相、l表示液相和g表示氣相；ρα，α相的本征密度；nα，α相占據(jù)的體積比例；vα，α相(組分)的絕對(duì)速率；ρα，α相的表觀密度；bα，單位質(zhì)量α相(組分)的體積力大小；σs，固相應(yīng)力張量；σl，液相應(yīng)力張量；σg，氣相應(yīng)力張量。

采用小變形假設(shè)和小速度假設(shè)，即vα≈0，整個(gè)過程看作準(zhǔn)靜態(tài)過程，式(1)可簡(jiǎn)化為式(2)。

契約是指“依照法律訂立的正式的證明、出賣、抵押、租賃等關(guān)系的文書”。1932年美國(guó)律師學(xué)會(huì)在《合同法重述》中所下的定義是：契約是“一個(gè)諾言或一系列諾言，法律對(duì)違反這種諾言給予救濟(jì)，或者在某種情況下，認(rèn)為履行這種諾言乃是一種義務(wù)”。市場(chǎng)經(jīng)濟(jì)其實(shí)就是契約經(jīng)濟(jì)。國(guó)有企業(yè)以及國(guó)有控股的混合所有制企業(yè)，有時(shí)會(huì)依賴自己占有的公共資源和壟斷地位對(duì)外漠視規(guī)則和契約精神，這種作風(fēng)延伸到自身的企業(yè)治理上由于不重視契約的管理帶來的管理漏洞也是比比皆是。很好地把契約化管理的思想應(yīng)用到企業(yè)管理過程中，一方面契約簽訂的過程可以促使科學(xué)決策，另一方面所有的經(jīng)濟(jì)業(yè)務(wù)行為的決策都有軌跡可跟蹤，利于監(jiān)督追責(zé)。

(2)

而對(duì)于飽和多孔介質(zhì)，不考慮氣相，式(2)可以寫為式(3)。

(3)

根據(jù)有效應(yīng)力原理，飽和多孔介質(zhì)的有效應(yīng)力(σ′)可以表示為式(4)。

σ′=σs+σl=σ-plδ-KβsTTδ-Kζcδ

(4)

聯(lián)立式(3)與式(4)并寫為時(shí)間導(dǎo)數(shù)的形式，得到飽和緩沖材料動(dòng)量守恒方程(式(5))。

(5)

式中：c，溶質(zhì)濃度；ζ，化學(xué)膨脹系數(shù)；T，溫度；σs，固相應(yīng)力；σl，孔隙水應(yīng)力；pl，孔隙水壓力；βsT，固相體積熱膨脹系數(shù)；K，土體的排水體積模量；δ，Kronecker Delta單位張量；σ，總應(yīng)力張量。

引入彈性本構(gòu)關(guān)系與Cauchy幾何方程(式(6))。

σ′=D∶ε

(6)

(7)

式中：D，土體的彈性張量；ε，應(yīng)變張量；u，位移向量。

(2) 飽和緩沖材料質(zhì)量守恒方程

在飽和狀態(tài)下，根據(jù)理論和實(shí)驗(yàn)研究，除水力梯度外，考慮溫度梯度和濃度梯度也會(huì)改變流體的運(yùn)動(dòng)，這種現(xiàn)象稱為Soret效應(yīng)和Dufour效應(yīng)。因此，考慮Soret效應(yīng)和Dufour效應(yīng)的廣義達(dá)西定律為式(8)。

(8)

同時(shí)，在達(dá)西定律中，速度其實(shí)是指流體相對(duì)于固體的速度，即為液體相對(duì)于固相骨架的平均流速(vr)，計(jì)算如式(9)。

(9)

式中：n，飽和多孔介質(zhì)的孔隙率；vf，飽和狀態(tài)下流體的實(shí)際流速。

依據(jù)非飽和緩沖材料質(zhì)量守恒方程中的水體質(zhì)量守恒方程(式(10))，因不考慮汽化現(xiàn)象，故去掉源匯項(xiàng)，可得流體的質(zhì)量守恒方程式(11)。

(10)

(11)

(12)

(13)

將式(9)帶入式(11)得式(14)。

(14)

將式(14)展開，可得式(15)。

(15)

(16)

(17)

將式(13)與式(16)相加，并聯(lián)立式(8)、(9)、(17)，可得整體的連續(xù)性方程式(18)。

(18)

而在飽和狀態(tài)下，對(duì)于多孔介質(zhì)密度變化較小的情形，多孔介質(zhì)密度可以表示為式(19)[48]。

(19)

并且可知式(20)。

(20)

式(19)對(duì)時(shí)間求導(dǎo)可寫為式(21)。

(21)

根據(jù)廣義Darcy定律，忽略慣性和黏性效應(yīng)，以及溫度梯度對(duì)液體流場(chǎng)的熱驅(qū)動(dòng)作用，則液相的相對(duì)表觀速度表示為式(22)。

(22)

由于α=1-K/Ks為Biot系數(shù), 將式(16)、(18)和式(22)帶入式(14)可得飽和緩沖材料質(zhì)量守恒方程式(23)。

(23)

式中：Ks，固體基質(zhì)的體積模量；βlT，液體的熱膨脹系數(shù)；ρs0，多孔介質(zhì)基質(zhì)初始密度；T0，多孔介質(zhì)內(nèi)初始溫度；pl，孔隙水壓力；clp，水體的壓縮性系數(shù)；εv，土體的體積應(yīng)變；σ′，有效應(yīng)力；σ′0，有效孔隙水應(yīng)力；μs，固相位移量；εij，應(yīng)變張量場(chǎng)；ui,i，位移場(chǎng)；δij，Kronecker符號(hào)。

(3) 飽和緩沖材料能量守恒方程

將非飽和緩沖材料能量守恒方程式(式(24))中的氣相耦合項(xiàng)去除，并將達(dá)西定律替換為式(8)，可得飽和狀態(tài)下的緩沖材料能量守恒方程式(25)。

{T[ρs(1-n)cs+ρlnScl+ρgn(1-S)cg]+

(24)

(25)

式中：cl，液相的比熱容；cs，固相的比熱容；cg，氣相的比熱容；λ，導(dǎo)熱系數(shù)；Kl，液相成分體積模量；Kg，氣相成分體積模量；kgr，氣相的相對(duì)滲透率；kgT，混合氣體熱驅(qū)動(dòng)系數(shù)；psv，飽和蒸汽壓力；pv，蒸汽壓力；ω，為液相遷移系數(shù)；μg為氣體的動(dòng)力黏滯系數(shù)；βl，液體的線熱膨脹系數(shù)；βs，固體的線熱膨脹系數(shù)；βg，氣體的線熱膨脹系數(shù)；Tg，氣體溫度。

(4) 飽和緩沖材料溶質(zhì)遷移方程

放射性核素在緩沖材料中的遷移擴(kuò)散受核素-緩沖材料相互作用與地下水水動(dòng)力學(xué)過程兩方面因素共同控制，因此，放射性核素遷移模式是由核素吸附模式和地下水水動(dòng)力學(xué)模式耦合而成的，在建立飽和緩沖材料中核素遷移方程時(shí)，需要考慮地下水的溶質(zhì)對(duì)流作用、水動(dòng)力彌散作用與核素的吸附作用[9,11]。

其中，溶質(zhì)對(duì)流是指當(dāng)?shù)叵滤诰彌_材料中滲流時(shí)，同時(shí)攜帶著核素離子以水的流動(dòng)速度在孔隙中遷移。這種隨水流運(yùn)動(dòng)而攜帶的溶質(zhì)的數(shù)量，稱為溶質(zhì)對(duì)流通量Jv。當(dāng)滲流為達(dá)西流時(shí)，隨之?dāng)y帶的溶質(zhì)對(duì)流通量與緩沖材料中核素的濃度c成正比[49](式(26))。

Jv=vfc

(26)

當(dāng)含放射性核素的地下水在緩沖材料中遷移時(shí)，由于孔隙存在使各點(diǎn)的流速向量和橫斷面平均流速向量的方向、大小各不相同，所以通過不同孔隙的核素離子在一定時(shí)間間隔所達(dá)到的位置也不同，這是由于含有核素的地下水具有黏滯性、孔隙大小不同及分布不均和固體顆粒的阻擋所以引起的。因此，在建立緩沖材料對(duì)核素的遷移阻滯模型時(shí)需要考慮水動(dòng)力彌散作用的影響。一般來說，水動(dòng)力彌散作用主要包括分子擴(kuò)散和機(jī)械彌散作用，它是一種不穩(wěn)定和不可逆的過程[50]。其中：分子擴(kuò)散是由于液體中所含核素濃度分布不均時(shí)，核素會(huì)在濃度梯度作用下從高濃度處向低濃度處遷移，它不受液體流動(dòng)本身的影響[51]。核素?cái)U(kuò)散通量可用Fick定律描述，即核素的擴(kuò)散通量與其濃度梯度成正比，擴(kuò)散方向與濃度梯度方向相反[52]，即Fick第一定律，如式(27)所示。

(27)

機(jī)械彌散也稱對(duì)流擴(kuò)散，指含核素地下水流體通過緩沖材料中孔隙的平均流速(即滲流速度)和濃度與緩沖材料斷面上的平均流速和濃度不一致所導(dǎo)致的分散現(xiàn)象[50]。在對(duì)流作用中考慮了通過緩沖材料斷面上平均流速產(chǎn)生的核素遷移問題。機(jī)械彌散實(shí)質(zhì)上是溶液通過緩沖材料中孔隙的平均流速和濃度的差值導(dǎo)致的核素遷移問題[53]。因此，由機(jī)械彌散作用引起的核素遷移通量也可用Fick定律來描述，如式(28)所示。

(28)

式中：Jd，機(jī)械彌散通量；Dd，機(jī)械彌散系數(shù)。

由于分子擴(kuò)散和機(jī)械彌散從引起核素遷移效果上看類似，所以，在實(shí)際應(yīng)用中把分子擴(kuò)散和機(jī)械彌散統(tǒng)稱為水動(dòng)力彌散，水動(dòng)力彌散系數(shù)D為分子擴(kuò)散系數(shù)Dm和機(jī)械彌散系數(shù)Dd之和，如式(29)所示。

D=Dm+Dd

(29)

以膨潤(rùn)土為基材的集成緩沖材料能夠作為人工屏障材料應(yīng)用于高放廢物處置庫(kù)中，其重要原因之一就是其對(duì)放射性核素具有良好的吸附作用。當(dāng)含有放射性核素的地下水在緩沖材料中發(fā)生滲流時(shí)，流體與緩沖材料顆粒的接觸面促使水中的放射性核素吸附在緩沖材料介質(zhì)表面上，進(jìn)而延遲其遷移。在等溫線性吸附條件下，緩沖材料對(duì)核素的吸附量與溶液中的放射性核素濃度(c)成正比，如式(30)所示。

F=Kdc

(30)

式中：F，緩沖材料對(duì)核素的平衡吸附容量，mg/g；Kd，吸附分配系數(shù)，mL/mg。

根據(jù)普遍守恒原理，利用式(26)—(30)，可得式(31)。

(31)

將式(31)兩邊同除n可得式(32)。

(32)

對(duì)于正壓流，液體的密度可以表示成壓力和溫度的函數(shù)ρl(pl,T)，因此水體密度的微分同水體密度之間的比值存在如下關(guān)系(式(33))。

(33)

(34)

(35)

式中：Rd，阻滯因子；D，水動(dòng)力彌散系數(shù)，當(dāng)水流速度較小時(shí)，D可近似為表觀擴(kuò)散系數(shù)Da。

2 飽和緩沖材料中核素遷移擴(kuò)散模型求解及長(zhǎng)期阻滯特性分析

模擬多因素耦合作用下核素在飽和緩沖材料中的遷移擴(kuò)散特征，實(shí)質(zhì)上就是對(duì)上述建立的多場(chǎng)耦合遷移擴(kuò)散控制方程進(jìn)行求解。這類耦合問題涉及到多過程、多參數(shù)的多重耦合，求解過程復(fù)雜且難度大。而COMSOL是一款高級(jí)數(shù)值仿真軟件，它是以有限元為基礎(chǔ)，通過求解偏微分方程(單場(chǎng))或偏微分方程組(多場(chǎng))來解決多場(chǎng)耦合問題[54]。它不僅自帶許多經(jīng)典的物理模塊，還支持與MATLAB、AutoCAD等相互導(dǎo)入、自定義PDE模塊，能更精準(zhǔn)地解決實(shí)際工程中的物理問題。因此，考慮到多物理場(chǎng)有限元軟件COMSOL Multiphysics具備直接全耦合求解的優(yōu)點(diǎn)，本次求解控制方程組主要采用內(nèi)置接口和添加上述守恒方程中的耦合項(xiàng)作為源項(xiàng)相結(jié)合方式，并利用課題組前期相關(guān)模型實(shí)驗(yàn)獲得的參數(shù)[11]，可實(shí)現(xiàn)多物理場(chǎng)核素在緩沖材料中遷移擴(kuò)散行為的直接耦合分析。

2.1 幾何模型處理

本課題組自主研制的多場(chǎng)耦合條件下緩沖材料長(zhǎng)期阻滯性能Mock-up實(shí)驗(yàn)系統(tǒng)主要由雙層中空罐體裝置系統(tǒng)、數(shù)據(jù)采集與處理系統(tǒng)、供水/供液控制系統(tǒng)構(gòu)成(圖2)。本工作的幾何模型以該實(shí)驗(yàn)系統(tǒng)中的多場(chǎng)耦合實(shí)驗(yàn)罐體裝置為參考，其內(nèi)腔直徑為750 mm，高度為750 mm，模擬固化體的加熱器殼體直徑為180 mm，高度為180 mm?；贑OMSOL多場(chǎng)耦合問題的計(jì)算模塊，取罐體內(nèi)腔簡(jiǎn)化截面作為水力邊界，長(zhǎng)寬均為750 mm；取加熱器外部截面為發(fā)熱邊界，長(zhǎng)寬均為180 mm(圖3(a))。對(duì)簡(jiǎn)化截面進(jìn)行網(wǎng)格剖分，采用COMSOL Multiphysics中的四邊形網(wǎng)格剖分(圖3(b))，并在邊界處添加邊界層，增加模型的收斂性。

圖2 緩沖材料多場(chǎng)耦合長(zhǎng)期阻滯性能Mock-up實(shí)驗(yàn)系統(tǒng)Fig.2 Mock-up experimental system for long-term retarding performance of buffering materials under multi-field coupling conditions

(a)——幾何模型及邊界， (b)——模型網(wǎng)格劃分圖3 數(shù)值求解幾何模型Fig.3 Numerical solution of geometric model

2.2 內(nèi)置模塊導(dǎo)入

(1) 固體力學(xué)接口

在固體力學(xué)接口中，添加線彈性材料，指定楊氏模量和泊松比組合，緩沖材料的模型下選擇“外部應(yīng)力”，選擇“應(yīng)力張量(空間)”，然后再添加源項(xiàng)作為耦合項(xiàng)。邊界選擇輥支撐，限制法向位移，并將1.2節(jié)構(gòu)建的飽和多孔介質(zhì)系統(tǒng)動(dòng)量守恒方程式(7)代入該接口。

(2) 多孔介質(zhì)傳熱接口

COMSOL中在添加物理場(chǎng)界面下，傳熱分支下的多孔介質(zhì)物理場(chǎng)中添加多孔介質(zhì)傳熱接口，用來模擬多孔介質(zhì)中的導(dǎo)熱和對(duì)流傳熱。多孔介質(zhì)中定義的溫度方程對(duì)應(yīng)于具有熱力學(xué)性質(zhì)平均模型的對(duì)流-擴(kuò)散方程，并同時(shí)考慮固體基質(zhì)和流體性質(zhì)。當(dāng)進(jìn)入多孔介質(zhì)的溫度和流體處于平衡狀態(tài)時(shí)，這個(gè)方程有效。如果不是，則使用傳熱界面中的局部熱非平衡接口代替。由COMSOL軟件給出的多孔介質(zhì)傳熱方程[55]與1.2節(jié)構(gòu)建的飽和狀態(tài)下的緩沖材料能量守恒方程式(25)聯(lián)合求解。

在“多孔介質(zhì)傳遞屬性”中，流體內(nèi)流速場(chǎng)設(shè)為dl·u/Rd，其余參數(shù)設(shè)為用戶定義，“多孔基體”內(nèi)所有參數(shù)設(shè)為用戶定義，方程中并未考慮體應(yīng)變對(duì)熱傳導(dǎo)的影響，因此在熱源項(xiàng)添加體應(yīng)變作用項(xiàng)，完成對(duì)飽和狀態(tài)下的緩沖材料能量守恒方程的導(dǎo)入。

(3) 達(dá)西定律

“達(dá)西定律”接口用來模擬多孔介質(zhì)孔隙中流體的流動(dòng)，以壓力梯度為主要驅(qū)動(dòng)力的前提下，可用于滲透率和孔隙率均很小且低速運(yùn)動(dòng)下的流體或介質(zhì)進(jìn)行建模[56]。由COMSOL軟件給出的達(dá)西方程[55]與1.2節(jié)構(gòu)建的飽和狀態(tài)下的緩沖材料質(zhì)量守恒方程式(23)聯(lián)合求解。

(4) 多孔介質(zhì)中的稀物質(zhì)傳遞

該接口專門用于模擬多孔介質(zhì)中的運(yùn)輸，包括固相和液相，其中化學(xué)物質(zhì)在多孔介質(zhì)內(nèi)有可能受到擴(kuò)散、對(duì)流、遷移、分散、吸附等的影響。該接口支持固相完全不變形或者氣相介質(zhì)不運(yùn)動(dòng)的情況。由COMSOL軟件給出的方程[55]與1.2節(jié)構(gòu)建的飽和緩沖材料溶質(zhì)遷移方程式(35)聯(lián)合求解。

2.3 邊界條件與初始值

(1) 溫度邊界條件

加熱器設(shè)定為溫度邊界，恒溫90 ℃，罐體外部設(shè)置為溫度邊界，恒溫25 ℃。

(2) 位移邊界條件

設(shè)置緩沖材料邊界為法向位移為0，切向可自由活動(dòng)。位移邊界均設(shè)置為第一類邊界條件(輥支撐邊界條件)。

(3) 水力邊界條件

水力邊界條件設(shè)置在加熱器邊界處，選擇壓力邊界，設(shè)為200 kPa，罐體外邊界為0 Pa。

(4) 溶質(zhì)濃度邊界

(5) 初始值

根據(jù)多場(chǎng)耦合條件下緩沖材料長(zhǎng)期阻滯性能Mock-up實(shí)驗(yàn)系統(tǒng)所需實(shí)驗(yàn)條件、緩沖材料的基本參數(shù)，以及實(shí)驗(yàn)測(cè)定的壓實(shí)膨潤(rùn)土的初始參數(shù)，并參考國(guó)內(nèi)外相關(guān)文獻(xiàn)[57-60]，給出耦合的初始值，結(jié)果列入表1。

表1 多場(chǎng)耦合模型初始值Table 1 Initial values of THMC coupling model

2.4 耦合模型參數(shù)

本模型采用的緩沖材料力學(xué)模型為線彈性模型，緩沖材料的性質(zhì)參數(shù)列入表2。

表2 緩沖材料性質(zhì)參數(shù)Table 2 Property parameters of buffer materials

此外，通過課題組前期相關(guān)動(dòng)態(tài)遷移柱實(shí)驗(yàn)獲得了鈾元素在膨潤(rùn)土基緩沖材料中遷移擴(kuò)散的相關(guān)參數(shù)：阻滯因子Rd=3.49、飽和滲透系數(shù)ks=1.00×10-12m/s和表觀擴(kuò)散系數(shù)Da=1.29×10-4m2/a[11]。

2.5 多因素耦合下核素遷移擴(kuò)散特征及預(yù)測(cè)分析

利用COMSOL軟件內(nèi)置瞬態(tài)求解器，在“直接”求解器中選擇“MUMPS”，求解因變量分別為濃度c、壓力p、溫度T以及(動(dòng)量方程對(duì)應(yīng)的)位移u，本次研究總計(jì)時(shí)長(zhǎng)為10萬年，設(shè)置初始步長(zhǎng)為0，時(shí)間步長(zhǎng)為10年，停止時(shí)間為10萬年，以此達(dá)到每10年輸出一次結(jié)果，直到10萬年時(shí)停止計(jì)算的目的，其他選項(xiàng)使用默認(rèn)設(shè)置。

為探討在地下水中通過浸潤(rùn)使緩沖層達(dá)到飽和后，緩沖材料中鈾在不同時(shí)間尺度內(nèi)的遷移范圍與遷移距離變化特征，依據(jù)上文中的多場(chǎng)耦合作用下核素遷移擴(kuò)散模型、初始條件、邊界條件和計(jì)算參數(shù)等進(jìn)行數(shù)值模擬。模擬分析結(jié)果表明(如圖4和圖5所示)，當(dāng)廢物包裝容器破裂并發(fā)生遷移的0～1 000 a時(shí)間尺度為鈾遷移擴(kuò)散的初期階段，由于初期階段緩沖材料對(duì)鈾的吸附容量大，加之溫度對(duì)緩沖材料吸附鈾的能力有正向影響，所以在多場(chǎng)耦合作用的初期階段時(shí)，核素在緩沖材料中遷移擴(kuò)散較緩慢，遷移距離隨時(shí)間增加的幅度均在1 m左右，表現(xiàn)為鈾濃度-遷移距離關(guān)系曲線比中后期階段(1萬～10萬年)的濃度-遷移距離關(guān)系曲線陡。

時(shí)間，a：(a)——100，(b)——500，(c)——1 000，(d)——10 000，(e)——50 000，(f)——100 000圖4 鈾在緩沖材料中遷移過程動(dòng)態(tài)變化Fig.4 Variation of radionuclide uranium migration with time in buffer materials

時(shí)間，a：1——100，2——500，3——1 000，4——10 000，5——50 000，6——100 000圖5 不同時(shí)間尺度下鈾在緩沖材料中的遷移距離Fig.5 Migration distances of U in buffer materials at different time

隨著核素遷移擴(kuò)散不斷進(jìn)行，在多場(chǎng)耦合作用的中后期階段(1萬～10萬年)，由于緩沖材料對(duì)鈾的吸附容量逐漸達(dá)到飽和并趨于穩(wěn)定后，其遷移距離較初期階段增加更明顯，表現(xiàn)為鈾濃度-遷移距離關(guān)系曲線變緩，遷移距離隨時(shí)間增幅均在3 m左右。

由此可知，若將緩沖材料的厚度分別設(shè)置成2、3 m和4 m，能夠分別確保在100、500 a和1 000 a內(nèi)將鈾阻滯在緩沖材料內(nèi)；若要保證10 000 a內(nèi)鈾不穿透緩沖材料而遷移到處置庫(kù)圍巖地下水中，則需要設(shè)置更大厚度的緩沖材料，其厚度至少應(yīng)為7 m。而Wang等[61]開展了單場(chǎng)環(huán)境下鈾元素在美國(guó)MX-80膨潤(rùn)土中遷移距離數(shù)值模擬，結(jié)果表明，鈾在100、1 000 a和10 000 a等不同時(shí)間尺度下分別遷移了0.66、2.1 m和6.6 m。此外，本課題組[11]也在源項(xiàng)釋放后100、1 000 a和10 000 a后，對(duì)鈾在B7AP 型集成緩沖材料(膨潤(rùn)土、凹凸棒石、黃鐵礦的質(zhì)量比為63∶27∶10)中的遷移距離進(jìn)行了數(shù)值模擬分析，結(jié)果表明，單場(chǎng)環(huán)境下鈾在上述不同時(shí)間尺度下分別遷移了0.25、0.8 m和3.0 m?？梢?，在處置庫(kù)的近場(chǎng)環(huán)境條件下，鈾在緩沖材料中的遷移擴(kuò)散過程受到了多物理場(chǎng)耦合作用，導(dǎo)致不同的緩沖材料對(duì)鈾的長(zhǎng)期阻滯效果受到了一定程度的削弱。因此，通過開展多場(chǎng)耦合作用下鈾在緩沖材料中的遷移特征模擬研究，能夠?yàn)楦叻艔U物處置庫(kù)中人工屏障設(shè)計(jì)及緩沖材料厚度設(shè)置等提供科學(xué)參考。

3 結(jié) 論

(1) 基于混合物理論、連續(xù)介質(zhì)理論、質(zhì)量守恒、動(dòng)量守恒、能量守恒及溶質(zhì)擴(kuò)散的Fick定律，推導(dǎo)飽和緩沖材料中核素遷移擴(kuò)散的熱-水-力-化耦合控制方程。

(2) 基于多物理場(chǎng)仿真軟件COMSOL Multiphysics具備直接全耦合求解的優(yōu)點(diǎn)，以自主研制的多場(chǎng)耦合條件下緩沖材料長(zhǎng)期阻滯性能Mock-up實(shí)驗(yàn)系統(tǒng)為幾何模型，采用內(nèi)置接口和添加熱-水-力-化耦合控制方程中的耦合項(xiàng)作為源項(xiàng)相結(jié)合方式，并在初始條件、邊界條件和前期相關(guān)模型實(shí)驗(yàn)獲得的參數(shù)的基礎(chǔ)上，利用COMSOL軟件內(nèi)置瞬態(tài)求解器，實(shí)現(xiàn)多物理場(chǎng)耦合作用下核素在緩沖材料中遷移擴(kuò)散行為的直接耦合分析。

(3) 多物理場(chǎng)仿真軟件COMSOL Multiphysics的多因素耦合作用下飽和緩沖材料中核素遷移擴(kuò)散行為及長(zhǎng)期阻滯特性分析表明，初期階段核素在緩沖材料中遷移擴(kuò)散較緩慢，遷移距離隨時(shí)間增幅均在1 m左右；中后期階段，由于緩沖材料對(duì)鈾的吸附容量逐漸達(dá)到飽和并趨于穩(wěn)定，其遷移距離較初期階段增加更明顯，遷移距離隨時(shí)間增幅均在3 m左右；若要保證1 000 a及10 000 a時(shí)間尺度內(nèi)鈾不穿透緩沖材料而遷移到處置庫(kù)圍巖地下水中，則需要分別設(shè)置緩沖厚度為4 m和7 m。通過開展多場(chǎng)耦合作用下，鈾在緩沖材料中的遷移特征數(shù)值模擬研究，能夠?yàn)槲覈?guó)高放廢物深地質(zhì)處置庫(kù)地下實(shí)驗(yàn)室開展1∶1工程尺度的工程屏障設(shè)計(jì)與安全性能評(píng)價(jià)提供技術(shù)方法和設(shè)計(jì)依據(jù)參考。