非均勻多孔介質(zhì)等效滲透率的普適表達式

張 東，劉曉麗，王恩志

(清華大學(xué)水沙科學(xué)與水利水電工程國家重點實驗室，北京 100084)

多孔介質(zhì)滲流問題由來已久，與生活實際和工業(yè)生產(chǎn)都密切相關(guān)。Henry Darcy在研究均勻砂柱問題時，提出了著名的達西定律，奠定了滲流學(xué)科的基礎(chǔ)。在此基礎(chǔ)之上，學(xué)者們針對地下水滲流問題相繼提出了地下水運動基本微分方程、Dupuit井穩(wěn)定滲流公式等，促進了滲流學(xué)科的迅速發(fā)展[1-2]。直到目前，多孔介質(zhì)滲流問題仍然是重要而廣泛的研究領(lǐng)域，包括孔隙模型和模擬、數(shù)值模擬方法改進、多孔介質(zhì)的非均勻性、非線性滲流、多相流等，整體而言，更趨于微觀尺度的非線性滲流性質(zhì)和滲流過程中的化學(xué)-物理機理研究[3-6]。微觀尺度考慮多孔介質(zhì)內(nèi)部的細觀結(jié)構(gòu)，認為介質(zhì)是非連續(xù)性的，存在主要的滲流和儲水空間——孔隙和吼道，可以形成明顯的滲流通道，計算時建立非連續(xù)介質(zhì)模型，如毛細管束模型、孔隙網(wǎng)絡(luò)模型等，采用LBM、逾滲理論等方法，根據(jù)研究問題可能還會考慮毛細作用、浸潤作用等，主要是巖石滲流問題，如油氣開采及地?zé)釂栴}等[7-11]；而宏觀尺度則忽略多孔介質(zhì)的微觀結(jié)構(gòu)，認為大于一定尺度(多孔介質(zhì)的表征單元體)，其微觀結(jié)構(gòu)的滲透率趨于穩(wěn)定，可以用滲透率來反映其滲透性能，認為多孔介質(zhì)為連續(xù)介質(zhì)，運動方程基于達西定律，多用于工程問題如地下工程、環(huán)境工程等。微觀尺度的關(guān)注并不代表宏觀大尺度計算不受重視，反而正如Thomas所總結(jié)的：“微觀尺度的物理特性決定宏觀尺度的模型”[12]，正是微觀尺度的機理性研究，為宏觀模型提供了改進的基礎(chǔ)，一系列理論和方法如體積平均方法、混合理論等建立并運用[13]。

縱觀微觀和宏觀尺度研究，多孔介質(zhì)的滲透性都是重要的研究指標(biāo)，對于宏觀實際工程問題，滲透率是數(shù)值計算的基礎(chǔ)，滲透率的準(zhǔn)確性直接決定模擬結(jié)果的精確度；而對于微觀尺度研究，滲透率是評價介質(zhì)特性的重要指標(biāo)，也是聯(lián)系宏觀問題的紐帶[14-16]。滲流問題的多尺度和尺度差異性是一個重大挑戰(zhàn)，尤其是工程大尺度模型所帶來的巨大計算量，提高計算效率勢在必行。

達西定律是線性滲流定律，其適用范圍有限，一般認為砂性土地下水滲流，在低雷諾數(shù)(Re<10)時服從達西定律[17-19]。但大部分實際工程問題仍然以達西定律作為基本假設(shè)進行處理簡化[20-22]，而目前的商用軟件如MODFLOW等許多模塊也是基于達西定律而建立的[23-24]。復(fù)雜的水文地質(zhì)條件決定了多孔介質(zhì)滲透性質(zhì)在空間的不均勻分布特性，J Koestel測量了某土壤層孔隙結(jié)構(gòu)，發(fā)現(xiàn)由于壓實和生物作用，1 m內(nèi)不同深度的土壤的滲透性差異很大[25]；水文地質(zhì)勘探常采用測井法確定地層的滲透性，但試驗范圍局限于百米，只是整個工程地質(zhì)模型的局部[26]。因此，研究不均勻介質(zhì)局部滲透率分布與整體等效滲透率的關(guān)系還是很有必要的。

本文基于連續(xù)介質(zhì)的假定，考慮介質(zhì)的不均勻性，推導(dǎo)多孔介質(zhì)滲透率空間分布與總體等效滲透率的關(guān)系，并進行數(shù)值模擬計算的比對。

