高放廢物處置庫緩沖層熱-水-力耦合數值模擬

王玉平，王 哲，易發成，曾志強

（1.宜賓學院國際應用技術學部，四川宜賓 644000；2.西南科技大學環境與資源學院，四川綿陽 621010；3.宜賓學院理學部，四川宜賓 644000）

核能的高效利用和核武器的發展，產生對生物圈有害的放射性廢物.特別是高水平放射性廢物（簡稱高放廢物,HLW），是對環境潛在危險最大的放射性廢物，其放射周期可達上萬年，是最難處置且最受關注的科學問題，到2050年，我國將產生8.3萬噸乏燃料[1].目前世界上主要有兩種處理方式：一是通過核嬗變法改變核素性質使之無毒化；二是深地質處置使之與生態圈完全隔離.深地質處置，是將HLW采用工程屏障系統和自然屏障系統深埋于地下500-1000米，使核廢物與生物圈永遠隔離.在HLW深地質處置的工程屏障體系中，緩沖材料一方面延滯廢物中核素向地下水遷移，另一方面阻止核素進入生物圈環境.對于HLW處置中熱-水-應力（THM）耦合過程的研究，許多學者從試驗和理論方面做了很多工作，包括室內試驗、原位試驗和數值模擬.因多場耦合試驗周期長、相互作用復雜，因此更多的是建立理論模型進行研究.

Biot[2]最先提出了流體飽和多孔介質的動力學問題，并推導了基本控制方程，在水文學、土力學等領域發揮了重要作用.Noorishad等[3]在Biot固結理論的基礎上推出飽和土體的THM耦合方程，Barton和Bandisoo[4]對考慮地下水的工程中應力場、滲流場以及溫度場之間的耦合關系做了更深一步的研究.1990年，Mu和Landanyi[5]首次提出THM三場耦合問題，并給出了簡化的模型，耦合問題開始成為研究的重點.國內外大多數多場耦合模型都是在連續介質理論基礎上建立的.Zhang等[6]對深地質處置庫中的Callovo和Opalinus黏土的THM行為進行了廣泛的研究，Yu等[7]研究了HLW釋放熱量導致孔隙壓力增加.Collin[8]提出一個THM模型解決粘土屏障中的耦合問題，并對壓實膨潤土進行濕熱試驗對比驗證.Kanno[9]對地質處置庫中緩沖材料的THM耦合過程進行了力學模型開發和數值模擬.Della[10]對壓實粘土的持水和應力應變，提出了一種考慮土體多相多尺度耦合的本構模型.Fernandez[11]通過試驗研究不同溫度對高壓實膨潤土滲透系數的影響.

本文以混合物理論為基礎，基于COMSOL有限元分析平臺建立軸對稱模型，對Mock-up試驗臺架中緩沖層的THM耦合進行有限元模擬，考察緩沖層中溫度場、飽和度和應力場的分布及變化規律，以期為高放廢物處置庫設計提供依據.

1 非飽和緩沖材料THM耦合模型

連續介質理論認為，在物質所占據的空間中，物質是無間隙地連續分布的.可忽略物質的微觀結構，用偏微分方程表達宏觀物理量.混合物理論是在現代連續體物理的精神下發展起來的，提出了一個新的構成公理“成分的等同性”，并用它來獲得與微觀對應相一致的宏觀構成方程.混合理論用于將連續介質力學原理推廣到幾個可相互滲透的連續介質來模擬多相系統.為構建彈性力學基礎上的多孔介質THM耦合模型，作如下假設：

（1）非飽和的介質被視為多相系統(固體、液體和氣體).孔隙一部分被液態水填充，一部分被氣體填充.

（2）多相介質是一種混合物，每一相都是連續的，每個空間點同時被一個物質點占據.

（3）氣相密度是溫度和壓力的函數，且滿足理想氣體狀態方程.

（4）流體滲流符合Darcy定律，且流體不與固體顆粒發生化學反應.飽和蒸氣壓力滿足Kelvin濕度定律，水蒸氣的擴散符合廣義Fick定律.

（5）多孔介質的變形滿足小變形假設，Terzaghi有效應力原理對非飽和介質適用.

（6）假定固、液、氣三相保持局部熱平衡.

多場耦合的場與場之間是密切相關的，物理場變量與其它物理場變量之間是相互影響的.高放廢物處置庫中的緩沖材料作為一種多孔介質，滿足多場耦合作用的理論體系，包括連續介質理論、混合物理論、滲流力學、傳熱學理論、普遍的守恒原理及Fick定律等.為了建立非飽和緩沖材料THM多場耦合數學模型，需結合混合物理論、連續介質理論、三大守恒定理（動量守恒、質量守恒、能量守恒）以及Fick定律推導非飽和緩沖材料的THM耦合控制方程.

