滲流井輻射管-含水層間水量交換數(shù)值模擬研究

暢俊斌，柯賢敏，王 瑋,葛佳亮，田國(guó)林

(1.陜西地礦九0八水文地質(zhì)工程地質(zhì)大隊(duì)，陜西 西安 710600；2.陜西省礦產(chǎn)資源勘與綜合利用重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710054；3.中煤能源研究院有限責(zé)任公司，陜西 西安 710054；4.長(zhǎng)安大學(xué)環(huán)境科學(xué)與工程學(xué)院，陜西 西安 710054；5.旱區(qū)地下水與生態(tài)效應(yīng)教育部重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710054)

滲流井取水工程是近年來(lái)興起的復(fù)雜取水構(gòu)筑物，利用河床表面覆蓋的砂石層凈化河水，將河水轉(zhuǎn)化成地下水，從而獲得水資源的取水工程[1]。滲流井取水工程因其低耗、高產(chǎn)、易于管理等優(yōu)點(diǎn)[2]被廣泛應(yīng)用。長(zhǎng)期以來(lái)，大多研究者將井管概化為等強(qiáng)度線匯的第二類(lèi)井孔邊界條件[3～6]，由于水流會(huì)產(chǎn)生水頭損失，故井管并不是等水頭邊界。陳崇希等[7]提出“滲流-管流耦合模型”，將輻射管內(nèi)外滲流管流很好地耦合起來(lái)，便于計(jì)算滲流井出水量。Eunhee Lee等[8]利用經(jīng)驗(yàn)系數(shù)刻畫(huà)水平井管并根據(jù)不同雷諾數(shù)下內(nèi)摩擦力的經(jīng)驗(yàn)公式，計(jì)算由于摩擦產(chǎn)生的水頭損失。劉基[9]針對(duì)以往模型的不足，建立改進(jìn)的滲流井計(jì)算模型，計(jì)算多個(gè)滲流井同時(shí)工作時(shí)的出水量。徐嘉璐[10]通過(guò)大型砂槽物理模擬試驗(yàn)，研究大降深條件下滲流井取水造成河水位與地下水位脫節(jié)現(xiàn)象，確定了河水位與地下水位脫節(jié)的臨界點(diǎn)及疏干單元的變化規(guī)律。王瑋等[11～12]針對(duì)滲流井取水問(wèn)題，在“等效滲透系數(shù)”和“滲流—管流耦合模型”等理論基礎(chǔ)上，建立滲流井取水計(jì)算模型，并成功將該計(jì)算模型用于悖牛川水源地滲流井取水工程地下水允許開(kāi)采量計(jì)算[2]。徐嘉璐等[13]以輻射管與含水層間交換水量為耦合點(diǎn)，進(jìn)一步建立了改進(jìn)滲流井取水計(jì)算模型，使得計(jì)算模型能夠更精細(xì)地刻畫(huà)輻射管特征，并實(shí)現(xiàn)多個(gè)滲流井同時(shí)取水的計(jì)算。之后改進(jìn)的滲流井計(jì)算模型成功應(yīng)用于陜北多個(gè)水源地計(jì)算中[14～15]。

以往滲流井計(jì)算模型中，輻射管-含水層間的水量交換系數(shù)大多通過(guò)反演得到，部分采用經(jīng)驗(yàn)系數(shù)確定，未有準(zhǔn)確的計(jì)算公式對(duì)其進(jìn)行定量計(jì)算。輻射管與含水層間的交換量是滲流井出水量的唯一來(lái)源，是滲流井取水工程的重要環(huán)節(jié)，為此研究輻射管與含水層間的交換量具有重要意義。本文通過(guò)大型砂槽物理模擬試驗(yàn)，對(duì)滲流井取水的影響因素進(jìn)行分析，推導(dǎo)出輻射管-含水層間水量交換系數(shù)的計(jì)算公式，并對(duì)部分試驗(yàn)方案建立數(shù)值模型，利用試驗(yàn)數(shù)據(jù)對(duì)數(shù)值模型可靠性及輻射管-含水層間水量交換系數(shù)計(jì)算公式的正確性進(jìn)行驗(yàn)證。

1 試驗(yàn)方法與過(guò)程

1.1 試驗(yàn)裝置