1 數(shù)學(xué)模型及驗證

1.1 數(shù)學(xué)模型

本文的滲流模型可以概述為圖1的物理模型，本研究僅以水平x方向的滲流特性為代表，左右邊界分別為入口、出口，且都為壓力邊界，即典型的達西壓降滲流問題。

壓降法是達西滲流試驗方法，試驗時以水頭表示，其本質(zhì)即水頭壓力差驅(qū)動滲流發(fā)生，其模型示意圖如圖1所示。

圖1 壓降滲流試驗示意圖Fig.1 Diagram of the percolation with pressure drop

數(shù)值試驗的目的在于得到多孔連續(xù)介質(zhì)的等效滲透率。數(shù)值試驗定義多孔介質(zhì)滲透率分布函數(shù)κ(x,y,z)是空間位置函數(shù)，在不同空間區(qū)域滲透率及其變化趨勢不同，即建立了非均勻多孔介質(zhì)模型，基于達西定律探討了滲透率空間分布與多孔介質(zhì)總體等效滲透率的關(guān)系，并與COMSOL數(shù)值模擬的計算值進行了比對分析。

根據(jù)壓降試驗條件建立相應(yīng)的數(shù)學(xué)模型(圖2)。模型建立所采用的假定為：

(1)多孔介質(zhì)為連續(xù)介質(zhì),Lx、Ly、Lz為滲流域的x、y、z方向的長度；

(2)該滲流服從線性達西滲流定律；

(3)假設(shè)為飽和滲流，忽略重力作用；

(4)邊界位置的速度很小，忽略邊界的速度水頭，故本研究僅僅考慮壓力水頭；

圖2 壓降滲流數(shù)學(xué)模型示意圖Fig.2 Diagram of the percolation mathematic model

根據(jù)上述假設(shè)，在多孔介質(zhì)域內(nèi)滲流滿足達西定律[27-28]，即：

v=-K

(1)

式中：K——滲透系數(shù)/(m·s-1)，由介質(zhì)和流體性質(zhì)共同決定；

κ——滲透率/m2，僅與多孔介質(zhì)的幾何結(jié)構(gòu)有關(guān)，本文中采用滲透率作為研究指標(biāo)；

v——斷面平均流速/(m·s-1);

μ——流體動力粘滯系數(shù)/(Pa·s)；

γ——流體重度/(N·m-3)。

(2)

(3)

基于連續(xù)介質(zhì)假定，則在滲流域內(nèi)任意點(x0,y0,z0)及其鄰域可定義滲透率κ(x0,y0,z0)，即滲流域內(nèi)可定義滲透率分布函數(shù)κ(x,y,z)，滲流域數(shù)學(xué)模型可表達如下：

滲流域(0≤x≤Lx,0≤y≤Ly,0≤z≤Lz)：

(4)

滲流邊界：

1)x=0：p=P_in

(5)

2)x=Lx：p=Po

(6)

(7)

式中：P_in——入口壓力/Pa；

Po——出口壓力/Pa；

n——邊界的外法向量；

ρ——流體密度/(kg·m-3)。

1.2 數(shù)值模擬驗證

對式(4)～(7)數(shù)值求解即可對滲流流域進行模擬，本文采用有限元COMSOL的達西滲流模塊，其控制方程和式(2)一致。按照圖2所示，設(shè)置相應(yīng)的邊界條件，采用自由四面體網(wǎng)格劃分方法和MUMPS算法進行求解。采用不同單元尺寸進行網(wǎng)格剖分，得到一致的計算結(jié)果，表明計算結(jié)果不隨網(wǎng)格尺寸變化。且對單次結(jié)果進行流量平衡分析，入口和出口的流量差值小于10-9m3，可認為計算是準(zhǔn)確可信的。

2 一維滲透率分布

考慮多孔連續(xù)介質(zhì)域內(nèi)κ(x,y,z)僅僅與x位置有關(guān)，而在y和z方向不發(fā)生改變，可退化為一維滲透率分布函數(shù)κ(x)。

將式(2)微分化，則在多孔連續(xù)介質(zhì)域內(nèi)任意一點(x0,y0,z0)及其鄰域滿足：

(8)

式中：vx0——滲流域內(nèi)任意一點(x0,y0,z0)的x方向滲流速率。

結(jié)合一維滲流域達西定律并考慮流動的連續(xù)性，則有：

(9)

