基于表皮效應(yīng)的單井循環(huán)地下?lián)Q熱系統(tǒng)地下水滲流方程的統(tǒng)一分析

宋 偉，王 浩，鄭昌金，劉遠(yuǎn)周

（北方工業(yè)大學(xué) 土木工程學(xué)院，北京100144）

0 前言

地下水源熱泵系統(tǒng)可以利用地?zé)豳Y源，其中單井循環(huán)地下?lián)Q熱系統(tǒng)是地下水源熱泵系統(tǒng)的一種形式。在取熱量相同的情況下，與地下埋管換熱器相比，單井循環(huán)地下?lián)Q熱系統(tǒng)所需的打井?dāng)?shù)量較少，初投資較低，從而得到廣泛應(yīng)用[1]，[2]。目前，單井循環(huán)地下?lián)Q熱系統(tǒng)主要有循環(huán)單井、抽灌同井和填礫抽灌同井[3]。循環(huán)單井是在基巖上打井，且為裸井，回水主要與井壁換熱，少部分回水流經(jīng)基巖含水層后直接進(jìn)入抽水區(qū)。抽灌同井和填礫抽灌同井的相同點(diǎn)：①均在井孔內(nèi)設(shè)置隔板，用于防止回水與抽水相互摻混，減少熱貫通效應(yīng)；②二者的回水均進(jìn)入含水層，與含水層中的原水和固體骨架換熱后進(jìn)入抽水區(qū)。不同點(diǎn)為填礫抽灌同井井壁周圍回填孔隙率較大的礫石，這層礫石區(qū)有利于回水進(jìn)入含水層[4]。

針對3種熱源井的地下水運(yùn)動(dòng)規(guī)律已有研究。倪龍建立了非穩(wěn)態(tài)、有越流的循環(huán)單井地下水滲流方程，武強(qiáng)在倪龍研究的基礎(chǔ)上，建立了含有蓄能顆粒與無儲(chǔ)能顆粒的循環(huán)單井地下水滲流數(shù)學(xué)模型[5]～[7]。李旻建立了無越流、非穩(wěn)態(tài)條件下抽灌同井的地下水滲流方程，武強(qiáng)在李旻研究的基礎(chǔ)上添加越流項(xiàng)。王玉林在武強(qiáng)研究的基礎(chǔ)上，分析了抽水與回水周期性變化和回灌率對地下水流動(dòng)的影響[8]～[11]。由于不同熱源井建立的地下水滲流方程不同，且邊界條件較為復(fù)雜，因而地下水滲流數(shù)學(xué)模型難以應(yīng)用于實(shí)際工程，因此，須要建立3種熱源井統(tǒng)一的理論模型。大多學(xué)者僅考慮滲透系數(shù)各向異性，未得到滲透系數(shù)各向異性對地下水流動(dòng)的影響。填礫抽灌同井的填礫區(qū)與表皮效應(yīng)理論中的表皮類似。目前，已有學(xué)者研究了表皮效應(yīng)對地下水流動(dòng)的影響，在實(shí)際生產(chǎn)中為了減小表皮效應(yīng)對開采地下水的影響，假設(shè)有效井孔半徑大于實(shí)際井孔半徑，并在此基礎(chǔ)上建立了數(shù)學(xué)模型[12]。王玉林發(fā)現(xiàn)井壁處的表皮效應(yīng)系數(shù)為非均勻分布，將含水層沿豎直方向分割成多個(gè)小含水層，每個(gè)小含水層的表皮效應(yīng)系數(shù)均勻分布，假設(shè)含水層的表皮效應(yīng)系數(shù)按某種規(guī)律分布，建立承壓含水層的穩(wěn)態(tài)滲流數(shù)學(xué)模型，通過模擬發(fā)現(xiàn)，表皮效應(yīng)系數(shù)越小，含水層的降深和抽水流量越大[13]，[14]。

