增強(qiáng)型地?zé)嵯到y(tǒng)熱流固耦合過程數(shù)值模擬研究

許 欣，李 宏

(大連理工大學(xué) 海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024)

地?zé)崮茏鳛橐环N清潔、可再生能源，越來越受到各國政府的重視[1]。傳統(tǒng)的增強(qiáng)型地?zé)嵯到y(tǒng)是一種基于石油鉆井技術(shù)的地?zé)衢_采方法，技術(shù)的關(guān)鍵是通過水力壓裂技術(shù)在鉆孔的底部巖層中產(chǎn)生貫通裂紋。1973年，美國國家實(shí)驗(yàn)室和能源部在Fenton Hill開始了干熱巖發(fā)電技術(shù)試驗(yàn)，井深3 000 m，利用花崗巖內(nèi)的裂縫連接兩口井，但正是由于這些裂縫的存在，注入水的回收率大大降低，嚴(yán)重影響水和熱的回收和抽取[2]。Gringarten等[3]提出了一種從斷裂熱干巖石中提取熱量的理論。冷水通過每個(gè)裂縫的底部進(jìn)入，熱通過水流從干熱巖傳遞到裂縫巖體。這一數(shù)學(xué)模型可以大大提高干熱巖地?zé)嵯到y(tǒng)的經(jīng)濟(jì)利用率。

在許多地?zé)崽镏?，有證據(jù)表明注入的流體可能會(huì)沿著“裂縫”這一優(yōu)先流動(dòng)路徑快速遷移，且也有大量現(xiàn)場證據(jù)表明，在裂縫性儲層中，注入井的水可以快速地遷移到生產(chǎn)井[4-7]。Shaw建立了多井系統(tǒng)中干熱巖水流吸熱的數(shù)學(xué)模型，研究了多井系統(tǒng)的井距和井徑尺寸，對于地?zé)醿拥臒岵尚实挠绊慬8]。不僅如此生產(chǎn)井和注入井的位置同樣會(huì)對人工地?zé)醿拥牟蔁岢潭纫约靶十a(chǎn)生影響。

以上研究可以看出，地?zé)醿拥拈_發(fā)具有“埋藏深、工程難、投資大”的難點(diǎn)，由于基于鉆井技術(shù)的增強(qiáng)型地?zé)崮Ｐ?EGS-D)難以建立大規(guī)模穩(wěn)定蓄熱，水流量小，易造成污染等缺點(diǎn)，Tang等[9]提出了一種基于開挖技術(shù)的增強(qiáng)型地?zé)崮Ｐ?EGS-E)開發(fā)干熱巖。它可以在一定程度上克服水力壓裂技術(shù)的弊端，避免大高程抽水二次耗能，結(jié)合深部礦產(chǎn)資源開采開發(fā)深地?zé)?，提高綜合經(jīng)濟(jì)效率。本文以EGS-E中鉆爆致裂區(qū)為研究對象，連接并順序執(zhí)行兩個(gè)計(jì)算機(jī)程序TOUGH2和FLAC3D，然后用一維固結(jié)問題的解析解驗(yàn)證耦合程序的正確性，最后將程序應(yīng)用于厚度、孔隙度和滲透率接近均勻的多孔介質(zhì)理想儲層模型，模擬人工地?zé)醿幼⒉删辛黧w的快速運(yùn)移及注水井和產(chǎn)能井的位置變化對模型"熱掃掠"范圍的影響。

1 TOUGH2和FLAC3D 力學(xué)耦合水熱模型程序流程

1.1 水熱過程數(shù)學(xué)模型

TOUGH2的控制方程除了考慮巖體變形過程中孔隙度和滲透率的變化，還考慮了孔隙壓力和流體溫度對變形的影響。流動(dòng)和熱對流的一般控制方程是[10]：

(1)

其中，Mκ是模型單元中每單位體積的質(zhì)量或能量；Fκ是質(zhì)量或者能量通量；qκ是控制單元中源/匯項(xiàng)。

對于多相流體的流動(dòng)有：

Mw=φ(SlρlXwl+SgρgXwg)

(2)

Fw=Xwlρlul+Xwgρgug

(3)

其中，Mw是水的質(zhì)量；φ是孔隙率；S是飽和度；ρ是密度；X是質(zhì)量分?jǐn)?shù)；u是達(dá)西滲流中的速度；下標(biāo)w是水；l代表液相； g代表氣相。

對于熱對流的傳導(dǎo)過程：

Mh=φ(SlρlUl+SgρgUg)+(1-φ)ρsUs

(4)

(5)

式中：Mh是總熱量；U是內(nèi)部能量；λ是導(dǎo)熱系數(shù)；T是溫度。下標(biāo)h是熱焓；s代表固相。

