牛赟凱,王全榮
(中國(guó)地質(zhì)大學(xué)(武漢)環(huán)境學(xué)院,湖北 武漢 430078)
管道是含水層中十分常見(jiàn)的水文地質(zhì)結(jié)構(gòu),包括自然條件下的巖溶管道、動(dòng)物的洞穴、蟲(chóng)孔,也包括人為的管道,如人工挖掘的水平井、隧道和探礦巷道等。由于管道與周?chē)慕橘|(zhì)在滲透性方面具有較大的差異,這些管道的存在使含水層呈現(xiàn)高度非均質(zhì)性,因此成為管道-孔隙雙重介質(zhì)系統(tǒng),其中高滲透的管道是一個(gè)系統(tǒng),低滲透的孔隙或裂隙基質(zhì)是一個(gè)系統(tǒng),導(dǎo)致地下水流場(chǎng)的模擬十分困難[1]。隨著工業(yè)化進(jìn)程的不斷發(fā)展,面對(duì)日益復(fù)雜的地下水環(huán)境變化,一個(gè)良好的管道流模型是研究地下水污染物遷移[2]、巷道涌水[3]等環(huán)境安全問(wèn)題的基礎(chǔ),對(duì)進(jìn)一步刻畫(huà)含水層復(fù)雜的地質(zhì)環(huán)境與滲流場(chǎng)具有十分重要的意義。
目前,針對(duì)含水層中不同類(lèi)型的管道流,已經(jīng)開(kāi)發(fā)了許多數(shù)值模型,包括等效滲透系數(shù)模型、多節(jié)點(diǎn)井流(multi-node well,MNW)模型和管道流(conduit flow process,CFP)模型3種類(lèi)型。陳崇希等[4-5]提出“滲流-管流模型”,該模型利用等效滲透系數(shù)法將管道處理為高滲透介質(zhì),并將沿水流方向的滲透系數(shù)設(shè)置得很大,其取值基于不同流動(dòng)狀態(tài)下管道流模型計(jì)算的等效滲透系數(shù);Halford等[6]和Konikow等[7]基于MODFLOW軟件平臺(tái)開(kāi)發(fā)了MNW模型,該模型為不同網(wǎng)格中的管道按序分配節(jié)點(diǎn),將所有節(jié)點(diǎn)處理為一個(gè)連通的管道,再利用三維地下水流動(dòng)方程來(lái)計(jì)算含水層介質(zhì)網(wǎng)格的水頭,則管道中水頭可利用含水層介質(zhì)網(wǎng)格水頭加管道中水頭損失項(xiàng)來(lái)進(jìn)行計(jì)算。多節(jié)點(diǎn)井流模型分為MNW1模型及其修訂版本MNW2模型,兩者的算法相同,其中MNW2模型較MNW1模型改善了對(duì)非垂直井與非完整井的模擬效果,并可以根據(jù)指定的泵性能(揚(yáng)程-容量)曲線對(duì)流量進(jìn)行調(diào)整,擴(kuò)充了MNW模型的適用范圍。如:Shoemaker等[8]在MODFLOW軟件中加入了管道流模塊CFPv1;Reimann等[9-11]和Kavousi等[12]后續(xù)擴(kuò)展了管道流(CFP)模型的應(yīng)用范圍,將其升級(jí)為CFPv2改進(jìn)程序,加入了管道抽水、溶質(zhì)與熱的轉(zhuǎn)輸?shù)裙δ?為含水層中地下水管道流的模擬研究提供了更豐富的場(chǎng)景。CFP模型采用管道-孔隙雙介質(zhì)系統(tǒng),利用有限差分法將多孔介質(zhì)含水層處理為矩形網(wǎng)格,將離散管道節(jié)點(diǎn)嵌入含水層網(wǎng)格模型中,使用Hagen-Poiseuille方程[13]或結(jié)合Colebrook-White公式的Darcy-Weisbach方程[14]計(jì)算含水層管道中的地下水流量,并使用三維地下水流動(dòng)方程計(jì)算含水層中的層流,最后通過(guò)管道邊界將交換流量耦合在一起。
MNW模型在進(jìn)行含水層中地下水管道流求解時(shí),采用管道節(jié)點(diǎn)所處網(wǎng)格的滲透系數(shù)進(jìn)行運(yùn)算[6],在模擬中一般假設(shè)存在管道的含水層介質(zhì)網(wǎng)格滲透系數(shù)為孔隙介質(zhì)滲透系數(shù),實(shí)際上,由于管道具有貫通性與高水力傳導(dǎo)性能,含水層的水力傳導(dǎo)能力將受到管道幾何形狀與分布狀態(tài)的影響[15-19]。管道占據(jù)有限差分網(wǎng)格一定空間,管道-含水層網(wǎng)格系統(tǒng)具有很強(qiáng)的非均質(zhì)性,使管道-含水層網(wǎng)格的整體滲透系數(shù)不能按照孔隙介質(zhì)的滲透系數(shù)處理,其取值應(yīng)大于孔隙介質(zhì)的滲透系數(shù)。等效滲透系數(shù)法可以定量計(jì)算并概化不同管道流動(dòng)狀態(tài)下含水層的平均等效滲透系數(shù),為CFP模型及在MODFLOW軟件中進(jìn)行含水層中地下水管道流與層流計(jì)算提供了滲透系數(shù)選取的理論依據(jù)。
因此,對(duì)于MNW模型中的管道網(wǎng)格滲透系數(shù)取值問(wèn)題,本研究采用MNW2版本的多節(jié)點(diǎn)井流模型和CFPv2版本的管道流模型,利用等效滲透系數(shù)法修正多節(jié)點(diǎn)井流模型,并基于室內(nèi)二維砂槽管道流試驗(yàn)對(duì)修正后的多節(jié)點(diǎn)井流模型的模擬結(jié)果與管道流(CFP)模型模擬流場(chǎng)、實(shí)測(cè)流場(chǎng)進(jìn)行了對(duì)比分析,為尋找更加準(zhǔn)確而高效的含水層中地下水管道流數(shù)值模擬方法提供參考。
由于管道的滲透性顯著高于含水層孔隙介質(zhì)的滲透性,因此需要為管道設(shè)置一個(gè)較大的滲透系數(shù),將管道系統(tǒng)等效為高滲透介質(zhì)系統(tǒng)。管道的滲透系數(shù)取值與管道中地下水的流動(dòng)狀態(tài)有關(guān)。
在層流狀態(tài)下,當(dāng)管道流的雷諾數(shù)Re<2 300時(shí),管道的等效滲透系數(shù)表達(dá)式如下:
(1)
式中:Kn為管道的等效滲透系數(shù)[L/T];d為管道直徑[L];γ為流體容重[M/L3];μ為流體動(dòng)力黏滯系數(shù)[M/LT]。
在紊流狀態(tài)下,管道的等效滲透系數(shù)表達(dá)式為:
(2)
式中:g為重力加速度[L/T2];f為摩擦系數(shù)[無(wú)量綱];v為流體流速[L/T]。
管道在含水層網(wǎng)格中的示意圖,如圖1所示。
圖1 管道在含水層網(wǎng)格中的示意圖Fig.1 Diagram of a conduit in aquifer grid
由于等效滲透系數(shù)法揭示的含水層中地下水管道流動(dòng)規(guī)律在形式上與達(dá)西定律一致,因此在管道截面積小于含水層孔隙介質(zhì)網(wǎng)格面積時(shí),可利用達(dá)西定律與質(zhì)量守恒定律計(jì)算管道網(wǎng)格的平均等效滲透系數(shù)[20-21]。在圖1中,Aw為管道截面面積[L2];Ap為含水層孔隙介質(zhì)網(wǎng)格面積[L2];Kw為管道等效滲透系數(shù)[L/T];Kp為含水層孔隙介質(zhì)滲透系數(shù)[L/T];Qw為管道流量[L3/T];Qp為孔隙介質(zhì)流量[L3/T]。
若假設(shè)管道網(wǎng)格的平均等效滲透系數(shù)為Kavg,則有:
(3)
若網(wǎng)格剖分與管道截面一致,則可忽略含水層孔隙介質(zhì)項(xiàng)KpAp,使
Kavg=Kw
(4)
該方法要求依據(jù)管道流不同流動(dòng)狀態(tài)來(lái)計(jì)算相應(yīng)的等效滲透系數(shù)并賦值,若由于資料問(wèn)題管道的等效滲透系數(shù)難以計(jì)算,可以初步賦予一個(gè)較大的滲透系數(shù),在后續(xù)模型驗(yàn)證的過(guò)程中將模擬值與實(shí)測(cè)值進(jìn)行對(duì)比校準(zhǔn),以調(diào)整管道的等效滲透系數(shù)。
MODFLOW-CFP模型的優(yōu)勢(shì)在于能夠較詳細(xì)地刻畫(huà)管道在含水層中的分布形態(tài)與幾何特征,是目前使用最廣泛的含水層中地下水管道流模擬程序。該模型采用管道-孔隙雙重介質(zhì)系統(tǒng),孔隙介質(zhì)含水層系統(tǒng)采用MODFLOW程序模擬,管道系統(tǒng)采用CFP模塊模擬,兩者之間的水量交換通過(guò)源匯項(xiàng)的方式加入各自的計(jì)算方程中。
根據(jù)質(zhì)量守恒定律建立管道內(nèi)的水量平衡方程,用于計(jì)算網(wǎng)格水頭與管道水頭的差值:
(5)
式中:Qip為第i段管道的流量[L3/T];Qex為交換流量[L3/T];QR為外源補(bǔ)給/排泄流量[L3/T];Qs為管道中的含水層與管道之間的儲(chǔ)量變化[L3/T],當(dāng)管道為承壓狀態(tài)時(shí),管道節(jié)點(diǎn)儲(chǔ)蓄量變化為0,可不計(jì)入公式。
管道流量公式采用Hagen-Poiseuille方程與結(jié)合Colebrook-White公式的Darcy-Weisbach方程。Hagen-Poiseuille方程中,管道兩端的水頭損失是由于管道中層流存在自身的流動(dòng)阻力造成的,流動(dòng)阻力與液體黏滯系數(shù)有關(guān);Darcy-Weisbach方程中考慮了管道壁的粗糙性。
含水層與管道之間的交換流量Qex可表示為:
Qex=αi,j,k(hin-hi,j,k)
(6)
式中:αi,j,k為i、j、k網(wǎng)格單元中含水層系統(tǒng)與管道系統(tǒng)的交換系數(shù)[L2/T];hin為第i段管道的水頭[L];hi,j,k為i、j、k網(wǎng)格單元中含水層的網(wǎng)格水頭[L]。
含水層與管道之間的水量交換通過(guò)其交換系數(shù)αi,j,k來(lái)控制,其計(jì)算公式為:
(7)
式中:(Kw)i,j,k為含水層i、j、k網(wǎng)格單元的滲透系數(shù)[L/T];Δlip為第i段管道長(zhǎng)度[L];dip為第i段管道直徑[L];τip為第i段管道粗糙度[無(wú)量綱];rip為第i段管道半徑[L]。
MNW模型在有限差分網(wǎng)格中根據(jù)管道分布按序分配節(jié)點(diǎn),將所有節(jié)點(diǎn)作為整體,處理為一個(gè)連通的管道,水流在通過(guò)管道運(yùn)動(dòng)的過(guò)程中會(huì)產(chǎn)生多項(xiàng)水頭損失[6],該模型考慮了水流在管道中運(yùn)動(dòng)的部分水頭損失,因此在對(duì)含水層中地下水管道流的模擬上能夠更加精確地刻畫(huà)管道水頭。
MNW模型采用一個(gè)通用的水頭損失方程來(lái)計(jì)算含水層網(wǎng)格水頭與管道水頭的差值:
(8)
式中:hwell為管道水頭[L];hn為含水層中地下水網(wǎng)格水頭[L];Qn為第n個(gè)管道節(jié)點(diǎn)的流量[L3/T];A為線性含水層水頭損失系數(shù)[T/L2];B為線性水頭損失系數(shù)[T/L2];C為非線性水頭損失系數(shù)[TP/L(3P-1)];P為非線性流量的冪。
公式(8)表明,管道的水頭等于管道所在含水層網(wǎng)格水頭(hn)加上幾個(gè)水頭損失項(xiàng)。在MODFLOW軟件中使用有限差分法求解地下水流動(dòng)偏微分方程,獲得hn[22]。
第一個(gè)水頭損失項(xiàng)(AQn)描述了由于管道的半徑小于井所在單元的尺寸而導(dǎo)致的含水層水頭損失。在忽略表皮效應(yīng)和局部紊流效應(yīng)的情況下,公式(8)中含水層水頭損失AQn可用Thiem穩(wěn)態(tài)流動(dòng)方程來(lái)表示:
(9)
式中:T=Kb為含水層的導(dǎo)水系數(shù)[L2/T](其中K為含水層的滲透系數(shù)[L/T];b為含水層的厚度[L]);ro為有限差分單元的有效半徑[L];rw為管道的實(shí)際半徑[L]。
公式(9)中計(jì)算含水層導(dǎo)水系數(shù)的含水層滲透系數(shù)為管道網(wǎng)格滲透系數(shù),系數(shù)A的計(jì)算公式為
(10)
第二個(gè)水頭損失項(xiàng)(BQn)是指在管道端點(diǎn)和篩管附近及內(nèi)部發(fā)生的水頭損失,系數(shù)B的計(jì)算公式為
(11)
式中:Skin為衡量表皮效應(yīng)大小的無(wú)量綱表皮系數(shù)。
修正后的水頭損失方程為
(12)
當(dāng)監(jiān)測(cè)數(shù)據(jù)不足時(shí),無(wú)法計(jì)算相關(guān)系數(shù),則可使用多節(jié)點(diǎn)井流模型中的“SPECIFYcwc”模式,利用MNW2模塊中封裝的網(wǎng)格到管道的水力傳導(dǎo)系數(shù)(CWCn),手動(dòng)賦予一個(gè)近似的初始值[6,23],并在后續(xù)模型驗(yàn)證的過(guò)程中將模擬值與實(shí)測(cè)值進(jìn)行對(duì)比校準(zhǔn)。
為了對(duì)比修正前后多節(jié)點(diǎn)井流模型的模擬效果,本文以文獻(xiàn)[24]、[25]中二維砂槽管道流試驗(yàn)為例,采用不同模型對(duì)二維砂槽管道流進(jìn)行了數(shù)值模擬,并將水頭模擬結(jié)果與實(shí)測(cè)結(jié)果進(jìn)行了對(duì)比。具體參數(shù)如下:砂槽整體尺寸為57.8 cm×2 cm×26 cm(長(zhǎng)×寬×高),管道直徑為2 cm,位于砂槽底部,砂槽兩側(cè)為長(zhǎng)5 cm的水箱,含水層規(guī)格為47.8 cm×2 cm×24 cm;承壓含水層兩側(cè)為定水頭邊界,管道進(jìn)水口為定水頭邊界;含水層為均質(zhì)各向同性,含水層滲透系數(shù)為0.061 9 cm/s,含水層儲(chǔ)水系數(shù)為2×10-5m-1,含水層初始水位為26 cm,左側(cè)進(jìn)水處定水頭值為30 cm,右側(cè)出水處定水頭值為23 cm,管道進(jìn)水口處定水頭值為34 cm(圖2)。整個(gè)試驗(yàn)時(shí)長(zhǎng)為240 s。
圖2 二維砂槽管道流試驗(yàn)裝置與實(shí)測(cè)水頭[24-25]Fig.2 Experimental device and measured head distribution for 2D sand tank pipeline flow[24-25]
根據(jù)砂槽形態(tài)構(gòu)建二維砂槽網(wǎng)格模型,將三維砂槽網(wǎng)格模型剖分為22列×13行,不同區(qū)域分別表示含水層、管道和定水頭裝置(圖3)。含水層網(wǎng)格從上至下12行,分為12層,從左至右分為20列,左右兩側(cè)各有一列寬5 cm的網(wǎng)格表示用于控制定水頭邊界條件的水槽,底部為直徑2 cm的管道。含水層中水文地質(zhì)參數(shù)和邊界條件與試驗(yàn)相同,初始水頭設(shè)置為26 cm。
圖3 二維砂槽網(wǎng)格模型示意圖Fig.3 Schematic diagram of the 2D sand tank grid model
采用MODFLOW-CFP模型對(duì)二維砂槽管道流進(jìn)行數(shù)值模擬時(shí),在砂槽網(wǎng)格模型第13層插入22個(gè)節(jié)點(diǎn)控制一條管道,管道網(wǎng)格采用孔隙介質(zhì)滲透系數(shù)。管道模型中,管道彎曲度為1,假設(shè)管道光滑,粗糙度設(shè)置為0.000 000 1;管道流的雷諾數(shù)Re采用經(jīng)驗(yàn)值,下界雷諾數(shù)設(shè)計(jì)為2 000,上界雷諾數(shù)設(shè)計(jì)為4 000;在管道系統(tǒng)中將進(jìn)水口與出水口設(shè)置為定水頭,管道系統(tǒng)與含水層網(wǎng)格系統(tǒng)的交換系數(shù)αi,j,k由MODFLOW程序通過(guò)管道形態(tài)與含水層滲透系數(shù)自動(dòng)計(jì)算。通過(guò)數(shù)值模擬后,可得到CFP模型的模擬流場(chǎng)。
采用MNW模型對(duì)二維砂槽管道流進(jìn)行數(shù)值模擬時(shí),在砂槽網(wǎng)格模型第13層插入22個(gè)節(jié)點(diǎn)控制一條管道,管道網(wǎng)格采用孔隙介質(zhì)滲透系數(shù),并將管道外源補(bǔ)排項(xiàng)的值設(shè)為0 m3/s,不考慮表皮效應(yīng)。經(jīng)調(diào)試發(fā)現(xiàn),調(diào)整系數(shù)C與P對(duì)模擬結(jié)果的影響并不顯著,嘗試在多節(jié)點(diǎn)井中額外加入外源流量注入或排泄后,參數(shù)的調(diào)節(jié)對(duì)流場(chǎng)有較大的影響,這表明在該定水頭邊界條件下,管道節(jié)點(diǎn)之間的流量較小,管道中的紊流水頭損失很小。
由于管道與附近含水層水頭差的變化較含水層中水頭差的變化大,一般插值方法如最小二乘法、三角剖分法等得到的流場(chǎng)線不平滑或波動(dòng)極大,插值標(biāo)準(zhǔn)誤差高,與流場(chǎng)的實(shí)際情況差距較大,而經(jīng)測(cè)試后發(fā)現(xiàn),克里金插值法得到的流場(chǎng)線自然平滑,且插值標(biāo)準(zhǔn)誤差低于0.05,優(yōu)于其他插值方法,準(zhǔn)確度較高,因此本文將二維砂槽管道流試驗(yàn)流場(chǎng)中的實(shí)測(cè)水頭值使用Surfer進(jìn)行克里金插值,重新繪制得到水槽中的整體流場(chǎng),并將CFP模型和MNW模型的模擬流場(chǎng)做同樣處理,得到原始MNW模型和CFP模型的模擬流場(chǎng)與實(shí)測(cè)流場(chǎng)的對(duì)比結(jié)果,如圖4所示。
圖4 原始MNW模型和CFP模型模擬流場(chǎng)與實(shí)測(cè)流場(chǎng) 的對(duì)比圖Fig.4 Comparison of simulated flow field using original MNW model and CFP model with experimental flow field
由圖4可知:原始MNW模型的模擬結(jié)果并未表現(xiàn)出含水層中存在雙重介質(zhì)系統(tǒng),管道中水頭從左至右降低,與實(shí)測(cè)流場(chǎng)情況不符;CFP模型中管道表現(xiàn)高滲透性,其模擬流場(chǎng)與實(shí)測(cè)流場(chǎng)趨勢(shì)相同,砂槽內(nèi)水頭中部高,向兩邊減小,而原始MNW模型的模擬結(jié)果與其擬合程度較差。
數(shù)值模型的模擬流場(chǎng)與實(shí)測(cè)流場(chǎng)的誤差,可采用均方根誤差(RMSE)進(jìn)行量化分析,其計(jì)算公式如下:
(13)
式中:M為實(shí)測(cè)水頭點(diǎn)位個(gè)數(shù);hi實(shí)測(cè)值為實(shí)測(cè)點(diǎn)i的水頭[L];hi模擬值為實(shí)測(cè)點(diǎn)i的數(shù)值模擬水頭[L]。
經(jīng)誤差分析可知,原始MNW模型模擬水頭與實(shí)測(cè)水頭的RMSE值為2.516 8 cm,整體誤差水平較大;CFP模型模擬水頭與實(shí)測(cè)水頭的RMSE值為0.925 8 cm,其模擬效果優(yōu)于原始MNW模型。因此,在原始MNW模型的基礎(chǔ)上,利用等效滲透系數(shù)對(duì)模型進(jìn)行修正,即將管道節(jié)點(diǎn)處有限差分網(wǎng)格的滲透系數(shù)由孔隙介質(zhì)滲透系數(shù)修正為等效滲透系數(shù)。在試驗(yàn)中,管道中流量較小,可視作層流,根據(jù)公式(3)計(jì)算得管道層流等效滲透系數(shù)Kavg為12 129 cm/s。采用修正MNW模型對(duì)二維砂槽管道流進(jìn)行了數(shù)值模擬,并將其模擬結(jié)果與CFP模型模擬流場(chǎng)和實(shí)測(cè)流場(chǎng)進(jìn)行了對(duì)比,其結(jié)果如圖5所示。
圖5 修正MNW模型模擬流場(chǎng)與CFP模型模擬流場(chǎng)、 實(shí)測(cè)流場(chǎng)對(duì)比圖Fig.5 Flow field comparison of modified MNW model with CFP model and experimental flow field
由圖5可知:原始MNW模型經(jīng)過(guò)等效滲透系數(shù)法修正后,修正MNW模型很好地表現(xiàn)出管道的存在使含水層產(chǎn)生了非均質(zhì)性,管道與含水層之間的水量交換增加,含水層水頭值抬升,模擬結(jié)果趨近實(shí)測(cè)流場(chǎng);修正MNW模型相較CFP模型的管道表現(xiàn)出更強(qiáng)的滲透性,由于定水頭邊界的影響,管道中水頭處處相同,處于承壓狀態(tài),采用修正MNW模型模擬得到的含水層整體水頭高于CFP模型模擬值。
經(jīng)誤差分析可知:修正MNW模型模擬水頭與實(shí)測(cè)水頭的RMSE值為0.591 4 cm,小于CFP模型模擬水頭與實(shí)測(cè)水頭的RMSE值0.925 8 cm,表明修正MNW模型的模擬效果更好;與未修正MNW模型對(duì)比,模擬水頭與實(shí)測(cè)水頭的RMSE值由2.516 8 cm減小為0.591 4 cm,模擬誤差降低了76.14%,說(shuō)明修正MNW模型的模擬效果得到了大幅度提升。
1) 對(duì)于含水層中的管道,在MODFLOW軟件中使用有限差分法進(jìn)行數(shù)值模擬時(shí),其存在會(huì)影響管道網(wǎng)格滲透系數(shù)的取值,進(jìn)而影響模擬精度。管道中流量較小的情況下,原始MNW模型并不能很好地展現(xiàn)管道的高滲透性質(zhì)。根據(jù)等效滲透系數(shù)法修正MNW模型,對(duì)比CFP模型模擬流場(chǎng),修正MNW模型的模擬結(jié)果優(yōu)于CFP模型;對(duì)比實(shí)測(cè)流場(chǎng),修正MNW模型的模擬結(jié)果趨近于實(shí)測(cè)流場(chǎng),相較原始MNW模型誤差降低了76.14%。
2) 修正MNW模型綜合考慮了高滲透管道的存在對(duì)有限差分網(wǎng)格滲透系數(shù)取值的影響,能夠更好地體現(xiàn)管道與含水層的雙重介質(zhì)性,經(jīng)等效滲透系數(shù)修正后,管道-含水層網(wǎng)格系統(tǒng)的滲透系數(shù)增加,水流在管道中運(yùn)動(dòng)的水頭損失降低,管道與含水層介質(zhì)中的水頭差增加,水流交互更加顯著。若管道內(nèi)水頭高于含水層中水頭,MNW模型在耦合等效滲透系數(shù)后會(huì)進(jìn)一步抬升含水層中水頭,反之則降低。
3) 本研究中管道內(nèi)的水體流態(tài)為層流,理論上MNW模型在刻畫(huà)紊流帶來(lái)的水頭損失上更具優(yōu)勢(shì),該修正MNW模型對(duì)管道流紊流狀態(tài)的模擬效果還需要進(jìn)一步驗(yàn)證。