本文利用表皮效應(yīng)統(tǒng)一3種熱源井近井壁處的滲透系數(shù)，得到滲透系數(shù)轉(zhuǎn)換關(guān)系式。本文建立填礫抽灌同井地下水滲流方程和邊界條件，獲得地下水滲流方程的解析解，考慮滲透系數(shù)各向異性對地下水流動(dòng)的影響；通過滲透系數(shù)轉(zhuǎn)換關(guān)系式，本文得到抽灌同井和循環(huán)單井地下水滲流方程的解析解。本文編制計(jì)算程序求解填礫抽灌同井地下水滲流方程的數(shù)值解。通過對填礫抽灌同井的地下水滲流方程的解析解與數(shù)值解進(jìn)行比較，驗(yàn)證了填礫抽灌同井的地下水滲流方程的數(shù)值解的合理性。

1 統(tǒng)一熱源井井壁處滲透系數(shù)

1.1 表皮效應(yīng)

在實(shí)際開發(fā)地下水資源或開采石油時(shí)，鉆井、完井和生產(chǎn)實(shí)踐活動(dòng)減小了井壁處的滲透系數(shù)。井壁處滲透系數(shù)改變的區(qū)域稱為“表皮”。圖1為表皮效應(yīng)示意圖。

圖1 表皮效應(yīng)示意圖Fig.1 Schematic diagram of skin effect

圖中曲線AB為降深曲線，在工程實(shí)踐中，井筒內(nèi)壓力點(diǎn)為E1，理論推導(dǎo)的壓力點(diǎn)為E，為解釋實(shí)際與理論的壓力差ΔP，本文引入表皮效應(yīng)系數(shù)。表皮效應(yīng)系數(shù)的計(jì)算式為[15]

式中：S為表皮效應(yīng)系數(shù)；K為含水層的滲透系數(shù)，m/s；Ks為表皮的滲透系數(shù)，m/s；rw為井壁半徑，m；rc為表皮層影響半徑，m。

當(dāng)表皮滲透系數(shù)小于含水層的滲透系數(shù)時(shí)，表皮效應(yīng)為正表皮效應(yīng)；當(dāng)表皮滲透系數(shù)大于含水層的滲透系數(shù)時(shí)，表皮效應(yīng)為負(fù)表皮效應(yīng)。

1.2 3種熱源井滲透系數(shù)

研究人員已對3種熱源井井壁處的滲透系數(shù)進(jìn)行了詳細(xì)地研究，并得到了3種熱源井井壁處的滲透系數(shù)數(shù)學(xué)模型，具體如下。

①抽灌同井

M K Hubbert在進(jìn)行達(dá)西實(shí)驗(yàn)的過程中發(fā)現(xiàn)，流體流過同一種多孔介質(zhì)時(shí)，抽灌同井條件下，井壁處的滲透系數(shù)與流體的運(yùn)動(dòng)粘滯系數(shù)成反比，與多孔介質(zhì)平均直徑的平方成正比，經(jīng)過數(shù)學(xué)推導(dǎo)得到抽灌同井滲透系數(shù)k1的表達(dá)式為[16]

式中：D1為抽灌同井條件下，井壁處顆粒的平均直徑，m；C為固體骨架結(jié)構(gòu)的比例常數(shù)，C與顆粒的大小、分布、球度、密實(shí)程度等有關(guān)；g為重力加速度，m/s2；v為液體的運(yùn)動(dòng)粘滯系數(shù)，m2/s。

②填礫抽灌同井

填礫抽灌同井的填礫區(qū)是由粒徑較大的礫石組成，填粒區(qū)是多孔介質(zhì)，多孔介質(zhì)的孔隙較多且孔隙的直徑較大，基于此假設(shè)多孔介質(zhì)中的孔隙為一束等長度的毛細(xì)管，利用N-S方程和Kozeny-Carman方程聯(lián)立求解，得到填礫抽灌同井條件下，井壁處滲透系數(shù)k2的表達(dá)式為[17]

式中：ρ為流體的密度，kg/m3；D2為礫石的平均直徑，m；μ為流體的動(dòng)力粘滯系數(shù)，N·s/m2；n為孔隙率。