1.2 力學(xué)過程數(shù)學(xué)模型

FLAC3D是由Itasca公司開發(fā)的商業(yè)軟件，可以對土壤、巖石和其他材料進(jìn)行熱流固耦合分析。 FLAC3D的力學(xué)模塊中的控制方程是[11]：

(6)

其中，ρm是巖體的平均密度；v是指定系統(tǒng)的速度。在一個(gè)時(shí)間步長內(nèi)的應(yīng)力應(yīng)變關(guān)系可寫為：

(7)

FLAC3D以節(jié)點(diǎn)為計(jì)算對象，力和質(zhì)量全都集中在節(jié)點(diǎn)上，然后通過運(yùn)動(dòng)方程使用時(shí)域有限差分法[12]和混合離散法[13]求解。節(jié)點(diǎn)運(yùn)動(dòng)方程可表示如下：

(8)

(9)

在FLAC3D中，為了確保精度和穩(wěn)定性，通常在運(yùn)動(dòng)方程和應(yīng)力應(yīng)變本構(gòu)方程之間設(shè)置足夠小的迭代步。在每個(gè)時(shí)間步長內(nèi)，首先根據(jù)初始速度和力進(jìn)行運(yùn)動(dòng)方程的計(jì)算，并且將計(jì)算的節(jié)點(diǎn)速度代入到方程(9)中以獲得新的應(yīng)力和應(yīng)變速率。

1.3 TOUGH2和FLAC3D的耦合

TOUGH2能夠?qū)崿F(xiàn)流體的流動(dòng)和熱量的傳輸，而不能反映變形與孔隙度和滲透率之間的關(guān)系，所以本文建立連接的TOUGH-FLAC3D模擬器來進(jìn)行熱流固耦合的計(jì)算。當(dāng)連接兩個(gè)代碼時(shí)，耦合方程不能同時(shí)被求解，而是順序求解。本研究通過連接外部耦合模塊，實(shí)現(xiàn)兩個(gè)代碼的順序執(zhí)行，并在實(shí)驗(yàn)室估計(jì)的經(jīng)驗(yàn)函數(shù)上，對各物理場的參數(shù)校準(zhǔn)，采用數(shù)學(xué)插值方法實(shí)現(xiàn)相互轉(zhuǎn)化。耦合函數(shù)通過連接兩個(gè)代碼，實(shí)現(xiàn)了最重要的耦合，包括巖石變形對孔隙度和滲透率的影響，以及孔隙壓力和溫度對巖石變形的影響。

在TOUGH2中，網(wǎng)格坐標(biāo)由單元的中心點(diǎn)定義。元素大小由相鄰單元格的中心點(diǎn)到連接平面的距離確定，而FLAC3D中的網(wǎng)格坐標(biāo)由角點(diǎn)坐標(biāo)來定義。由于TOUGH2和FLAC3D分別是以單元中心點(diǎn)和節(jié)點(diǎn)為計(jì)算單位的，為了實(shí)現(xiàn)TOUGH2和FLAC3D的耦合計(jì)算，需要通過插值對主要變量進(jìn)行由單元到單元節(jié)點(diǎn)轉(zhuǎn)換。耦合流程如圖1所示，可歸為以下四步：(1)在TOUGH2中進(jìn)行多相滲流計(jì)算；(2)將單元中心點(diǎn)的溫度和孔隙壓力通過插值的方法傳至節(jié)點(diǎn)；(3)在FLAC3D中設(shè)置集合收斂和偏差，進(jìn)行應(yīng)力應(yīng)變分析；(4)用得到的平均應(yīng)力修正孔隙度和滲透率。然后將得到的新的孔隙度和新的滲透率放在TOUGH2的輸入文件中進(jìn)行多相滲流計(jì)算，循環(huán)直至到達(dá)指定時(shí)間。

2 一維固結(jié)沉降模型解析解和數(shù)值解對比

為了驗(yàn)證 TOUGH2-FLAC3D程序的準(zhǔn)確性，本文選取存在解析解的Terzaghi一維固結(jié)沉降問題進(jìn)行驗(yàn)證。該模型描述了地面施加載荷后土體排水固結(jié)過程(見圖2)，究其本質(zhì)為滲流與力的耦合。固結(jié)模型的具體參數(shù)如表1所示。

圖1 一維固結(jié)沉降的彈簧活塞模型[14]

表1 一維固結(jié)沉降參數(shù)設(shè)置

2.1 太沙基單向固結(jié)微分方程及其解析解

太沙基單向固結(jié)微分方程可表示為如下形式:

(10)

其中,Cv稱為土的固結(jié)系數(shù)，cm2/s，其值為：

(11)

