(北京交通大學(xué)土木建筑工程學(xué)院，北京，100044)

對于高放廢物的處置，建造地下處置設(shè)施進行隔離是目前國際上普遍接受且技術(shù)上可行的方案[1]。處置系統(tǒng)的安全性取決于工程屏障系統(tǒng)(廢物體、廢物罐、緩沖材料和回填材料等子系統(tǒng))和天然屏障系統(tǒng)(近場、遠(yuǎn)場、地下水、生物圈和環(huán)境等子系統(tǒng))的綜合阻隔能力。其中，地下水子系統(tǒng)由巖體中的孔隙水和裂隙水組成，是核污染物運移的主要載體。在處置系統(tǒng)運行期間，高放廢物會在較長時間內(nèi)發(fā)生衰變并產(chǎn)生很高的熱量，使得近場圍巖溫度升高，改變其物理、化學(xué)性質(zhì)，可能形成傳熱(T)、滲流(H)、應(yīng)力(M)及化學(xué)反應(yīng)(C)的多場耦合作用，降低多重屏障系統(tǒng)的隔離功能。作為中國高放廢物處置庫首預(yù)選區(qū)，甘肅北山的主巖為致密的花崗巖，其孔隙極不發(fā)育，但有稀疏裂隙分布[2]，因此，裂隙水流場成為影響北山處置庫阻滯能力的關(guān)鍵因素之一。要評價處置主巖的性能，需要研究滲流-傳熱耦合條件下近場裂隙巖體的響應(yīng)。一般而言，裂隙巖體滲流-傳熱過程涉及非線性的數(shù)學(xué)問題，通常采用數(shù)值方法求解。例如，對于單一孔隙介質(zhì)，張玉軍[3]建立了適合飽和與非飽和孔隙介質(zhì)的二維THM耦合方程并給出了相應(yīng)的有限元求解方法，對高放廢物處置庫近場圍巖TMH 耦合過程進行了初步模擬研究，分析了溫度、飽和度和孔隙水壓力隨時間的變化趨勢，并把二維分析拓展到三維分析[4]；對于雙重孔隙介質(zhì)，張玉軍等[5]建立了二維THM 耦合方程及其有限元計算方式。但是，上述計算方式不能很好地處理巖石滲透率極低而巖體卻有稀疏裂隙分布的滲流-傳熱問題。鑒于此，基于二維等效孔隙介質(zhì)模型，柴軍瑞[6]建立了考慮溫度和滲流雙向耦合作用的有限元數(shù)學(xué)模型，但由于溫度與滲透系數(shù)的關(guān)系是通過經(jīng)驗公式來表述，因而，該方法對于裂隙網(wǎng)絡(luò)介質(zhì)的求解存在一定的局限性；劉學(xué)艷等[7]利用積分有限差分程序TOUGH2模擬了規(guī)則裂隙米尺度巖體的加熱/注水試驗，進一步分析了恒定熱源作用下巖體中TH耦合特征的演化規(guī)律，但沒有考慮熱源強度隨時間的變化。對于離散裂隙網(wǎng)絡(luò)模型，王如賓等[8-9]假定裂隙水流溫度沿流動方向呈線性變化，采用解析法和FEPG有限元軟件分別計算了單裂隙水流對巖體溫度場的影響，實質(zhì)上施加了裂隙水流溫度已知的邊界條件。此外，由于需要對整個物理模型進行離散，現(xiàn)有的商用軟件被用于高放廢物處置庫的多場耦合分析時對計算機的存儲能力和計算性能要求很高。為了更加有效地描述多場耦合過程，學(xué)者們采用了半解析方法對高放廢物處置裂隙巖體多場耦合的簡化模型進行了求解。例如，項彥勇等[10]考慮巖石中的二維熱傳導(dǎo)，建立了單裂隙巖體二維滲流-傳熱模型，提出了巖石和裂隙水溫度的拉氏變換半解析方法，通過數(shù)值拉氏逆變換得到時間域上的解答，之后把二維分析拓展到了三維分析[11-12]。拉氏變換半解析方法雖然避免了對時間積分，但存在計算過程復(fù)雜和計算精度較低等問題。為此，本文作者將裂隙巖體中的水-巖熱交換集度視作虛擬熱源，利用瞬時熱源作用下巖石溫度的基本解，建立關(guān)于時間和空間的邊界積分方程，通過卷積積分將其簡化為空間域上的積分方程，應(yīng)用有限元方法離散裂隙水的熱對流方程，建立線性代數(shù)方程組，經(jīng)過迭代計算直接得到任意時刻裂隙水和巖石的溫度，稱之為時域半解析計算方法。本文先計算無內(nèi)熱源時裂隙水的溫度，并與拉氏變換半解析解進行對比，最后分析衰變熱源作用下裂隙巖體的滲流-傳熱過程。