根據(jù)式(8)～(9)可以得到多孔連續(xù)介質(zhì)一維滲透率分布與等效滲透率的關(guān)系表達式：

(10)

由不同滲透率的介質(zhì)串聯(lián)拼接形成的孔隙介質(zhì)稱為復(fù)合多孔介質(zhì)。由2種不同滲透率的介質(zhì)串聯(lián)形成的復(fù)合介質(zhì)以數(shù)學(xué)語言表達，即：0≤x≤αLx區(qū)域的介質(zhì)滲透率為κ1，在αLx

(11)

采用COMSOL有限元法進行滲流數(shù)值模擬計算，計算幾何模型如圖2所示，其中Lx=10 m、Ly=5 m、Lz=5 m、P_in=1 000 Pa，Po=0。此外設(shè)定不同區(qū)域的滲透率參數(shù)：κ1=4×10-12m2,κ2=8×10-12m2(下同)，以對比驗證基于式(10)所得到的推論式(11)。

考慮兩種不同滲透率的介質(zhì)串聯(lián)形成的復(fù)合多孔介質(zhì)，改變?nèi)≈挡⒎謩e利用COMSOL和式(11)計算得到等效滲透率，計算結(jié)果如表1所示，相應(yīng)的曲線如圖3所示。

表1 滲透率分段分布函數(shù)計算參數(shù)及計算結(jié)果

圖3 數(shù)值模擬計算值與理論值對比曲線Fig.3 Numerical and analytical solutions

