井內嵌入同軸多孔與固體圓柱地下水滲流流型分析

朱 琳，雷海燕，馬 非，戴傳山

（天津大學中低溫熱能高效利用教育部重點實驗室，天津 300350）

垂直于圓管外的繞流現象是自然界較普遍存在的流體力學現象之一，如卡門渦街。在多孔介質滲流層內垂直于圓管外的繞流問題也是石油鉆探、地下水開采[1]、地熱能開發、流化床等許多工業領域經常遇到的物理現象，如地熱井筒與熱儲層的耦合傳熱問題[2]，甚至在臨床醫學研究中也有相關的微血管外的間質滲流研究[3]。在飽和多孔介質內嵌入物體，被嵌物體外的繞流、傳熱傳質問題在很多工業應用領域中有廣泛的應用。針對多孔介質內嵌有簡單幾何形狀的物體結構，若能獲得其物體外部的繞流與熱質擴散的問題的簡單表達形式的理論解析解，對指導類似結構下的工業設計有重要的作用。

Pop等[4]首次獲得了考慮流體黏性耗散Brinkman模型圓柱外繞流的理論解析解。相對Pop等[4]研究的有滲流通過的多孔介質內嵌入圓柱的流場結構而言，多孔介質內開挖有一圓形孔洞，且孔洞內充滿純流體的空間結構下的流場理論解析解獲得的比較早。Drost等[5]采用忽略慣性力僅考慮流體黏性力的Stokes流模型獲得了無限大多孔介質內有一圓形孔洞的流場理論解析解，得出孔內平均流速與地層遠端流速相比（即平均收斂因子）約等于2。Sano[6]進一步分析了圓形孔洞的存在對周圍達西流的影響，得出在多孔介質滲透系數較低情況下，孔洞圓心處速度約是遠端地層滲流速度的3倍等重要結論。這一結論隨后得到了Momii等[7]在實驗室內基于激光多普勒方法進行的實驗驗證。Raja等[8?10]進一步推導了有均勻滲流通過的不同滲透系數下多孔介質層內有一規則或微變形的圓形孔洞或球形孔洞結構下的流場內流體速度（即收斂因子α）分布理論解。

Li等[11]針對一種相對前人研究更為復雜的有滲流通過的無限大多孔介質內嵌一固體圓柱，且圓柱外包裹有水環的幾何結構下的各區域流型結構進行了研究，其研究目的是獲取多孔介質內嵌入細長吸收碳纖維捕捉流體內膠體顆粒物的機理。Sekhar[12]將流體環內固體替換為多孔介質，得出該幾何結構下的流動問題理論解，分析了內、外多孔介質在不同滲流系數情況下，整體多孔介質表現出的宏觀平均水文特征、局部速度分布以及流動阻力特性。

前人研究結果表明，有均勻滲流通過的多孔介質內嵌入固體圓柱、多孔層、有水環或無水環等復雜幾何結構下流場的理論解析研究文獻并不多見，應用領域或研究目的卻比較廣泛，如復合材料封井、井內流速觀測[13]、控制污染物在地下水和土壤中的擴散[14]等。值得注意的是，前人最多只對3層結構，即無限大多孔介質-水環-環狀多孔層/圓形固體進行了研究。本文借鑒以往的研究方法，對有滲流通過的無限大多孔介質流場內嵌入固體圓柱，且圓柱外部順序包裹有多孔介質環、純流體環的復雜幾何結構下的流場速度解析解進行了理論分析，著重研究了不同滲透系數下的內外多孔介質、水環間隙情況等對圓柱外繞流流型以及縱橫界面上的速度分布的影響規律。由于地下水滲流速度較小，以往研究中常不考慮由彌散效應引起的黏度系數的變化，即多孔介質內黏度系數與純流體內的黏度系數視為相同[15]。此外，在井筒附近的地下水滲流可能會產生顆粒的沉積堵塞[16]，但本研究忽略這一影響，認為局部滲透系數不變。本文研究結果對類似結構下的地埋管套管換熱器、地下水污染物吸收裝置設計、地下水測速裝置研發等有一定的理論指導意義。

1 物理模型的建立與求解方法

