基于表皮效應的單井循環地下換熱系統地下水滲流方程的統一分析

宋 偉，王 浩，鄭昌金，劉遠周

（北方工業大學 土木工程學院，北京100144）

0 前言

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

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

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

1 統一熱源井井壁處滲透系數

1.1 表皮效應

在實際開發地下水資源或開采石油時，鉆井、完井和生產實踐活動減小了井壁處的滲透系數。井壁處滲透系數改變的區域稱為“表皮”。圖1為表皮效應示意圖。

圖1 表皮效應示意圖Fig.1 Schematic diagram of skin effect

圖中曲線AB為降深曲線，在工程實踐中，井筒內壓力點為E1，理論推導的壓力點為E，為解釋實際與理論的壓力差ΔP，本文引入表皮效應系數。表皮效應系數的計算式為[15]

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

當表皮滲透系數小于含水層的滲透系數時，表皮效應為正表皮效應；當表皮滲透系數大于含水層的滲透系數時，表皮效應為負表皮效應。

1.2 3種熱源井滲透系數

研究人員已對3種熱源井井壁處的滲透系數進行了詳細地研究，并得到了3種熱源井井壁處的滲透系數數學模型，具體如下。

①抽灌同井

M K Hubbert在進行達西實驗的過程中發現，流體流過同一種多孔介質時，抽灌同井條件下，井壁處的滲透系數與流體的運動粘滯系數成反比，與多孔介質平均直徑的平方成正比，經過數學推導得到抽灌同井滲透系數k1的表達式為[16]

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

②填礫抽灌同井

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

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

③循環單井

陳崇希通過研究發現，水平井的回水區分為兩部分：一部分為圓形管道，另一部分為含水層。圓形管道的滲透系數與管道內流體的流動狀態有關，當圓形管道內流體的流動狀態為層流時，圓形管道的滲透系數可以由文獻[18]進行推導。

圓形管道的水頭損失Δh的計算式為

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

f的計算式為

式中：Re為雷諾數。

Re的計算式為

v的計算式為

圓形管道的滲流速度V的計算式為

當介質為水時，設定孔隙率n為1，則V=u，并帶入式（4）～（7）化簡整理得到：

綜上可得，循環單井的滲透系數k3為

1.3 統一滲透系數

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

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

由于3種熱源井的結構不同，導致3種熱源井井壁處的滲透系數不同，本文將3種熱源井井壁周圍區域放大，得到3種熱源井表皮效應，具體如圖3所示。

圖3 3種熱源井表皮效應Fig.3 Skin effect of three thermalwells

假設抽灌同井、循環單井存在與填礫區類似的區域（回填區），抽灌同井的回填物與含水層的物理性質相同，循環單井的回填物為水。根據滲流力學中的表皮效應，將3種熱源井的回填區類比為表皮，由此對3種熱源井的滲透系數進行統一。基于表皮效應理論，將抽灌同井條件下井壁的滲透系數等效于填礫抽灌同井和循環單井條件下井壁的滲透系數。由于抽灌同井回填區材料與含水層的物理性質相同，即：

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

2 地下水滲流方程的求解

2.1 地下水滲流方程

填礫抽灌同井結構如圖4所示。填礫抽灌同井具有井軸對稱的特點，填礫抽灌同井的井軸為縱軸，含水層的底部為橫軸，含水層上半部分為回水區，含水層下半部分為抽水區，回水經回水濾網進入回水區，含水層原水從抽水區進入抽水濾網，抽水與回水在含水層同一徑向不同深度處發生。

圖4 填礫抽灌同井結構圖Fig.4 Structure diagram of FECSCW

含水層應滿足的假設條件、滲流方程及邊界條件如下。

①假設條件

假設含水層為無越流承壓含水層，當含水層的壓力下降時，含水層立刻彈性釋水；含水層為均質、各向異性的物質，主滲透方向為水平方向和豎直方向，含水層厚度始終不變，含水層的側向無限延伸；流體在含水層中的流動符合達西定律；回灌和抽水流量相等，且抽、回水量在抽、回水濾網處均勻分布；當系統穩定運行時，地下水流動處于穩態滲流。

②滲流方程及邊界條件

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

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

邊界條件的表達式為

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

2.2 解析解

邊界條件的表達式為

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

式中：K0（nπr/ρM），K1（nπr1/ρM）為第二類虛宗量貝塞爾函數；w為填礫抽灌同井的結構計算變量。

通過3種熱源井井壁處的滲透系數關系式，本文得到抽灌同井與循環單井的降深方程。

填礫抽灌同井的結構計算變量的表達式為

2.3 數值解

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

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

網格的計算程序步驟：①設置網格的初始值及迭代計算的收斂條件（2次計算結果的差值小于10-6）；②進行網格的迭代計算，當網格的數值滿足收斂條件后停止計算。

3 模型驗證

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

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

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

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

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

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

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

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

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

由于圖6中抽水區與回水區之間存在地下水壓力差，迫使圖7中含水層的原水和回灌水從回水區流動到抽水區，從而發生流貫通，因此，填礫抽灌同井回水濾網與抽水濾網之間應設置擋板。回水區z=0.70 處水流的速度方向為向右，z=0.45處水流的速度方向為向下，水流速度方向的改變，消耗了大量動能，導致z=0.45 處水流的速度很小。由圖7可知，抽、回水區水流的速度隨含水層半徑的增大而減小。

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

圖8 在同一深度下不同點處含水層的滲透系數之比對地下水壓力的影響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 處，地下水壓力隨著不同點處含水層滲透系數之比的增大而增大，在z=0.23 處，地下水壓力隨著不同點處含水層滲透系數之比的增大而減小；抽水區與回水區之間的地下水壓力差隨著不同點處含水層滲透系數之比的增大而增大。通過分析發現：地下水壓力隨著不同點處含水層滲透系數之比的增大而增大；含水層中地下水流動方向與抽、回水口水流方向不一致會導致地下水消耗大量動能；地下水動能來源于抽、回水水泵的揚程；對于滲透系數之比大于1的含水層應選用揚程大的水泵。

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

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

4 結論

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

①通過表皮效應建立循環單井、抽灌同井和填礫抽灌同井壁面處滲透系數的轉換關系式。

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

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