滲流砂槽長(zhǎng)3.63 m，寬0.5 m，高1.17 m，砂槽頂端放置水箱，用于試驗(yàn)前期飽水及后期向砂槽內(nèi)供水，砂槽左側(cè)設(shè)置豎井，用于控制滲流井出水量。砂槽前后壁共設(shè)置18×10個(gè)水位觀測(cè)孔，并插入不同長(zhǎng)度的銅管(銅管另一端與測(cè)壓板相連)，用于觀測(cè)試驗(yàn)過(guò)程中砂槽內(nèi)部不同位置的水位。試驗(yàn)裝置示意圖見(jiàn)圖1。

圖1 試驗(yàn)裝置示意圖Fig.1 Design of the test equipment

1.2 試驗(yàn)方案設(shè)計(jì)

本次設(shè)計(jì)輻射管仰角為：0°、10°、13°(設(shè)計(jì)仰角為15°，由于工藝原因?qū)嶋H為13°)和30°，輻射管長(zhǎng)度為：1 m、2 m、3 m。設(shè)計(jì)彩條布(由于彩條布隔水效果過(guò)強(qiáng)，本次沿x方向每隔40 cm，設(shè)置一列直徑為0.5 cm的圓孔，圓孔間距為10 cm)和紗網(wǎng)(本次試驗(yàn)采用篩孔尺寸為1 mm的紗網(wǎng))分別代表滲漏性能弱和強(qiáng)的河床類(lèi)型。試驗(yàn)前期對(duì)上述部分方案進(jìn)行了小、中、大3次豎井降深試驗(yàn)，后期經(jīng)過(guò)改進(jìn)對(duì)上述大多數(shù)方案進(jìn)行自小而大5次豎井降深的試驗(yàn)。上述各條件互相組合，確定為設(shè)計(jì)試驗(yàn)方案。

1.3 試驗(yàn)過(guò)程

為了使輻射管內(nèi)外呈現(xiàn)多種流態(tài)及方便試驗(yàn)裝置的更換，本次砂槽物理模擬試驗(yàn)選擇中粗砂(粒徑為0.25～1.0 mm)為含水層介質(zhì)。為了減少試驗(yàn)誤差，對(duì)試驗(yàn)用砂進(jìn)行水洗，去除雜質(zhì)，采用水沉固結(jié)，分層裝砂，以滿足含水層均質(zhì)的要求。試驗(yàn)結(jié)束后沿x軸(0～363 cm)進(jìn)行6組豎管試驗(yàn)，測(cè)得豎直滲透系數(shù)依次為：28.81，37.81，38.89，40.07，46.73，33.50 m/d，該滲透系數(shù)可以為模型滲透系數(shù)(Kz)提供參考依據(jù)。

為了充分排除砂槽內(nèi)空氣，在砂槽底部裝入一層卵石層，厚度約為8 cm，以便在飽水過(guò)程中，水能從砂槽底部自下而上進(jìn)入含水層。裝砂達(dá)到一定高度時(shí)，將一端用紗布包裹的銅管(長(zhǎng)度：0.05 m、0.15 m、0.25 m)、滲流井微縮模型及測(cè)壓探頭進(jìn)行安裝，最終裝砂總高度為0.97 m。裝砂結(jié)束后，在砂層上布設(shè)彩條布(或尼龍紗網(wǎng))，并在彩條布(或尼龍紗網(wǎng))上鋪一層均勻的卵石層，模擬河床介質(zhì)。

飽水完成后，用洗耳球排除測(cè)壓孔中的氣泡，打開(kāi)供水水箱的開(kāi)關(guān)，使水流入“河流”中，當(dāng)“河流”水位達(dá)到溢流高度時(shí)，記錄測(cè)壓孔、豎井、三角堰的初始水位。然后打開(kāi)與豎井相連水龍頭，調(diào)整與水龍頭相連的橡膠管至第一次降深試驗(yàn)的豎井高度，記錄此時(shí)的時(shí)間，即為試驗(yàn)初始時(shí)刻，待壓力傳感器上數(shù)據(jù)達(dá)到穩(wěn)定后，記錄三角堰的堰高、豎井的水位、并測(cè)量三角堰流量。記錄頻率為10 min/次，共記錄5次，直到三角堰水位和豎井水位均穩(wěn)定不變。重復(fù)上述實(shí)驗(yàn)過(guò)程，完成其他降深試驗(yàn)。

