趙 軍，劉泉聲，張程遠

（1. 安徽理工大學 土木建筑學院，安徽 淮南 232001；2. 中國科學院武漢巖土力學研究所 巖土力學與工程國家重點實驗室，武漢 430071）

1 研究背景

地下水源熱泵技術是一種采集淺層低溫地能，同時滿足供暖和制冷的需求，并且實現(xiàn)零污染排放的能源利用方式。2005年該技術被中國建設部列為建筑業(yè)十項新技術之一，在建筑物中的推廣應用是國家列為節(jié)約資源節(jié)約工作重點之一，同時，許多地方都把發(fā)展地下水源熱泵作為發(fā)展本地經(jīng)濟的一個契機，對“節(jié)能減排”和“兩型社會”建設具有建設性的意義。在地下水源豐富的地區(qū)，如大江大河流域，地下水源熱泵是可以普遍采用的一種地源熱泵形式。地下水源熱泵利用了地下水，必須涉及取用水的回灌，不回灌可以避免地下熱積累對熱泵系統(tǒng)效能的影響，但是，只取水不進行有效回灌或回灌不慎，都會造成地面沉降和已有地下管線的破壞。當前水資源稅的開征，表明國家正著力改變長期以來實行資源和環(huán)境無價制度導致資源和環(huán)境惡化的現(xiàn)狀，通過征收資源與環(huán)境稅，體現(xiàn)資源和環(huán)境的價值。

關于單井系統(tǒng)（見圖 1），國內外不少專家和學者[1－8]有相當多的貢獻。在滲流場研究方面，李旻等[9]對單井回灌的含水層滲流場給出了解析解，為數(shù)值模擬分析奠定了基礎。何滿潮等[10－12]針對地熱水，對井回灌滲流場中的滲透系數(shù)進行了研究，得出受地下水溫度的影響，井體周圍產(chǎn)生的物理堵塞，使得井體周圍的滲透系數(shù)減小，增加了回灌的難度，并得出 K=K0e-λt（K0為初始孔隙系數(shù)；t為時間變量；λ為常數(shù)）這一衰減方程。雖然國內的不少學者在地源熱泵井的回灌問題上做了大量的工作，但關于地下水源熱泵系統(tǒng)熱-溫度-力全耦合研究文獻鮮見，僅有少量以將熱泵井群與含水層作為一個整體系統(tǒng)研究其整體儲熱性能這樣的理念為基礎的研究。

圖1 單井回灌系統(tǒng)簡示意圖Fig.1 Stetch of standing well system

本文從能量守恒定理、質量守恒定理等原理出發(fā)，分別建立熱傳遞方程、滲流方程、滲流方程、連續(xù)性方程以及 THM 耦合本構方程，將水源熱泵單井和地下水環(huán)境作為整個系統(tǒng)來進行多場的耦合分析研究，利用數(shù)值軟件對 THM 耦合模型給出定性的評價方法。

基本假定：①巖（砂）體是均質的，各向同性材料；②地下水流動服從Brinkman方程，Darcy定律；③巖體的變形為小變形，連續(xù)。

2 控制方程的建立

2.1 熱傳遞方程的建立

單井熱泵系統(tǒng)是利用地下水和地表溫差提取能量。在夏季，由于地下水水溫低于地表水，可以利用它來降溫，而冬季則恰好相反。被利用過的冷（熱）水往往都是要回灌到地下，它與原有地下水一般來說存在著溫度差異，將引起地下水形成溫度峰面，勢必對再次抽取地下水將產(chǎn)生影響。

引出參數(shù)有效熱傳導率keff用以定性的描述熱傳遞過程：

式中：n為巖體的孔隙率；kl、ks分別為流體和固體的熱傳導率?？刂企w的能量平衡方程為

式中：ρl為流體的密度；Cl為流體的比熱容；Cs為固體的比熱容；T為溫度； vx、vy、vz分別為流體在x、y、z方向的流體速度。