③循環(huán)單井

陳崇希通過研究發(fā)現(xiàn)，水平井的回水區(qū)分為兩部分：一部分為圓形管道，另一部分為含水層。圓形管道的滲透系數(shù)與管道內(nèi)流體的流動(dòng)狀態(tài)有關(guān)，當(dāng)圓形管道內(nèi)流體的流動(dòng)狀態(tài)為層流時(shí)，圓形管道的滲透系數(shù)可以由文獻(xiàn)[18]進(jìn)行推導(dǎo)。

圓形管道的水頭損失Δh的計(jì)算式為

式中：f為管壁的摩擦系數(shù)；l為管道長度，m；D3為管道直徑，m；u為管道內(nèi)流體的平均流速，m/s。

f的計(jì)算式為

式中：Re為雷諾數(shù)。

Re的計(jì)算式為

v的計(jì)算式為

圓形管道的滲流速度V的計(jì)算式為

當(dāng)介質(zhì)為水時(shí)，設(shè)定孔隙率n為1，則V=u，并帶入式（4）～（7）化簡整理得到：

綜上可得，循環(huán)單井的滲透系數(shù)k3為

1.3 統(tǒng)一滲透系數(shù)

3種熱源井示意圖如圖2所示。

圖2 3種熱源井示意圖Fig.2 Schematic diagram of three thermalwells

由于3種熱源井的結(jié)構(gòu)不同，導(dǎo)致3種熱源井井壁處的滲透系數(shù)不同，本文將3種熱源井井壁周圍區(qū)域放大，得到3種熱源井表皮效應(yīng)，具體如圖3所示。

圖3 3種熱源井表皮效應(yīng)Fig.3 Skin effect of three thermalwells

假設(shè)抽灌同井、循環(huán)單井存在與填礫區(qū)類似的區(qū)域（回填區(qū)），抽灌同井的回填物與含水層的物理性質(zhì)相同，循環(huán)單井的回填物為水。根據(jù)滲流力學(xué)中的表皮效應(yīng)，將3種熱源井的回填區(qū)類比為表皮，由此對3種熱源井的滲透系數(shù)進(jìn)行統(tǒng)一。基于表皮效應(yīng)理論，將抽灌同井條件下井壁的滲透系數(shù)等效于填礫抽灌同井和循環(huán)單井條件下井壁的滲透系數(shù)。由于抽灌同井回填區(qū)材料與含水層的物理性質(zhì)相同，即：

將式（11）帶入式（1）中，得到S=0。由于填礫抽灌同井回填區(qū)材料為粒徑較大的礫石，因此，回填區(qū)的滲透系數(shù)大于含水層的滲透系數(shù)，即S＜0，rw＜rc。將S＜0，rw＜rc帶入式（1）得到k2/K＞1，將k2/K＞1帶入式（2），（3），（11）整理得到：

2 地下水滲流方程的求解

2.1 地下水滲流方程

填礫抽灌同井結(jié)構(gòu)如圖4所示。填礫抽灌同井具有井軸對稱的特點(diǎn)，填礫抽灌同井的井軸為縱軸，含水層的底部為橫軸，含水層上半部分為回水區(qū)，含水層下半部分為抽水區(qū)，回水經(jīng)回水濾網(wǎng)進(jìn)入回水區(qū)，含水層原水從抽水區(qū)進(jìn)入抽水濾網(wǎng)，抽水與回水在含水層同一徑向不同深度處發(fā)生。

圖4 填礫抽灌同井結(jié)構(gòu)圖Fig.4 Structure diagram of FECSCW

含水層應(yīng)滿足的假設(shè)條件、滲流方程及邊界條件如下。

①假設(shè)條件

假設(shè)含水層為無越流承壓含水層，當(dāng)含水層的壓力下降時(shí)，含水層立刻彈性釋水；含水層為均質(zhì)、各向異性的物質(zhì)，主滲透方向?yàn)樗椒较蚝拓Q直方向，含水層厚度始終不變，含水層的側(cè)向無限延伸；流體在含水層中的流動(dòng)符合達(dá)西定律；回灌和抽水流量相等，且抽、回水量在抽、回水濾網(wǎng)處均勻分布；當(dāng)系統(tǒng)穩(wěn)定運(yùn)行時(shí)，地下水流動(dòng)處于穩(wěn)態(tài)滲流。