1 概念模型與控制微分方程

本研究基于以下假定條件[9]：

1)視巖體為無限大介質(zhì)，忽略巖石本體的滲透性，地下水僅在裂隙內(nèi)流動；

2)將裂隙視為平行板狀窄縫，裂隙長度遠(yuǎn)大于隙寬；

3)巖石和裂隙的熱物理性質(zhì)都不隨溫度變化而改變，其中裂隙內(nèi)水流為穩(wěn)定不可壓層流。

單裂隙巖體三維滲流-傳熱模型如圖1所示。圖1中，裂隙面A的邊界可以為任意形狀，裂隙開度為w；裂隙水的單寬流量為q，質(zhì)量密度為ρw，比熱容為cw、導(dǎo)熱系數(shù)為λw；裂隙兩側(cè)巖石Ω的質(zhì)量密度為ρr，比熱容為cr、導(dǎo)熱系數(shù)為λr；巖石內(nèi)部存在線熱源Γ，其強度為H~；裂隙水入流截面Aw上的溫度為T*；裂隙水和巖石的初始溫度、巖體無窮遠(yuǎn)處的溫度都為T0。

假設(shè)裂隙平行于x-y平面。由于裂隙開度很小，取z= 0平面為水-巖交界面。采用水-巖局部熱平衡假定，忽略熱存儲、熱傳導(dǎo)和熱流散效應(yīng)，則裂隙水的熱能守恒方程為[13]

圖1 單裂隙巖體三維滲流-傳熱模型Fig.1 Three dimensional water flow and heat transfer model of single planar fractured rocks

式中：(x,y)∈A；?2為二維哈密爾頓算子；h為水-巖熱交換集度；t為時間；T為溫度。

考慮三維熱傳導(dǎo)，則巖石的熱能守恒方程為[14]

式中：?23為三維拉普拉斯算子。

模型的初始條件和邊界條件為

2 積分方程

2.1 時空域邊界積分方程

根據(jù)三維熱傳導(dǎo)理論，t'時刻位于巖石基點X'=(x',y',z')處的單位點熱源在t＞t'時刻巖石場點X=(x,y,z)處產(chǎn)生的溫增θih(X-X',t-t')為[14]

將水-巖熱交換集度視為虛擬熱源，并考慮線熱源放出的熱量，則點X的溫度為

2.2 時間域離散

把0 ～t離散為N個間隔，認(rèn)為水-巖熱交換集度和線熱源強度在每個間隔期間保持不變，則

式(5)可表示為

式中：tn=nΔt，Δt=t N；h(X,tn)和H~(X~,tn)分別為tn-1～tn期間的水-巖熱交換集度和線熱源強度。

對式(6)中的t'運用卷積積分，則有

式中：θchN-n+1(X-X')為θih(X-X',t-t')在tN-n～tN-n+1上的積分，即為余誤差函數(shù)；αr=λr/(ρrcr)為巖石的熱擴散系數(shù)。

2.3 空間域離散

2.3.1 關(guān)于裂隙面的積分

裂隙面網(wǎng)格劃分如圖2所示，把裂隙面A離散為M個四邊形等參單元，共包含K個節(jié)點。

圖2 裂隙面網(wǎng)格劃分Fig.2 Meshing of the fracture planar

根據(jù)雙線性插值理論，單元內(nèi)部的物理量可以由節(jié)點值近似表示為[15]

式中：下標(biāo)m和l分別為單元編號和單元節(jié)點編號；Nl，T1m，hml和qml分別為l節(jié)點對應(yīng)的插值形函數(shù)、裂隙水溫度、水-巖熱交換集度和裂隙單寬水流量。

插值形函數(shù)取Lagrange多項式，即[15]

將式(8)和式(9)代入式(7)，則有：

式中：Am表示第m個單元的區(qū)域；N=[N1,N2,N3,N4]， 為插值形函數(shù)橫向量；hm=[hm1,hm2,hm3,hm4]T，為單元節(jié)點水-巖熱交換集度列向量。

令

式中：h=[h1,…,hk,…,hK]T，為所有節(jié)點水-巖熱交換集度列向量；Xk為第k個節(jié)點的坐標(biāo)。由于矩陣C1和C2中的核函數(shù)含有1/r，當(dāng)r→0時積分出現(xiàn)奇異，需要用特殊方法進行處理。

1)不包含奇點的單元積分。坐標(biāo)變換如圖3所示。當(dāng)計算點(x,y)位于m單元外時，利用式(13)[15]把全局坐標(biāo)系下的Am區(qū)域轉(zhuǎn)換成局部坐標(biāo)系下的Aˉm區(qū)域。

圖3 坐標(biāo)變換示意圖Fig.3 Diagram of coordinate transformation

采用高斯數(shù)值法[15]，則矩陣C1和C2中關(guān)于不包含奇點單元的積分為

式中：Δxm=xm2-xm1；Δym=ym2-ym1；Ψ=(ξ,η)，為X的局部坐標(biāo)；Ψij=(ξi,ηj)為高斯點的坐標(biāo)；φi和φj為高斯點Ψij的權(quán)重；Nξ和Nη分別為ξ和η的不同取值數(shù)。

2)包含奇點的單元積分。若(x,y)=(xm1,ym1)為計算點，則該點為奇點，以此為例說明矩陣C1和C2中的奇異積分。