2 試驗(yàn)結(jié)果分析

由于試驗(yàn)方案和影響因素均較多，故采用逐步優(yōu)選的方法對(duì)各影響因素進(jìn)行分析，最終優(yōu)選出取水效率最高的試驗(yàn)方案。

2.1 輻射管仰角對(duì)滲流井出水量影響

不同輻射管仰角條件下，滲流井出水量隨豎井降深變化見(jiàn)圖2。由圖2可以看出，不論在強(qiáng)滲漏補(bǔ)給能力河流條件下還是在弱滲漏補(bǔ)給能力河流條件下，試驗(yàn)出水量隨著豎井降深的增大均不斷增大。當(dāng)豎井降深增大時(shí)，導(dǎo)致輻射管內(nèi)外水頭差增大，含水層中水不斷向輻射管流動(dòng)，從而出水量增大。當(dāng)其他取水條件相同時(shí)，輻射管仰角大的試驗(yàn)方案出水量大于輻射管仰角小的試驗(yàn)方案出水量。這是由于隨著輻射管仰角的增大，相同長(zhǎng)度輻射管在水平方向上投影長(zhǎng)度相差不大，但在豎直方向上的投影長(zhǎng)度越大，輻射管影響面積增大。同時(shí)，輻射管仰角越大，輻射管末端距河流越近，容易接受河流滲漏補(bǔ)給，從而滲流井出水量大。

圖2 不同輻射管仰角條件下滲流井出水量隨豎井降深變化圖Fig.2 Water yield of the seepage well with the drawdown of vertical shaft under different elevation angles of the radiant tube

鑒于試驗(yàn)裝置限制，輻射管仰角為30°方案只進(jìn)行了1 m試驗(yàn)。故對(duì)取水效率較高的輻射管仰角為13°的方案進(jìn)行分析，研究該方案下輻射管長(zhǎng)度及河床類(lèi)型對(duì)滲流井出水量的影響。

2.2 輻射管長(zhǎng)度對(duì)滲流井出水量影響

在出水效果較好的13°方案基礎(chǔ)上，分析輻射管長(zhǎng)度對(duì)滲流井取水效果的影響，結(jié)果見(jiàn)圖3。

由圖3可知，隨著豎井降深增大，試驗(yàn)出水量也不斷增大。不同滲漏補(bǔ)給作用下，當(dāng)豎井降深相同時(shí)，隨著輻射管長(zhǎng)度的增加，滲流井出水量也不斷增大。輻射管長(zhǎng)度越長(zhǎng)，輻射管與含水層接觸面積就越大，且在輻射管仰角存在的情況下，輻射管末端距河更近，更容易接受河流的滲漏補(bǔ)給，在豎直方向也增大了與含水層的接觸面積，從而出水量較大?？梢?jiàn)本次試驗(yàn)方案中3 m 13°方案優(yōu)于其他長(zhǎng)度取水方案。

2.3 河床類(lèi)型對(duì)取水構(gòu)筑物取水效果的影響分析

在輻射管仰角及輻射管長(zhǎng)度研究的基礎(chǔ)上，選擇出水效果較好的3 m 13°方案，進(jìn)行河床類(lèi)型對(duì)滲流井取水效果的影響分析，結(jié)果見(jiàn)圖4。

圖4 不同河流類(lèi)型作用下試驗(yàn)出水量隨豎井降深變化圖Fig.4 Water yield of the seepage well with the drawdown of vertical shaft in different rivers

由圖4可以看出，在輻射管仰角和輻射管長(zhǎng)度相同條件下，隨著豎井降深的增大，試驗(yàn)出水量不斷增大，且強(qiáng)滲漏補(bǔ)給能力河流方案的試驗(yàn)出水量大于弱滲漏補(bǔ)給能力河流方案的試驗(yàn)出水量。兩種試驗(yàn)方案下，僅因?yàn)楹哟步橘|(zhì)不同即河床滲透系數(shù)不同，獲得河流補(bǔ)給量不同。河流滲漏補(bǔ)給量為：

(1)

式中：Q——河流滲漏補(bǔ)給量；