考慮無限大2D均質多孔介質（Pe）內有單向滲流，垂直井筒內另嵌入同軸均質多孔介質（Pi）與固體圓柱（S），在內外多孔介質之間有開放水環（Wen）的情況，假設多孔介質均為各向同性材料：滲透系數為ke的無限大多孔介質內有黏性流體穩定通過，繞過井中一包裹有滲透系數為ki的多孔層且半徑為r0的圓形固體（圖1）。水環區域（λ1r0

圖1 物理模型示意圖Fig.1 Schematic diagram of the physical model.

水環內、外多孔介質區由Brinkman方程[18]和連續性方程控制：

式中：μeあ—多孔介質內流體有效黏度系數；

μen—水環內純流體黏度系數，μeあ=μen=μ；

en—水環；

l=i、e—內部、外無限大多孔介質區。

速度矢量（V）、橫縱坐標（x、y）、各點到固體圓柱圓心半徑（r）和壓力（P），進行無量綱化：

式中：U∞—無限遠處的流速大小。

水環區控制式（1）（2）無量綱化：

多孔區控制式（3）～（4）無量綱化：

式中：χ2—滲透性能，χ2=r20/k。

由速度分量和剪切應力（T）的連續性[19?20]，得到邊界條件（每個邊界條件都可以分解為切向和法向，切向剪切應力包含壓力P）：

引入極坐標系下的流函數（ψ）：

式中：θ—速度分量Vθ與橫坐標軸x夾角。

式（6）—（8）簡化為：

結合無窮遠處均勻流動條件ψe∞→?rsinθ，各區域式（13）（14）的解為：

式中：I1—第一類修正一階貝塞爾函數；

K1—第二類修正一階貝塞爾函數。

其中系數A、B、C、D、E、F、G、H、M、N可由式（9）—（11）確定。

值得注意的是，若r0→0，解的形式與Sekhar[12]一致：

此外，若r0→0,λ1→0，解的形式為：

笛卡爾坐標系下的流速U：

2 方法驗證

為了驗證計算方法的正確性，本文以Raja等[8]的計算模型為例，進行了程序正確性驗證，該模型中包括純流體孔洞部分的Stokes流模型以及外部多孔介質部分流動的Brinkman模型。輸入參數與文獻[8]一致，結果見圖2。

圖2 χe變化下縱向中心剖面（Y=0），橫向中心剖面（X=0）無量綱速度分布圖Fig.2 Dimensionless velocity distribution map of the longitudinal center profile（Y=0） and the lateral center profile（X=0） under the changes of χe

該計算模型對應于無限大多孔介質內有一孔洞的情況，即r0→0,λ1→0。其中，橫縱坐標無量綱化：χe變化下經過圓心的縱、橫剖面速度分布見圖2。值得注意的是，無量綱化特征長度為孔洞半徑r0c=λ2r0。隨著χe的增加（滲透系數變小），孔洞的存在感變強，中心流速與遠端流速比值αc變大到3，若χe為0（滲透系數極大），則αc→1。本文結果與Raja等[8]研究結果符合非常好。

3 多層結構計算結果

在方法驗證的基礎上，對圖1所示的多層結構模型進行理論求解。由于影響參數比較多，僅對不同半徑、以及內外滲透系數變化情況下的流型結構求解。計算工況見表1。值得注意的是，為了更清晰地表示外部尺寸固定（即λ1/λ2為定值）情況下，內部圓柱尺寸的變化對于流場的影響，此處給出的r0是一組以r0=r02為標準的假設值（流場的變化只與λ1、λ2、χi、χe有關）。同理，整個流場均如式（5）所述，計算過程中使用相應的r0進行參數的無量綱化為了更加直觀地看出不同參數的影響，在計算結果的整理中，橫縱中心剖面坐標軸的無量綱化均以r02半徑為準：

表1 不同工況的計算參數Table 1 Summary of calculation conditions under different parameters.

3.1 固體圓柱半徑大小r0對流型結構的影響

圖3（a）（b）可以看出，固體半徑r0的變化對水環外的多孔區幾乎沒有影響。但隨著r0增大，水環內流速峰值變高，這是由于流體經過包裹有多孔層的固體結構時流動阻力變大（圖4）。此外，r0變大，固體圓柱在多孔介質層中對流體的影響半徑增大，逐漸呈繞流趨勢，如圖5（a）（b）。更進一步地，水環內的多孔區縱向剖面的平均流速值變小，橫向剖面上的速度變大。這是因為該區域內部阻力變大，流入的流量越來越少，但繞流增加了橫向剖面上的流速。另外，在r0極大的情況下（工況3），多孔層內流速急速下降到0，這是由于極薄多孔層對流場沒有太大影響，如圖6（a）。

