表皮效應(yīng)對(duì)井筒內(nèi)地下水流速影響數(shù)值模擬研究

汪啟航，李 旭,*，郭 強(qiáng)，蒯大偉

(1.煤炭開采水資源保護(hù)與利用全國重點(diǎn)實(shí)驗(yàn)室，北京 100011；2.安徽理工大學(xué) 地球與環(huán)境學(xué)院，安徽 淮南 232001；3.安徽省煤田地質(zhì)局第二勘探隊(duì)，安徽 蕪湖 241000)

近年來，我國地下水資源與環(huán)境問題日益嚴(yán)峻，對(duì)工農(nóng)業(yè)生產(chǎn)用水、居民日常用水安全以及礦山安全開采等都造成了嚴(yán)重的影響[1-4]，解決地下水資源與環(huán)境問題的有效途徑是對(duì)地下水資源進(jìn)行科學(xué)管理和環(huán)境評(píng)價(jià)[5-7]。而準(zhǔn)確獲得地下水流速，對(duì)于區(qū)域水資源評(píng)價(jià)、污染物運(yùn)移預(yù)測(cè)以及礦山水害治理等都具有重要的現(xiàn)實(shí)意義[8-9]。目前，比較常見的測(cè)定地下水流速的方法有單孔示蹤試驗(yàn)法[10-11]、測(cè)溫法[12]以及Aqua Vision流速流向儀測(cè)定法[13]等，其中Aqua Vision 流速流向儀測(cè)定地下水流速是將設(shè)備放入井筒，采用視頻顯微影像方法，實(shí)時(shí)拍攝溶于水中的膠體顆粒物以測(cè)量地下水的流速，同時(shí)融合了高分辨率磁通量羅盤和高放大率膠質(zhì)顆粒追蹤影像系統(tǒng)，來測(cè)定地下水的流向。例如，郭緒磊等[13]利用Aqua Vision 流速流向儀對(duì)位于湖北興山香溪河流域的ZK08 水文地質(zhì)鉆孔不同深度的含水層進(jìn)行地下水流速流向的測(cè)量。

實(shí)際上，井筒內(nèi)的地下水流速與含水層實(shí)際地下水流速存在一定差異，誤差的產(chǎn)生是受到了井筒附近表皮效應(yīng)的影響，導(dǎo)致井筒附近的地下水流場(chǎng)發(fā)生變化，造成井筒內(nèi)的流速與含水層中的流速之間存在偏差[14]。在鉆井施工工藝過程中受到井壁和鉆井泥漿侵入、濾料填充以及生物化學(xué)堵塞等因素的影響，導(dǎo)致抽水井區(qū)域滲透系數(shù)、孔隙率等水力性質(zhì)發(fā)生變化，這種現(xiàn)象被稱為表皮效應(yīng)(skin effect)[15-18]。根據(jù)表皮與含水層滲透系數(shù)之間差異，可將表皮分為正、負(fù)兩種類型。若表皮滲透系數(shù)小于含水層滲透系數(shù)，則定義為正表皮，否則為負(fù)表皮[19-22]。此外，表皮區(qū)的厚度通常從幾毫米到幾米不等，其厚度的差異也會(huì)對(duì)井筒附近的地下水流場(chǎng)產(chǎn)生較大影響[23]。為了解決井筒附近表皮結(jié)構(gòu)對(duì)井筒地下水流速的影響，M.Schubert 等[24]提出了利用井筒內(nèi)的地下水流速換算含水層地下水達(dá)西流速經(jīng)驗(yàn)方程(后文簡(jiǎn)稱經(jīng)驗(yàn)公式)。由于表皮結(jié)構(gòu)較為復(fù)雜性，且目前對(duì)經(jīng)驗(yàn)公式的研究還不夠完善，尤其是表皮厚度、滲透性等特征的變化對(duì)經(jīng)驗(yàn)公式計(jì)算精度的影響，還有待進(jìn)一步驗(yàn)證。因此，有必要開展井筒附近表皮特征(滲透性、厚度等)對(duì)井筒內(nèi)地下水流速影響機(jī)制的研究，進(jìn)一步評(píng)估M.Schubert 等提出的經(jīng)驗(yàn)公式的可靠性。