以M和N來分別替代等式（2）左端系數(shù)，即

式中：Q為匯源項。

利用坐標矢量關系后可進一步簡化為

在單井系統(tǒng)中，往往更多的時候是關心徑向的熱傳遞，由于考慮了巖體的各向同性性，式（5）可進一步簡化為

式中：ρs為固體骨架的密度；V為流體的速度。

式（6）就是單井系統(tǒng)中熱傳遞的能量表達式。

2.2 流動方程

2.2.1 區(qū)域流動方程

在利用水源熱泵井的過程中，抽取的地下水的流動速度的快慢對利用能量的過程起到關鍵的作用，接近井體壁處，流體速度較大，而遠離井體處流體的速度相對較小，因此，對井壁處和遠離井壁處分別建立不同的流動方程加以描述：

式中：?為拉普拉斯算子；η為黏度系數(shù)（kg/(ms)）；u為流體對固體顆粒的相對速度矢量（m/s）；ui,j、uj,i分別為I、j向的位移分量；k為滲透率（m2）；p為水壓力（kPa）。

2.2.2 連續(xù)性方程

地下水運動的連續(xù)性方程可以從質量守恒原理出發(fā)來考慮，即滲流場中水在某一單元體內的增減速率等于進出該單元流量速率之差，分別得到流體和固體的質量守恒方程[13－14]為

式中：s為飽和度；ρf為流體密度；vf為流體速度，其他符號意義同前。

2.3 應力本構方程的建立

2.3.1 有效應力原理

巖體中飽和-非飽和巖體的有效應力原理[8]可表示為式中：εp為孔隙水壓力引起的應變；εw為巖體裂隙吸水后膨脹引起的應變；βs為固相線熱膨脹系數(shù)；βw為吸水線膨脹系數(shù)；p為孔隙水壓力和空氣壓力（kPa）；C為巖體的彈性剛度張量；m為法向應力單位矩陣；Ks為固相的壓縮系數(shù)。

2.3.2 幾何方程

根據(jù)假定③可得幾何方程和體積應變：

式中：u為巖體骨架的位移量；εv為體應變。

2.3.3 巖/砂體中力學平衡方程

單井回灌系統(tǒng)中，由于冷（熱）水的抽取和回灌都將對巖體的內部應力（ Fi）狀態(tài)產(chǎn)生影響，而對其主要貢獻的是水頭壓力的變化和溫度的變化。

將式（10）中第二、三式代入式（12）可得

因孔隙水壓力是時間的函數(shù)，對其求偏微分可得

對式（14）取對時間的偏微分，忽略飽和度對時間變化的影響：

將式（14）代入式（15），可得

式中：h為水頭高度；γ為流體重度；c為常變量；pk為孔隙水壓力盒空氣壓力（kPa）。

對式（17）求關于時間t的偏微分：

根據(jù)假定③可得

式中：β1為熱膨脹系數(shù)。

將式（17）代入式（19），可得

再將式（20）代入式（16），經(jīng)過組合整理后，可得

式（21）中只含有溫度T、水頭壓力h和位移分量u以及內力Fi的本構方程。

為便于說明數(shù)值模擬的需求，下面給出單井系統(tǒng)巖土體破壞準則，考慮到實際工程，特給出三維狀態(tài)下的庫侖破壞準則：

式（23）中N和Q可分別用下式來表達：

式中：so、φ分別為黏聚力和內摩擦角。

通過式（23）可以看出，破壞準則是通過 fail來表達的，具體的判斷準則可表述為：當 fail = 0時，單井周圍的巖體開始破壞；當fail＞0時，單井周圍巖體處于穩(wěn)定的狀態(tài)；當fail＜0時，單井周圍巖體處于破壞狀態(tài)。

3 工程實例

