裂隙水流運動的入滲試驗及數(shù)值模擬研究

程勤波，陳 喜，張志才，張潤潤，丘 寧，黃日超，蔡煉彬(1.河海大學(xué)水文水資源學(xué)院, 江蘇 南京 210098；2.河海大學(xué)水文水資源與水利工程科學(xué)國家重點實驗室, 江蘇 南京 210098)

由于構(gòu)造、風(fēng)化、卸荷等作用，巖體裂隙在自然界中廣泛存在。由于巖體裂隙的滲透性遠大于基質(zhì)巖體，裂隙是巖體的主要導(dǎo)水介質(zhì)。在可溶性石灰?guī)r地區(qū)，裂隙控制水流運動，因此，研究巖體裂隙中流體運動規(guī)律對于了解喀斯特流域水動力過程以及溶質(zhì)運移具有重要意義[1]。

Poiseuille在均質(zhì)恒定運動條件下,導(dǎo)出了單片光滑平行板裂隙中水流運動的理論公式，即立方定律[2]。Lomize[3]和Louis[4]各自采用平行玻璃板模擬裂隙壁，證明了在層流時立方定律的有效性。然而，天然裂隙面多粗糙不光滑，并且裂隙寬度多存在小范圍的不平整和大范圍的起伏，不符合立方定律假設(shè)條件。Djik和Berkowitz運用核磁共振技術(shù)研究飽和粗糙裂隙中水流場分布，結(jié)果顯示裂隙橫斷面水流速度呈類似拋物線，但不完全對稱[5]。鑒于立方定律的局限性，眾多學(xué)者基于室內(nèi)外實驗和數(shù)值模擬結(jié)果提出了各種形式的修正公式，其核心思想是修正裂隙寬度[2～6]。

隙寬等裂隙特征是控制裂隙導(dǎo)水能力及水流運動的主要因素，目前主要根據(jù)自然界巖體出露面測量裂隙隙寬，但由于自然巖體裂隙中存在大量粗糙面，當(dāng)前技術(shù)還無法對自然界裂隙壁的粗糙程度及有效隙寬精確測算。對于非充填裂隙，Navier-Stokes(N-S)方程組描述裂隙水流運動比連續(xù)性方程具有更高的精度，但數(shù)值求解N-S方程組非常困難，因而將N-S方法應(yīng)用于野外原位實驗研究裂隙水流運動規(guī)律鮮有報道[6]。目前大多數(shù)軟件基于立方定律的連續(xù)性方程數(shù)值解模擬裂隙水流運動，主流裂隙水模擬軟件有：TOUGH，RockFlow，F(xiàn)racMan/MAFIC等[7]。

為研究自然巖體裂隙對入滲水流運動作用，本文選取具有垂向與橫向裂隙交錯的典型巖體，利用現(xiàn)場單環(huán)注水入滲試驗，分別基于N-S方程組、立方定律的連續(xù)性方程建立裂隙入滲水流運動數(shù)值模型，以實測單環(huán)水位變化為目標(biāo)，采用試算法，推求裂隙水力隙寬及入滲水流在不同裂隙中的比例。對比兩種模型計算結(jié)果，分析裂隙水?dāng)?shù)值模擬中立方定律適用性。

1 野外裂隙入滲實驗

試驗地位于貴州省普定縣陳旗村董家山陰坡巖石開挖剖面，該處巖石屬于三疊系關(guān)嶺組灰?guī)r。試驗巖塊有兩條較為發(fā)育的縱向裂隙面(圖1中①、③號)，其中，①號裂隙面實測隙寬為1 mm，沿巖石面垂向(圖2中Z方向)裂隙長56 cm，垂直于巖石面(Y方向)裂隙長約60 cm；③號裂隙面實測隙寬為2 mm，沿Z方向、Y方向裂隙面裂隙長度都約為60 cm。在橫向(X方向)，有條較發(fā)育的巖層解理面(②號裂隙)，裂隙面隙寬為1 mm，傾向164°，傾角8°，沿X、Y方向裂隙面裂隙長度分別為60 cm和47 cm(圖1)。