為此，筆者通過COMSOL Multiphysics 軟件建立考慮表皮效應(yīng)影響下的地下水流數(shù)值模型，包括不同表皮區(qū)滲透系數(shù)、表皮厚度以及井筒半徑等因素，分析上述因素對(duì)井筒內(nèi)地下水流速的影響，主要開展以下3 個(gè)方面研究：(1)井筒半徑對(duì)井筒內(nèi)地下水流速的影響；(2)表皮區(qū)域范圍和滲透系數(shù)對(duì)井筒內(nèi)地下水流速的影響；(3)根據(jù)模擬結(jié)果進(jìn)一步驗(yàn)證利用井筒內(nèi)地下水流速換算含水層地下水達(dá)西流速數(shù)學(xué)模型的可靠性，以期通過定量分析表皮區(qū)滲透系數(shù)、厚度等因素對(duì)經(jīng)驗(yàn)公式計(jì)算誤差的影響，為實(shí)際含水層流速精準(zhǔn)計(jì)算提供理論依據(jù)。

1 數(shù)學(xué)模型的建立與求解

1.1 數(shù)學(xué)模型建立

為了簡(jiǎn)化數(shù)學(xué)模型，本文提出如下假設(shè)條件：(1)承壓含水層為等厚且水平無限延伸，其厚度為B；(2)地下水流為穩(wěn)定流；(3)以水井中心為原點(diǎn)建立坐標(biāo)系，水井為完整井且半徑為rw；(4)水井附近的表皮和含水層為均質(zhì)且各向同性；(5)表皮為均勻表皮，表皮厚度rh=rs-rw，rs為表皮區(qū)域的半徑，表皮的外部區(qū)域?yàn)槌袎汉畬?，如圖1 所示；(6)模型中含水層地下水為三維流動(dòng)，故地下水滲流方程可以表示為：

圖1 表皮結(jié)構(gòu)對(duì)地下水流場(chǎng)影響Fig.1 Influence of skin structure on groundwater flow field

表皮區(qū)域地下水滲流方程可以表示為：

為建立一個(gè)穩(wěn)定的地下水流場(chǎng)，模型中承壓含水層的左邊界和右邊界均為定水頭邊界，分別為：

1.2 數(shù)值模型建立

本文研究表皮區(qū)域范圍和滲透系數(shù)對(duì)地下水流速的影響規(guī)律，為準(zhǔn)確獲取地下水流速提供必要的理論和模型支撐，通過COMSOL Multiphysics 軟件建立考慮表皮區(qū)域的三維承壓含水層結(jié)構(gòu)，由于井結(jié)構(gòu)為完整井，因此，模型中地下水流為水平二維流動(dòng)。利用該模型來研究表皮區(qū)范圍和滲透系數(shù)對(duì)地下水流速的影響。在該模型中井筒取無窮大的滲透系數(shù)來代替井筒內(nèi)滲透系數(shù)[25]，本數(shù)值模型的范圍是40 m×40 m×5 m，模型網(wǎng)格采用自由四面體網(wǎng)格，網(wǎng)格包含392 231 個(gè)域單元、20 534 個(gè)邊界元和772 個(gè)邊單元，構(gòu)建模型三維網(wǎng)格如圖2 所示。

圖2 模型三維網(wǎng)格剖分(細(xì)化)Fig.2 3D mesh generation (refinement) of the model

1.3 井筒內(nèi)地下水流速計(jì)算方法

圖3a 為井筒附近表皮結(jié)構(gòu)示意圖，中心白色區(qū)域?yàn)榫?，藍(lán)色區(qū)域?yàn)楸砥^(qū)域，K1、K2分別為含水層和表皮的滲透系數(shù)。在實(shí)際井筒地下水流速測(cè)定過程中，主要測(cè)定的是井筒內(nèi)的平均流速。由于該數(shù)值模型中地下水流場(chǎng)為穩(wěn)定流，所以含水層中地下水達(dá)西流速可以用達(dá)西定律求得；井筒內(nèi)地下水流速vj選擇深度為2.5 m時(shí)，從井筒中心點(diǎn)向四周等距離取值后求平均獲得，取井筒內(nèi)具有代表性的13 個(gè)地下水流速數(shù)值，求取平均值代替井筒內(nèi)地下水流速，取點(diǎn)如圖3b 所示，則井筒內(nèi)地下水流速vj可以由下式計(jì)算得出。