以武漢某小區(qū)作為工程為例，施工前對場地的地質狀況進行了解，特別注意是否有地下管線及其準確位置。對地面進行清理，鏟除地面雜草、雜物和浮土，平整地面，確定鉆孔的具體位置，圖2、3分別為現(xiàn)場鉆井圖和測試井管道連接示意圖。整理現(xiàn)場監(jiān)測到的數(shù)據(jù)后，得到地下水溫與時間的變化關系，圖4為理論值與實測值的關系。從圖4可以得出，理論值和實測值的誤差控制在10%內，說明文中所提出的傳熱模型對預測水源熱泵井的井周溫度場的變化是有效的。

圖2 鉆孔Fig.2 Drill hole

圖3 測試管道連接示意圖Fig.3 Sketch of test tube

圖4 理論值與實測值的比較Fig.4 Comparison between theoretical values and actual measurement values

4 數(shù)值模擬

數(shù)值模擬前，先對地質資料圖進行適當?shù)暮喕?，?jīng)過簡化后的地址柱狀圖依次為①雜填土層，厚6 m；②黏土層，厚24 m；③細砂層，厚45 m；④粗砂礫層，厚25 m；地下穩(wěn)定水位位于6 m水頭處。建立在第四系地層上，呈砂、砂礫石層的交互疊置，為了便于模型的模擬，取圖5的模型，厚4 m，長5 m，寬1 m，為完全飽和砂土，各向同性彈性材料，井體貫穿地層100 m，井體半徑為0.25 m，材料的力學參數(shù)見表1。

圖5 井體模型（單位：m）Fig.5 The model of well (unit: m)

表1 井體材料參數(shù)Table1 Parameters of well

經(jīng)過參數(shù)取值后，首先得到井體的水壓力梯度同軸向距離的關系圖，如圖6所示。從圖中可以看出，當井體半徑r ＞ 3時，曲線較平緩；當0.5 ＜ r ＜1.5時，曲線將發(fā)生較大的變化，此時可以用Darcy-Brinkman方程來描述；當0 ＜ r ＜ 0.5時，曲線呈“陡峭”型，可以用N-S方程來描述；水壓力隨著井體半徑的增加而增加，直至達到一個穩(wěn)定的值。因此，本文所提到的流動耦合方程是適合本模型的。

圖6 水壓力與井體半徑關系圖Fig.6 Relationship between r and P

在實際工程中必須要考慮到井體成形后井體的整體變形圖，為了便于說明問題，得到了單井的整體位移圖，如圖7所示，圖中x、y、z分別為模擬整體位移圖，單位為m。

圖7 單井整體位移及變形圖Fig.7 The whole displacement and deformation diagram of single-well

模型在計算過程中只考慮了半對稱結構，從圖7可以看出，井體的最大位移發(fā)生在井體的邊緣處，最大變形量達到3.309×10-3m，在井壁處變形量接近為0。

現(xiàn)在對井體的破壞準則進行討論，利用式（23）作為理論依據(jù)，對fail取不同的值后得到結果如圖8所示，為無量綱。從圖中可以看出，考慮到模型的對稱性，因此只得出井體的兩個角處產(chǎn)生了熱位移，隨著 fail值越來越小，破壞影響范圍將越大，施工前對巖體的性質了解對井體是否失效起到關鍵作用。

圖8 井體破壞準則模型圖（單位：m）Fig.8 The model of the failure criterion for well (unit: m)

5 結 論

（1）通過井體的破壞準則可以得出，井體隨著破壞值的變化而變化，在實際工程中應充分了解巖體的性質后才能施工。

（2）單井系統(tǒng)中水壓力隨著井體半徑的增加而增加，直至達到一個穩(wěn)定的值。

（3）現(xiàn)場溫度實測值和理論值具有較好的吻合性，本文的模型對預測水源熱泵井溫度的變化是有效的。