Kz——河床垂向滲透系數(shù)；

M——河床介質(zhì)厚度；

A——河流發(fā)生滲漏的面積；

ΔH——河流與含水層間的水頭差。

由式(1)可以看出，兩種河流試驗(yàn)方案中，除了河床的滲透系數(shù)不同外，其他各項(xiàng)均相等?？梢?jiàn)河床滲漏補(bǔ)給量受河床的滲透系數(shù)影響，河床滲透系數(shù)越大，河流滲漏補(bǔ)給量越大，滲流井出水量越大。

3 數(shù)值模擬

3.1 輻射管-含水層交換系數(shù)推導(dǎo)

輻射管與含水層的交換量計(jì)算公式可用達(dá)西定律表示，水頭可采用迭代法求解，故而問(wèn)題轉(zhuǎn)化為確定輻射管-含水層間的交換系數(shù)C：

Qex=Cj,i,k(hin-hj,i,k)

(2)

式中:Qex——輻射孔與含水層間交換量；

hin——節(jié)點(diǎn)in的水頭；

hj,i,k——單元格(j,i,k)的水頭；

Cj,i,k——輻射管-含水層間的水量交換系數(shù)。

Shoemaker[16]提出均質(zhì)各向同性介質(zhì)條件下輻射管與含水層間水量交換系數(shù)計(jì)算公式：

(3)

式中：Kj,i,k——輻射管所在單元格(j,i,k)的滲透系數(shù)；

dip——與節(jié)點(diǎn)in相連的管道ip的直徑；

Δlip——與管道ip相連的兩個(gè)管道節(jié)點(diǎn)的直線距離；

τip——管道ip相鄰節(jié)點(diǎn)in的曲率；

rip——輻射管的半徑。

式(3)中Cj,i,k值計(jì)算公式只適用于均質(zhì)各向同性介質(zhì)且缺少試驗(yàn)驗(yàn)證，其正確性有待考證。實(shí)際工程中，含水層介質(zhì)大多不是均質(zhì)各向同性，故均質(zhì)各向異性條件下的C1值的計(jì)算對(duì)指導(dǎo)滲流井的設(shè)計(jì)、模型的建立和計(jì)算具有重要意義。本文在式(3)的基礎(chǔ)上，先分別根據(jù)Kx、Ky和Kz計(jì)算出Cx、Cy、Cz。然后得到Cx、Cy、Cz垂直輻射管的分量，然后將各分量相乘開(kāi)方：

(4)

式中：α——輻射管仰角；

n——與輻射管節(jié)點(diǎn)i相連的節(jié)點(diǎn)總個(gè)數(shù)。

以往的輻射管-含水層間的交換系數(shù)均是通過(guò)試算法、經(jīng)驗(yàn)公式及參數(shù)反演得出，未有定量的計(jì)算公式。本文推導(dǎo)出C1值計(jì)算公式，并建立數(shù)值模型，將計(jì)算公式應(yīng)用于模型中，通過(guò)試驗(yàn)數(shù)據(jù)對(duì)模型及計(jì)算公式進(jìn)行驗(yàn)證。

3.2 數(shù)值模型

3.2.1概念模型

本次以整個(gè)砂槽為模擬區(qū)，即X：0～3.63 m；Y：0～0.5 m；Z：0～1.13 m。將含水層介質(zhì)概化成均質(zhì)軸對(duì)稱(chēng)各向異性多孔介質(zhì)；由于強(qiáng)滲漏補(bǔ)給能力河流(紗網(wǎng)方案)對(duì)水流下滲阻礙作用很小，將其概化為定水頭邊界；弱滲漏補(bǔ)給能力河流(彩條布方案)對(duì)水流下滲阻礙作用比較明顯，將其概化為河流邊界。砂槽四壁及底部為鋼板，將其概化為隔水邊界。試驗(yàn)過(guò)程中，穩(wěn)定時(shí)刻的豎井水位保持不變，將其概化為定水頭，由于豎井位于砂槽外部，故將豎井與平巷相連處單元格設(shè)置為定水頭。試驗(yàn)開(kāi)始后，砂槽內(nèi)含水層只接受河流補(bǔ)給，豎井水位下降后，地下水從各個(gè)方向向輻射管匯集，通過(guò)與豎井相連的平巷流入三角堰，整個(gè)過(guò)程中，地下水流動(dòng)形態(tài)呈現(xiàn)顯著的三維流特征，水流運(yùn)動(dòng)規(guī)律服從達(dá)西定律。