圖3 表皮結(jié)構(gòu)及井筒內(nèi)代表性流速點(diǎn)Fig.3 Skin structure and representative flow velocity points in wellbore

2 含水層地下水流速計(jì)算方法

為了能夠確定含水層的地下水達(dá)西流速與井筒內(nèi)地下水流速之間的關(guān)系，M.Schubert 等提出了在表皮結(jié)構(gòu)影響下井筒內(nèi)地下水流速和含水層地下水達(dá)西流速之間的定量相關(guān)關(guān)系，相關(guān)公式如下：

3 井筒內(nèi)地下水流速影響因素分析

本節(jié)利用構(gòu)建的數(shù)值模型分析井徑、表皮區(qū)域的滲透系數(shù)以及表皮厚度對(duì)井筒內(nèi)地下水流速的影響，同時(shí)對(duì)比經(jīng)驗(yàn)公式計(jì)算結(jié)果，進(jìn)一步驗(yàn)證經(jīng)驗(yàn)公式計(jì)算精度。對(duì)于以下分析，默認(rèn)的模型參數(shù)見表1。

表1 數(shù)值模型參數(shù)Table 1 Numerical model parameter

3.1 井 徑

為了驗(yàn)證經(jīng)驗(yàn)公式在無表皮影響下的可靠性，分析不同井徑對(duì)計(jì)算精度的影響，則設(shè)置模型參數(shù)為：K1=1.157 4×10-5m/s、rw=0.10～0.60 m(每隔0.05 m 取一個(gè)值)，rh=0 m。利用式(7)計(jì)算含水層地下水達(dá)西流速，即vw=K1J=2.893 5×10-7m/s。由式(8)-式(9)可知，在無表皮狀態(tài)下井筒內(nèi)地下水流速為含水層地下水達(dá)西流速的2 倍，故利用式(9)可以反算井筒內(nèi)地下水流速為5.787 0×10-7m/s。圖4 為利用經(jīng)驗(yàn)公式和數(shù)值模擬計(jì)算的不同井徑條件下井筒內(nèi)地下水流速。由圖4可知，數(shù)值模擬計(jì)算的井筒內(nèi)地下水流速與經(jīng)驗(yàn)公式計(jì)算的結(jié)果相差不大，說明在不考慮表皮效應(yīng)時(shí)，經(jīng)驗(yàn)公式是可靠的，且井筒內(nèi)地下水流速近似于含水層地下水達(dá)西流速的2 倍。

圖4 不同井徑條件下井筒內(nèi)地下水流速經(jīng)驗(yàn)公式和數(shù)值模擬計(jì)算結(jié)果對(duì)比Fig.4 Comparison of calculated groundwater flow velocity in the wellbore with different diameters by empirical formula and numerical simulation

同時(shí)，為了進(jìn)一步評(píng)估計(jì)算誤差，本文定義一個(gè)相對(duì)誤差參數(shù)Ev=|vgw-vw|/vw×100%，來評(píng)估M.Schubert等提出的經(jīng)驗(yàn)公式的可靠性，誤差計(jì)算結(jié)果見表2。由表2 可知，不同井徑條件下含水層地下水達(dá)西流速與換算后地下水達(dá)西流速誤差很小，其相對(duì)誤差范圍為0.01%～0.47%，說明經(jīng)驗(yàn)公式可以較好地刻畫井筒內(nèi)地下水流速與含水層地下水達(dá)西流速之間的關(guān)系。僅在井徑較小時(shí)，誤差有所增加，最大相對(duì)誤差為0.47%，因此，在實(shí)際應(yīng)用中由于井徑變化帶來的影響較小?？傊?，在無表皮影響下，式(8)-式(9)計(jì)算精度較高。