上述固結(jié)微分方程表示了超靜孔隙水壓力u與位置z及時(shí)間t的關(guān)系，然后根據(jù)土層滲流固結(jié)的起始條件與邊界條件，求出其解析解，當(dāng)附加應(yīng)力P沿土層均勻分布時(shí)孔隙水壓力u(z,t)的解答如下：

(12)

其中，m為奇正整數(shù)(1，3，5，……)；

(13)

其中，Es=κ(1+v)/[3(1-v)]為壓縮模量。

2.2 結(jié)果分析

圖2(a)為在t=25 s、409 s、1 638 s、6 553 s和26 214 s時(shí)，孔隙壓力隨深度變化的分布。

在外荷載P的作用下，隨著固結(jié)時(shí)間的增大，土中孔隙水逐漸排出，超靜孔隙水壓力逐步消散，土骨架的有效應(yīng)力逐步增大，直至超靜孔隙水壓力為0，水完全排出。在固結(jié)過程中，隨著孔隙水的排出，土體產(chǎn)生壓縮，使得土體更加密實(shí)，強(qiáng)度也逐漸增大。

從圖2可以看出，隨著時(shí)間的推移，孔隙水被排出，土體產(chǎn)生壓縮；頂部單元的豎向位移隨著孔隙水的排出逐漸增大，在孔隙水全部排出之后，豎向位移逐漸趨于穩(wěn)定。圖2(a)通過對數(shù)值解與解析解進(jìn)行對比，發(fā)現(xiàn)整體趨勢基本一致，誤差較小，這個(gè)誤差可能由模型網(wǎng)格稀疏以及滲流和力學(xué)計(jì)算過程中數(shù)據(jù)插值誤差引起的；圖2(b)為模型頂部單元的位移與太沙基一維固結(jié)問題的解析解對比分析，可以看出，TOUGH2和FLAC3D對一維固結(jié)問題模擬的數(shù)值模擬結(jié)果，與其解析解比較基本吻合，印證了上述力學(xué)水熱耦合程序的準(zhǔn)確性。

圖2 一維固結(jié)沉降解析解和TOUGH2-FLAC3D數(shù)值解對比結(jié)果

3 人工地?zé)醿拥哪M

3.1 模型建立

以EGS-D難以建立大規(guī)模穩(wěn)定蓄熱，水流量小，易造成污染等缺點(diǎn)，唐春安等人提出了一種新的EGS-E開發(fā)深部干熱巖儲層：該方法主要由豎井、上下水管道、換熱池、冷- 熱水流通道和爆破致裂區(qū)組成(如圖3所示)。其工作原理：通過開鑿豎井、換熱池、橫巷，形成地下空間。在橫巷中實(shí)施鑿巖爆破致裂，形成水-巖換熱接觸面積。地下空間施工完成后，一次性注水，并通過爆破致裂區(qū)換熱，形成高溫高壓地下熱湖。

圖3 基于開挖技術(shù)的增強(qiáng)型地?zé)嵯到y(tǒng)(EGS-E)[15]

以1 km3的花崗巖儲層地?zé)崮荛_發(fā)為例，以崩落區(qū)的十分之一為目標(biāo)儲層，此工程結(jié)構(gòu)包含25個(gè)分支“熱湖”構(gòu)件，則每個(gè)熱儲的容積為0.04 m3，若硐室崩落長1 m，則正方形截面邊長為200 m，模擬該平面問題。假設(shè)此200 m2的截面內(nèi)，絕大部分是崩落的巖體，孔隙度為2%。人工地?zé)醿幽Ｐ鸵圆豢蓾B透巖石的半無限半空間為界，整個(gè)初始溫度為300℃。在100℃溫度下的水以4 kg/s的恒定速率注入儲層的一側(cè)，而另一側(cè)則在特定的井筒壓力下進(jìn)行生產(chǎn)，供熱產(chǎn)能。在表2中給出了人工地?zé)醿幼⒉赡Ｐ蛦栴}的參數(shù)，并且低溫水由左側(cè)點(diǎn)I處注入，在右側(cè)點(diǎn)P處產(chǎn)能，模型的詳細(xì)信息如圖4所示。

圖4 人工地?zé)醿幼⒉赡Ｐ褪疽鈭D

表2 人工地?zé)醿訁?shù)設(shè)置

3.2 優(yōu)先流動(dòng)路徑分析

在經(jīng)過37個(gè)時(shí)間步后達(dá)到規(guī)定的最大時(shí)間1.577 88×108s，此時(shí)生產(chǎn)速度為4.991 4 kg/s，熱焓為872.461 kJ/kg； 生產(chǎn)元件的溫度為204.73℃。用軟件Tec plot 360對模擬結(jié)果進(jìn)行整理，如圖5所示為人工地?zé)醿幽Ｐ驮?.5 a、9 a、10 a、13 a、13.6 a和16 a的溫度分布云圖。