圖3 縱向、橫向中心剖面無量綱速度分布圖Fig.3 Dimensionless velocity distribution map of the longitudinal center profile and the lateral center profile

圖4 典型工況下的壓力云圖（工況2）Fig.4 Pressure contours under typical working condition 2

圖5 典型工況下的流線圖Fig.5 Streamline diagram under the typical working conditions

以工況2為對比工況，分別研究參數λ1、λ2、χi、χe的變化對流型結構的影響。此外，由于繞流，流體在水環內橫向剖面上的流速變化更為明顯，因此水環內 流型變化主要用該剖面上的流速峰值分析。

3.2 λ1變化對流型結構的影響

圖3（c）（d）為λ1變化下的中心剖面無量綱速度分布圖。λ1的變化對水環外的多孔區（下文簡稱外多孔區）幾乎沒有影響。但是，隨著λ1變小，水環內流速峰值先變高后變低，這是由于縫隙變大，進而剖面上速度分布展寬變大。即進入水環內多孔區（下文簡稱內多孔區）的流體變少，而水環內流量增加流速變大，直至流體幾乎不通過內多孔區（工況5），見圖6（b），且過程中隨著展寬變大，水環內流量分配空間變大導致其流速峰值逐漸降低。類似地，由于縫隙變大，內多孔區阻力相對水環阻力變大，該區域流量降低流速變小，因此流場也更容易平穩。同時，λ1的減小會使內多孔層中的繞流效果相對增強，因此該區域流型平穩程度要耦合區域內流量大小和圓柱外繞流強度兩種因素。

圖6 典型工況下的速度云圖Fig.6 Speed contours under the typical working conditions

3.3 λ2變化對流型結構的影響

圖3（e）（f）為λ2變化下的中心剖面無量綱速度分布圖。當λ2很大時，縫隙很大，極少流體流入內多孔區（工況6）。隨著λ2的增大，水環內流速峰值先升高后降低，這是由于縫隙變大流速分布展寬變大，流體越來越少地進入內多孔層且水環內流量分配空間變大導致的。由于內多孔區阻力相對水環阻力增大，該區域平均流速成倍降低且流型變化趨勢更平穩。λ2減小，意味著整個障礙物尺寸相對減小，因此，外部多孔介質區的流型會更加平緩。

3.4 χi變化對流型結構的影響

圖3（g）（h）為χi變化下的中心剖面無量綱速度分布圖。χi的變化對外多孔區的流場幾乎沒有影響。但是，當χi極大時，流體幾乎不進入內多孔區。隨著χi的減小，水環內流速峰值降低。這是因為內多孔區滲透性加強，內多孔層阻力相對水環阻力變小，內多孔區流量變大水環內流量減小。同理，內多孔區橫向和縱向剖面流速均相應變大，且變化明顯。當內多孔區滲透系數極大時（工況9），見圖6（c），水環與內多孔層之間的速度階躍消失，出現“穿透”現象，這是因為內多孔層阻力與水環阻力同量級。

3.5 χe變化對流型結構的影響

圖3（i）（j）為χe變化下的中心剖面無量綱速度分布圖。隨著χe變小水環峰值降低，內多孔區流速越低流型越平緩，但外多孔區流速增大。這是由于內外多孔區滲透系數差變小，整體結構流型過渡更為平緩，水環內的加速效應相應減弱，內多孔區的流速相應降低。χe≤10時外多孔區流速變化趨勢與χe較大時相反，這是由于與流體向水環內的滲透作用相比，水環外的繞流作用占了主導地位，見圖5（c），流體極少進入內多孔區。

4 結論

（1）水環內部參數如r0、λ1、χi的變化對外層多孔介質區影響很小；水環直徑的減小會使外部流型更為平穩；水環外部流型主要受控于外部滲透系數，但由于障礙物的存在，滲透系數對外部流型的影響是非單調的函數。

（2）橫向剖面上水環內的速度變化較大，固體圓柱r0減小，水環內峰值增加。內部與外部滲透系數同時增大，水環內加速效應減弱。水環間隙對水環內流速峰值影響較大。

（3）研究發現內部滲透系數增加到某一臨界值情況下會出現“穿透”現象，橫向速度剖面從階梯形變為拋物線形。