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

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

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

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

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

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

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

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

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

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

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

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

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

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

1.1 動(dòng)量守恒方程

式中：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為水體的相對(duì)滲透率，μl是水的動(dòng)力黏滯系數(shù)，pl為液體壓力，S為飽和度，clp為壓縮性系數(shù)，csp為飽和度的吸力影響系數(shù)，βsT為飽和度的溫度影響系數(shù)，βlT是熱膨脹系數(shù)，ω為液相遷移系數(shù)，psv為飽和蒸汽壓力，pv為蒸汽壓力.

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

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

1.3 能量守恒方程

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

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

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

2.1 高放廢物處置概念庫(kù)

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

圖1 試驗(yàn)臺(tái)架的數(shù)值計(jì)算模型

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

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

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

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

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

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

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

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

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

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

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

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

(2)水 體 動(dòng) 力 粘 滯 系 數(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ū)動(dòng)系數(shù)[12]：kgT=6×10-12(1-S).

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

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

3.1 溫度場(chǎng)變化過程

緩沖材料的溫度場(chǎng)變化趨勢(shì)如圖2所示，可見緩沖材料的溫度隨時(shí)間不斷升高.

圖2 緩沖材料溫度場(chǎng)分布變化

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

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

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

圖4 溫度隨時(shí)間的變化

3.2 飽和度變化過程

圖5 所示的緩沖材料的飽和度變化趨勢(shì)表明，在水平方向上，飽和度值從外向內(nèi)逐漸增大，在豎直方向變化規(guī)律一致，沒有出現(xiàn)因?yàn)榧铀绞胶退闹亓ψ饔茫瑢?dǎo)致緩沖層上、下部含水量差異較大的情況.在豎直方向上出現(xiàn)飽和度值中間大，兩頭小的輕微差異，這是因?yàn)槟Ｐ蛯由舷滤吔鐥l件為不透水邊界，上下部分的水體只能向中間滲透.

圖5 緩沖層飽和度隨時(shí)間變化

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

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

圖7 飽和度隨時(shí)間的變化

3.3 膨脹力變化過程

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

圖8 緩沖層膨脹力隨時(shí)間變化

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

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

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

圖10 膨脹力隨時(shí)間的變化

4 結(jié)語(yǔ)

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