②滲流方程及邊界條件

基于上述假設(shè)條件，得到承壓含水層中滲流微分方程為

式中：H為地下水的壓力，kPa；krr為含水層水平方向的滲透系數(shù)，m/s；kzz為含水層豎直方向的滲透系數(shù)，m/s。

邊界條件的表達(dá)式為

式中：r1為井中心與井壁之間的距離，m；r為井中心距含水層內(nèi)任一點(diǎn)的半徑，m；kh為填礫抽灌同井壁面處的滲透系數(shù)，m/s；Q1為回水流量，m3/s；Q2為抽水流量，m3/s；M為含水層的厚度，m；d1為抽水濾網(wǎng)下邊緣與含水層底部之間的距離，m；l1為抽水濾網(wǎng)長度，m；l2為回水濾網(wǎng)長度，m；d2為回水濾網(wǎng)下邊緣與抽水濾網(wǎng)上邊緣之間的距離，m；d3為回水濾網(wǎng)上邊緣與含水層頂部之間的距離，m。

2.2 解析解

邊界條件的表達(dá)式為

利用分離變數(shù)法對式（15）進(jìn)行處理，并將式（16）按余弦級數(shù)展開，得到填礫抽灌同井的降深方程。

式中：K0（nπr/ρM），K1（nπr1/ρM）為第二類虛宗量貝塞爾函數(shù)；w為填礫抽灌同井的結(jié)構(gòu)計(jì)算變量。

通過3種熱源井井壁處的滲透系數(shù)關(guān)系式，本文得到抽灌同井與循環(huán)單井的降深方程。

填礫抽灌同井的結(jié)構(gòu)計(jì)算變量的表達(dá)式為

2.3 數(shù)值解

圖4中填礫抽灌同井關(guān)于z軸對稱，因此只取填礫抽灌同井和含水層的一半進(jìn)行研究，并利用外節(jié)點(diǎn)法對井壁外含水層區(qū)域進(jìn)行網(wǎng)格劃分。為了與地下水滲流方程的解析解進(jìn)行比較，本文假設(shè)滲透系數(shù)各項(xiàng)同性，即krr=kzz，ρ=1，Z=z，通過有限差分法對式（15）進(jìn)行離散，得到網(wǎng)格中心點(diǎn)、回水濾網(wǎng)處、抽水濾網(wǎng)處、上頂板、下頂板處的地下水壓力計(jì)算式分別為

式中：i為網(wǎng)格的縱坐標(biāo)；i+1為i上方相鄰縱坐標(biāo)；i-1為i下方相鄰縱坐標(biāo)；j為網(wǎng)格的橫坐標(biāo)；j+1為j上方相鄰橫坐標(biāo)；j-1為j下方相鄰橫坐標(biāo)；H（i，j）為網(wǎng)格的地下水壓力，kPa；ri為網(wǎng)格的半徑，m；Δr為網(wǎng)格的寬度，m；Δz為網(wǎng)格的高度，m。

網(wǎng)格的計(jì)算程序步驟：①設(shè)置網(wǎng)格的初始值及迭代計(jì)算的收斂條件（2次計(jì)算結(jié)果的差值小于10-6）；②進(jìn)行網(wǎng)格的迭代計(jì)算，當(dāng)網(wǎng)格的數(shù)值滿足收斂條件后停止計(jì)算。

3 模型驗(yàn)證

