多因素耦合作用下緩沖材料對鈾的長期阻滯效應數(shù)值模擬

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

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

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

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

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

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

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

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

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

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

1.1 模型假設

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(1)

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

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

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

(2)

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

(3)

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

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

(4)

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

(5)

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

引入彈性本構關系與Cauchy幾何方程(式(6))。

σ′=D∶ε

(6)

(7)

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

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

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

(8)

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

(9)

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

依據(jù)非飽和緩沖材料質(zhì)量守恒方程中的水體質(zhì)量守恒方程(式(10))，因不考慮汽化現(xiàn)象，故去掉源匯項，可得流體的質(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)下，對于多孔介質(zhì)密度變化較小的情形，多孔介質(zhì)密度可以表示為式(19)[48]。

(19)

并且可知式(20)。

(20)

式(19)對時間求導可寫為式(21)。

(21)

根據(jù)廣義Darcy定律，忽略慣性和黏性效應，以及溫度梯度對液體流場的熱驅(qū)動作用，則液相的相對表觀速度表示為式(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，土體的體積應變；σ′，有效應力；σ′0，有效孔隙水應力；μs，固相位移量；εij，應變張量場；ui,i，位移場；δij，Kronecker符號。

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

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

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

(24)

(25)

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

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

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

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

Jv=vfc

(26)

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

(27)

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

(28)

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

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

D=Dm+Dd

(29)

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

F=Kdc

(30)

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

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

(31)

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

(32)

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

(33)

(34)

(35)

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

2 飽和緩沖材料中核素遷移擴散模型求解及長期阻滯特性分析

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

2.1 幾何模型處理

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

圖2 緩沖材料多場耦合長期阻滯性能Mock-up實驗系統(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)置模塊導入

(1) 固體力學接口

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

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

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

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

(3) 達西定律

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

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

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

2.3 邊界條件與初始值

(1) 溫度邊界條件

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

(2) 位移邊界條件

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

(3) 水力邊界條件

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

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

(5) 初始值

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

表1 多場耦合模型初始值Table 1 Initial values of THMC coupling model

2.4 耦合模型參數(shù)

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

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

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

2.5 多因素耦合下核素遷移擴散特征及預測分析

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

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

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

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

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

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

3 結 論

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

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

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