圖5 不同時(shí)刻人工地?zé)醿訙囟确植荚茍D(℃)

從圖5中我們可以看出，在模型左側(cè)用100 ℃的水以7.5 kg/s的速率注入后，且在右側(cè)以7.5 kg/s的速率抽取。在1.5 a時(shí)，在注入井周圍由于低溫水的注入，巖體與流體之間產(chǎn)生強(qiáng)烈的熱交換，此時(shí)生產(chǎn)井周圍單元的溫度急劇下降；在產(chǎn)能井周圍，由于液體被抽離，即溫度也會(huì)有一定程度的下降。隨著生產(chǎn)井模擬時(shí)間的增大，整個(gè)模型的溫度逐漸下降。注入模型的流體沿注入井和生產(chǎn)井之間的優(yōu)先流動(dòng)路徑運(yùn)移，在換熱效率足夠高的情況下，巖石與流體之間的換熱效果最好，即在這個(gè)連線上的溫度梯度變化最快。在9 a之后，模型注入井周圍的溫度都散落在注入流體的溫度附近，由于生產(chǎn)井中抽出流體的溫度在下降，因此，模型的產(chǎn)熱效率逐漸的下降。

由生產(chǎn)的流體溫度隨時(shí)間變化的關(guān)系曲線(見圖6)，可以看到隨著抽取時(shí)間的增大，生產(chǎn)井單元的溫度在逐漸下降，產(chǎn)生流體的溫度變化的斜率也在逐漸減小，最后趨于恒定，說明隨著時(shí)間的增大，模型中的熱量被逐漸抽離，且整個(gè)人工地?zé)醿赢a(chǎn)熱的效率在逐漸下降。

圖6 生產(chǎn)的流體溫度隨時(shí)間變化的關(guān)系

3.3 井位置變化對產(chǎn)能的影響

影響人工地?zé)醿娱_發(fā)的效率除了儲層的滲透率、孔隙度、注入井中流體的速率和溫度、以及產(chǎn)能井抽取流體的速率和井底壓力外，注入井和產(chǎn)能井的位置也會(huì)在很大程度上對其產(chǎn)生影響。為此，本小節(jié)通過對同一理想的人工地?zé)醿幽Ｐ?，僅改變注入井和生產(chǎn)井的位置來探究其對人工地?zé)醿赢a(chǎn)熱的效能的影響。通過改變注入井和生產(chǎn)井的位置，如圖7所示，對稱布置的①、②和③處，分別在左側(cè)設(shè)置注入井，在右側(cè)設(shè)置產(chǎn)能井。

圖7 井位置分布

圖8顯示了注入井和生產(chǎn)井分別布置于①、②和③處時(shí)，人工地?zé)醿釉?3.6 a時(shí)的溫度云圖。從圖可以看出，未設(shè)置井的兩角的區(qū)域面積(圖中標(biāo)注)，隨著注入井和生產(chǎn)井在相對位置距離的變近在逐漸的增大，即在注入的流體在沿優(yōu)先流動(dòng)路徑遷移時(shí)，掃掠過的面積在逐漸減??；隨著注入井和生產(chǎn)井在相對位置距離的變近，各模型中產(chǎn)能井單元溫度的變化在逐漸的增大；此外，在同樣時(shí)間刻度，按對角布置井的位置，儲層的熱抽取更快。結(jié)果表明，按對角布置井的位置可以將儲層的熱在最大程度上且最快的時(shí)間內(nèi)產(chǎn)出。

圖8 在13.6 a不同模型的溫度云圖(℃)

4 結(jié) 論

本文在TOUGH2軟件框架的基礎(chǔ)上，引入了FLAC3D的應(yīng)力應(yīng)變分析，結(jié)合TOUGH2共同建立通用的三維力學(xué)水熱耦合模型，并用解析解驗(yàn)證程序的耦合功能的準(zhǔn)確性，得出結(jié)論如下：

(1) 建立人工地?zé)醿幼⒉赡Ｐ?，?yàn)證了人工地?zé)醿又辛黧w由注入井到生產(chǎn)井這一“優(yōu)先流動(dòng)路徑”的快速遷移，與現(xiàn)場地?zé)崽镒C據(jù)相吻合[4-7]。同時(shí)也詳細(xì)監(jiān)測了在非等溫注入下的生產(chǎn)溫度隨時(shí)間的變化曲線。

(2) 通過改變注入井和生產(chǎn)井的位置，研究了人工地?zé)醿幼⒉赡Ｐ碗S井位置的變化所引起采熱效率的不同，得出井沿對角這一布置方式，注入流體在優(yōu)先流動(dòng)路徑下掃掠的面積最大，即儲層的生產(chǎn)效率高，開發(fā)也較為完全。