喬 蘭，曹勝飛，，陳 亮，劉月妙，謝敬禮

（1.北京科技大學土木與環(huán)境工程學院，北京 100083；2.核工業(yè)北京地質(zhì)研究院，北京 100029）

高放射性廢物是一種放射性強、毒性大、核素半衰期長并且發(fā)熱的特殊廢物，需要把它們與人類生存環(huán)境長期、可靠地永久隔離，使其在進入生物圈之前衰變到人類可接受的水平。對于高放廢物的最終處置，目前公認的安全可靠，且技術(shù)上可行的方法為深地質(zhì)處置，即把高放廢物埋在距離地表500～1000m的地質(zhì)體中，使之永久與人類的生存環(huán)境隔離，埋藏高放廢物的地下工程即稱為“高放廢物處置庫”［1－2］。高放廢物處置庫采用的是“多重屏障系統(tǒng)”設計思路，即把廢物貯存在廢物罐中，外面包裹緩沖材料，再向外為圍巖。一般把廢物體、廢物罐和緩沖回填材料稱為工程屏障，把周圍的地質(zhì)體稱為天然屏障。由于地下結(jié)構(gòu)開挖擾動、高放廢物衰變放熱、地下水滲流等作用的影響，緩沖材料在地質(zhì)處置庫中處于典型的熱－水－力多場耦合環(huán)境。

為了揭示處置庫環(huán)境熱－水－力－化學多場耦合作用下緩沖回填材料行為變化的規(guī)律，國內(nèi)外學者進行了大量的試驗和理論研究。ZHANG Chunliang［3］利用瑞士的 Mont Terri地下實驗室對Opalinus硬泥巖現(xiàn)場加熱試驗，對其受熱所產(chǎn)生的熱－水－力全程進行了模擬，深入闡述了硬泥巖在THM耦合作用下的變化特征。LI Xiangling［4］在HADES的眾多現(xiàn)場試驗中，對Boom Clay和工程屏障體系進行熱－水－力性狀的研究。西班牙ENRESA負責在Grimsel地下試驗室實施的FEBEX試驗［5－8］，研究了高放廢物處置庫近場所發(fā)生的THM過程，獲得了大量的試驗數(shù)據(jù)。劉月妙等［9］在考慮地應力、地溫梯度、水力、高放廢物釋熱和緩沖回填材料膨脹力條件下，分別模擬高放廢物處置庫熱－力耦合、水－力耦合和熱－水－力耦合條件下的溫度場、應力場合應變場的變化規(guī)律，預測高放廢物處置庫的長期穩(wěn)定性。葉為民等［10］研究了我國高廟子膨潤土作為緩沖回填材料在多場耦合條件下的行為特征，建立了適合描述緩沖回填材料在熱－吸力等多場耦合條件下的本構(gòu)關(guān)系。

基于熱－水－力多場耦合理論，針對緩沖材料試驗臺架，本文采用LAGAMINE數(shù)值軟件對高壓實高廟子膨潤土在多場耦合條件下的長期性能進行了數(shù)值模擬分析，研究了溫度、飽和度、吸力和膨脹力隨時間的變化規(guī)律，為高放廢物地質(zhì)處置庫的建設提供參數(shù)和借鑒。

1 熱－水－力多場耦合數(shù)學模型

為了開展膨潤土在多場耦合條件下的長期性能研究，本研究建立了系統(tǒng)的孔隙介質(zhì)多場耦合數(shù)值模型。該數(shù)學模型考慮了THM條件下的耦合現(xiàn)象，包括熱的傳播（熱傳導和熱對流）、水流場、蒸氣擴散和冷凝、膨潤土在非飽和條件下的力學特性等。

1.1 水、熱擴散模型

膨潤土由固體顆粒、液態(tài)水和氣體三相組成，在高放廢物深地質(zhì)處置條件下，膨潤土的整體平衡是一個復雜的多場耦合現(xiàn)象。在本數(shù)學模型中，假設不同物質(zhì)質(zhì)量守恒，因此在總體平衡方程中不存在物質(zhì)之間的相互轉(zhuǎn)化，描述系統(tǒng)流動體系的變量包括液態(tài)水壓、氣壓和溫度。

1.1.1 水的擴散

在緩沖材料中，水分為液態(tài)水和水蒸氣，其中達西定律被用于描述液態(tài)水的流動。