表2 地下水流速計(jì)算數(shù)據(jù)Table 2 Calculation data of groundwater flow velocity

3.2 表皮滲透系數(shù)

由于表皮效應(yīng)的存在，井筒附近區(qū)域水文地質(zhì)參數(shù)發(fā)生變化，必然對(duì)井筒內(nèi)的地下水流速產(chǎn)生影響。因此，有必要進(jìn)一步研究滲透系數(shù)對(duì)井筒內(nèi)地下水流速的影響機(jī)制。為了方便討論表皮滲透系數(shù)的影響，定義表皮指數(shù)δ=K2/K1，反映表皮與含水層之間滲透性差異[25]。主要參數(shù)設(shè)置如下：rw=0.1 m，rs=0.3 m，K1、K2取值見表1。因此，δ的取值范圍為0.5～1.4，若δ＜ 1 則表示為正表皮，δ＞ 1 則表示為負(fù)表皮，當(dāng)δ=1 表示為無表皮[26]。圖5 表示表皮區(qū)域滲透系數(shù)與井筒內(nèi)地下水流速之間的關(guān)系。結(jié)果表明，隨著表皮滲透系數(shù)的增大，井筒內(nèi)地下水流速逐漸增大。根據(jù)達(dá)西定律可知，水力梯度不變的條件下，滲透系數(shù)增大，達(dá)西流速增大。因此，表皮區(qū)域滲透系數(shù)增大，則井筒附近的高滲區(qū)域增大，從而導(dǎo)致井筒內(nèi)地下水流速增大；相反，如果表皮區(qū)域滲透系數(shù)減小，那么相當(dāng)于在井筒附近形成低滲區(qū)域，從而導(dǎo)致井筒內(nèi)地下水流速減小。

圖5 井筒內(nèi)流速隨表皮指數(shù)變化Fig.5 Variation of flow velocity in wellbore with skin index

3.3 表皮厚度

圖6 為不同表皮條件下井筒內(nèi)流速隨表皮厚度變化。在正表皮狀態(tài)下，δ取 0.6，負(fù)表皮狀態(tài)下，δ取1.3；其他參數(shù)相同，rw=0.2 m，rh=0.1、0.2、0.3、0.4、0.5 m。由圖6 可知，在正表皮狀態(tài)下，隨著表皮厚度的增大，井筒內(nèi)地下水流速逐漸減小。由達(dá)西定律可知，在水力梯度不變的條件下，滲透系數(shù)減小，則達(dá)西流速減小。因此，處于正表皮狀態(tài)時(shí)，隨著表皮區(qū)域厚度增大，則低滲區(qū)域增大，井筒附近滲透系數(shù)減小，從而導(dǎo)致井筒內(nèi)地下水流速減小。處于負(fù)表皮狀態(tài)時(shí)，恰好與正表皮情況相反，隨著負(fù)表皮厚度的增大，井筒內(nèi)地下水流速逐漸增大。這主要是由于當(dāng)δ＞1 時(shí)，表示負(fù)表皮區(qū)域滲透系數(shù)大于含水層滲透系數(shù)，表皮區(qū)域相當(dāng)于一個(gè)高滲區(qū)域，隨著高滲區(qū)域增大，井筒附近滲透系數(shù)增大，從而引起井筒內(nèi)地下水流速不斷增大。

圖6 不同表皮條件下井筒內(nèi)流速變化Fig.6 Variation of flow velocity in wellbore under different skin conditions

4 經(jīng)驗(yàn)公式計(jì)算誤差分析

4.1 表皮滲透性