圖1 巖體剖面照片及主要裂隙Fig.1 Photo of the rock profile and fractures

裂隙面均無明顯填充物，實測裂隙寬是裂隙出露面處的平均隙寬。為便于數(shù)值模型中網(wǎng)格劃分及提高數(shù)值解的收斂性，裂隙面均概化為平面，即裂隙寬度沿程不變。

為研究各裂隙對入滲水流影響程度，按上述三條裂隙組合分為三種模式(圖2)：(a)單條垂向裂隙(①號裂隙)；(b)橫向與垂向裂隙組合(①②號裂隙)；(c)一條橫向與兩條垂向裂隙組合(①②③號裂隙)。

單環(huán)裂隙入滲試驗步驟：(1)將單環(huán)罩于上表面裂隙上、水位計置于單環(huán)內(nèi)(圖1，單環(huán)半徑為15 cm)；(2)瞬時將單環(huán)水位加注到H0(=0.0743 m)；(3)自記水位計自動記錄單環(huán)水位下降過程。

圖2 3種巖石裂隙概化模式Fig.2 Three cases of the conceptual rock fractures

根據(jù)雷諾數(shù)(Re)判別裂隙水流的流態(tài)，用以選擇合適的數(shù)學(xué)模型[2]。根據(jù)實測單環(huán)水位下降最大速度得裂隙水流最大雷諾數(shù)為68，遠小于臨界值1 150[2]，所以本次試驗裂隙水流流態(tài)為層流。

2 數(shù)值模擬方法

2.1 Navier-Stokes方程組

非飽和無填充裂隙水流入滲問題可視為水流驅(qū)替空氣問題，該問題屬于流體力學(xué)中的水和空氣的兩相流問題。描述水、空氣運動的Navier-Stokes方程組[8]為(簡稱N-S方程)：

連續(xù)性方程：

(1)

動量守恒方程：

(2)

式中：ρ——流體密度/(kg·m-3)；

ν——流體速度/(m·s-1)；

t——時間/s；

W——源匯項/(kg·(m3·s)-1)；

p——靜壓/Pa；

G——重力體積力/N；

f——其它體積力(如兩相流體間的作用)/N；

?！獞?yīng)力張量/N。

對于可壓縮流體，應(yīng)力張量計算公式為：

對于不可壓縮流，應(yīng)力張量計算公式為：

Γ=μ(·ν)

式中：μ——動力黏度/(kg·(m·s)-1)；

μν——第二動力黏度/(kg·(m·s)-1)；

對可壓流體(空氣)，采用理想氣體方程描述密度與壓強的關(guān)系[8]：

(3)

式中：R——普適氣體常數(shù)/(m2·(T·s2)-1)；

T——絕對溫度/K。

對于水、空氣兩相流問題可采用容積比率模型(簡稱VOF)模擬：VOF模型假設(shè)水與空氣不相容，并且兩流體不會互相穿插，在整個流場內(nèi)采用統(tǒng)一N-S方程組求解，但在不同相中方程組參數(shù)不同(如流體密度、動力粘度、壓縮公式)。在流體交界面處，計算單元格各相的容積比率(容積比率之和為1)，根據(jù)容積比率計算方程組在該單元格內(nèi)參數(shù)的平均值[8]，單元格平均密度ρ可采用如下公式計算：

ρ=α1ρ1+α2ρ2

(4)

式中：α1，α2——水和空氣所占容積比；

ρ1，ρ2——水和空氣密度/(kg·m-3)。

2.2 單環(huán)水位變化模擬

相對于裂隙，由于單環(huán)容積大，水位變化緩，單環(huán)內(nèi)水動力現(xiàn)象不突出，單環(huán)水位變化控制方程采用如下水量平衡方程：

(5)

