流體慣性效應和裂隙幾何屬性對非達西系數的影響

汪慧香, 馬 雷, 邢 坤, 張昊明

(合肥工業大學 資源與環境工程學院,安徽 合肥 230009)

裂隙介質中流體流動的非達西系數β的確定,對研究裂隙水污染運移、放射性廢物地下儲存和石油工程中的流體非達西流動具有重要意義[1-4]。對于流體在地質多孔介質中流動的β的影響因素和定量計算方法,相關研究成果較豐富[5-7]。文獻[8]通過將介質和流體的性質與Forchheimer的系數聯系起來,改進了用于描述多孔介質中非達西流動的Forchheimer方程,同時,首次將β作為湍流因子引入Forchheimer方程中。改進的Forchheimer方程為:

(1)

其中:P為壓力;x為流動路徑的長度;ρ為流體密度;g為流體重力加速度;J為總壓降;Jv為黏滯力導致的壓降;Ji為慣性力導致的壓降;a1為黏性項系數;b1為慣性項系數;v為流體平均速度;μ為流體的動態黏度;k為介質的滲透率;β為流體流動時所受到的慣性阻力,也稱為非達西系數或Forchheimer系數。

然而,文獻[4,9]已經證明裂隙中的非達西流是由足夠大的流體慣性效應引起的。參數β被廣泛接受為“非達西系數”或“慣性阻力系數”,β作為直接決定Forchheimer方程中非線性系數大小的參數,對于量化和表征多孔介質中非達西流動具有重要意義。

基于實驗室實驗,文獻[10-11]研究表明,多孔介質中的β僅取決于多孔介質的幾何屬性,與流體性質無關。文獻[12]將多孔介質中β與滲透率之間的關系統一歸納為負冪律關系,即

(2)

其中:A為比例因子;b為無量綱冪指數(b>0);ko為僅取決于介質幾何屬性的固有滲透率。

相關研究結果表明,多孔介質中β的經驗表達式通常與多孔介質的滲透率和孔隙度有關[13]。

近年來,含有ko和β的改進的Forchheimer方程被廣泛應用于單裂隙的研究中[14-15]。為得到裂隙介質中β的量化模型,需要研究β與裂隙幾何屬性之間的關系。文獻[16]研究得到圍壓條件下β與水力開度eh之間的負冪律關系;文獻[6]研究單裂隙的幾何屬性對裂隙中非達西流動的影響,發現粗糙度對β有顯著影響;文獻[17]通過分析圍壓下不同裂隙的幾何形狀對β的影響,發現連通孔隙區域的不均勻分布對β會產生影響;文獻[18]研究表明,eh僅取決于單裂隙的幾何特性,且與流體慣性效應無關;文獻[5]通過使用表觀滲透率ka取代固有滲透率ko,將式(2)應用于單裂隙介質中。單裂隙的表觀滲透率ka計算公式為:

(3)

eh=[12μQ/(Pw)]1/3

(4)

其中:Q為流量;w為裂隙寬度。

上述研究結果表明,單裂隙的幾何屬性(如粗糙度、開度等)決定β的大小,β與eh之間的數值關系也表明其與單裂隙幾何屬性的聯系密切。然而,需要注意的是,無論是在多孔介質還是裂隙介質中,已有研究結果大多是在忽略流體慣性效應條件下的相關實驗或數值模擬結果。由于流體慣性效應對β的影響幾乎可以被忽略,β只表現出與裂隙幾何屬性密切相關的特性。一旦流體流動不再遵循達西定律,或者流體慣性效應不能被忽略,β是否仍然僅取決于介質的幾何性質,目前尚不清楚。

理論上,含有β的Forchheimer方程的二次項系數代表流體因慣性力而產生的能量耗散,因此流體慣性效應與β可能密切相關。關于流體慣性效應對β的影響,相關研究中有一些間接證據。文獻[2]通過用砂和玻璃珠填充的柱實驗發現,β的變化可能與流體流態的轉變相關;文獻[19]認為β可能與流體性質有關;文獻[20]通過數值模擬發現,單裂隙的ka會隨著雷諾數的增加,即流體慣性效應的增加而降低;文獻[21]進一步證明粗糙單裂隙的ka減小與足夠大的流體慣性效應有關。由于ka與β之間的密切關系,流體慣性效應對ka的影響將進一步影響β。