[1]MILLARD A,REJIB A,CHIJIMATSU M,et al.Numerrical study of the THM effects on the near-field safety of a hypothetical nuclear waste repository—BMT1 of the DECOVALEX Ⅲ project. Part 2: Effects of THM coupling in continuous and homogeneous rocks[J]. Int. J.Rock Mech. Min. Sci.,2005,42: 731－744.

[2]MASSEI N,LACROIX M,WANG H Q,et al. Transport of particulate material and dissolved tracer in a highly permeable porous medium:comparision of the transfer parameters[J]. Journal of Contaminant Hydrology,2002,57: 21－39.

[3]STEPHEN E S. Particle transport through twodimensional,saturated porous media:influence of physical structure of the medium[J]. Journal of Hydrology,1995,167: 79－98.

[4]STEPHEN E S. The importance of the third dimension on transport through saturated porous media:case study based on transport of particles[J]. Journal of Hydrology,1996,179: 181－195.

[5]AHFIR N D,BENAMAR A,ALEMA,et al. Influence of internal structure and medium length on transport and deposition of suspended particles:a laboratary study[J].Transport Porous Media,2009,76: 589－307.

[6]BARTELDS G A,BRUINING J.,AND MOLENAAR J.The modeling of velocity enhancement in polymer flooding[J]. Transport Porous Media,1997,26: 75－88.

[7]NORIO T,KASUMI Y,GEORGE Z. Model study of the thermal storage system[J]. Geothermics,2003,32: 603－607.

[8]MASSEI N,LACROIX M,WANG H Q,et al. Transport of particulate material and dissolved tracer in a highly permeable porous medium: Comparison of the transfer parameters[J]. Journal of Contaminat Hydrology,2002,57: 21－39.

[9]李旻,刁乃仁,方肇洪. 單井回灌地源熱泵承壓含水層滲流解析解[J]. 山東建筑工程學院學報,2006,21(1): 1－5.LI Min,DIAO Nai-ren,FANG Zhao-hong. Analytical solution of seepage flow in a confined aquifer with a standing column well[J]. Journal of Shandong Institute of Architecture and Engineering,2006,21(1): 1－5.

[10]何滿潮,李啟民. 地熱資源在移民小區(qū)可持續(xù)發(fā)展應用研究[J]. 太陽能學報,2005,25(2): 223－226.HE Man-chao,LI Qi-ming. Reserch of sustainable development of geothermal resources in immigrant communities[J]. Acta Energiae Solaris Sinica,2005,25(2): 223－226

[11]何滿潮,劉斌,姚磊華,等. 地熱水對井回灌滲流場理論研究[J]. 中國礦業(yè)大學學報,2004,33(3): 245－248.HE Man-chao,LIU Bin,YAO Lei-hua,et al. Study on theory of seepage field a round geothemal productionreinfection doublets wells[J]. Journal of China University of Mining &Technology,2004,33(3): 245－248.

[12]何滿潮,劉斌,姚磊華,等. 地下熱水回灌過程中滲透系數(shù)研究[J]. 吉林大學學報(自然科學版),2002,32(4):374－377.HE Man-chao,LIU Bin,YAO Lei-hua,et al. Study on hydraulic conductivity during geothermal reinjection[J].Journal of Jilin University(Earth Science Edition),2002,32(4): 374－377.

[13]孔祥言,李道倫,徐獻芝,等. 熱-流-固耦合滲流的數(shù)學模型研究[J]. 水動力學研究與進展,2005,20(2): 270－275.KONG Xiang-yan,LI Dao-lun,XU Xiang-zhi,et al.Study on the mathematical models of coupled therma hydrological mechanical process[J]. Journal of Hydrodynamics,2005,20(2): 270－275.

[14]張強林,王媛. 巖體 THM 耦合模型控制方程建立[J].西安石油大學學報[J]. 2007,22(2): 139－145.ZHANG Qiang-lin,WANG Yuan. Establishment of THM coupling governing equations in rock masses[J]. Journal of Xi′an Shiyou University,2007,22(3): 139－145.