1.1 動量守恒方程

式中：D是彈性張量，u為位移向量，T為溫度，χ是Bishop系數，Pl是液體壓力，Pg是氣體壓力，K是土體的排水體積模量，βsT是固相體積熱膨脹系數，εsw是土體濕化膨脹體積應變，δ是Kronecke單位張量，F為體積力.

1.2 質量守恒方程

（1）固相骨架質量守恒方程

式中：n為孔隙率，εv為土體的體應變，孔隙水壓為pl，孔隙氣壓為pg.

（2）水體質量守恒方程

式中：ρl為水體的密度，k是本征滲透率張量，krl為水體的相對滲透率，μl是水的動力黏滯系數，pl為液體壓力，S為飽和度，clp為壓縮性系數，csp為飽和度的吸力影響系數，βsT為飽和度的溫度影響系數，βlT是熱膨脹系數，ω為液相遷移系數，psv為飽和蒸汽壓力，pv為蒸汽壓力.

（3）氣體質量守恒方程

式中：vgr為氣相相對于固相的本征流速，klr為液相相對滲透率，kgr為氣相相對滲透率，kgT為混合氣體熱驅動系數，μg為氣體的動力粘滯系數，ρg為混合氣體的本征密度，cgp為壓縮系數，βgT為混合氣體熱膨脹系數，cap是干空氣壓縮系數，βaT是干空氣熱膨脹系數，H為Henry溶解系數，ρl、ρa分別為水、干空氣的密度.

1.3 能量守恒方程

給定任意單元體，在單位時間內α相總內能的變化率是一定等于該相穿過單元體的熱流量加上熱源，三相混合的能量守恒方程可以寫為：

式中：vgr為氣相相對于固相的本征流速，kgr為氣相的相對滲透卒，kgT為混合氣體熱驅動系數，μg為氣體的動力粘滯系數，Ks、Kl和Kg分別為固、液和氣相體積應變模量，βs、βl和βg分別為固、液和氣相體積熱膨脹系數，cs、cl和cg分別為固、液和氣相比熱容，λ為熱傳導系數.

2 高放廢物處置概念庫及耦合模型參數

2.1 高放廢物處置概念庫

本文的多場耦合緩沖材料長期阻滯性能Mock-up試驗腔體直徑750 mm，高750 mm，內部加熱器高180 mm，外徑200 mm；罐體為密封狀態，內部的空氣、水蒸氣不與外界交換，故氣體邊界均設為零通量邊界.基于COMSOL多場耦合問題的計算模塊，試驗臺架的數值計算模型如圖1所示.取沿加熱器邊緣中點水平方向為截面A，圖1（c）為截面A上11個計算結果輸出點的位置，其中A點（x=0.09 m）位于加熱器外緣中部，B點（x=0.23 m）位于緩沖層中間，C點（x=0.37 m）位于靠近圍巖的緩沖層.COMSOL中理查茲方程接口可用于描述變飽和多孔介質流體運動，理查茲方程模型中飽和液體體積分數設定為1，殘余液體體積分數設為0，滲透率模型選擇“滲透率”，儲水模型選擇“用戶定義”，保留模型選擇Van Genuchten(VG模型)，VG模型的參數分別為a=0.009，b=1.3386，l=2.46×106J/kg.

圖1 試驗臺架的數值計算模型

2.2 緩沖材料的相關參數和本構方程

緩沖材料力學模型為線彈性模型，其性質參數如表1所示.

表1 膨潤土性質參數[12-14]

本構關系如下：

（1）膨潤土的水體熱驅動[15]系數：klT=2.5×10(n×S-15.5).

（2）通過試驗得到膨潤土的導熱系數：λ=2.1545ω+1.0316ρd-1.9875n+0.6181S-0.5508.

2.3 流體的相關參數和本構方程

采用的VG模型可使飽和土的吸力為0，符合吸濕過程中土的吸力變化特點[16]，相關參數見表2.

表2 流體的相關參數

流體的相關本構關系如下：

(1)水體相對滲透率[17]

其中n為孔隙率，n0為初始孔隙率.

(2)水 體 動 力 粘 滯 系 數[18]μl=1.984×10-6×

(3)水體的熱膨脹系數[12]βl=1.7769×10-10T3-5.7556×10-8T2+1.1383×10-5T+3.024×10-6.