3.2.2數(shù)學(xué)模型

在上述概念模型分析的基礎(chǔ)上，建立本次砂槽物理模擬試驗(yàn)地下水三維穩(wěn)定流數(shù)學(xué)模型：

式中：H——地下水位標(biāo)高/m；

x、y、z——坐標(biāo)變量/m；

K——含水層滲透系數(shù)/(m·d-1)；

D——模擬區(qū)的范圍；

Γ2——第二類(lèi)邊界；

n——隔水邊界的外法線方向；

np——潛水面邊界的內(nèi)法線方向；

C1——輻射管與含水層間交換系數(shù)/(m2·d-1)；

Qe——輻射孔-含水層間的交換量/(m3·d-1)；

Qp——輻射管內(nèi)出水量/(m3·d-1)；

Hr——河水位/m；

ΔH——水頭損失/m；

Hp——輻射管內(nèi)水位/m；

ν——水流運(yùn)動(dòng)黏滯系數(shù)/(m2·s-1)；

d——輻射管的直徑/m；

v——滲流速度/(m·s-1)；

Kr——河床介質(zhì)滲透系數(shù)；

Mr——河床介質(zhì)厚度/m；

qr——河流滲漏補(bǔ)給量/(m3·d-1)。

管流計(jì)算公式中包含有水頭、流速、流量等多個(gè)變量，可采用迭代法進(jìn)行求解。可通過(guò)把含水層與滲流井間的交換量作為匯項(xiàng)把該方程耦合進(jìn)MODFLOW。計(jì)算時(shí)，先把各管道節(jié)點(diǎn)的頂部標(biāo)高作為滲流井井管內(nèi)水頭迭代初值，并假設(shè)初始流態(tài)為層流，在此基礎(chǔ)上，結(jié)合MODFLOW中含水層水頭的迭代初值，計(jì)算出含水層與滲流井間的交換量，進(jìn)而可計(jì)算出井管內(nèi)各段的流量變化。根據(jù)這個(gè)流量及井管直徑，可計(jì)算井管內(nèi)水流的雷諾數(shù)，據(jù)此選擇不同的管流公式可以計(jì)算出新的管內(nèi)水頭分布，再利用新的管內(nèi)水頭值重復(fù)上述計(jì)算，直到連續(xù)兩次迭代的井管內(nèi)水頭值絕對(duì)誤差滿足精度為止。此時(shí)含水層與水平井間交換量可作為匯項(xiàng)加入到滲流方程中，采用MODFLOW中的迭代可求解出含水層水頭，再將迭代出的含水層水頭作為MODFLOW中含水層水頭的迭代初值，重復(fù)上述迭代計(jì)算過(guò)程，直至連續(xù)兩次迭代的含水層水頭值絕對(duì)誤差滿足精度要求為止。再利用輻射管與含水層間水量交換系數(shù)、最終確定的含水層水頭及井管內(nèi)的水頭得出滲流井出水量。其中，MODFLOW軟件通過(guò)CFP文件將模型與離散的管道耦合起來(lái)，井管由大量管道及節(jié)點(diǎn)組成，通過(guò)設(shè)置每個(gè)管道節(jié)點(diǎn)的高程及管道參數(shù)(輻射管與含水層間水量交換系數(shù)、彎曲度、粗糙度、層流與紊流相互轉(zhuǎn)換的兩個(gè)臨界雷諾數(shù)、管徑等)并將其編號(hào)來(lái)模擬實(shí)際管道。

3.2.3數(shù)值模型

采用有限差分法，對(duì)模擬區(qū)進(jìn)行網(wǎng)格剖分，垂向上砂槽頂部的“河流”占一層、底部的卵石層占2層，河流與卵石層厚度均為8 cm，中間含水層共分28層。整個(gè)模擬區(qū)沿X軸方向被剖分為92列，沿Y軸方向被剖分為12行，沿Z軸方向被剖分為31層(第25層為平巷所在層)，單層活動(dòng)單元共1 104個(gè)，總活動(dòng)單元數(shù)為34 224個(gè)。在改進(jìn)的滲流井計(jì)算模型基礎(chǔ)上，以輻射管-含水層間交換量為耦合點(diǎn)，建立“滲流-管流耦合模型”，利用穩(wěn)定狀態(tài)下測(cè)得的滲流井出水量、豎井降深及觀測(cè)孔實(shí)測(cè)降深與模型計(jì)算結(jié)果進(jìn)行擬合，驗(yàn)證模型可靠性。

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

