張太平,劉 洋,丁若偉,凡倩瑩
(1.山東省地質科學研究院,山東 濟南 250013;2.自然資源部金礦成礦過程與資源利用重點實驗室,山東 濟南 250013;3.山東省金屬礦產成礦地質過程與資源利用重點實驗室,山東 濟南 250013;4.中國地質大學(武漢)環境學院,湖北 武漢 430078)
地熱能是一種清潔、可持續和可再生的能源,由于其具有儲量巨大、產量穩定且基本不受季節、氣候影響等優勢,其開發利用與存儲問題越來越受到人們的廣泛關注。實際上,更多的熱能貯存于地下的巖石中,也被稱為“干熱巖型地熱資源”。在干熱巖型地熱資源的開發利用過程中,裂隙-基巖滲流傳熱問題的基礎研究是認識和揭示相關地熱能開發利用的熱學理論基礎,也是解決能源與環境領域諸多問題的關鍵環節之一,具有較大的工程應用價值。注抽試驗廣泛用于地熱能的開發利用,即先將冷水注入干熱巖中獲取熱量,然后將熱水抽出來使用,有時完整的干熱巖的儲水能力有限,必須采用水壓致裂方法提高地熱能的開采效率。通常注抽試驗分為多井注抽試驗和單井注抽(Single-Well Push-Pull,SWPP)試驗,其中多井注抽試驗至少需要兩口井,開展多井注抽試驗時,首先將冷水通過一口或者多口井注入到裂隙中,然后通過觀測井將熱水抽出;而單井注抽試驗只需要一口井,其操作是首先將冷水通過一口井注入裂隙介質中,經過一段時間后再將熱水從該井中抽出。相比之下,單井注抽試驗具有成本低,耗時短和操作簡單等優點。
到目前為止,國內外學者已經建立了大量的模型用于研究熱在裂隙-基巖含水層介質中物質和能量的遷移轉化機理,如Neretnieks、Tang等、Mahmoudzadeh等、Zhu等、Zhou等、Zhu等、Trinchero等。在裂隙-基巖含水層系統中的SWPP熱質運移解析解方面,Kocabas采用雙重拉普拉斯變換法構建了兩階段的SWPP熱質運移試驗模型;Jung等采用拉普拉斯變換法和傅立葉變換法構建了三階段的SWPP熱質運移試驗解析解;Wang等考慮混合效應物理過程,采用拉普拉斯變換法和格林函數法構建了三階段的SWPP熱質運移試驗解析解。混合效應是指注入的冷水與井筒中的熱水混合的過程。現有的SWPP試驗模型包含很多的假定條件,實際很難滿足,例如現有的模型忽略了裂隙-基巖含水層介質中的熱彌散和橫向熱擴散的影響,導致現有的模型難以推廣。在數值解方面,目前有很多商業的數值模擬軟件(FEFLOW,GMS,COMSOL,TOUGH2等)可以構建SWPP熱質運移試驗模型,如魏永霞等采用GMS構建了韶關溫泉三維地熱數值模型,開展了地熱資源量評價;劉東林等采用TOUGH2軟件開展了東麗湖地區霧迷山組地熱資源可開采潛力評價研究;單丹丹等采用COMSOL軟件研究了地熱深井生產過程中井筒熱能損失。然而,Wang等指出現有的商業模擬軟件在刻畫混合效應過程時假定井筒的水位等于含水層的厚度,因此現有的商業模擬軟件難以準確地開展地熱資源量評價、地熱相關參數的敏感性分析和參數反演方面的研究。
為了更好地理解裂隙-基巖含水層系統中的熱質運移機理,本文采用對流-擴散方程來構建裂隙-基巖含水層系統中SWPP熱質運移試驗的數學模型,但由于數學模型考慮了裂隙和基巖中的橫向和縱向熱傳導和熱擴散以及井筒中的混合效應等多個過程,很難導出該數學模型的解析解,因此本文采用有限差分法構建SWPP熱質運移試驗的數值模型來研究這些過程對試驗結果的影響,數值模型采用MATLAB軟件中的ODE15s求解器求解。
x
方向,垂直裂隙方向為z
方向。SWPP熱質運移試驗模型的控制方程中包括熱對流、熱擴散和熱傳導過程,SWPP熱質運移試驗模型的控制方程為
圖1 SWPP試驗概念模型示意圖Fig.1 Schematic diagram of the SWPP test conceptual model

