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

(合肥工業(yè)大學(xué) 資源與環(huán)境工程學(xué)院,安徽 合肥 230009)

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

(1)

其中:P為壓力;x為流動(dòng)路徑的長度;ρ為流體密度;g為流體重力加速度;J為總壓降;Jv為黏滯力導(dǎo)致的壓降;Ji為慣性力導(dǎo)致的壓降;a1為黏性項(xiàng)系數(shù);b1為慣性項(xiàng)系數(shù);v為流體平均速度;μ為流體的動(dòng)態(tài)黏度;k為介質(zhì)的滲透率;β為流體流動(dòng)時(shí)所受到的慣性阻力,也稱為非達(dá)西系數(shù)或Forchheimer系數(shù)。

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

基于實(shí)驗(yàn)室實(shí)驗(yàn),文獻(xiàn)[10-11]研究表明,多孔介質(zhì)中的β僅取決于多孔介質(zhì)的幾何屬性,與流體性質(zhì)無關(guān)。文獻(xiàn)[12]將多孔介質(zhì)中β與滲透率之間的關(guān)系統(tǒng)一歸納為負(fù)冪律關(guān)系,即

(2)

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

相關(guān)研究結(jié)果表明,多孔介質(zhì)中β的經(jīng)驗(yàn)表達(dá)式通常與多孔介質(zhì)的滲透率和孔隙度有關(guān)[13]。

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

(3)

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

(4)

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

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

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

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

1 數(shù)值模擬方法

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

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

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

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

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

表1中:em為機(jī)械開度;Rp為峰值粗糙度;RSD為相對(duì)粗糙度,用相對(duì)標(biāo)準(zhǔn)偏差(relative standard deviation)表示;Rrms為均方根粗糙度;bmax為最大開度。

1.2 裂隙流數(shù)值模擬

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

ρ(u)u=μ2u-P

(5)

u=0

(6)

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

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

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

COMSOL Multiphysics軟件的計(jì)算流體力學(xué)模塊是基于有限元方法對(duì)Navier-Stokes方程進(jìn)行求解的,因此數(shù)值模型求解前需要對(duì)物理模型進(jìn)行網(wǎng)格剖分。為保證數(shù)值模擬結(jié)果的準(zhǔn)確性,本研究均采用極細(xì)化網(wǎng)格剖分方法。

1.3 非達(dá)西效應(yīng)量化方法

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

(7)

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

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

(8)

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

(9)

E可反映流體在流動(dòng)過程中慣性力耗散導(dǎo)致的壓降占總壓降的比例。

當(dāng)Jv>Ji時(shí),黏滯力引起的壓降占主導(dǎo)地位,流動(dòng)為弱慣性流,對(duì)應(yīng)于0Jv時(shí),慣性效應(yīng)引起的壓降占主導(dǎo)地位,流動(dòng)為強(qiáng)慣性流,對(duì)應(yīng)于0.5

2 數(shù)值模擬結(jié)果與分析

2.1 流體慣性效應(yīng)對(duì)β的影響

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

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

流體慣性效應(yīng)對(duì)單裂隙ka的影響如圖2所示。從圖2可以看出,ka也隨著雷諾數(shù)的增加而降低[28]。圖1、圖2結(jié)果表明,粗糙單裂縫中的流體流動(dòng)出現(xiàn)顯著的非達(dá)西效應(yīng),即形成了非達(dá)西流。

圖2 流體慣性效應(yīng)對(duì)粗糙單裂隙滲透率的影響曲線

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

圖3 單裂隙粗糙度和流體慣性效應(yīng)對(duì)β的影響

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

2.2 流體慣性效應(yīng)影響下β的量化

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

(10)

其中,C為標(biāo)度參數(shù)。

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

圖4 流體慣性效應(yīng)對(duì)冪指數(shù)-2b、比例因子A及參數(shù)m的影響

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

-2b=28.74Rrms(E-0.637)

(11)

A為確定β大小的因素之一,如果將A假定為常數(shù),那么冪指數(shù)-2b增加將導(dǎo)致β增加,這標(biāo)志著非達(dá)西效應(yīng)和β之間的正相關(guān)關(guān)系。

由圖4b可知,即使在同一個(gè)裂隙中,A也會(huì)隨著非達(dá)西效應(yīng)的不同程度而變化。因此,需要通過進(jìn)一步確定A的變化趨勢來解釋2.1節(jié)中發(fā)現(xiàn)的β與流體慣性效應(yīng)之間的負(fù)相關(guān)關(guān)系。

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

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

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

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

m=lRrms

(12)

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

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

r1=-14.37Rrms(E-0.637),

r2=28.74Rrms(E-0.637)

(13)

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

3 結(jié) 論

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

2) 本文著重量化分析0

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