本次對(duì)3 m 13°兩種河流方案進(jìn)行數(shù)值模擬，模擬結(jié)果見(jiàn)表1。

由表1可以看出，滲流井出水量的實(shí)測(cè)值與模型計(jì)算值近似一致，強(qiáng)滲漏補(bǔ)給河流方案實(shí)測(cè)流量與模型計(jì)算流量相對(duì)誤差差值范圍在-6.36%～4.88%。弱滲漏補(bǔ)給河流方案實(shí)測(cè)流量與模型計(jì)算流量的相對(duì)誤差值范圍在-4.56%～2.42%。由圖5可知，模型計(jì)算結(jié)果與實(shí)測(cè)結(jié)果的“流量-降深”曲線基本重合，擬合良好。且隨著豎井降深的增加滲流井出水量也增大，這與試驗(yàn)得到規(guī)律一致。

表1 流量擬合結(jié)果表Table 1 Fitting results of the flow

注：表1中D1-D5分別表示第一次—第五次降深試驗(yàn)。

圖5 強(qiáng)滲漏和弱滲漏補(bǔ)給河流出水量擬合曲線Fig.5 The fitting curve of strong and weak leakage capacity river

為進(jìn)一步驗(yàn)證所建模型的合理性，本次對(duì)測(cè)壓板觀測(cè)孔水位降深與模型計(jì)算降深也進(jìn)行擬合(以第一、三、五次降深為例繪圖)，擬合結(jié)果見(jiàn)圖6。

由圖6可以看出，除個(gè)別觀測(cè)點(diǎn)外，各觀測(cè)孔計(jì)算水位降深與實(shí)測(cè)降深組合的點(diǎn)均位于y=x+0.03與y=x-0.03之間，說(shuō)明模型計(jì)算降深值與實(shí)測(cè)降深值的誤差在-0.03～0.03 m之間。對(duì)于小降深條件下，誤差僅為0.005 m。且從圖中還可看出，隨著豎井降深的增大，砂槽內(nèi)各觀測(cè)點(diǎn)的降深也變大，可見(jiàn)當(dāng)豎井降深增大時(shí)，含水層中的水均向輻射管流動(dòng)，從而導(dǎo)致各觀測(cè)點(diǎn)水位降深增大，當(dāng)豎井降深增大時(shí)，輻射管的影響面積也變大。

圖6 強(qiáng)滲漏和弱滲漏補(bǔ)給河流測(cè)壓點(diǎn)水位降深擬合結(jié)果Fig.6 Fitting results of the pressure observation drawdown in strong and weak leakage capacity rivers

由表1及圖5可以看出實(shí)測(cè)流量與模型計(jì)算流量近于一致，由圖6可以看出，各觀測(cè)孔模型計(jì)算降深與試驗(yàn)實(shí)測(cè)降深誤差很小，滿足模型精度要求。由此可知，試驗(yàn)實(shí)測(cè)結(jié)果與模型計(jì)算結(jié)果擬合良好，說(shuō)明本次所建模型是合理的，輻射管-含水層間水量交換系數(shù)計(jì)算公式是可行的。

4 結(jié)論

(1)其他條件不變，隨著豎井降深的增大，輻射管仰角的增大，輻射管長(zhǎng)度的增加，河流滲漏補(bǔ)給能力的增強(qiáng)，滲流井出水量均不斷增大。

(2)在本次物理模擬試驗(yàn)中，3 m 13°強(qiáng)滲漏補(bǔ)給能力河流方案為取水效果最佳方案。

(3)在以往研究的基礎(chǔ)上，推導(dǎo)出了適用于均質(zhì)軸對(duì)稱(chēng)各項(xiàng)異性介質(zhì)條件下的輻射管-含水層水量交換系數(shù)的計(jì)算公式。