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

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

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

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

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

1 野外裂隙入滲實(shí)驗(yàn)

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

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

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

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

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

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

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

2 數(shù)值模擬方法

2.1 Navier-Stokes方程組

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

連續(xù)性方程：

(1)

動(dòng)量守恒方程：

(2)

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

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

t——時(shí)間/s；

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

p——靜壓/Pa；

G——重力體積力/N；

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

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

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

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

Γ=μ(·ν)

式中：μ——?jiǎng)恿︷ざ?(kg·(m·s)-1)；

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

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

(3)

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

T——絕對(duì)溫度/K。

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

ρ=α1ρ1+α2ρ2

(4)

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

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

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

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

(5)

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

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

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

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

單環(huán)水位對(duì)時(shí)間的離散形式可寫(xiě)為：

(6)

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

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

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

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

(·

(7)

式中：h——水勢(shì)/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]。針對(duì)裂隙入滲實(shí)驗(yàn)，使用FLUENT軟件建模，其中，裂隙中水流驅(qū)替空氣過(guò)程模擬選用FLUENT中的兩相流模型(VOF模型)；針對(duì)入滲單環(huán)水位變化的模擬，由于FLUENT提供的現(xiàn)有邊界模式中，沒(méi)有變水頭入滲邊界模式，本文采用FLUENT提供的用戶接口(UFD)，遵照C語(yǔ)言書(shū)寫(xiě)規(guī)范，編制了模擬單環(huán)水位變化的UFD代碼模塊，作為裂隙水流運(yùn)動(dòng)邊界條件。

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

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

3.1.2裂隙水力寬度推求

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

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

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

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

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

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

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

圖5 FLUENT軟件模擬末期裂隙中壓強(qiáng)、含水量、水流速度分布圖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ò)程需要考慮水氣兩相流運(yùn)動(dòng)，難度較大，而且?guī)r體裂隙寬度較小，裂隙容水量少，水分會(huì)瞬時(shí)充滿裂隙，運(yùn)用立方定律可忽略裂隙水的浸潤(rùn)過(guò)程，僅考慮飽和裂隙水流問(wèn)題。SUTRA是美國(guó)地質(zhì)調(diào)查局(USGS)開(kāi)發(fā)的飽和、非飽和帶水流計(jì)算程序，廣泛應(yīng)用于土壤水及地下水?dāng)?shù)值模擬[19]。

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

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

(8)

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

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

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

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

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

當(dāng)α→∞時(shí)，K→0，同時(shí)θ→θr。這說(shuō)明當(dāng)α取較大值時(shí)，SUTRA可模擬出裂隙中的非飽和節(jié)點(diǎn)特征?？紤]到算法的穩(wěn)定性，經(jīng)過(guò)多次測(cè)試，本研究采用如下非飽和水動(dòng)力參數(shù)值：α=0.002/Pa，n=2，θ=1，θr=0。依據(jù)FLUENT求得的最佳水力裂隙寬度，及立方定律可求得①②③三條裂隙的飽和滲透系數(shù)分別為：0.051, 0.051和0.204 m/s。由于N-S方程在模擬裂隙水流運(yùn)動(dòng)中存在誤差(圖3)，因此N-S方法推求的裂隙寬度可能與實(shí)際情況存在差異，從而導(dǎo)致SUTRA模型中裂隙飽和滲透系數(shù)估計(jì)不準(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é)點(diǎn)編號(hào)及坐標(biāo)，再組織成SUTRA所需的格式即完成數(shù)值計(jì)算所需的網(wǎng)格剖分。

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

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

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

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

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

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

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

3.3 兩種模擬方法比較

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

兩模型單環(huán)水位模擬結(jié)果(圖3、圖6)對(duì)比表明：SUTRA模擬的單環(huán)水位下降速率小于FLUENT計(jì)算結(jié)果。原因可能是SUTRA模型(基于立方定律)計(jì)算的是裂隙斷面平均流速，而N-S模型能刻畫(huà)裂隙橫斷面水流流速的拋物線分布(即“優(yōu)勢(shì)流”)，其實(shí)際流速要快于SUTRA計(jì)算的平均流速。同時(shí)，F(xiàn)LUENT模擬的單環(huán)水位下降速率的變化率(圖6)也大于SUTRA(圖3)，這反映出FLUENT模型對(duì)水位變化的敏感性比SUTRA模型大，因?yàn)槠湓诹严稊嗝嫔狭魉俚牟町愲S著單環(huán)水位(壓強(qiáng))減小而減小。

兩模型模擬的監(jiān)測(cè)斷面出流過(guò)程差異較大(圖4、圖7)，總體說(shuō)來(lái)FLUENT模擬的水流過(guò)程更加細(xì)致，因?yàn)槠淠M了水、氣驅(qū)替過(guò)程，而SUTRA模擬結(jié)果則相對(duì)平順。

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

4 結(jié)論

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

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

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

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