z
≤b
,t
>0)(1a)

(1b)


(2a)

(2b)

(2c)
式中:k
和k
分別為裂隙和基巖的熱傳導率[W/(L·℃)];β
和β
分別為裂隙的縱向和橫向熱彌散度。在許多前人的研究中,裂隙和基巖中的初始溫度分布都假定為常數。實際上,在很多含水層中,溫度的分布在空間上并不是均勻的。類似于Chen等對于初始溫度的處理方式,本研究假設裂縫和基巖中的初始溫度分布是一個任意函數:
T
(x
,z
,t
)|=0=T
0(x
,z
)(3a)
T
(x
,z
,t
)|=0=T
0(x
,z
)(3b)
式中:T
0(x
,z
)和T
0(x
,z
)為分布在裂隙和基巖上的溫度的任意函數,在交界面上T
0(x
,z
=b
)=T
0(x
,z
=b
)。與Wang等的工作類似,本研究的數學模型的內邊界條件考慮混合效應,在SWPP熱質運移試驗的注入階段,考慮混合效應的內邊界條件為


(4a)
T
(x
,t
)|=0=T
0(r
)(4b)

在SWPP熱質運移試驗的抽取階段,考慮混合效應的內邊界條件為

(5)

其外邊界條件為

(6)
裂隙-基巖交界面的溫度和溫度通量都是連續的,因此交界面的邊界條件為
T
(x
,z
,t
)|==T
(x
,z
,t
)|=(7a)

(7b)
式中:φ
為基巖的孔隙度。x
,x
]離散為N
-1個區域,即N
個節點,這里采用x
表示x
方向上外邊界的長度,則x
方向上的空間離散為
(8)
其中,x
+12的計算公式如下:
(9)
裂隙介質區域[-b
,0]和基巖區域[0,Z
]在z
方向上分別離散成M
-1(M
個節點)個子區域和M
-1(M
個節點)個子區域,這里采用Z
表示z
方向上外邊界的長度,則z
方向上的空間離散為
(10a)

(10b)
其中,z
+12和z
+12的計算公式分別如下:
(11a)

(11b)
因此,本研究構建的SWPP熱質運移試驗數學模型的空間離散為


i
=0,1,…,N
;j
=0,1,…,M
)(12a)


(12b)

(12c)

t
>t
)(12d)

(12e)
T
,(x
,z
,t
)|==T
(x
,z
,t
)|=(12f)