式中：H——單環(huán)水位/m；

Q——裂隙入流率/(m·s-1)；

νf——裂隙入滲面的流速/(m·s-1)；

S——裂隙入滲面的面積/m2。

單環(huán)水位對時間的離散形式可寫為：

(6)

式中：Qtn——t時刻，n次迭代期裂隙面入流率/(m·s-1)；

Qt-Δt——t-Δt時刻裂隙面入流率/(m·s-1)。

2.3 基于立方定律的連續(xù)性方程

對于具有充填物的裂隙或裂隙隙寬小的非填充裂隙，水流速度一般較小，多屬符合達西水流規(guī)律的層流，可采用基于達西定律和立方定律的連續(xù)性方程研究裂隙水流運動[7,9]，其非飽和流微分方程形式為：

(·

(7)

式中：h——水勢/m；

K——滲透系數(shù)/(m·s-1)；

C——立方定律修正系數(shù)；

d——裂隙隙寬/m；

Θ——飽和度；

Ss——貯水率/m-1；

c——比水容量/m-1。

3 數(shù)值模擬及結(jié)果分析

3.1 基于N-S方程組的數(shù)值模型構(gòu)建及模擬結(jié)果

3.1.1數(shù)值模型構(gòu)建

本文采用基于有限體積法求解N-S方程組的FLUENT軟件建模求解[8]。針對裂隙入滲實驗，使用FLUENT軟件建模，其中，裂隙中水流驅(qū)替空氣過程模擬選用FLUENT中的兩相流模型(VOF模型)；針對入滲單環(huán)水位變化的模擬，由于FLUENT提供的現(xiàn)有邊界模式中，沒有變水頭入滲邊界模式，本文采用FLUENT提供的用戶接口(UFD)，遵照C語言書寫規(guī)范，編制了模擬單環(huán)水位變化的UFD代碼模塊，作為裂隙水流運動邊界條件。

裂隙水流數(shù)值模型構(gòu)建中非常重要、也較為困難的是網(wǎng)格剖分。本文先采用Gambit軟件[8]將裂隙面剖分成1 cm×1 cm的四邊形網(wǎng)格，再沿裂隙壁將裂隙剖分成6層，構(gòu)建出裂隙空間的三維網(wǎng)格系統(tǒng)。

初始狀態(tài)設(shè)置為空氣，初始含水量設(shè)置為零，無殘余含水量。N-S模型參數(shù)(如動力粘度(μ)，流體密度(ρ)和重力加速度(g))均采用FLUENT軟件默認值，這些值均源于大量實驗測定結(jié)果[8]。

3.1.2裂隙水力寬度推求

由于天然裂隙壁粗糙不光滑，裂隙寬度沿程變化，因此實測裂隙寬度與裂隙的水力寬度一般不一致，研究表明水力寬度遠小于實際寬度[6]。本文以數(shù)值模型計算的單環(huán)水位與實測水位之差最小為目標(biāo)，采用試錯法推求實際剖面中三條裂隙組合情形下(圖2c)各裂隙水力寬度。由于每次調(diào)整裂隙寬度都需要重新做網(wǎng)格剖分，比較繁瑣，因此為減少計算量，三條裂隙水力寬度按裂隙實測寬度同比調(diào)整。推求的①②③號裂隙水力寬度分別為：0.25 mm、0.25 mm、0.50 mm。數(shù)值模型計算的單環(huán)水頭與實測水頭變化過程見圖3，兩者較為接近。圖3中模擬前期，模擬值明顯偏小于實際值，這主要是由于在水流浸潤過程中裂隙的水動力寬度比水流穩(wěn)定之后的水動力寬度大，而試錯法擬合得到的隙寬僅為該裂隙在時間、空間上的綜合(平均)水動力寬度。

3.1.3不同裂隙組合情形下入滲水流模擬

三種裂隙組合情形下(圖2)分別模擬的單環(huán)水位變化過程見圖3?？梢钥闯?，水頭模擬誤差隨著巖體裂隙概化合理性增大而減小，僅一條垂向裂隙和兩條裂隙組合情形比三條裂隙情形所模擬的單環(huán)水位高。在①、②兩條裂隙間設(shè)定監(jiān)測斷面I、在②、③兩條裂隙間設(shè)定監(jiān)測斷面II(圖2)，計算通過斷面I、II水流速度過程見圖4，入滲結(jié)束時裂隙中壓強、含水量、水流速度分布見圖5。

根據(jù)圖4中設(shè)定的監(jiān)測斷面水流速度模擬結(jié)果，得到通過斷面I、II的入滲水量：兩條裂隙情形下(圖2b)，83.3%入滲水沿垂向裂隙(①號)下滲，剩余16.7%通過水平向斷面I入滲；三條裂隙模擬情形下(圖2c)，沿垂向裂隙(①號)入滲水流占75.9%，通過斷面I的水平流量占總?cè)霛B量的24.1%，該部分水流中，通過斷面II的流量占總?cè)霛B量的21.8%。因此，三條裂隙情形下通過水平向I斷面的入滲水量比兩條裂隙情形大7.4%，說明入滲水運動不僅取決于垂向、水平向裂隙發(fā)育程度，而且取決于裂隙連通程度。但總體而言，與單環(huán)入滲面直接連接的垂向裂隙入滲水量占主要部分。