式中：pw是液態(tài)水壓，g是重力加速度，μw是液態(tài)水的動力粘滯系數(shù)，kint是固有滲透系數(shù)，kr，w是相對滲透系數(shù)。

水蒸氣的流動傳播服從菲克定律，見式（2）。

式中：Datm是摩爾擴散系數(shù)，τv是孔隙介質(zhì)的折曲率，Sr，g是氣體的飽和度。

1.1.2 熱的傳播

在本研究模型中，采用唯一的溫度變量，假設非飽和土中所有的物質(zhì)均處于等溫狀態(tài)。熱的傳播主要與三種物理現(xiàn)象有關(guān)：熱傳導、熱對流、蒸氣化現(xiàn)象，如式（3）所示。

式中：Γ是介質(zhì)的熱傳導系數(shù)；cp，α是α介質(zhì)的比熱；L是潛熱。在本研究中，為了簡化，不考慮溫度對熱傳導系數(shù)的影響。

1.2 力學模型

在緩沖材料中，飽和度和吸力變化明顯，大量的研究表明，吸力對土體的力學特性有顯著地影響。本文采用巴塞羅那基本模型（BBM 模型）［11］，該模型可以較好的反映吸力對力學特性的影響。

1.2.1 屈服面

BBM模型是在劍橋模型的基礎(chǔ)上發(fā)展起來的，為了考慮飽和度對土體力學特性的影響，吸力被作為一個獨立的變量引入到模型中。因此，模型的塑性屈服面可以在（p，q，s）三維空間中進行描述，其主要包括以下內(nèi)容。

在（p，q）平面內(nèi)，對于給定的吸力值，塑性屈服面可以寫成如下的形式：

式中：p是平均凈應力，q是偏應力，M是臨界狀態(tài)線斜率，ps是土體的張拉強度，p0是土體在吸力為s時的預固結(jié)壓力。

式中：s是吸力，s0是SI屈服面的硬化參數(shù)，代表歷史上土體承受的最大吸力值。

1.2.2 硬化準則

該模型中，塑性屈服面的硬化演變?nèi)Q于總體的塑性體積應變增量，由硬化參數(shù)和s0控制，如式（6）、式（7）所示。

式中：e是土體的孔隙比，Pat是大氣壓，λs和κs分別是彈性和塑性剛度參數(shù)。

2 緩沖材料熱－水－力耦合試驗臺架模擬和數(shù)值分析

2.1 緩沖材料試驗臺架系統(tǒng)

試驗臺架系統(tǒng)由進水系統(tǒng)、試驗腔體、加熱系統(tǒng)和數(shù)據(jù)采集系統(tǒng)組成。試驗材料為內(nèi)蒙古高廟子膨潤土，傳感器包括溫度、溫濕度、土壓力、孔隙水壓力傳感器。簡化如圖1所示。試驗中，加熱器由室溫逐漸加熱至90℃并恒定運行，進水壓力維持在0.05MPa。

圖1 試驗結(jié)構(gòu)剖面圖

2.2 邊界條件

考慮到試驗臺架是標準的軸對稱結(jié)構(gòu)，作為近似簡化，取縱向截面的一半進行數(shù)值模擬，計算域尺寸如圖2所示，該模型有1447個單元，4161個節(jié)點。

系統(tǒng)初始溫度為20℃室溫，膨潤土初始飽和度為48%，孔隙比為0.57。根據(jù)高廟子膨潤土的土水特征曲線，初始吸力確定為80MPa。為了得到更好的數(shù)值收斂性，氣壓在計算中被設定為常數(shù)。

關(guān)于邊界條件，計算域的周邊位移固定，試驗臺架與周圍環(huán)境的熱對流現(xiàn)象，由邊界熱交換單元來模擬。模擬過程分兩部分進行。

2.1 兩組血漿NT-proBNP水平比較 病例組患兒的NT-proBNP水平在治療3 d、治療7d 、治療14 d時均顯著低于組內(nèi)治療前(F=176.405，P<0.05)，病例組患兒的NT-proBNP水平在治療前、治療3 d、治療7 d、治療14 d時均顯著高于對照組(F=286.557，P<0.05)。見表1。

1）第一步，在10小時內(nèi)，計算域與加熱器接觸邊界溫度從室溫（20℃）上升至預定加熱溫度90℃，同時外側(cè)水壓從初始吸力值上升到0.05MPa。

2）隨后3年內(nèi)，邊界條件保持穩(wěn)定，模擬膨潤土內(nèi)部的物理力學變化過程。