(4)水體的比熱容Cl=6872.8945-23.09T+0.069T2-7.35×10-5T3.

(5)液相遷移系數[19]ω=其中C為調節系數，取2×10-5；R=8.3144 J/(mol·K).

2.4 氣體的相關參數和本構方程

氣體相關參數如表3所示.

表3 氣體的相關參數[20,21]

氣體的相關本構關系如下：

(1)氣體密度.氣體包含水蒸氣和干空氣，N2和O2組成干空氣（體積比4∶1）.則水蒸氣密度、空氣密度和干空氣密度分別為：

式中：Ml為水分子摩爾質量（0.018 kg/mol），pv為蒸汽壓，MO為O2摩爾質量（0.032 kg/mol），MN為N2摩爾質量（0.028 kg/mol）.

(2)氣體熱驅動系數[12]：kgT=6×10-12(1-S).

(5)干空氣溶解系數.干空氣主要由氧和氮組成，Henry溶解系數H[12]可以表示為：

3 THM耦合數值模擬結果及分析

3.1 溫度場變化過程

緩沖材料的溫度場變化趨勢如圖2所示，可見緩沖材料的溫度隨時間不斷升高.

圖2 緩沖材料溫度場分布變化

第10天，緩沖材料溫度變化不大；第100天，將加熱器設定為30℃并保持恒定，并在緩沖材料外邊界處施加0.1 MPa水頭壓力，隨著熱量向往擴散，加熱器周圍溫度場已有明顯小幅度上升，而緩沖層外邊界溫度場未發生明顯變化.這是因為膨潤土的滲透速率較慢，滲透系數較低.當試驗進行到150天時，加熱器升溫到40℃，而緩沖材料外邊界處水頭壓力約為0.5 MPa，緩沖材料由內向外溫度梯度逐漸增大，加熱器周圍緩沖材料溫度上升較快，緩沖材料外邊界約為30℃，在351天時，將溫度升到90℃并保持恒定，此時緩沖材料外邊界處水頭壓力為1 MPa，相對于360天而言，1000天時加熱器周圍的緩沖材料溫度場分布穩定，緩沖材料外邊界溫度值較低，且變化較小，因為水頭壓力加速了水的滲流，帶走一部分熱量.500天和1000天相比變化不大，說明在500天時，緩沖材料就可以把廢物罐產生的熱量傳到HLW處置庫外，防止處置庫內部熱量聚集，有利于HLW處置庫的安全和穩定.

圖3為溫度隨徑向距離的變化趨勢.距加熱器越近，溫度越高，溫度從靠近加熱器到外邊界逐漸降低，呈線性分布.第10天緩沖層內溫度保持在30℃，加熱后第100天溫度迅速上升.第1000天在靠近圍巖的緩沖層處溫度低于30℃，說明熱傳遞和地下水滲透已經貫穿整個緩沖層.在500天和1 000天的溫度分布變化很小，說明膨潤土的導熱性能較好，溫度趨于穩定.圖4顯示了截面A不同位置的溫度隨時間的演變.可以看出：各點的溫度都隨時間的增加而增大，且距離加熱器越近，曲線斜率越陡，溫度增長速率越大.靠近加熱器處的緩沖層，溫度很快達到80℃，隨后溫度增速減緩趨于穩定.緩沖層中間的溫度值很快達到50℃穩定下來，反映出緩沖層具有較好的熱穩定性.

圖3 溫度隨徑向距離的變化

圖4 溫度隨時間的變化

3.2 飽和度變化過程

圖5 所示的緩沖材料的飽和度變化趨勢表明，在水平方向上，飽和度值從外向內逐漸增大，在豎直方向變化規律一致，沒有出現因為加水方式和水的重力作用，導致緩沖層上、下部含水量差異較大的情況.在豎直方向上出現飽和度值中間大，兩頭小的輕微差異，這是因為模型層上下水力邊界條件為不透水邊界，上下部分的水體只能向中間滲透.

圖5 緩沖層飽和度隨時間變化