包含奇點的單元再劃分如圖4所示。以(xm1,ym1)為圓心、rc為半徑，將m單元劃分為扇形區(qū)域和不規(guī)則區(qū)域A2m這2個部分。

矩陣C1和C2中關(guān)于奇點單元的積分為

圖4 包含奇點的單元再劃分Fig.4 Remeshing an element including a singular point

在極坐標(biāo)系下，式(15)中關(guān)于A1m的積分為

根據(jù)式(13)，區(qū)域A2m的y取值范圍Y(ξ)在局坐標(biāo)系下為

對式(15)關(guān)于A2m的積分進行高斯數(shù)值計算[15]：

對于矩陣C1和C2在其他包含奇點單元的積分，可以采用上述類似的方法計算。

2.3.2 關(guān)于線熱源的積分

把線熱源Γ進行F等分，用{0,…,f,…,F}對節(jié)點依次編號。根據(jù)Simpson 積分，各節(jié)點的權(quán)重為[10]

把式(19)代入式(7)，則有：

式中：下標(biāo)f為0,1,…,F，為曲線節(jié)點編號；為節(jié)點f的坐標(biāo)；為節(jié)點f放熱強度。

令

將式(7)應(yīng)用于裂隙面上的所有節(jié)點，則有

式中：T(X,t)=[T(X1,t),…,T(Xk,t),…,T(XK,t)]T和E=[T(X1,0),…,T(Xk,0),…,T(XK,0)]T分別為節(jié)點t時刻溫度列向量和節(jié)點初始溫度列向量。

3 數(shù)值求解

3.1 有限元方程

令

式中：qm=[qm1,qm2,qm3,qm4]T為單元節(jié)點單寬裂隙水流量列向量；-NT=[-Nx,-Ny]T為修正系數(shù)，用以減小對流占優(yōu)問題存在的數(shù)值振蕩。

修正系數(shù)為[16]

式中：Δxk為單元沿xk方向的邊長；為單元中心xk方向的單寬裂隙水流量；kw為水力擴散率。

對式(1)進行有限元離散（裂隙面網(wǎng)格劃分同2.3.1節(jié)），則有

3.2 迭代計算

將式(23)代入式(27)可得

采用高斯賽德爾迭代計算式(28)（取容差為1%），得到不同時刻的水-巖熱交換集度，再代入式(23)，即可求出任意時刻裂隙水和巖石的溫度。對于稀疏裂隙巖體系統(tǒng)，亦可采用同樣的思路進行求解。

4 算例分析

4.1 方法驗證

單一矩形裂隙巖體滲流-傳熱模型如圖5所示，長L、寬W的矩形裂隙沿x-y平面延伸；裂隙水沿x軸方向均勻流動，其入流溫度恒為20°C；裂隙水和巖石的初始溫度均為90°C；其他參數(shù)見表1。當(dāng)無線熱源時，本文采用時域半解析法計算裂隙水溫度，與一維解析解[17]和三維拉氏變換半解析解[11]作對比以驗證該方法的可靠性。選擇圖5中線段lf(x=800 m，y=40 m，z= 0 m)為觀察對象，定義量綱一溫差為(T0-T)/(T0-Tin)。3 種方法的計算結(jié)果如圖6所示。