(12g)
上述方程(12a)~(12g)即為考慮混合效應、橫向熱傳導和熱彌散作用的二維SWPP熱質運移試驗數值模型,數值模型采用MATLAB軟件中的ODE15s求解器求解。
D
=D
=0和β
=0以及t
=0時,本研究構建的數值模型將還原為Wang等的模型。本研究所選取的參數源于Wang等和Jung等研究中,參數設置為:k
=2.1 W/(m·℃),k
=2.1 W/(m·℃),ρ
=2 650 kg/m,c
=1 000 J/(kg·℃),ρ
=1 000 kg/m,c
=4 200 J/(kg·℃),φ
=0.5,φ
=0.5,b
=0.1 m,B
=5 m,Q
=-Q
=1.44 m/h,r
=0.5 m,t
=250 min,t
=400 min,h
,inj=50 m,h
,res=47.5 m,h
,ext=45 m,D
=D
=0,β
=β
=0,T
=1℃,初始條件設置為0。圖2為本研究構建的數值解在特定條件下與Wang等的數值解的對比。
圖2 兩模型在井壁上溫度穿透曲線的對比Fig.2 Comparison of BTCs at the wellbore computed by two models
由圖2可見,修改后兩個模型在井壁上的溫度穿透曲線能夠很好地擬合,從而驗證了本研究構建的數值模型的正確性。
D
值下井壁上的溫度穿透曲線,其計算結果見圖3。其中,D
=0表示不考慮橫向熱擴散作用,其他參數值的設置同圖2。
圖3 不同有效熱擴散系數Dfz值下井壁上溫度穿透 曲線的對比Fig.3 Comparison of BTCs at the wellbore for the different values of Dfz
由圖3可見,隨著D
值的增加,注入階段井壁上的溫度降低,但是在抽取階段井壁上的溫度升高,拖尾現象嚴重。以上分析表明:裂隙介質中的橫向熱擴散作用會導致溫度穿透曲線發生明顯的變化,說明橫向熱擴散作用不容忽視。b
)的影響,通常來說,當裂隙介質的寬度較小時,橫向熱擴散作用較小。因此,為了研究裂隙寬度對SWPP熱質運移試驗結果的影響,采用本研究建立的數值模型,計算了t
=240 min時不同裂隙寬度下在x
=0.5 m位置處裂隙-基巖中的溫度分布曲線,其計算結果見圖4。其中,D
=0.063 1 m/d,其他參數設置同圖2。
圖4 不同裂隙寬度b下在x=0.5 m處裂隙-基巖中的 溫度分布曲線Fig.4 Distribution curves of temperature at x=0.5 m in fracture-matrix for the different values of fracture width (b)
由圖4可見,在橫向擴散系數一定的條件下,裂隙寬度為1.00 m裂隙和基巖介質中的溫度值最大。但值得注意的是,當z
<b
時,裂隙介質中的溫度在橫向熱擴散作用下緩慢降低。以上分析表明:裂隙介質的寬度會對SWPP熱質運移試驗造成明顯的影響,裂隙介質的寬度越大,裂隙和基巖介質中溫度值越高。x
方向上的熱交換量,其計算結果見圖5。其中,t
=0.5 d和t
=1.0 d表示注入階段裂隙-基巖交界面上熱交換量的分布,t
=20 d和t
=25 d表示抽取階段裂隙-基巖交界面上熱交換量的分布。參數設置如下:k
=2.1 W/(m·℃),k
=2.1 W/(m·℃),ρ
=2 650 kg/m,c
=1 800 J/(kg·℃),c
=1 000 J/(kg·℃),ρ
=1 000 kg/m,c
=4 200 J/(kg·℃),φ
=0.5,b
=0.1 m,B
=5 m,Q
=-Q
=1.44 m/h,r
=0.5 m,t
=1 d,t
=35 d,h
,inj=50 m,h
,res=47.5 m,h
,ext=45 m。
圖5 不同時間下裂隙-基巖交界面在x方向上的熱交換量Fig.5 Distribution of the heat flux across the interface of fracture-matrix along the x-direction at different times
由圖5可見:在注入階段,t
=0.5 d和t
=1.0 d時裂隙-基巖交界面上的熱交換量為負值,這說明在注入階段熱量由裂隙介質往基巖介質中運移,但值得注意的是,t
=0.5 d時裂隙-基巖交界面上的熱交換量大于t
=1.0 d時的熱交換量,這主要是因為在注入階段的初期,裂隙-基巖中的溫度梯度較大,隨著注入時間的增加,裂隙-基巖中的溫度梯度逐漸減小,進而造成裂隙-基巖交界面上的熱交換量減小;相反地,在抽取階段,裂隙-基巖交界面上的熱交換量為正值,這說明在抽取階段熱量由基巖介質往裂隙介質中運移,隨著抽取時間的增加,裂隙-基巖交界面上的熱交換量增加,這主要是因為隨著抽取時間的增加,裂隙中的熱量逐漸被抽出,造成裂隙-基巖中的溫度梯度逐漸增加,從而造成基巖介質中的熱量往裂隙介質中運移。本文采用考慮橫向熱擴散作用下的對流擴散方程建立了SWPP熱質運移試驗數學模型,數學模型的內邊界條件考慮混合效應,數學模型的求解采用有限差分法構建數值模型并用MATLAB軟件中的ODE15s求解器求解,通過與現有的解析解模型進行對比,驗證了數值模型的穩定性和精度,并分析了橫向熱擴散作用對SWPP熱質運移試驗結果的影響,得到以下結論:
(1) 在SWPP熱質運移試驗中,橫向熱擴散作用對熱質運移的影響不能忽視,裂隙介質中的橫向熱擴散作用會導致注入階段井壁上的溫度穿透曲線值降低,而抽取階段井壁上的溫度穿透曲線值則升高,拖尾現象明顯。
(2) 裂隙介質的寬度對SWPP熱質運移試驗結果的影響明顯,裂隙寬度越大,裂隙和基巖介質中溫度值越高。
(3) 在SWPP熱質運移試驗的注入階段,隨著注入時間的增加,裂隙介質中的熱量往基巖介質中的運移量逐漸降低,在抽取階段,隨著抽取時間的增加,基巖介質中的熱量往裂隙介質中的運移量逐漸增加。