引用參考文獻(xiàn)[22]中填礫抽灌同井的數(shù)據(jù)，表1為填礫抽灌同井參數(shù)取值。本文通過地下水滲流方程的解析解驗(yàn)證地下水滲流方程的數(shù)值解，具體步驟：①地下水滲流方程的解析解與數(shù)值解的地下水壓力分布趨勢是否相同；②地下水滲流方程的解析解選取3個(gè)不同位置的地下水壓力值，地下水滲流方程的數(shù)值解同樣選取相同位置的地下水壓力值，解析解的地下水壓力值與數(shù)值解的地下水壓力值之間的差值是否不超過5%。同時(shí)滿足上述兩個(gè)條件，則地下水滲流方程的數(shù)值解正確。

表1 填礫抽灌同井參數(shù)取值Table 1 Parameter selection of FECSCW

圖5為地下水滲流方程的數(shù)值解。當(dāng)滲透系數(shù)各向同性時(shí)，krr=kzz，ρ=1，Z=z；當(dāng)滲透系數(shù)各向異性時(shí)，krr/kzz=0.5，krr/kzz=2。由式（17），（18）得到地下水壓力和地下水流動(dòng)速度。地下水滲流方程的解析解如圖6所示。

圖5 地下水滲流方程的數(shù)值解Fig.5 Numerical solution of groundwater seepage equation

圖6 地下水滲流方程的解析解Fig.6 Analytic solution of groundwater seepage equation

由圖5，6可知，地下水壓力分布關(guān)于z=0.48反對稱，抽水區(qū)的地下水壓力小于含水層的初始地下水壓力時(shí)，含水層中的原水進(jìn)入抽水濾網(wǎng)；回水區(qū)的地下水壓力大于含水層的初始地下水壓力時(shí)，回水回灌到含水層。由圖5，6還可以看出，含水層上、下頂板的地下水壓力沿含水層豎直方向沒有變化，與含水層上、下頂板的邊界條件相對應(yīng)。在圖5中選取3個(gè)不同點(diǎn)的地下水壓力值，這3個(gè) 點(diǎn)的 坐 標(biāo) 分 別 為H（0.4 ，0.6 ），H（0.1 ，0.7 ），H（0.5 ，0.8 ），相應(yīng)的地下水壓力值分別為118.88 ，129.72 ，118.71 kPa；在圖6中選取3個(gè)不同點(diǎn)的地下水壓力值，這3個(gè)點(diǎn)坐標(biāo)分別為H（0.4 ，0.6 ），H（0.1 ，0.7 ），H（0.5 ，0.8 ），相應(yīng)的地下水壓力值分別為119.68 ，127.72 ，114.21 kPa。

3組誤差的平均值為2.0%，因此，利用地下水滲流方程的解析解可以驗(yàn)證數(shù)值解，并且地下水滲流方程的數(shù)值解滿足上文的2個(gè)條件，因此，地下水滲流方程的數(shù)值解正確。地下水滲流方程的解析解與數(shù)值解之間的偏差由數(shù)學(xué)化簡所致，數(shù)學(xué)簡化可分為兩步，第一步為假設(shè)地下水滲流方程中2個(gè)自變量互不影響，地下水滲流方程由偏微分方程化簡為微分方程；第二步為采用有限差分法將微分方程化簡為差分方程。提高數(shù)值解精度的方法為通過加密網(wǎng)格，使差分方程無限逼近微分方程。

含水層的地下水速度矢量圖如圖7所示。

圖7 含水層的地下水速度矢量圖Fig.7 Groundwater velocity vector of aquifer

由于圖6中抽水區(qū)與回水區(qū)之間存在地下水壓力差，迫使圖7中含水層的原水和回灌水從回水區(qū)流動(dòng)到抽水區(qū)，從而發(fā)生流貫通，因此，填礫抽灌同井回水濾網(wǎng)與抽水濾網(wǎng)之間應(yīng)設(shè)置擋板?；厮畢^(qū)z=0.70 處水流的速度方向?yàn)橄蛴?，z=0.45處水流的速度方向?yàn)橄蛳?，水流速度方向的改變，消耗了大量?dòng)能，導(dǎo)致z=0.45 處水流的速度很小。由圖7可知，抽、回水區(qū)水流的速度隨含水層半徑的增大而減小。

