一種在滲流傳熱耦合作用下巖體溫度的計算方法

一種在滲流傳熱耦合作用下巖體溫度的計算方法
【專利摘要】本發明涉及一種在滲流傳熱耦合作用下巖體溫度的計算方法，包括以下步驟：建立裂隙一維流動傳熱模型；確定需要計算溫度的時間和位置；計算s值；對應s建立關于拉氏變換域內裂隙水溫度的包含n個未知量的n維線性方程組；根據向量，計算拉氏變換域內巖層指定點的溫度；判斷是否每個s都計算完畢；根據拉氏變換域內巖層指定點的溫度，通過數值拉氏逆變換，計算時間域內裂隙水溫度和巖層指定點溫度；判斷是否每個時間和位置都計算完畢；輸出時間域內裂隙水溫度和巖層內指定點溫度。
【專利說明】
一種在滲流傳熱耦合作用下巖體溫度的計算方法
技術領域
[0001 ]本發明涉及一種在滲流傳熱耦合作用下巖體溫度的計算方法。
【背景技術】
[0002] 巖體(特別是裂隙巖體)是處于一定的地質環境中的，地下水、地應力和溫度是地 質環境中的三個主要的因素，且這三者之間相互聯系，相互作用，相互制約，形成巖體滲流 場，應力場，溫度場三場耦合效應。在研究三場耦合作用機制時一般先進行等效化處理，將 研究對象看作是連續介質或等效連續介質，在此基礎上先進行兩場相互作用的研究，進一 步發展到三場耦合。
[0003] 目前很多學者己經就此方面做出了研究：
[0004] 國外方面，Parts( 1966)用解析方法探討了在多孔介質中自由對流的規律，同時討 論了水平流對自由熱對流的影響;Harlan，R.L(1973)提出的土體凍結過程中的熱質迀移數 學模型，第一次將滲流場與溫度場耦合起來分析，將凍脹理論推向一個新的階段;Hodgki nson(1980)建立了核廢料圍巖水熱對流的數學模型，采用級數方法探討其運動規律，討論 水平流對自由熱對流的影響;Morrow(1981)進行了模擬核廢料環境條件的室內滲流試驗， 測定熱源溫度在200-300°C，軸壓在10、20MPa，圍壓為30、60MPa條件下的滲透性，發現滲透 率隨溫度而下降，原因在于花崗巖礦物成分在高溫下溶解度增加，遇低溫區溶解度下降及 熱膨脹所致;Wicken S(1984)用有限元法模擬核廢料圍巖地下水流動模式，并同時考慮了潛 水面及熱對流2種驅動力的作用;Johansen等（1988)針對在阿拉斯加費爾班克附近地區的 多年凍土區的試驗隧道，研究了多年凍土蠕變特性和多年凍區地下洞的開挖方法及特殊支 護方法;L 〇Well(1990)探討了熱彈性對平行裂隙及裂隙滲流的影響，認為因熱膨脹作用，裂 隙產生自閉，可驅動流體流動;Abdallah(1995)采用離散單元法研究了在不連續面由于巖 體與水體的熱傳導性不同，存在著巖體與水體有能量交換條件下的熱對流問題。
[0005] 國內方面，宋曉晨、徐衛亞對裂隙巖體進行了等效化處理，在一定假設和簡化的基 礎上將巖體視為飽和多孔介質，利用混合物理論，推導出裂隙巖體等效連續介質溫度場、滲 流場和應力場三場耦合的全耦合數學模型及其控制方程。
[0006] 張樹光基于傳熱學和滲流理論建立了深埋巷道圍巖的熱擴散數學模型，結合邊界 條件采用MATLAB對其進行了數值求解，獲得了在風流和滲流耦合作用下圍巖的溫度場和溫 度矢量呈對稱分布。研究結果表明，無滲流狀態下溫度場和溫度矢量呈對稱分布，風流速度 對溫度分布有明顯的影響，但不改變其對稱分布的狀態。滲流改變了溫度場和溫度矢量原 有的對稱分布的狀態，熱交換平衡區隨著滲流速度的增加，將向順滲流的方向移動。
[0007] 張玉軍等通過對一些巖體滲透性能與熱物理性能的等效連續化處理，初步建立了 巖體溫度場、滲流場、應力場耦合作用數學模型，得出裂隙中的滲流通過與巖體的熱量傳遞 影響著溫度的分布。

