裂隙巖體滲流-傳熱的時(shí)域半解析計(jì)算方法

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

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

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

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

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

2)將裂隙視為平行板狀窄縫，裂隙長(zhǎng)度遠(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)部存在線熱源Γ，其強(qiáng)度為H~；裂隙水入流截面Aw上的溫度為T*；裂隙水和巖石的初始溫度、巖體無(wú)窮遠(yuǎn)處的溫度都為T0。

假設(shè)裂隙平行于x-y平面。由于裂隙開度很小，取z= 0平面為水-巖交界面。采用水-巖局部熱平衡假定，忽略熱存儲(chǔ)、熱傳導(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為時(shí)間；T為溫度。

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

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

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

2 積分方程

2.1 時(shí)空域邊界積分方程

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

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

2.2 時(shí)間域離散

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

式(5)可表示為

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

對(duì)式(6)中的t'運(yùn)用卷積積分，則有

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

2.3 空間域離散

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

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

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

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

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

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

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

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

令

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

1)不包含奇點(diǎn)的單元積分。坐標(biāo)變換如圖3所示。當(dāng)計(jì)算點(diǎn)(x,y)位于m單元外時(shí)，利用式(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)于不包含奇點(diǎn)單元的積分為

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

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

包含奇點(diǎn)的單元再劃分如圖4所示。以(xm1,ym1)為圓心、rc為半徑，將m單元?jiǎng)澐譃樯刃螀^(qū)域和不規(guī)則區(qū)域A2m這2個(gè)部分。

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

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

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

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

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

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

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

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

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

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

令

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

式中：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é)點(diǎn)t時(shí)刻溫度列向量和節(jié)點(diǎn)初始溫度列向量。

3 數(shù)值求解

3.1 有限元方程

令

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

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

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

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

3.2 迭代計(jì)算

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

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

4 算例分析

4.1 方法驗(yàn)證

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

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

4.2 核廢處置庫(kù)近場(chǎng)溫度

圖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)均勻分布，其放熱強(qiáng)度為[19]

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

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

4.3 增強(qiáng)型地?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

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

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

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

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

為了分析裂隙空間分布對(duì)裂隙水流溫度的影響，設(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)為對(duì)稱問題，任一工況裂隙中的水流-傳熱響應(yīng)相同。在同步回灌和抽采下，系統(tǒng)運(yùn)行至10 a時(shí)3種工況的裂隙水流溫度分布如圖10所示。

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

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

5 結(jié)論

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

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