圖3 三種裂隙情形中模擬的單環(huán)水位變化與實測值對比圖Fig.3 Comparison of the measured and simulated water levels for three rock fracture cases

圖4 FLUENT中兩監(jiān)測斷面流速變化圖Fig.4 Simulated flow velocity of two monitoring sections with the FLUENT model

圖5壓強分布表明，F(xiàn)LUENT軟件計算出的壓強具有較好的連續(xù)性；含水量分布中三種情形求得的裂隙①的含水量分布差異較大，表明FLUENT中兩相流模型(VOF模型)在復(fù)雜水流情況下會出現(xiàn)異常結(jié)果，但這并沒有影響壓強及速度計算結(jié)果；流速分布表明，滲出面附近流速較大，水流多從溢流面底部滲出。

圖5 FLUENT軟件模擬末期裂隙中壓強、含水量、水流速度分布圖Fig.5 Pressure, water content and velocity distribution at the end of simulation with FLUENT

3.2 基于立方定律的數(shù)值模型構(gòu)建及模擬結(jié)果

3.2.1數(shù)值模型構(gòu)建

非填充裂隙中水分入滲過程需要考慮水氣兩相流運動，難度較大，而且?guī)r體裂隙寬度較小，裂隙容水量少，水分會瞬時充滿裂隙，運用立方定律可忽略裂隙水的浸潤過程，僅考慮飽和裂隙水流問題。SUTRA是美國地質(zhì)調(diào)查局(USGS)開發(fā)的飽和、非飽和帶水流計算程序，廣泛應(yīng)用于土壤水及地下水?dāng)?shù)值模擬[19]。

利用SUTRA模擬裂隙入滲水流實驗時，有兩類邊界需要特殊處理：(1)單環(huán)入滲邊界。根據(jù)單環(huán)水位變化公式編寫單環(huán)水位模擬程序，將計算的單環(huán)內(nèi)水壓設(shè)置為SUTRA中的變水頭邊界。(2)裂隙溢流邊界?？刹捎靡缌髅嬖O(shè)置方法，將溢流面處計算節(jié)點設(shè)置為SUTRA中的變水頭邊界，根據(jù)臨近節(jié)點壓強計算該處節(jié)點壓強，當(dāng)臨近結(jié)點壓強大于0(飽和)時，則該節(jié)點壓強設(shè)置為0，反之，該節(jié)點壓強設(shè)置與臨近節(jié)點壓強相等。該方法可有效防止回流，確保滲流面出流通暢?？紤]到裂隙的儲水能力小，裂隙含水量對裂隙水量影響較小，因此將裂隙初始壓強值設(shè)置為0，確保裂隙滲流面出流通暢。