圖2 幾何和網(wǎng)格劃分圖

2.3 參數(shù)確定

水、熱特征參數(shù)基于CHEN Bao的試驗成果［12］，通過對試驗數(shù)據(jù)擬合，高廟子膨潤土吸力與飽和度的關(guān)系可以由式（8）進行描述。

式中：Sr，u是土體的最大飽和度，Sr，res代表在吸力控制下土體可以達到的最小飽和度。根據(jù)試驗結(jié)果，Sr，u和Sr，res數(shù)值確定為1.0和0.1，通過對試驗結(jié)果的擬合，式中參數(shù)a1＝3.5×10－6Pa－1，a2＝0.8，a3＝90。

根據(jù)YE Weimin的試驗結(jié)果［13］，高廟子膨潤土的固有滲透系數(shù)kint＝2.0×10－21m2，相對滲透系數(shù)曲線如式（9）所示。

在熱傳導試驗中，對不同飽和度的膨潤土試樣進行熱傳導系數(shù)的測量，根據(jù)試驗結(jié)果，可采用如下線性關(guān)系來反映熱傳導系數(shù)隨飽和度的變化。

土體其他水、熱參數(shù)如表1所示。

Cui開展了不同吸力和溫度條件下高廟子膨潤土的相關(guān)力學試驗［15］。根據(jù)λ（s）隨吸力的變化規(guī)律，參數(shù)γ和β得到了確定。表2匯總了高廟子膨潤土的主要力學參數(shù)。

2.4 數(shù)值模擬結(jié)果及分析

依據(jù)構(gòu)建的熱－水－力多場耦合數(shù)學模型，采用LAGAMINE多場耦合軟件，完成了緩沖材料試驗臺架熱－水－力多場耦合數(shù)值模擬。

試驗臺架中部水平方向溫度隨時間的變化趨勢如圖3所示。由圖中可以發(fā)現(xiàn)，在試驗過程中溫度不斷升高，并且距離加熱器越遠，加熱器的干燥作用越小，升溫幅度趨于平緩。另外，數(shù)值結(jié)果表明熱邊界單元可以很好的模擬在熱對流作用下試驗設備外部溫度不斷升高的現(xiàn)象。

表1 高廟子膨潤土基本物理參數(shù)［14］

表2 高廟子膨潤土力學參數(shù)［14］

圖3 水平向溫度分布隨時間變化趨勢圖

在水平方向上，吸力分布隨時間的變化趨勢如圖4所示。由圖4可以看出，隨著膨潤土的飽和，吸力值逐漸降低，且距離加熱器越近，降低的趨勢越明顯，特別是在加熱器附近區(qū)域，在試驗前期階段，吸力值大于初始吸力值，意味著在這一區(qū)域存在著先干燥再飽和的過程。

這一現(xiàn)象可以通過圖5更直觀的反映出來。如圖5所示，在加熱器附近的A點，吸力值首先增加，然后逐漸減小。此干燥－飽和過程可以通過水蒸氣蒸發(fā)引起的干燥作用和水的滲透導致的飽和作用的相互關(guān)系來解釋。在試驗初期，隨著溫度的升高，靠近加熱器的水分蒸發(fā)，并通過水蒸氣擴散向外側(cè)傳播。而由于極低的滲透系數(shù)，水的滲透作用相對較弱，因此整體上表現(xiàn)出干燥現(xiàn)象。隨著試驗的進行，裂隙水通過滲透作用逐步抵達這一區(qū)域并起主導作用，因此隨后表現(xiàn)出飽和現(xiàn)象。

圖4 水平向吸力分布隨時間變化趨勢圖

圖5 特征點吸力分布隨時間變化趨勢圖

水平方向上A、B、C三點處的膨脹力隨時間的變化趨勢如圖6所示。由圖6可以看出，在試驗初期膨脹力增加很快，隨后趨勢逐漸減弱。且靠近進水管C點處，由于飽和度較大，產(chǎn)生的膨脹力最大，三年后，膨脹力大概為2.5MPa左右?？紤]到飽和過程尚未結(jié)束，膨脹力上升曲線并未變得平緩，故該數(shù)值尚未達到膨脹力的最大值。

圖6 特征點膨脹力隨時間變化趨勢圖

3 結(jié)論

通過對緩沖材料中型試驗臺架模擬結(jié)果的定性分析，可以認為本數(shù)值計算方法可以很好地描述高廟子膨潤土的熱－水－力多場耦合特性。模擬結(jié)果得出以下結(jié)論。