【發明內容】

[0008] 以上所述現有技術無針對飽和裂隙巖體中高放核廢料地下處置庫裂隙巖體滲流-傳熱耦合的計算方法。為解決上述問題，本發明的目的在于提供一種在滲流傳熱耦合作用 下巖體溫度的計算方法，本發明所述方法以飽和裂隙巖體中高放核廢料地下處置庫方案規 劃為背景，采用格林函數法，建立了半解析法計算流程，可以準確的描述裂隙巖體滲流-傳 熱耦合的真實情況，并且便于考慮如邊界熱源等與高放核廢料地下處置庫概念模型相關的 情況，為高放核廢料地下處置庫在飽和裂隙巖體中的規劃建設服務。
[0009] 本發明的目的是通過以下技術方案實現的：
[0010] -種在滲流傳熱耦合作用下巖體溫度的計算方法，包括以下步驟：
[0011] 步驟1，建立裂隙一維流動傳熱模型，所述裂隙一維流動傳熱模型具有以下定義： 兩側隔水巖層不滲透，水僅在裂隙內部沿裂隙方向做一維流動，且流動狀態為穩定層流;巖 層中既有垂直于裂隙方向的熱傳導也有平行于裂隙方向的熱傳導；不考慮水的相變，以及 溫度變化所帶來的流體、巖層的物性變化，且流體不可壓縮，裂隙開度恒定；以平行于裂隙 方向的方向為X軸，以垂直于裂隙方向的方向為y軸；
[0012] 定義以下參數:裂隙水和巖層的初始溫度1'(1，7,0)=1'()，裂隙計算起點的水溫丁 (0,0,1:)=1' ：1：^(1，7,1:)為相對溫度差，711)為1〇(1，7,1:)的拉氏變換；
[0013] 步驟2,確定需要計算溫度的時間t和位置(x，y);
[0014] 步驟3,根據s = ns*ln2/t，計算s值，其中，s為拉普拉斯變換系數，ns = l，2,. . .，N，N 為離散拉氏變換參數的個數；
[0015] 步驟4,對應s建立關于拉氏變換域內裂隙水溫度的包含η個未知量；Τα的η維線性 方程組，如下：
[0016] [Α]{η = !，)丨
[0017] 其中，[Α]表示η維線性方程組中的系數矩陣，{b}表示η維線性方程組的右端項向 量，通過該方程組求解出.ΙΓ丨，即為向量；
[0018] 步驟5:根據向量，計算拉氏變換域內巖層指定點的溫度；
[0019] 步驟6:判斷是否每個s都計算完畢，是則進行下一步，否則返回步驟3;
[0020] 步驟7:根據拉氏變換域內巖層指定點的溫度，通過數值拉氏逆變換，計算時間域 內裂隙水溫度和巖層指定點溫度；
[0021]
[0022]
[0023] 步驟8:判斷是否每個時間和位置都計算完畢，是則進行下一步，否則返回步驟2;
[0024] 步驟9:輸出時間域內裂隙水溫度和巖層內指定點溫度。
[0025] 進一步的，所述步驟4中的η維線性方程組具體為：
[0028] pw為水的密度，(^為水的比熱容，pr為兩側隔水巖層的密度，cr為兩側隔水巖層的 比熱容，Ar為兩側隔水巖層的熱傳導系數，uw為裂隙水流的平均速度，b為裂隙半開度，1\為 裂隙長度，η表示積分區間劃分的單元數，△ x = lx/n表示單元間隔的大小，xi = i △ x(i=0~ η)和Xj = j Δ χ( j = 〇~η)表示橫坐標；
[0029] Κο表示零階修正的第二類貝塞爾函數，即：
[0030]
， 為歐拉常數，I〇(x)為零階修正的第一類貝塞爾函數
[0031] 心表示一階修正的第二類貝塞爾函數，dKoUVdxi-Kdx);
[0032] ^表示積分權重，利用辛普森3/8公式的積分權重為：
[0033] Wj = 3/8,7/6,23/24，l，l，· · ·，l，l，23/24,7/6,3/8(j = 0~η)。
[0034] 進一步的，步驟5中拉氏變換域內巖層指定點的溫度的計算根據下式得到：
[0035]
〇
[0036] 進一步的，步驟7中所述數值拉氏逆變換根據Stehfest方法，函數f (t)的拉氏變換 i.：QQ ^ /(i)e'W，根據F(s)的值求得f(t)在t = T的值。
[0037] 本發明的有益效果為：
[0038] 通過本發明所述方法，可計算得到裂隙巖體滲流-傳熱耦合作用下任意點在時間 域內的溫度值。本發明是針對核廢料處置庫滲漏的滲流場與溫度場耦合簡化模型的格林函 數計算方法。
【附圖說明】
[0039]圖1為巖層二維熱傳導的滲流傳熱耦合計算模型；
[0040] 圖2為裂隙巖體滲流傳熱耦合模型計算的格林函數法時間域內的分布熱源圖；
[0041] 圖3為裂隙巖體滲流傳熱耦合模型計算的格林函數法拉氏變換域內的分布熱源 圖；
[0042] 圖4為巖層內沿y方向上各點溫度的變化圖；
[0043]圖5為巖層內y方向上三點溫度隨時間的變化圖；
[0044] 圖6為巖層內X方向上三點溫度隨時間的變化圖。
【具體實施方式】
[0045] 為了使本發明的目的、技術方案及優點更加清楚明白，以下結合【具體實施方式】和 附圖，對本發明進行進一步詳細說明。應當理解，此處所描述的具體實施例僅僅用以解釋本 發明，并不用于限定本發明。
[0046] -種在滲流傳熱耦合作用下巖體溫度的計算方法，包括以下步驟：
[0047] 1.1概念模型
[0048] 本發明中要計算的裂隙巖體滲流-傳熱耦合模型是考慮巖層二維熱傳導的情況， 即在巖層中既有垂直于裂隙方向的熱傳導也有平行于裂隙方向的熱傳導，如圖1所示：
[0049] 根據考慮巖層二維熱傳導的滲流-傳熱模型計算需要，現作如下假設：
[0050] 1、巖層不可滲透，流體僅在裂隙（或含水層）內部沿裂隙方向做一維流動，且流動 狀態為穩定層流；
[0051 ] 2、由于溫度場分布變化相對較小，不考慮水的相變，以及溫度變化所帶來的流體、 巖層的微小的物性變化(求解過程中視為常量）；裂隙開度(或含水層厚度)恒定;且流體不 可壓縮；
[0052] 3、系統內無內熱源，忽略粘性耗散和輻射傳熱，不計位能、動能的變化；
[0053] 4、建立控制方程時考慮熱彌散作用（熱彌散系數SDlx，從后面的實際分析中，可以 得知在層流條件下，熱彌散作用可以忽略不計）；
[0054] 5、若進一步忽略裂隙水溫度控制方程中隨時間的熱存儲項，則所得結 果可以進一步簡化。
[0055] 1.2控制方程及定解條件
[0056] 1、控制方程
[0057]對于圖1所示的裂隙一維流動傳熱模型，在控制方程方面，主要考慮兩側隔水巖層 中的二維熱傳導。所以在裂隙熱流的控制方面，忽略裂隙的熱存儲和熱擴散，則裂隙的熱量 守恒方程可以表示為：
[0058]
( 1 )
[0059 ]其中，pw為水的密度;cw為水的比熱;Ar為巖層的熱傳導系數;U W為裂隙水流的平均 速度;2b為裂隙開度;T(x，y，t)為溫度。
[0060]兩側隔水巖層的熱量守恒方程還考慮了與巖體裂隙相平行方向上的熱傳導，即：
[0061 ]
C2)
[0062] Pr為兩側隔水巖層的密度，Cr為兩側隔水巖層的比熱。
[0063] 2、定解條件
[0064] 考慮巖層二維熱傳導的滲流-耦合模型的初始條件中，裂隙水和巖層的初始溫度 被假定為連續的，即T (X，y，0) = Το。在t = 0+時水注入，裂隙計算起點的水溫T (0，0，t) = f。 隨著水在裂隙中流動，裂隙水溫度會隨位置變化，因此，裂隙水溫也是最終結果的很重要一 部分。
[0065] 上述初始條件及邊界條件可總結為：
[0066]
^3)
[0067] 1.3格林函數法與數值積分求解
[0068] 為了方便求解式(1)和(2)，首先引入相對溫度差TD(x，y，t)，即：
[0069]
(4)
[0070] 將相對溫度差TD(x，y，t)代入微分方程（1)、（2)及定解條件，并進行拉普拉斯變 換;變換后的微分方程(1 )、（2)為：
[0073] 式中，S--拉普拉斯變換系數；[0074] To--拉普拉斯變換后的TD(x，y，t)。[0075] 變換后的微分方程(5)、（6)要受以下裂隙水入口處的邊界條件控制：
[0071] (5)
[0072] 16)
[0076]
(7)
[0077] 此處使用到了拉普拉斯變換基本概念和微分性質，即：
[0078]
(.8)
[0079]
(9)
[0080] 式(5)、式(6)和式(7)是定義在二維坐標內的方程，要想求解它們，可以利用格林 函數法將它們轉化為一維積分形式。
[0081] 如圖2、3所示，在無限平面和拉氏變換域內，點P'（x'，y'）處的單位集度熱源（即時 間域內單位集度熱源的拉氏變換)在點P(x，y)處所產生的溫度匕（即時間域內相對溫度差 的拉氏變換)為：
[0082]
(1:〇)
[0083]式中，Ko表示零階修正的第二類貝塞爾函數;即：
[0084]
，其中，γ =0.5772157 稱為歐拉常數，I〇(x)為零階修正的第一類貝塞爾函數：
表示 熱源點(X'，y'）與任意點(x，y)的距離。
[0085]設有沿Γ分布的熱源，集度為q(X'，y'，t)|(x，, y，)eΓ;在拉氏變換域內，該分布熱源 在點P(x，y)處所產生的溫度可以表示為：
[0086]
[0087] 其中4?1，>) = 1{?(1',.1，',?)}為分布熱源集度9&'，7'，0的拉氏變換。
[0088] 此處利用上述格林函數理論，巖層中任意點的溫度就是熱源點所產生的場沿裂隙 長度lx的積分，即：
[0089]
(12)
[0090]式中，?為拉普拉斯變換后的裂隙-巖層界面處(忽略裂隙開度，取裂隙-巖層界面 為y=〇)的熱流密度。
[0091] 熱流密度^應該與裂隙邊緣的熱傳導相等，故：
[0092]
(13)
[0093]利用式(5)，式（13)可變為：
[0094](14)
[0095] 將式（14)代入式（12)可得：
[0096]
[0097] 對式（15)進行分部積分可得：
[0098]
[0099]
[0100]式中，心表示一階修正的第二類貝塞爾函數(dKoUVdxz-KKx));前述中提到的 邊界條件式(7)也在此處用到。
[0101] 將式(17)用到裂隙與巖層交界面處，即(y = 0,0彡X彡lx)可得：
[0102]
(18)
[0103]這是一個關于Γ0(λ·,0,4的積分方程（菲德海姆(Fredholm)第二類積分方程），可以 采用數值積分求解。但是，方程中的心含有一階（Ι/r)奇點，需要先正則化。為此，用 替代方程(18)中奇異積分的，解析積分得：
[0104
(丨 9)
[0105]用方程(19)減去方程(18)得： (20)
[0106] ' '
[0107] 可以看到，這個積分方程沒有奇點，[Γ"χ0,d- Γ0沁0,4消除了原來積分的奇點。
[0108] 為了求解方程(20)需要引入一種求積方法，一般來說，高斯求積法是最為有效的， 但為了使結果在x = 〇，lx時不用差值且足夠精確，此處應該使用辛普森求積公式，那么方程 (20)將變為：
[0110] 式中，η表示積分區間劃分的單元數；Δ χ=1χ/η表示單元間隔的大小;xk = kA x(k =0~η)表不橫坐標;wj表不積分權重，利用辛普森3/8公式的積分權重為：[0111] wj = 3/8,7/6,23/24，l，l，· · ·，l，l，23/24,7/6,3/8(j = 0~n) (22)[0112] 求積公式(21)中的一些節點（如Xl = Xj)時需要特殊處理，雖然-，消除了原來 積分的奇點，但是當xi = xj時，該積分成為〇/〇不定式，可取其極限值： (21)
[0109] i = \- a
[0113]
(23)
[0114] 式(23)可利用差分近似的方法進行簡化得：
[0115](2.4): Λ.-χ
[0116] 但當i = 1和i =η時不能使用中心差分而要使用向后差分和向前差分。
[0117] 根據以上分析式(21)可化為：
[0118]
[0119]
[0120] 式(25)可化為：

[0121]
[0122]
[0123] 注意到，式(26)是包含η個未知量維線性方程組，即：
[0124] [Α] |Γ； - {Λ； (21)
[0125] 其中，[Α]表示方程組中的系數矩陣，{b}表示方程組的右端項向量。該方程組不含 任何未知參量，所以可以直接通過該方程組求解出丨7'丨，CT)·即為向量，并使用拉式近似 逆變換Stehf est方法來求出時間域內的裂隙各點的溫度值。
[0126] 得到= 后，對式(17)進行數值積分，計算巖層內的b，得：
[0127]
(28)
[0128] 利用式(27)和式(28)得到的：ΓΛ;和都是拉普拉斯變換域內的值，要利用拉普拉 斯逆變換來得到對應的時間域內的值，此處使用的方法是一種拉普拉斯近似逆變換的方 法--Stehfest 方法。
[0129] Stehfest方法的基本思路是:假設函數f(t)的拉氏變換為F(s)，
[0130]
(29)
[0131] 倘若能夠計算F(s)，那么函數f(t)在t = T的值可由下式算得：
[0132] (30)
[0133]
[0134] (31)
[0135] 這樣，利用上述中的拉普拉斯近似逆變換方法即可得到裂隙巖體滲流-傳熱耦合 作用下任意點在時間域內的溫度值。
[0136] 實施例
[0137] 對于計算方法所要研究的內容，作為核廢料地下處置庫的圍巖研究，除了關注裂 隙中流體溫度隨位置和時間的變化外，更重要的是要了解巖層溫度隨位置和時間的變化， 圖4~6顯示的就是巖層溫度的計算結果。
[0138] 1、巖層溫度隨位置的變化
[0139] 此處研究的是巖層溫度在垂直于裂隙各點上的變化，此處選取的是x = 30m，即距 離入水口 30米且垂直于裂隙上的各點，分別計算出在1年、5年、10年時的溫度，計算結果如 圖4所示。
[0140] 從圖4可以看出，溫度沿y方向正方向不斷降低，直至巖層的初始溫度，越接近裂 隙，溫度變化越劇烈；時間越久，相同位置的溫度越高；在距離裂隙80米左右開始巖層溫度 變化基本可以忽略。
[0141] 2、巖層溫度隨時間的變化
[0142] 如圖5和圖6所示，研究的是巖層溫度隨時間的變化。對于圖5,此處選取了距離入 水口 50米且垂直于裂隙方向上的3個點，分別是(50,50)、（50,100)和(50,150)，分別計算出 它們的溫度隨時間的變化，時間變化的區間為100年。
[0143] 由圖5可見，巖層溫度隨時間的增加會有不同程度的升高，理論上直至升高到與裂 隙流體一致。由圖5可見，時間變化越短巖層溫度變化的幅度越明顯，距離裂隙越遠，巖層溫 度的變化越不明顯，從圖5可以看出，（50,150)位置的巖層溫度變化幅度最小并且在100年 范圍內的變化在3%以內。
[0144] 對于圖6,此處選取了距離裂隙50米且平行裂隙方向上的3個點，分別是（50,50)、 (100，50)和（150，50)。分別計算出它們的溫度隨時間的變化，同樣，時間變化的區間為100 年。
[0145] 由圖6可見，巖層溫度隨時間的增加會有不同程度的升高，理論上直至升高到與裂 隙流體一致。由圖6可見，時間變化越短巖層溫度變化的幅度越明顯，距離入水口越遠，巖層 溫度的變化越平緩，從圖上可以看出，（150,50)位置的巖層溫度變化幅度最小并且在100年 范圍內的變化在3%以內。
[0146]以上所述僅為本發明的較佳實施例而已，并不用以限制本發明，凡在本發明的精 神和原則之內所作的任何修改、等同替換和改進等，均應包含在本發明的保護范圍之內。
【主權項】
1. 一種在滲流傳熱禪合作用下巖體溫度的計算方法，其特征在于，包括W下步驟： 步驟1，建立裂隙一維流動傳熱模型，所述裂隙一維流動傳熱模型具有W下定義:兩側 隔水巖層不滲透，水僅在裂隙內部沿裂隙方向做一維流動，且流動狀態為穩定層流;巖層中 既有垂直于裂隙方向的熱傳導也有平行于裂隙方向的熱傳導;裂隙開度恒定；W平行于裂 隙方向的方向為X軸，W垂直于裂隙方向的方向為y軸； 定義W下參數:裂隙水和巖層的初始溫度1'^，7,0)=1'〇，裂隙計算起點的水溫1'(0,0，*) =1'*^〇山7,〇為相對溫度差，托，為1〇山7，*)的拉氏變換； 步驟2,確定需要計算溫度的時間t和位置(x，y); 步驟3,根據s = ns*ln2/t，計算S值，其中，S為拉普拉斯變換系數，ns = l，2, . . .，N，N為離 散拉氏變換參數的個數； 步驟4,對應S建立關于拉氏變換域內裂隙水溫度的包含n個未知量;T。,的n維線性方程 組，如下：其中，[A]表示n維線性方程組中的系數矩陣，{M表示n維線性方程組的右端項向量，通 過該方程組求解出171，即為T。,向量； 步驟5:根據向量，計算拉氏變換域內巖層指定點的溫度； 步驟6:判斷是否每個S都計算完畢，是則進行下一步，否則返回步驟3; 步驟7:根據拉氏變換域內巖層指定點的溫度，通過數值拉氏逆變換，計算時間域內裂 隙水溫度和巖層指定點溫度； 其中，步驟8:判斷是否每個時間和位置都計算完畢，是則進行下一步，否則返回步驟2; 步驟9:輸出時間域內裂隙水溫度和巖層內指定點溫度。2. 根據權利要求1所述的在滲流傳熱禪合作用下巖體溫度的計算方法，其特征在于，所 述步驟4中的n維線性方程組具體為：其中Pw為水的密度，Cw為水的比熱容，化為兩側隔水巖層的密度，(:r為兩側隔水巖層的比熱 容，、為兩側隔水巖層的熱傳導系數，Uw為裂隙水流的平均速度，b為裂隙半開度，Ix為裂隙 長度，n表示積分區間劃分的單元數，A X = Ix/n表示單元間隔的大小，Xi = i A X (i = O~n)和 Xj = j A x( j = 0~n)表示橫坐標； Ko表示零階修正的第二類貝塞爾函數，即：為歐拉常數，Io(X)為零階修正的第一類貝塞爾函數：Ki表示一階修正的第二類貝塞爾函數，dKo(x) Ak =-Ki(X); Wj表示積分權重，利用辛普森3/8公式的積分權重為： wj = 3/8,7/6,23/24，l，l,. . .，l，l，23/24,7/6,3/8(j = 0~n)。3. 根據權利要求2所述的在滲流傳熱禪合作用下巖體溫度的計算方法，其特征在于，步 驟5中拉氏變換域內巖層指定點的溫度的計算根據下式得到：O4. 根據權利要求1所述的在滲流傳熱禪合作用下巖體溫度的計算方法，其特征在于，步 驟7中所述數值拉氏逆變換根據Stehfest方法，函數f(t)的拉氏變換為F(S)，根據F( S)的值求得f (t)在t = T的值。
【文檔編號】G06F17/50GK106033485SQ201610454727
【公開日】2016年10月19日
【申請日】2016年6月21日
【發明人】郭家奇, 趙允濤, 程峰, 李鵬, 馬春莉, 毛輝, 王寬, 宋紅兵, 肖志均, 王小寶, 崔超, 龔劍鋒
【申請人】中國大唐集團科技工程有限公司