對于n(n>2)種或以上串聯(lián)所形成的復(fù)合介質(zhì)亦可根據(jù)式(10)求解，當(dāng)n趨近無窮時，κ(x)即為連續(xù)函數(shù)。本文中設(shè)置κ(x)為不同形式的連續(xù)函數(shù)，分別利用COMSOL和式(10)計算得到等效滲透率結(jié)果如表2所示(其中部分滲透率分布函數(shù)難以直接得到數(shù)學(xué)理論解，則利用MATLAB進行數(shù)值積分求解(下同)。一維滲透率分布函數(shù)條件下等效滲透率的理論值和模擬計算值的誤差小于2%，考慮數(shù)值計算的誤差，可以認為兩者的結(jié)果相當(dāng)吻合，理論推導(dǎo)式(10)是準(zhǔn)確的。

表2 滲透率連續(xù)分布函數(shù)計算參數(shù)及計算結(jié)果

3 二維滲透率分布

考慮多孔連續(xù)介質(zhì)域內(nèi)κ(x,y,z)僅僅與x和y位置有關(guān)，而在z方向不發(fā)生改變，可退化為二維滲透率分布函數(shù)κ(x,y)。則在點(x0,y0,z0)的x方向滲透速率滿足：

(12)

則在x=x0(?x0∈[0,Lx])斷面平均流量為：

(13)

對多孔介質(zhì)滲透域整體使用達西定律，則有：

(14)

根據(jù)流動連續(xù)性，則得到二維滲透率分布的等效滲透率滿足：

(15)

顯然式(15)是基于x=x0的一系列方程，但壓力分布函數(shù)p=p(x,y)與κ(x,y)密切相關(guān)，難以事先確定，因此式(15)難以直接求得理論解。因此，此處利用上一節(jié)的推導(dǎo)結(jié)論來簡化推導(dǎo)二維滲流率分布函數(shù)與等效滲透率的關(guān)系。

(16)

對該層使用達西定律，則得到該層的斷面平均流量：

(17)

則對于整個多孔連續(xù)介質(zhì)區(qū)域的斷面平均流量，即各分層滲流平均流量的總和，有：

(18)

根據(jù)流動連續(xù)性則可得到二維滲透率分布函數(shù)與等效滲透率的關(guān)系式：

(19)

對于不同滲透率的介質(zhì)并聯(lián)拼接形成的介質(zhì)可以用數(shù)學(xué)語言表達，即0≤x≤αLx區(qū)域的介質(zhì)滲透率為κ1，在αLx

(20)

考慮兩種不同滲透特性的多孔介質(zhì)并聯(lián)模型，改變α的取值并分別利用COMSOL滲流模擬和式(20)計算得到等效滲透率，計算結(jié)果如表3所示，相應(yīng)曲線見圖4。

表3 滲透率二元分段分布函數(shù)計算參數(shù)及計算結(jié)果

圖4 數(shù)值模擬計算值與理論值對比曲線Fig.4 Numerical and analytical solutions

對于n(n>2)種或以上并聯(lián)所形成的復(fù)合介質(zhì)亦可根據(jù)式(19)求解，當(dāng)n趨近無窮時即為連續(xù)函數(shù)。本文中設(shè)置為不同形式的二元連續(xù)函數(shù)，分別利用COMSOL和式(19)計算得到等效滲透率結(jié)果(表4)。

二維滲透率分布函數(shù)條件下等效滲透率的理論值和模擬計算值的誤差小于5%，考慮數(shù)值計算的誤差，可以認為兩者的結(jié)果相當(dāng)吻合，理論推導(dǎo)式(19)是準(zhǔn)確的。

表4 滲透率二元連續(xù)分布函數(shù)計算參數(shù)及計算結(jié)果

4 三維滲透率分布

考慮多孔連續(xù)介質(zhì)域內(nèi)κ(x,y,z)在x,y,z方向都具有差異性，則在點(x0,y0,z0)的x方向滲透速率滿足：

(21)

則在x=x0(?x0∈[0,Lx])斷面平均流量為：

(22)

對多孔介質(zhì)整體使用達西定律，則有：

(23)

根據(jù)流動連續(xù)性則得到三維滲透率分布的等效滲透率滿足：

(24)

同上，式(24)難以直接求得理論解。因此，此處利用上一節(jié)的推導(dǎo)結(jié)論來簡化推導(dǎo)三維滲流率分布函數(shù)與等效滲透率的關(guān)系。

(25)

對該層使用達西定律，則得到該層的斷面平均流量為：

(26)

則對于整個多孔連續(xù)介質(zhì)區(qū)域的斷面平均流量即各分層滲流平均流量的總和為：

(27)

根據(jù)流動連續(xù)性則可得到三維滲透率分布函數(shù)與等效滲透率的關(guān)系式:

(28)

構(gòu)建三維非均勻多孔介質(zhì)，其滲透率分布函數(shù)如式(29)所示。

(29)

根據(jù)式(28)，可以得到理論解：

(1-α)(1-β)κ4

(30)

根據(jù)建立的非均勻多孔介質(zhì)模型，改變α,β的取值并分別利用COMSOL和式(30)計算得到等效滲透率，計算結(jié)果如圖5所示，理論值和模擬計算值幾乎完全重合，誤差極小。

圖5 數(shù)值模擬計算值與理論值對比曲線Fig.5 Numerical and analytical solutions

本文中設(shè)置κ(x,y,z)為不同形式的三元連續(xù)函數(shù)，分別利用COMSOL和式(28)計算得到等效滲透率結(jié)果(表5)。

三維滲透率分布函數(shù)條件下等效滲透率的理論值和模擬計算值的誤差小于5%，考慮數(shù)值計算的誤差，可以認為兩者的結(jié)果相當(dāng)吻合，理論推導(dǎo)公式(28)是準(zhǔn)確的。

表5 滲透率三元連續(xù)分布函數(shù)計算參數(shù)及計算結(jié)果

5 結(jié)論

(1)針對一維、二維、三維滲流模型，基于連續(xù)介質(zhì)假設(shè)和達西定律，通過理論推導(dǎo)得到了等效滲透率的理論表達式；

(2)針對不同區(qū)域簡單常數(shù)κ(x,y,z)分布，即不同滲透特性的多孔介質(zhì)的串并組合模型，發(fā)現(xiàn)理論表達式和數(shù)值解誤差小于5%，充分表明理論推導(dǎo)的準(zhǔn)確性；

(3)進一步設(shè)定了非線性κ(x,y,z)分布，對比分析了理論值和數(shù)值模擬解，相對誤差也很小，證明了理論表達式的準(zhǔn)確性，但隨著非線性的增加，誤差逐漸增大，這是由于數(shù)值計算是通過線性求解所導(dǎo)致的。

本文所構(gòu)建的非均勻多孔介質(zhì)等效滲透率理論表達式，為簡化復(fù)雜地質(zhì)模型，減少計算量，為快速計算模擬提供了理論支持，還有利于模擬計算結(jié)果的合理性快速評估。此外，本文所設(shè)定的滲透率分布函數(shù)都比較簡單，針對實際情況構(gòu)建復(fù)雜的滲透率分布函數(shù)以及其對于等效滲透率的影響還有待更進一步的研究。