圖6 顯示了飽和度隨徑向距離的變化趨勢，在試驗前100天，靠近圍巖的飽和度值迅速上升，靠近加熱器的飽和度值有所下降，中間緩沖層變化很小；當試驗進行到500天，靠近加熱器一側緩沖層飽和度值達到了85%；靠近圍巖的緩沖層已經飽和；當試驗進行到1000天，靠近加熱器一側飽和度值達到89%.隨著試驗的繼續，整個緩沖層飽和度值分布逐漸趨于均勻.圖7顯示了截面A上飽和度值隨時間的變化曲線，每個點的飽和度值都隨時間的增加而增大，且距離注水邊界越近，曲線斜率越陡，飽和度值增長速率越大；在加熱器附近，其飽和度值先減小再達到穩定，隨后緩慢上升.整個試驗過程飽和度值以不同的增長速率進行，這可能是由于水沿著傳感器電纜滲透和膨潤土混合物的不均勻性造成的.靠近加熱器到緩沖層一半的飽和度值，在試驗前400天呈現波浪形態，這是因為模型考慮了水蒸氣的蒸發和孔隙氣壓導致飽和度值暫時減小.A、B兩點的飽和度值在開始階段變化不大，分別為57.61%，55.35%，隨后慢慢增大，大約在500天分別為84.78%，88.17%，這是因為在試驗初期，加熱器溫度為30℃，靠近熱源一側緩沖層的溫度迅速上升，濕度逐漸減小，熱效應大于濕效應.隨著高壓水頭作用，飽和度值升高，上升速度為前期慢后期快，而C點的飽和度很快上升到100%，因為該點受到地下水的充分滲透.

圖6 飽和度隨徑向距離的變化

圖7 飽和度隨時間的變化

3.3 膨脹力變化過程

緩沖材料的膨脹力變化趨勢如圖8所示.試驗初期，進水壓力小，緩沖層滲透性較低，膨脹力增長緩慢.當應力傳感器與緩沖材料砌塊完全接觸后，隨著進水壓力的增大，水力邊界處的膨脹力快速增大，隨后其他緩沖層的膨脹力也不同程度的增大.第0天，緩沖材料處于非飽和狀態，整個緩沖層的膨脹力為負值；第360天，注水壓力恒定為1 MPa，水力邊界處緩沖層迅速膨脹，膨脹力增大；第500天，注水壓力升到1.5 MPa并保持恒定值，加熱器溫度為90℃，緩沖層進一步向里滲透，膨脹力繼續不斷增大.第1000天，隨著試驗的繼續，膨脹力由外向內逐漸變大并達到峰值.

圖8 緩沖層膨脹力隨時間變化

圖9 是膨脹力隨徑向距離的變化趨勢，說明膨脹力從靠近圍巖邊界開始增大，在靠近加熱器附近，由于溫度升高導致膨脹力有所下降.

圖9 膨脹力隨徑向距離的變化

圖10 為截面A上各點膨脹力隨時間的變化曲線.總的來看，每個點的膨脹力都隨時間的增加而增大，膨脹力時程曲線經歷了快速膨脹、緩慢膨脹并逐漸趨于穩定.距離注水邊界越近，曲線斜率越陡，膨脹力增長速率越大；距離加熱器越近，曲線斜率越緩，膨脹力增長速率越小.含水量從加熱器向緩沖層逐漸下降，水分重新分布導致局部膨脹力增加.注水壓力一直保持恒定的增長，由于不斷向緩沖層注水，使得緩沖層的飽和度值不斷升高，并在其周圍形成孔隙壓力，緩沖層的孔隙壓力和飽和度不斷重新分布調整，因此膨脹力曲線出現小波動變化.

圖10 膨脹力隨時間的變化

4 結語

本文通過COMSOL有限元數值模擬，探討了高放廢物處置庫緩沖層非飽和介質中的熱-水-力多場耦合現象，（1）基于混合物理論，將緩沖材料視為多組分連續介質，針對耦合概念模型提出了基本假設，為建立數學模型提供理論依據和前提.考慮蒸發潛熱的影響，基于質量守恒方程、能量守恒方程和動量守恒方程，建立了非飽和緩沖材料的THM耦合數學模型，得到關于水壓、溫度等數學物理方程.（2）發現緩沖層中溫度隨時間先增大后趨于穩定，且距離加熱器越近，曲線斜率越陡，溫度值增長速率越大，說明緩沖層具有良好的熱傳導性.溫度對緩沖層飽和度的影響較大，特別是在靠近加熱器周圍，飽和度在水平方向上，從外向內逐漸增大，在豎直方向變化規律一致，沒有出現因為加水方式和水的重力作用，使得下部水量增多的情況.緩沖層中膨脹力受到水壓和溫度的共同影響.試驗初期，進水壓力小，膨脹力不斷調整，膨脹力增長緩慢.隨著進水壓力的增大，緩沖層的膨脹力不同程度的增大.靠近水力邊界的緩沖層膨脹力受高壓水頭作用，膨脹力較大，本文建立的模型能夠體現溫度-滲流-膨脹力耦合現象，可為高放廢物地質處置庫設計提供依據.