由圖6可見，時域半解析解與解析解的吻合程度比拉氏變換半解析解的高。原因在于這二種方法的計算格式不同：拉氏變換半解析解需要通過數(shù)值逆變換得到時間域上的解答，但數(shù)值逆變換存在一定的計算精度問題[18]；而時域半解析解是時間域上的解答，避免了數(shù)值拉氏逆變換過程。因此，時域半解析方法的計算結(jié)果更精確，計算過程也更簡單。

4.2 核廢處置庫近場溫度

圖5 單一矩形裂隙巖體滲流-傳熱模型[11]Fig.5 Water flow and heat transfer model in rocks with a rectangular-shaped fracture[11]

圖6 裂隙水的量綱一溫差分布Fig.6 Normalized water temperature difference distribution in the fracture

在圖5中，設(shè)廢物罐線熱源沿Γ(x=50～100 m，y=40 m，z=5 m)均勻分布，其放熱強度為[19]

表1 計算參數(shù)Table1 Calculation parameters

在衰變熱源作用下，觀測線lf的裂隙水溫度和lr(x=0～800 m，y=40 m，z=10 m)的巖石溫度分別如圖7和圖8所示，圖中熱源溫差表示有熱源和無熱源2 種工況時溫度的差值。由圖7和圖8可見，裂隙水和巖石的溫度曲線都往裂隙水流方向傾斜，即上游熱量經(jīng)裂隙水傳輸至下游，加快了核廢處置庫近場熱量向遠(yuǎn)場傳遞。由于廢物罐的放熱功率不斷減小，裂隙水和巖石的溫度增幅都隨時間的推進而降低。

4.3 增強型地?zé)嵯到y(tǒng)溫度

圖7 裂隙水溫度隨衰變熱源的變化Fig.7 Change of water temperature with varying heat source

圖8 巖石溫度隨衰變熱源的變化Fig.8 Change of rock temperature with varying heat source

將單裂隙對井增強型地?zé)岣拍钅Ｐ蚚20]擴展至圖9所示的多裂隙系統(tǒng)，假設(shè)巖體中含有M'個完全相同的平行圓狀裂隙，回灌井和生產(chǎn)井都穿過每個裂隙的同1條直徑，兩井關(guān)于裂隙圓心對稱。

圖9 平行圓狀裂隙對井地?zé)嵯到y(tǒng)模型Fig.9 Dipole-well geothermal model with parallel discshaped fractures

根據(jù)單一圓狀裂隙的水流方程[20]，圖9中每個裂隙的單寬水流量為

式中：qx(X)和qy(X)分別為x方向和y方向的單寬裂隙水流量；Qi為作業(yè)井的注水量，以注水為正抽水為負(fù)，其中i為作業(yè)井?dāng)?shù)；(xi,yi)為i作業(yè)井圓心的平面坐標(biāo)；rw為作業(yè)井半徑。

為了分析裂隙空間分布對裂隙水流溫度的影響，設(shè)注水量Qin和抽水量Qex都為0.05 m3/s，回灌井與生產(chǎn)井的半徑都為0.1 m，兩井相距60 m，裂隙數(shù)為2 條，裂隙半徑都為50 m，裂隙開度都為0.001 m，3種工況的裂隙間距分別為20，35和50 m，其他參數(shù)同表1。由于該系統(tǒng)為對稱問題，任一工況裂隙中的水流-傳熱響應(yīng)相同。在同步回灌和抽采下，系統(tǒng)運行至10 a時3種工況的裂隙水流溫度分布如圖10所示。

由圖10可見：隨著裂隙間距增大，裂隙水流低溫區(qū)范圍逐漸縮小并且收縮速度逐漸下降，這是因為多裂隙內(nèi)的水-巖熱交換過程存在著相互影響，且影響程度與裂隙間距呈負(fù)相關(guān)。

圖10 運行至10 a時3種工況的裂隙水溫度分布Fig.10 Temperature distributions in different fracture planes at the operation time of 10 years

5 結(jié)論

1)相比于拉氏變換半解析方法，時域半解析方法避免了數(shù)值拉氏逆變換及該變換過程存在的誤差，具有計算結(jié)果更精確和計算過程更簡單的特點，還能用于分析含有平行稀疏裂隙巖體的水流-傳熱，其適用范圍更廣泛。

2)裂隙水將上游核廢物釋放的熱量傳輸至下游，加快了處置庫近場熱量向遠(yuǎn)場傳遞，有利于降低近場溫度。因此，對核廢處置庫近場溫度進行預(yù)測時，有必要考慮裂隙水流的影響。