為仿真模擬裂隙水流運動，在SUTRA中將所有計算節(jié)點區(qū)分為：飽和節(jié)點和非飽和節(jié)點。飽和節(jié)點含水量為1，滲透系數(shù)為飽和滲透率，而非飽和節(jié)點含水量接近0，滲透系數(shù)也應(yīng)該近似為0。SUTRA采用VGM模型描述非飽和水動力特征：

(8)

式中：θ——土壤飽和含水量；

θr——土壤殘余含水量；

K——土壤非飽和滲透系數(shù)/(m·s-1)；

Ks——土壤飽和滲透系數(shù)/(m·s-1)；

α，n，m——VGM模型經(jīng)驗參數(shù)，其中m=1-1/n。

當(dāng)α→∞時，K→0，同時θ→θr。這說明當(dāng)α取較大值時，SUTRA可模擬出裂隙中的非飽和節(jié)點特征。考慮到算法的穩(wěn)定性，經(jīng)過多次測試，本研究采用如下非飽和水動力參數(shù)值：α=0.002/Pa，n=2，θ=1，θr=0。依據(jù)FLUENT求得的最佳水力裂隙寬度，及立方定律可求得①②③三條裂隙的飽和滲透系數(shù)分別為：0.051, 0.051和0.204 m/s。由于N-S方程在模擬裂隙水流運動中存在誤差(圖3)，因此N-S方法推求的裂隙寬度可能與實際情況存在差異，從而導(dǎo)致SUTRA模型中裂隙飽和滲透系數(shù)估計不準(zhǔn)。所以此處將集中討論SUTRA與FLUENT模擬結(jié)果的差異，并據(jù)此探討裂隙水?dāng)?shù)值模擬中立方定律的適用性，而忽略SUTRA模型模擬單環(huán)水頭的準(zhǔn)確度。

數(shù)值模型SUTRA網(wǎng)格劃分同樣采用Gambit網(wǎng)格剖分軟件：先將裂隙面剖分成1 cm×1 cm網(wǎng)格，再沿裂隙壁將裂隙剖分3層。網(wǎng)格劃分好后，可輸出網(wǎng)格節(jié)點編號及坐標(biāo)，再組織成SUTRA所需的格式即完成數(shù)值計算所需的網(wǎng)格剖分。

3.2.2不同裂隙組合情形下入滲水流模擬

根據(jù)裂隙入滲實驗的概化結(jié)果，利用SUTRA建立三種裂隙組合情形的數(shù)值模型，模擬的單環(huán)水位變化見圖6，通過監(jiān)測斷面I、II的流速見圖7，入滲結(jié)束時裂隙中壓強、含水量、水流速度分布見圖8。

圖6 三種裂隙分布情形中單環(huán)水位模擬與實測值對比圖Fig.6 Comparison of the measured and simulated water levels for three rock fracture cases

圖6顯示隨著裂隙條數(shù)的增加，SUTRA模擬的單環(huán)水位下降速度增大。圖7中，兩條裂隙模擬情形下通過I的流量占總?cè)霛B量的19.5%，三條裂隙模擬情形下通過I的流量占總?cè)霛B量的33.8%，通過II的流量占總?cè)霛B量的31.0%；與FLUENT相比，沿水平側(cè)滲量要偏大，但總體而言，沿主滲流面(①號裂隙面)入滲量占主導(dǎo)，三條裂隙情形下，所占比例為66.2%。

圖8表明，SUTRA計算出的壓強和含水量分布具有較好的連續(xù)性，尤其是含水量分布并未出現(xiàn)不合理現(xiàn)象，表明SUTRA模型具有較好的穩(wěn)定性；流速分布顯示滲出面底部水流速度較大，表明水流多從溢流面底部滲出。