1）溫度場隨時間變化逐漸升高，且隨著距離加熱器越遠，升溫趨勢逐漸變緩。

2）隨著膨潤土的飽和，土中吸力不斷降低，但是在試驗前期階段，由于加熱水蒸氣蒸發(fā)引起的干燥作用和水的滲透導致的飽和作用的相互影響，靠近加熱器的區(qū)域存在一個先干燥后飽和的過程。

3）在水的滲透作用下，高廟子膨潤土隨時間逐漸飽和。不過由于其極低的滲透系數(shù)，三年后膨潤土塊依然沒有完全飽和。三年后，膨脹力大約達到2.5MPa，膨脹力上升曲線并未變得平緩，故該數(shù)值尚未達到膨脹力的最大值。

這里需要說明的是，由于缺乏足夠的試驗數(shù)據(jù)，本研究僅進行了定性分析。因此，本數(shù)值模擬分析只可以被作為對提出的數(shù)值計算模型的初步驗證。隨著試驗的進行，將會展開進一步的數(shù)值研究。

［1］王駒，陳偉明，蘇銳，等.高放廢物地質(zhì)處置及其若干科學問題［J］.巖石力學與工程學報，2006，25（4）：801-812.

［2］潘自強，錢七虎.高放廢物地質(zhì)處置戰(zhàn)略研究［M］.北京：原子能出版社，2009.

［3］ZHANG Chunliang，Rothfuchs Tilmann，Wieczorek Klaus，et al.Monitoring and modeling of responese of the Opalinus clay to heating［J］.Chinese Journal of Rock Mechanics and Engineering，2006，25（4）：659-669.

［4］Li Xiangling，Bernier Frederic，Bel Johan.The Belgian HLW repository design and associated R＆D on the THM behavior of the host rock and EBS［J］.Chinese Journal of Rock Mechanics and Engineering，2006，25（4）：681-692.

［5］A.Lloret，M.V.Villar.Advances on the knowledge of the thermo-hydro-mechanical behaviour of heavily compacted“FEBEX”bentonite［J］.Physics and Chemistry of the Earth，2007，32，701-715.

［6］J.Samper，Liange Zheng，Luis Montenegro，et al.Coupled thermo-hydro-mechanical models of compacted bentonite after FEBEX in situ test［J］.Applied Geochemistry，2008，23，1186-1201.

［7］E.E.Alonso，J.Alcoverro，F(xiàn).Coste，et al.The FEBEX benchmark test：case definition and comparison of modeling approaches［J］.International Journal of Rock Mechanics ＆Mining Sciences，2005，42，611-638.

［8］M.Victoria Villar，Antonio Lloret.Dismantling of the first section of the FEBEX in situ test：THM laboratory tests on the bentonite blocks retrieved［J］.Physics and Chemistry of the Earth，2007，32，716-729.

［9］劉月妙，王駒，蔡美峰，等.高放廢物處置庫近場熱－水－力耦合長期穩(wěn)定性分析［C］.第二屆廢物地下處置學術(shù)研討會論文集，2008.

［10］葉為民，萬敏，陳寶，等.非飽和膨潤土水－力學性狀的溫度效應研究進展［C］.第二屆廢物地下處置學術(shù)研討會論文集，2008.

［11］Alonso.E.E，Gens.A，Josa.A.A constitutive model for partially saturated soils［J］.Geotechnique，1990，40（3）：405-430.

［12］CHEN Bao，QIAN Lixin，YE Weimin，et al.Soil-water characteristic curve of Gaomiaozi bentonite［J］.Chinese Journal of Rock Mechanics and Engineering，2006，25（4）：788-793.

［13］Ye W M，Cui Y J，Qian L X，Chen Bao.An experimental study of the water transfer through confined compacted GMZ bentonite［J］.Engineering Geology，2009b，108（3/4）：169-176.

［14］Liang Chen，Ju Wang，Yuemiao Liu，et al.Numerical thermohydro-mechanical modeling of compacted bentonite in Chinamock-up test for deep geological disposal［J］.Journal of Rock Mechanics and Geotechnical Engineering，2012，4（2）：183-192.

［15］Cui yujun，Tang Anhminh，Qian Lixin，et al.Thermalmechanical behavior of compacted GMZ bentonite［J］.Soil and Foundations，2011，51（6）：1065-1974.