基于溫度場-滲流場-應(yīng)力場耦合的雙重介質(zhì)流固熱耦合分析

1裂隙巖體耦合對巖體（尤其是斷裂巖體）的溫度場耦合研究具有重要的理論意義，并具有廣泛的實踐背景。高溫干熱巖體的開發(fā)、核廢料貯存庫的安全評估、環(huán)境工程中污染源的遷移及深部開采等都迫切需要研究巖體的溫度場–滲流場–應(yīng)力場三場耦合作用。從20世紀80年代以來,對該問題的研究取得了重要進展。J.Noorisha等提出了飽和裂隙巖體的THM耦合基本方程式,N.Barton和S.C.Bandis針對工程巖體對三場耦合進行初步探討;L.Jing等給出了相對系統(tǒng)的三場耦合模型,但對模型簡化和實用方面的研究還不夠深入。近十幾年來,提出了許多耦合模型,如單純的熱傳導(dǎo)模型,水流動模型、水熱耦合模型、固流熱耦合模型;考慮巖體介質(zhì)性態(tài),有均質(zhì)模型、裂隙介質(zhì)模型等。相比之下較全面的模型是T.Kohl等提出的固流熱耦合數(shù)學(xué)模型,T.A.Sprecker等提出并完善了裂隙巖體的溫度場、滲流場與應(yīng)力場耦合的有限元解法。當前,隨著能源儲庫、核廢料貯存庫及干熱巖體地熱開發(fā)等大型巖土工程的出現(xiàn),三場耦合研究成果已逐步服務(wù)于工程,用于評價工程的長期穩(wěn)定性、安全性[9~12]。國內(nèi)學(xué)者李寧等建立了一些耦合模型,但大多是基于多孔連續(xù)介質(zhì)或離散介質(zhì)假設(shè)。本文采用雙重介質(zhì)(擬連續(xù)介質(zhì)+裂隙介質(zhì))提出研究裂隙巖體溫度場–滲流場–應(yīng)力場耦合數(shù)學(xué)模型,開發(fā)了三場耦合數(shù)值模擬軟件,通過與滲流場–應(yīng)力場耦合,應(yīng)力場–溫度場耦合進行數(shù)值比較可知,進行雙重介質(zhì)溫度場–滲流場–應(yīng)力場耦合分析是十分必要的。2裂隙巖體膜滲流規(guī)律模型裂隙巖體是巖土工程中廣泛遇到的一類復(fù)雜巖體,是由隨機分布的裂隙和被裂隙切割的巖塊組成的不連續(xù)介質(zhì)。裂隙巖體溫度場–滲流場–應(yīng)力場耦合作用是一個相對復(fù)雜的研究問題,主要表現(xiàn)在裂隙巖體賦存地質(zhì)環(huán)境的各個組成部分,即滲流場、應(yīng)力場與溫度場都隨時間、空間而發(fā)生變化;與此同時,以上各個組成部分之間的耦合作用處于一種復(fù)雜動態(tài)變化過程中。傳統(tǒng)的離散介質(zhì)假設(shè)(忽略巖塊孔隙系統(tǒng)的透水性)和連續(xù)介質(zhì)假設(shè)(將裂隙中的滲流平均到巖體的滲流中去)往往很難真實地反映裂隙巖體的溫度場—滲流場–應(yīng)力場耦合作用,而雙重介質(zhì)(擬連續(xù)介質(zhì)+裂隙介質(zhì))三場耦合模型在一定程度上能真實地模擬裂隙巖體滲流、傳熱過程,對大型裂隙(斷層、節(jié)理、人工裂縫等)可以查明,按裂隙介質(zhì)對待;對被這些大中型裂隙切割而成的含眾多低序次的小裂隙介質(zhì)則不可能查明,認為是不確定的,采用擬連續(xù)介質(zhì)模型模擬。然后根據(jù)兩類介質(zhì)接觸處水頭、溫度、位移相等來建立裂隙巖體的三場耦合模型。本文建立的三場耦合模型基于如下假設(shè):(1)巖體是由擬連續(xù)介質(zhì)、裂隙介質(zhì)組成的結(jié)構(gòu)體,擬連續(xù)介質(zhì)為均質(zhì)同性的彈塑性體,裂縫介質(zhì)服從節(jié)理單元模型。(2)滲流規(guī)律服從非線性達西定律。(3)不考慮水的汽化,巖體被單相水飽和。(4)熱質(zhì)在固相、液相介質(zhì)中傳遞方式以傳導(dǎo)和強迫對流為主,各介質(zhì)的比熱容及熱傳導(dǎo)系數(shù)不隨溫度和壓力而變化。(5)介質(zhì)的變形以孔隙變形為主。2.1裂隙介質(zhì)擬連續(xù)介質(zhì)及kfi式中:e為體積變形;βT,βP分別為流體熱膨脹系數(shù)和壓縮系數(shù)。對于擬連續(xù)介質(zhì):對于裂隙介質(zhì):式(3)~(7)中,kij,kfi分別擬連續(xù)介質(zhì)和裂隙介質(zhì)滲透系數(shù);Θe為有效體積應(yīng)力;b為裂隙寬度;v(T)為運動黏滯系數(shù);φ為孔隙率;qi為比流量;p為水壓;si為切向坐標。2.2耦合的能量平衡方程根據(jù)熱力學(xué)第一定律:外界傳入系統(tǒng)的能量與內(nèi)部熱源產(chǎn)生的能量之和等于物質(zhì)內(nèi)能的增量與對外做功之和,建立分別考慮巖體(固相)與溫度場、液相與溫度場耦合的能量平衡方程,該方程考慮了液體熱對流(強迫對流為主)和固相熱傳導(dǎo)而導(dǎo)致的能量遷移。對于擬連續(xù)介質(zhì),有裂縫中主要介質(zhì)為水,能量的遷移是通過裂隙水流的熱傳導(dǎo)、熱對流(強迫對流為主)及巖體裂隙表面和裂隙水流之間的熱傳導(dǎo)的雙過程實現(xiàn)的,即式(8),(9)中:λr,λw分別為巖體和水的熱傳導(dǎo)系數(shù);ρr,ρw分別為巖體和水的密度;Cpr,Cpw分別為巖體和水的熱容系數(shù);Tr,Tw分別為巖體和水的溫度;Trb為巖體裂縫邊緣的溫度。2.3溫度場–滲流場—靜力平衡控制方程對于擬連續(xù)介質(zhì),靜力平衡控制方程為對于裂隙介質(zhì),靜力平衡控制方程為式(10)~(12)中:λ,μ均為拉梅常數(shù);uj,ji,ui,jj均為巖體位移;F為體力;β為巖體熱膨脹系數(shù);σn′,σs′分別為裂縫法向和切向有效應(yīng)力;kn,ks分別為裂縫法向與切向剛度;εn,εs分別為裂縫法向與切向變形;α為等效孔隙壓系數(shù)。式(10)~(12)這3組控制方程,再輔以初始和邊界條件,就構(gòu)成了基于雙重介質(zhì)假設(shè)的裂隙巖體溫度場–滲流場–應(yīng)力場耦合數(shù)學(xué)模型。這一模型具有如下特點:流體滲流方程中考慮了巖體的變形和溫度對滲透系數(shù)與水份遷移的耦合影響。能量守恒控制方程(式(8))中左邊第1項為溫度變化引起的能量變化;第2項為熱傳導(dǎo)引起的能量變化,第3項為熱對流引起的能量變化,第4,5項分別表示熱–液耦合和熱–固耦合項。固體變形方程式(10)中(αp,)i,βT,i描述滲流和溫度對變形的耦合作用,如果是等溫過程,能量方程不必要,在滲流方程和應(yīng)力本構(gòu)方程中略去含溫度T的各項,則方程組退化為流–固耦合滲流的控制方程。如果只考慮熱–固耦合,則在靜力平衡方程中沒有(αp,)i項。在能量守恒控制方程式(8)中去掉左邊第3項熱對流引起的能量變化,第4項熱–液耦合項即可。圖1給出了溫度場–滲流場–應(yīng)力場耦合模型(以擬連續(xù)介質(zhì)為例忽略變形生熱)。3雙介質(zhì)三場耦合模型的數(shù)值分析3.1高溫巖體計算模型雙重介質(zhì)滲流場–應(yīng)力場–溫度場耦合非常復(fù)雜,基本求解策略是將不同介質(zhì)的固體變形、能量方程和滲流看作獨立的子系統(tǒng),根據(jù)兩類介質(zhì)接觸處水頭、溫度、位移相等,溫度、應(yīng)力、水頭在同一介質(zhì)內(nèi)及不同介質(zhì)間隨時間傳遞,而將各子系統(tǒng)耦合求解。對時間采用變時步差分方法,對空間采用有限元方法離散求解,對擬連續(xù)介質(zhì)采用空間等參單元,裂縫采用八節(jié)點等參節(jié)理單元,開發(fā)雙重介質(zhì)三場耦合模型的三維有限元程序。下面用本文開發(fā)的三維有限元計算程序?qū)Ω邷馗蔁釒r體地熱開發(fā)系統(tǒng)進行多場耦合數(shù)值模擬分析,探討多物理場之間的耦合效應(yīng)。圖2給出了地下2200m的高溫巖體的計算模型(人工控制致裂形成眾多近似平行的人工裂縫構(gòu)成了水滲流的裂隙網(wǎng)絡(luò)通道,模型底部為固定約束,四周取不透熱邊界)及有限元網(wǎng)格圖,(x,y)坐標系下注入井坐標為(37.5m,37.5m);生產(chǎn)井坐標為(162.5m,162.5m),計算參數(shù)見表1。本文具體模擬了3種工況:工況1:溫度場–滲流場–應(yīng)力場耦合;工況2:滲流場–應(yīng)力場耦合(不考慮熱應(yīng)力、對流項對溫度的影響和溫度變化引起水壓力的變化);工況3:溫度場–應(yīng)力場耦合(不考慮流體滲流)。選取裂隙面、過注入井和生產(chǎn)井的縱剖對角面及裂隙面上靠近注入井的5428#節(jié)點、靠近生產(chǎn)井的5598#節(jié)點,研究3a內(nèi)巖體應(yīng)力、溫度、水壓分布和裂隙單元的張開度變化及介質(zhì)滲透系數(shù)的變化。3.2圍巖應(yīng)力分析三場耦合作用下,水壓和溫度變化梯度產(chǎn)生附加節(jié)點載荷,注入的低溫水使得注入井附近圍巖溫度下降,圍巖發(fā)生收縮變形,引起巖體應(yīng)力上升,為協(xié)調(diào)圍巖收縮變形,裂隙單元應(yīng)力降低;水壓變化梯度產(chǎn)生附加節(jié)點載荷,導(dǎo)致研究區(qū)域的有效應(yīng)力低于原巖應(yīng)力。三場耦合作用下,巖體應(yīng)力總體上要小于原巖應(yīng)力,并隨滲流方向應(yīng)力逐步上升。圖3,4分別給出了三場耦合下,縱剖對角面和裂隙面上σz分布圖,擬連續(xù)介質(zhì)內(nèi)應(yīng)力為40.5~52.5MPa,而注入井附近裂隙單元應(yīng)力只有23.0MPa左右。(2)裂隙水土流失區(qū)的黃土壤壓射在三場耦合作用下,水在基巖(擬連續(xù)介質(zhì))中滲流速度遠不如裂隙介質(zhì)滲流快,主干裂隙系統(tǒng)的強導(dǎo)水性使水壓沿此通道迅速上升,在非穩(wěn)定滲流階段,裂隙水壓要高于擬連續(xù)介質(zhì)內(nèi)的水壓(見圖5),主干裂隙顯著地控制著研究區(qū)域的水壓分布。由于固定井口壓力的影響,使得井口附近區(qū)域裂隙水壓梯度大,中間區(qū)域水壓梯度少,表現(xiàn)為圖6中井口附近區(qū)域水壓等值線相對密集,中間區(qū)域水壓等值線相對稀疏。(3)儲層流變質(zhì)巖滲流場的耦合圖7,8分別給出了三場耦合下縱剖對角面和裂隙面上溫度分布圖。流體滲流過程也是能量的傳遞過程,能量在滲流過程中通過熱對流和熱傳導(dǎo),從巖體中傳到了流體中,能量被注入水提取,這體現(xiàn)了能量方程中的φCpwρwkiip,iiTw,i,λrTr,ii(擬連續(xù)巖石介質(zhì))和(ρwCpwTwkfip,i)i,2λr(Trb-Tw)/b(裂隙介質(zhì))項,在擬連續(xù)巖石介質(zhì)中由于滲透系數(shù)遠低于裂隙介質(zhì)的滲透系數(shù),孔隙水流的熱對流速度遠不如裂隙水流的熱對流速度,使裂隙介質(zhì)的溫度比它毗鄰的擬連續(xù)巖石介質(zhì)的溫度要低(見圖7)。這體現(xiàn)了滲流場對溫度場的耦合。由于主干裂隙的存在使水頭與溫度分布呈現(xiàn)出各向異性,因而這是單純連續(xù)介質(zhì)耦合模型無法描述的。4雙介質(zhì)三場耦合效應(yīng)分析(1)滲流場及其應(yīng)力場耦合分析圖9~12為不同耦合工況下圍巖應(yīng)力分布情況。在工況3中,溫度場只對注入井附近圍巖應(yīng)力產(chǎn)生顯著影響,該區(qū)域溫度降低梯度大,產(chǎn)生較大的熱應(yīng)力,注入井附近80m區(qū)域內(nèi)圍巖應(yīng)力明顯上升,高于原巖應(yīng)力0%~12%,裂隙介質(zhì)為協(xié)調(diào)巖體的收縮變形,出現(xiàn)了應(yīng)力降低區(qū),在注入井附近裂隙單元應(yīng)力僅為45MPa(見圖9,10);考慮工況2時,滲流場對圍巖應(yīng)力的作用是全局性的,滲流場對應(yīng)力場的耦合效應(yīng)通過修正太沙基有效應(yīng)力來實現(xiàn),巖體有效應(yīng)力隨滲流方向逐步上升(見圖11,12)。在三場耦合分析下,存在復(fù)雜的耦合路徑如滲流場和溫度場之間存在熱對流耦合作用(見圖1)。因此,三場耦合分析下圍巖應(yīng)力不是圖9,11的簡單疊加。在注入井附近區(qū)域內(nèi),溫度場、滲流場及應(yīng)力場的三場耦合作用是顯著的,而靠近生產(chǎn)井區(qū)域,滲流場–應(yīng)力場耦合作用控制著此區(qū)域內(nèi)應(yīng)力分布。由此可見,一定空間、時間域內(nèi),往往存在起主導(dǎo)控制作用的雙場耦合系統(tǒng)。對于裂隙介質(zhì)來說滲流場、溫度場對應(yīng)力場的耦合作用可歸結(jié)為裂隙水壓、溫度梯度產(chǎn)生附加節(jié)點載荷導(dǎo)致裂隙張開度的變化。從圖13中可以看出:裂縫張開度隨時間的延長總體呈增加趨勢,在t=1a內(nèi),裂縫張開度急驟增加,注入井附近5428#裂隙單元三場耦合分析,裂隙張開度變化最大?b=0.0356cm,其次是滲流場–應(yīng)力場耦合分析,?b=0.0288cm,t=1a后裂隙張開度趨于穩(wěn)定;而應(yīng)力場–溫度場耦合分析時,t=5100h該裂隙單元張開度變化僅為?b=0.0085cm,在5100h后該裂隙單元張開度趨于穩(wěn)定。這顯示了滲流場對裂隙應(yīng)力場的強耦合效應(yīng)。三場耦合分析下,生產(chǎn)井附近5598#裂隙單元水壓、溫度梯度產(chǎn)生的附加載荷相對較小,在t=1a內(nèi),該裂隙單元張開度變化量?b=0.0192cm,1a后,隨著圍巖應(yīng)力的調(diào)整,該區(qū)域內(nèi)裂隙單元的應(yīng)力略有增加致使單元張開度稍有下降,t=3a時,?b=0.0185cm。(2)裂隙單元分布應(yīng)力場對滲流場的耦合作用主要體現(xiàn)在擬連續(xù)巖石介質(zhì)的滲透系數(shù)是有效體積應(yīng)力的負指數(shù)函數(shù);裂隙介質(zhì)的滲透系數(shù)是裂隙變形的函數(shù)。取z=60m的水平巖體面,在耦合工況1下,該面上的滲透系數(shù)分布如圖14所示,滲透系數(shù)從低應(yīng)力區(qū)的3.4×10-4m/s到高應(yīng)力區(qū)的0.7×10-4m/s變化;應(yīng)力場、溫度場變化導(dǎo)致裂隙張開度和流體黏性的改變,在耦合工況1,2下裂隙滲透系數(shù)總體上呈時間相依的增加趨勢。在t=10079h時對5428#裂隙單元進行三場耦合分析,裂隙滲透系數(shù)為0.87m/s,滲流場–應(yīng)力場耦合分析裂隙滲透系數(shù)僅為0.65m/s,為三場耦合分析所得裂隙滲透系數(shù)的74.7%(見圖15)。這是由于未考慮溫度場耦合作用所致,可見裂隙滲流場對溫度場的變化較為敏感;這顯示出了溫度場對裂隙滲流場的較強耦合作用。5流固熱三種耦合作用機理(1)雙重介質(zhì)溫度場–滲流場–應(yīng)力場全耦合數(shù)學(xué)模型包含了眾多耦合項的作用,并定義了不同介質(zhì)的本構(gòu)關(guān)系及耦合變量;該模型考慮了不同介質(zhì)間的水量和熱量交換,體現(xiàn)了主干裂隙的優(yōu)先滲流和強熱對流,具有較好的擬真性和操作性。(2)雙重介質(zhì)三場耦合路徑非常復(fù)雜,本文將固體變形、能量方程和滲流看作獨立的子系統(tǒng),根據(jù)兩類介質(zhì)接觸處水頭、溫度、節(jié)點位移相等的原則,采用有限元和有限差分法相結(jié)合,將具有耦合效應(yīng)的參數(shù)值在各子系統(tǒng)相互傳遞耦合求解。(3)建立了高溫干熱巖體地熱