圖7 SUTRA模擬的斷面流量變化圖Fig.7 Simulated flow velocity of two sections with SUTRA

圖8 SUTRA模擬末期裂隙中壓強、含水量和水流速度分布圖Fig.8 Pressure, water content and velocity distribution at the end of simulation with SUTRA

3.3 兩種模擬方法比較

理論上，N-S方程可以推導(dǎo)出達西定律與立方定律[3]，換而言之，基于立方定理的連續(xù)方程水流模型僅為N-S方程在裂隙介質(zhì)中的簡化形式。然而N-S方程求解過于復(fù)雜，因此在實際應(yīng)用中常采用基于立方定理的連續(xù)性方程模型。

兩模型單環(huán)水位模擬結(jié)果(圖3、圖6)對比表明：SUTRA模擬的單環(huán)水位下降速率小于FLUENT計算結(jié)果。原因可能是SUTRA模型(基于立方定律)計算的是裂隙斷面平均流速，而N-S模型能刻畫裂隙橫斷面水流流速的拋物線分布(即“優(yōu)勢流”)，其實際流速要快于SUTRA計算的平均流速。同時，F(xiàn)LUENT模擬的單環(huán)水位下降速率的變化率(圖6)也大于SUTRA(圖3)，這反映出FLUENT模型對水位變化的敏感性比SUTRA模型大，因為其在裂隙斷面上流速的差異隨著單環(huán)水位(壓強)減小而減小。

兩模型模擬的監(jiān)測斷面出流過程差異較大(圖4、圖7)，總體說來FLUENT模擬的水流過程更加細致，因為其模擬了水、氣驅(qū)替過程，而SUTRA模擬結(jié)果則相對平順。

SUTRA模擬的裂隙壓強和流速分布(圖8)與FLUENT模擬結(jié)果(圖5)近似，但含水量分布與FLUENT有明顯差異，這反映出SUTRA中用于概化兩相流運動的方法存在誤差，主要表現(xiàn)為非飽和節(jié)點含水量消退較慢。同時，SUTRA模擬的流速分布(尤其在滲流面附近)比FLUENT模擬結(jié)果平緩，其主要原因可能是SUTRA模擬的單元格出流能力(即平均流速)小于FLUENT模擬結(jié)果(即可描述“優(yōu)勢流”的拋物線型流速)。該因素同時也導(dǎo)致SUTRA模擬的溢流高度大于FLUENT結(jié)果?？傮w而言，基于Darcy定律和立方定律的連續(xù)性方程模型可用于近似模擬裂隙水流運動。

4 結(jié)論

(1)以實測單環(huán)水位變化為目標(biāo)，F(xiàn)LUENT推求的裂隙等效水力寬度僅為實測裂隙寬度的1/4，遠小于實測隙寬，這與以往研究成果相符。

(2)FLUENT計算結(jié)果表明，巖體裂隙概化越符合實際情況，其單環(huán)水頭模擬誤差越??；裂隙水流大部分沿垂向裂隙下滲，但也有部分(約20%)水流沿橫向裂隙側(cè)向運動。

(3)SUTRA模擬結(jié)果表明，隨著裂隙條數(shù)的增加，SUTRA模擬的單環(huán)水位下降速度增大，但單環(huán)水位下降速率的變化率小于FLUENT計算結(jié)果，這是由于N-S方程能刻畫裂隙橫斷面水流速度的拋物線分布(即“優(yōu)勢流”)，而立方定律則只能反映裂隙斷面的平均流速；SUTRA模擬的壓強、流速分布與FLUENT模擬結(jié)果接近，因此基于達西定律和立方定律的連續(xù)性方程模型可用于近似模擬裂隙水流過程。

(4)在實際應(yīng)用中，由于N-S方程求解復(fù)雜繁瑣，不便于區(qū)域數(shù)值模擬計算，因此采用連續(xù)性方程求解裂隙水流問題更為方便有效。