本文通過揭示流體慣性效應對粗糙單裂隙中β的影響,提出在考慮流體慣性效應條件下β的量化模型。

1 數值模擬方法

1.1 二維單裂隙物理模型的建立

本文采用數值模擬的方法來探究流體慣性效應對單裂隙β的影響。雖然三維數值模擬在維度上更接近于真實裂隙,但是三維數值模擬方法需要大量的計算機資源和計算時間[22-23];相對于三維數值模擬,二維數值模擬雖然在維度上與真實裂隙維度有一定的差距,但大量研究表明,其結果也是可靠的[20],并且所需的計算機資源和計算時間也大大縮短,因此本文采用二維數值模擬方法。

在COMSOL Multiphysics軟件中建立二維粗糙單裂隙的物理模型,物理模型參考文獻[24]。為了實現單裂隙的不同幾何屬性,通過插入多邊形粗糙元(梯形、矩形和三角形)、復制變換和布爾分割獲得粗糙的裂隙上壁面,將下壁面設置為光滑壁。文獻[24]研究表明,在類似具有規則粗糙元的單裂隙數值模擬中,粗糙元的密度D(即兩個粗糙元之間的間隔與粗糙元的凸起度Δ之比)也會影響流體流動。為了消除D的影響,本文將具有不同凸起度(Δ分別為1.0、0.8、0.4 mm)的3種不同形狀(矩形、梯形和三角形)的粗糙元密度設置為一個定值(D=6)。為獲得足夠的數據量,通過設置單裂隙的水力開度(eh)、粗糙元的形狀和Δ,共獲得45條具有不同幾何屬性的二維粗糙單裂隙。

為了量化單裂隙的幾何屬性,使用文獻[5]提出的裂隙幾何屬性量化方法。二維粗糙單裂隙幾何屬性取值見表1所列。

表1 3種形狀粗糙元二維粗糙單裂隙幾何屬性取值

表1中:em為機械開度;Rp為峰值粗糙度;RSD為相對粗糙度,用相對標準偏差(relative standard deviation)表示;Rrms為均方根粗糙度;bmax為最大開度。

1.2 裂隙流數值模擬

為符合裂隙流動的實際條件,選擇層流模型作為整個數值模擬的模型框架,此模型將單裂隙中的流體設定為不可壓縮牛頓流體,由流體連續性方程和Navier-Stokes方程共同控制,控制方程分別為:

ρ(u)u=μ2u-P

(5)

u=0

(6)

其中:u為速度向量,u=[uxuy];為微分算子。

數值模擬中使用的材料(水,默認為不可壓縮的牛頓流體)需要在構建物理模型和選擇流體模型后進行選擇,這將直接決定流體的密度和黏度等物理參數。

物理模型的最左端設置為流體入口,最右端設置為流體出口。將物理模型的邊界設置為不可滲透的邊界,不發生滑動。入口和出口的邊界條件設定為壓力。入口壓力定義為參數P-in,出口壓力的大小定義為0,入口壓力與出口壓力之差即為壓降(壓降的大小直接決定著單裂隙中流體慣性效應的大小)。

COMSOL Multiphysics軟件的計算流體力學模塊是基于有限元方法對Navier-Stokes方程進行求解的,因此數值模型求解前需要對物理模型進行網格剖分。為保證數值模擬結果的準確性,本研究均采用極細化網格剖分方法。

1.3 非達西效應量化方法

流體慣性效應是β的重要影響因子,需要對其進行精確的量化。作為量化慣性效應的傳統參數,雷諾數(Re)已被廣泛用于多孔和裂隙介質的研究,其計算公式為:

(7)

其中,d為裂隙的特征長度。

然而使用經典的雷諾數來量化單裂隙中的流體慣性效應存在以下問題:式(7)沒有考慮可能導致強慣性效應的因素,如單裂隙的曲折程度和粗糙度;d的確定具有歧義,文獻[25-26]使用單裂隙的機械開度作為特征長度來計算雷諾數,文獻[4]使用水力開度。為解決特征長度的選擇問題,文獻[27]定義一種新類型的雷諾數,即Forchheimer數(Fo)來更準確地量化非達西效應大小,Fo計算公式為:

(8)

從式(8)可以看出,Fo為Forchheimer方程中慣性項和黏性項的比值,在計算時不需要確定特征長度。對Fo進行歸一化處理,可得非達西效應因子E,即

(9)

E可反映流體在流動過程中慣性力耗散導致的壓降占總壓降的比例。

當Jv>Ji時,黏滯力引起的壓降占主導地位,流動為弱慣性流,對應于0Jv時,慣性效應引起的壓降占主導地位,流動為強慣性流,對應于0.5

2 數值模擬結果與分析

2.1 流體慣性效應對β的影響

作為Forchheimer方程慣性項系數的一個重要參數,β問題的探究必須基于非達西流。單裂隙中出現非達西流動的原因有很多,如復雜的幾何屬性[4]、局部渦流的形成[21,28]和足夠大的接觸面[29],其中普遍公認的是當流體慣性效應或流速足夠大時,非達西效應便會出現[9,13]。本文通過增加單裂隙入口和出口之間的壓力梯度(P-in),來獲得足夠的流體慣性效應,以實現單裂隙中的非達西流這一基本流態條件。以5個具有矩形粗糙元且Δ=1.0 mm的單裂隙為例,數值模擬結果如圖1所示。從圖1可以看出,裂隙中的水流偏離達西定律,但可以用非線性的Forchheimer方程很好地描述。

圖1 粗糙單裂隙中流體流動的流速-水力梯度曲線

流體慣性效應對單裂隙ka的影響如圖2所示。從圖2可以看出,ka也隨著雷諾數的增加而降低[28]。圖1、圖2結果表明,粗糙單裂縫中的流體流動出現顯著的非達西效應,即形成了非達西流。

圖2 流體慣性效應對粗糙單裂隙滲透率的影響曲線

基于Forchheimer方程可以得到不同單裂隙中不同慣性效應(用E量化)下的β。單裂隙粗糙度和流體慣性效應對β的影響如圖3所示。從圖3可以看出,單裂隙的β與相對粗糙度RSD成正比,與E的大小成反比。

圖3 單裂隙粗糙度和流體慣性效應對β的影響

傳統上,根據β的大小來判斷流體非達西效應的強度[13],然而,從圖3可以看出,對于具有恒定相對粗糙度的單裂隙,β會隨著流體慣性效應增加而減小。根據式(9),E的大小由式(1)中的黏性項系數a1和慣性項系數b1確定,然而,根據文獻[30]中a1與b1之間的量化關系,變化的β也會引起a1的變化,因此直接使用β的大小來判斷非達西效應的程度并不準確。

2.2 流體慣性效應影響下β的量化

為了量化單裂隙中的β,將β與多孔介質中滲透率k之間的關系擴展到裂隙介質中[5,12],即將式(3)帶入式(2)中,可得:

(10)

其中,C為標度參數。

根據式(10),除eh外,β的大小由比例因子A和冪指數-2b決定。因此,研究流體慣性效應影響下β量化模型的關鍵,在于如何量化慣性效應對上述2個參數的影響。粗糙度作為單裂隙最重要的幾何屬性之一,其大小會直接影響流體流動中的非達西效應[31-33],這也表明流體慣性效應與裂隙粗糙度對β的影響是相互耦合的,下面分析這些復雜的耦合效應。流體慣性效應(用E量化)對冪指數-2b、比例因子A及參數m的影響如圖4所示。圖4中,m為取決于單裂隙幾何屬性的參數。

圖4 流體慣性效應對冪指數-2b、比例因子A及參數m的影響

1) 冪指數(-2b)的量化分析。從圖4a可以看出:在弱慣性流下,冪指數-2b隨著E增加而增加;在強慣性流下,冪指數-2b沒有規律性的變化趨勢,此現象也表明,當流體的慣性力占主導地位時,流體的流動狀態會變得更加復雜;在弱慣性流下,具有相同Δ、不同形狀粗糙元的單裂隙,其冪指數增長率(斜率)相似,冪指數-2b與E之間正比例關系的斜率隨著粗糙元Δ(或相對粗糙度)的增加而增加,這也驗證了粗糙度和非達西效應之間的正相關性[34]。根據粗糙度和流體慣性效應對冪指數-2b的影響進行無量綱分析,得到單裂隙粗糙度Rrms和E耦合影響下的冪指數量化模型為:

-2b=28.74Rrms(E-0.637)

(11)

A為確定β大小的因素之一,如果將A假定為常數,那么冪指數-2b增加將導致β增加,這標志著非達西效應和β之間的正相關關系。

由圖4b可知,即使在同一個裂隙中,A也會隨著非達西效應的不同程度而變化。因此,需要通過進一步確定A的變化趨勢來解釋2.1節中發現的β與流體慣性效應之間的負相關關系。

2) 比例因子A的量化分析。無論是強慣性流還是弱慣性流,圖4b結果與文獻[7,24-25]的數據類似,即都表明A會隨著E的增大而減小。然而,在圍壓或剪切條件下,A會隨著E增加而增加[5-6,35]。應力過程將對巖石裂隙的幾何屬性產生較大影響,如在圍壓或剪切過程中,裂隙開度將進一步閉合,導致裂隙滲透率、粗糙度的變化和水體與固體相互作用加強[16],這可能是導致A與E之間異常關系的原因。本文未向裂隙施加應力,因此得到A與E之間的負相關關系,并與未施加應力條件的文獻[7,24-25]的數據相吻合。

由式(10)可知,冪指數-2b通過乘數12b影響A。由2.2.1節分析可知,在弱慣性流下,E增加會導致冪指數-2b的增加,而12b則會降低。這解釋了圖4b中的弱慣性流條件下,A隨著E增加而減小。受冪指數控制的乘數12b的存在會導致在分析E對A影響的過程中產生歧義。為了避免上述現象,A可以分解為式(10)所示的12b與參數C的乘積形式,然后對參數C進行初步分析。但在去除E對冪指數(或乘數12b)的影響后,參數C隨著E的增加沒有明顯規律。為了解釋這一現象,對參數C做進一步的數值分析。

本文定義一個新的參數m,其僅取決于單裂隙的幾何屬性,不會受到慣性效應的影響;C=mEn可以表示C受到流體慣性效應和裂隙幾何屬性的耦合影響,其中n為E的無量綱階數。從圖4c可以看出:對于Δ=1.0 mm的裂隙,當n=1時,m在層流(弱慣性流)下保持穩定;對于Δ=0.8 mm的裂隙,n=0;對于Δ=0.4 mm的裂隙,n=-1;在強慣性流下,m呈下降趨勢。這表明E與C之間存在多階反比關系。

此外,m的大小將受到單裂隙幾何屬性的影響。從圖4c還可以看出:m隨著粗糙度增加而增加;具有更大Δ(或相對粗糙度)且形狀相同的粗糙元會產生更大的m;具有相同Δ、不同形狀粗糙元的裂隙,m也不同;三角形的m大于梯形和矩形的m,這證實m也由幾何形狀決定。進一步分析表明,m與Rrms之間存在線性關系,即

m=lRrms

(12)

其中,l為與粗糙元形狀相關的無量綱參數,l矩形=7.437,l梯形=12.361,l三角形=33.540。

最后,結合式(10)～式(12),得到β、E和m之間的量化模型為:

r1=-14.37Rrms(E-0.637),

r2=28.74Rrms(E-0.637)

(13)

式(13)為β量化模型,該模型考慮了流體慣性效應和裂隙幾何形狀對β的耦合影響。

3 結 論

1) 本文為了表征和量化流體慣性效應對非達西系數β的影響,對具有規則形狀粗糙元的粗糙單裂隙進行一系列數值模擬。數值模擬結果表明,流體慣性效應(用非達西效應因子E表征)對β的影響十分顯著,β會隨流體慣性效應的增大而減小,并隨裂隙相對粗糙度的增大而增大。

2) 本文著重量化分析0

3) 基于無量綱分析得到一個同時考慮流體慣性效應和裂隙幾何屬性的β量化模型。