為了進(jìn)一步評(píng)估表皮區(qū)域滲透性對(duì)經(jīng)驗(yàn)公式計(jì)算精度的影響，本節(jié)通過改變表皮區(qū)域滲透系數(shù)，利用式(8)-式(9)計(jì)算并分析其計(jì)算結(jié)果的準(zhǔn)確性。模型的主要參數(shù)設(shè)置為：rw=0.1 m，rs=0.3 m，K1、K2取值見表1。因此，δ分別為0.5、0.6、0.7、0.8、0.9、1.0、1.1、1.2、1.3、1.4，則利用經(jīng)驗(yàn)公式可以換算得到含水層地下水達(dá)西流速，其計(jì)算的相對(duì)誤差范圍在0.14%～0.23%。圖7 表示相對(duì)誤差隨滲透系數(shù)變化關(guān)系，結(jié)果表明，當(dāng)表皮區(qū)域的滲透系數(shù)增大時(shí)，經(jīng)驗(yàn)公式計(jì)算結(jié)果的相對(duì)誤差越高，這說明表皮區(qū)域滲透系數(shù)會(huì)影響經(jīng)驗(yàn)公式的計(jì)算精度，在實(shí)際操作中，應(yīng)關(guān)注表皮滲透系數(shù)對(duì)計(jì)算結(jié)果的影響。

圖7 相對(duì)誤差隨滲透系數(shù)變化Fig.7 Variation of relative error with hydraulic conductivity

4.2 表皮厚度

圖8 表示正表皮條件下相對(duì)誤差隨表皮厚度變化的關(guān)系。由于井筒附近表皮區(qū)域水力性質(zhì)變化，必然影響經(jīng)驗(yàn)公式的計(jì)算精度。模型主要參數(shù)分別為：rw=0.2 m，rh=0.1、0.2、0.3、0.4、0.5 m，δ=0.6，則利用經(jīng)驗(yàn)公式可以換算得到含水層地下水達(dá)西流速，其計(jì)算的相對(duì)誤差范圍為0%～0.08%，表明在正表皮狀態(tài)下，由經(jīng)驗(yàn)公式換算的值與實(shí)際值差異較小。由圖8 可知，隨著表皮厚度的增大，相對(duì)誤差呈逐漸增大的趨勢(shì)，相關(guān)系數(shù)R2=0.988 1，表明趨勢(shì)線擬合情況很好，說明經(jīng)驗(yàn)公式精度較高，在一定程度上可以較好地刻畫井筒內(nèi)地下水流速與含水層地下水達(dá)西流速之間的關(guān)系。

圖8 正表皮狀態(tài)相對(duì)誤差變化趨勢(shì)Fig.8 Trend of relative error of positive skin state

圖9 表示負(fù)表皮條件下相對(duì)誤差隨表皮厚度變化的關(guān)系。模型主要參數(shù)分別為：rw=0.2 m，rh=0.1、0.2、0.3、0.4、0.5 m，δ=1.3。則由經(jīng)驗(yàn)公式計(jì)算的相對(duì)誤差范圍為0.03%～0.07%，表明由式(8)-式(9)計(jì)算得出的換算后含水層地下水達(dá)西流速與含水層地下水達(dá)西流速較為接近。此外，由圖9 可知，相對(duì)誤差呈現(xiàn)逐漸減小的趨勢(shì)，但相關(guān)系數(shù)R2=0.052 9，表明此時(shí)擬合效果較差，且誤差波動(dòng)性較大，說明經(jīng)驗(yàn)公式在負(fù)表皮條件下計(jì)算含水層流速的穩(wěn)定性較差。

圖9 負(fù)表皮狀態(tài)相對(duì)誤差變化趨勢(shì)Fig.9 Trend of relative error of negative skin state

5 極端條件下收斂因子分析

在實(shí)際工程應(yīng)用中表皮厚度通常在0.001～1 m 范圍內(nèi)[27]，本文取最大表皮厚度rh=1 m，井筒半徑rw=0.1 m，則rs=1.1 m。此外，幾種典型的介質(zhì)滲透性如下：細(xì)砂、礫石、黏土的滲透系數(shù)范圍分別為0.20×10-6～1.89×10-4、4.75×10-4～0.995、1.0×10-12～9.9×10-10m/s[28]。