圖8為在同一深度下不同點(diǎn)處含水層的滲透系數(shù)之比對地下水壓力的影響

圖8 在同一深度下不同點(diǎn)處含水層的滲透系數(shù)之比對地下水壓力的影響Fig.8 The influence of the ratio of permeability coefficients of aquifers at different points at the same depth on groundwater pressure

由圖8可知：在z=0.69 處，地下水壓力隨著不同點(diǎn)處含水層滲透系數(shù)之比的增大而增大，在z=0.23 處，地下水壓力隨著不同點(diǎn)處含水層滲透系數(shù)之比的增大而減??；抽水區(qū)與回水區(qū)之間的地下水壓力差隨著不同點(diǎn)處含水層滲透系數(shù)之比的增大而增大。通過分析發(fā)現(xiàn)：地下水壓力隨著不同點(diǎn)處含水層滲透系數(shù)之比的增大而增大；含水層中地下水流動(dòng)方向與抽、回水口水流方向不一致會(huì)導(dǎo)致地下水消耗大量動(dòng)能；地下水動(dòng)能來源于抽、回水水泵的揚(yáng)程；對于滲透系數(shù)之比大于1的含水層應(yīng)選用揚(yáng)程大的水泵。

圖9為在同一半徑下不同點(diǎn)處含水層的滲透系數(shù)之比對地下水壓力的影響。

由圖9可知，當(dāng)krr/kzz=2時(shí)，r=0.07 ～0.23 m處，地下水壓力變化幅度劇烈；當(dāng)krr/kzz=0.5 時(shí)，r=0.07 ～0.23 m處，地下水壓力變化幅度平緩，抽、回水的影響范圍隨著含水層滲透系數(shù)之比的減小而減小。此外，圖9中地下水壓力曲線相交于一點(diǎn)，該點(diǎn)坐標(biāo)為z=0.48 m，這與抽、回水區(qū)地下水壓力分布關(guān)于z=0.48 m反對稱相對應(yīng)。

4 結(jié)論

本文統(tǒng)一3種熱源井壁面處的滲透系數(shù)，建立填礫抽灌同井地下水滲流方程，使用分離變數(shù)法、傅里葉變換和有限差分法，分別得到地下水滲流方程的解析解和數(shù)值解，并對地下水滲流方程的解析解和數(shù)值解進(jìn)行對比得到如下結(jié)論。

①通過表皮效應(yīng)建立循環(huán)單井、抽灌同井和填礫抽灌同井壁面處滲透系數(shù)的轉(zhuǎn)換關(guān)系式。

②以填礫抽灌同井為例，本文建立穩(wěn)態(tài)承壓含水層地下水滲流方程及邊界條件，獲得填礫抽灌同井的降深方程。通過滲透系數(shù)轉(zhuǎn)換關(guān)系式，本文得到抽灌同井和循環(huán)單井的降深方程。本文利用有限差分法對地下水滲流方程進(jìn)行離散，并編寫計(jì)算程序，求得填礫抽灌同井的地下水滲流方程數(shù)值解。計(jì)算結(jié)果表明，地下水滲流方程的數(shù)值解和解析解的平均誤差為2.0%，誤差來源于數(shù)學(xué)化簡，減少誤差的方法為加密網(wǎng)格，加密網(wǎng)格使得地下水滲流方程的數(shù)值解無限逼近解析解，填礫抽灌同井的地下水滲流方程可用于解決該熱源井影響下的穩(wěn)態(tài)無越流承壓含水層的地下水流動(dòng)問題。

③抽、回水區(qū)的壓力分布關(guān)于z=0.48 m反對稱。填礫抽灌同井出現(xiàn)流貫通現(xiàn)象。抽、回水的影響范圍隨著不同點(diǎn)處含水層滲透系數(shù)之比的增大而增大。滲透系數(shù)之比大于1的含水層應(yīng)選擇揚(yáng)程大的水泵；在填礫抽灌同井回水濾網(wǎng)與抽水濾網(wǎng)之間的含水層設(shè)置擋板后，能夠有效減緩流貫通的趨勢。