U型地源熱泵系統(tǒng)熱-流-力學(xué)耦合模型研究

趙 軍，劉泉聲，張程遠(yuǎn)

（1. 安徽理工大學(xué) 土木建筑學(xué)院，安徽 淮南 232001；2. 中國(guó)科學(xué)院武漢巖土力學(xué)研究所 巖土力學(xué)與工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430071）

1 研究背景

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

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

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

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

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

2 控制方程的建立

2.1 熱傳遞方程的建立

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

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

式中：n為巖體的孔隙率；kl、ks分別為流體和固體的熱傳導(dǎo)率。控制體的能量平衡方程為

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

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

式中：Q為匯源項(xiàng)。

利用坐標(biāo)矢量關(guān)系后可進(jìn)一步簡(jiǎn)化為

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

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

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

2.2 流動(dòng)方程

2.2.1 區(qū)域流動(dòng)方程

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

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

2.2.2 連續(xù)性方程

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

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

2.3 應(yīng)力本構(gòu)方程的建立

2.3.1 有效應(yīng)力原理

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

2.3.2 幾何方程

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

式中：u為巖體骨架的位移量；εv為體應(yīng)變。

2.3.3 巖/砂體中力學(xué)平衡方程

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

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

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

對(duì)式（14）取對(duì)時(shí)間的偏微分，忽略飽和度對(duì)時(shí)間變化的影響：

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

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

對(duì)式（17）求關(guān)于時(shí)間t的偏微分：

根據(jù)假定③可得

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

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

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

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

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

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

式中：so、φ分別為黏聚力和內(nèi)摩擦角。

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

3 工程實(shí)例

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

圖2 鉆孔Fig.2 Drill hole

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

圖4 理論值與實(shí)測(cè)值的比較Fig.4 Comparison between theoretical values and actual measurement values

4 數(shù)值模擬

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

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

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

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

圖6 水壓力與井體半徑關(guān)系圖Fig.6 Relationship between r and P

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

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

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

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

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

5 結(jié) 論

（1）通過井體的破壞準(zhǔn)則可以得出，井體隨著破壞值的變化而變化，在實(shí)際工程中應(yīng)充分了解巖體的性質(zhì)后才能施工。

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

（3）現(xiàn)場(chǎng)溫度實(shí)測(cè)值和理論值具有較好的吻合性，本文的模型對(duì)預(yù)測(cè)水源熱泵井溫度的變化是有效的。

[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]. 山東建筑工程學(xué)院學(xué)報(bào),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]何滿潮,李啟民. 地?zé)豳Y源在移民小區(qū)可持續(xù)發(fā)展應(yīng)用研究[J]. 太陽(yáng)能學(xué)報(bào),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]何滿潮,劉斌,姚磊華,等. 地?zé)崴畬?duì)井回灌滲流場(chǎng)理論研究[J]. 中國(guó)礦業(yè)大學(xué)學(xué)報(bào),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]. 吉林大學(xué)學(xué)報(bào)(自然科學(xué)版),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]孔祥言,李道倫,徐獻(xiàn)芝,等. 熱-流-固耦合滲流的數(shù)學(xué)模型研究[J]. 水動(dòng)力學(xué)研究與進(jìn)展,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]張強(qiáng)林,王媛. 巖體 THM 耦合模型控制方程建立[J].西安石油大學(xué)學(xué)報(bào)[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.