在實(shí)際地層中，成井條件往往非常復(fù)雜，受到鉆井泥漿侵入、生物化學(xué)堵塞等因素的影響，抽水井區(qū)域滲透系數(shù)減小[14-17]?？紤]極端正表皮條件下，此時(shí)由于受到上述因素的影響，表皮區(qū)域的滲透系數(shù)趨近于黏土，則滲透系數(shù)取1.0×10-12m/s。假定含水層介質(zhì)為細(xì)砂，則取含水層滲透系數(shù)為1.0×10-5m/s。將相關(guān)參數(shù)代入式(8)，計(jì)算極端正表皮條件下的水平收斂因子為4.03×10-7，通過數(shù)值模擬計(jì)算分析表明，在該條件下經(jīng)驗(yàn)公式計(jì)算誤差為2.23%。

在濾料充填的情況下，考慮極端負(fù)表皮狀態(tài)，此時(shí)表皮區(qū)域的滲透系數(shù)趨近于礫石，其滲透系數(shù)為0.995 m/s。同樣，假定含水層為細(xì)砂含水層，則含水層滲透系數(shù)為1.0×10-5m/s。將相關(guān)參數(shù)代入式(8)，計(jì)算極端負(fù)表皮條件下的水平收斂因子為3.694，其計(jì)算誤差為98%。綜合以上分析，極端條件下水平收斂因子的取值范圍為4.03×10-7～3.694。計(jì)算結(jié)果也表明，正表皮的影響范圍較大，因此，在實(shí)際計(jì)算過程中要關(guān)注泥漿侵入、生物化學(xué)堵塞等因素的影響。

6 結(jié)論

a.通過COMSOL Multiphysics 軟件建立考慮表皮效應(yīng)的試驗(yàn)數(shù)值模型，分析無表皮狀態(tài)下井筒井徑、表皮滲透系數(shù)以及表皮厚度對(duì)井筒內(nèi)地下水流速的影響，并驗(yàn)證M.Schubert 提出經(jīng)驗(yàn)公式的可靠性。

b.不考慮表皮效應(yīng)時(shí)，隨著井徑增加，僅在井徑較小時(shí)相對(duì)誤差有所增大，相對(duì)誤差均小于0.47%，說明經(jīng)驗(yàn)公式是可靠的，且孔內(nèi)地下水流速近似于含水層地下水達(dá)西流速的2 倍。

c.固定井筒井徑大小與表皮區(qū)域范圍，隨著表皮滲透系數(shù)增大，井筒內(nèi)地下水流速逐漸增大，相對(duì)誤差增大。通過計(jì)算結(jié)果對(duì)比，說明負(fù)表皮狀態(tài)下經(jīng)驗(yàn)公式的計(jì)算結(jié)果更接近實(shí)際。

d.在正表皮狀態(tài)下，隨著表皮厚度的增大，井筒內(nèi)地下水流速逐漸減小，相對(duì)誤差逐漸增大；在負(fù)表皮狀態(tài)下，隨著表皮厚度的增大，井筒內(nèi)地下水流速逐漸增大，計(jì)算相對(duì)誤差總體上呈現(xiàn)減小的趨勢(shì)，但其誤差波動(dòng)性較大，因此在該條件下經(jīng)驗(yàn)公式應(yīng)用過程中穩(wěn)定性較差。

符號(hào)注釋：

H為水頭，m；H1為邊界S1的水頭，m；H2為邊界S2的水頭，m；J為含水層水力梯度；K1為含水層滲透系數(shù)，m/s；K2為表皮區(qū)滲透系數(shù)，m/s；Kax、Kay、Kaz分別為含水層沿x、y、z方向的滲透系數(shù)，m/s；Ksx、Ksy、Ksz分別為表皮區(qū)域沿x、y、z方向上的滲透系數(shù)，m/s；r為與井筒中心的距離，m；rw為井筒半徑，m；rs為表皮區(qū)半徑，m；S1、S2為定水頭邊界；S3、S4為通量邊界；n為邊界S3、S4的外法線方向；vj為井筒內(nèi)地下水流速，m/s；vi為井筒內(nèi)具有代表性的地下水流速，m/s，i取13；vw為含水層地下水達(dá)西流速，m/s；vgw為換算后含水層地下水達(dá)西流速，m/s；α為水平收斂因子。