王玉平，王 哲，易發(fā)成，曾志強

（1.宜賓學(xué)院國際應(yīng)用技術(shù)學(xué)部，四川宜賓 644000；2.西南科技大學(xué)環(huán)境與資源學(xué)院，四川綿陽 621010；3.宜賓學(xué)院理學(xué)部，四川宜賓 644000）

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

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

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

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

連續(xù)介質(zhì)理論認為，在物質(zhì)所占據(jù)的空間中，物質(zhì)是無間隙地連續(xù)分布的.可忽略物質(zhì)的微觀結(jié)構(gòu)，用偏微分方程表達宏觀物理量.混合物理論是在現(xiàn)代連續(xù)體物理的精神下發(fā)展起來的，提出了一個新的構(gòu)成公理“成分的等同性”，并用它來獲得與微觀對應(yīng)相一致的宏觀構(gòu)成方程.混合理論用于將連續(xù)介質(zhì)力學(xué)原理推廣到幾個可相互滲透的連續(xù)介質(zhì)來模擬多相系統(tǒng).為構(gòu)建彈性力學(xué)基礎(chǔ)上的多孔介質(zhì)THM耦合模型，作如下假設(shè)：

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

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

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

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

（5）多孔介質(zhì)的變形滿足小變形假設(shè)，Terzaghi有效應(yīng)力原理對非飽和介質(zhì)適用.

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

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

1.1 動量守恒方程

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

1.2 質(zhì)量守恒方程

（1）固相骨架質(zhì)量守恒方程

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

（2）水體質(zhì)量守恒方程

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

（3）氣體質(zhì)量守恒方程

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

1.3 能量守恒方程

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

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

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

2.1 高放廢物處置概念庫

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

圖1 試驗臺架的數(shù)值計算模型

2.2 緩沖材料的相關(guān)參數(shù)和本構(gòu)方程

緩沖材料力學(xué)模型為線彈性模型，其性質(zhì)參數(shù)如表1所示.

表1 膨潤土性質(zhì)參數(shù)[12-14]

本構(gòu)關(guān)系如下：

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

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

2.3 流體的相關(guān)參數(shù)和本構(gòu)方程

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

表2 流體的相關(guān)參數(shù)

流體的相關(guān)本構(gòu)關(guān)系如下：

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

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

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

(3)水體的熱膨脹系數(shù)[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)液相遷移系數(shù)[19]ω=其中C為調(diào)節(jié)系數(shù)，取2×10-5；R=8.3144 J/(mol·K).

2.4 氣體的相關(guān)參數(shù)和本構(gòu)方程

氣體相關(guān)參數(shù)如表3所示.

表3 氣體的相關(guān)參數(shù)[20,21]

氣體的相關(guān)本構(gòu)關(guān)系如下：

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

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

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

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

3 THM耦合數(shù)值模擬結(jié)果及分析

3.1 溫度場變化過程

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

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

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

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

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

圖4 溫度隨時間的變化

3.2 飽和度變化過程

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

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

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

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

圖7 飽和度隨時間的變化

3.3 膨脹力變化過程

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

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

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

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

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

圖10 膨脹力隨時間的變化

4 結(jié)語

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