《南京工業(yè)大學(xué)學(xué)報(bào)》--承壓含水層地下水穩(wěn)定流的有限層分析-諸宏博

1、承壓含水層地下水穩(wěn)定流的有限層分析諸宏博,王旭東,宰金珉( 南京工業(yè)大學(xué)土木工程學(xué)院，江蘇 南京 )摘 要：有限層法是一種對(duì)空間某一方向進(jìn)行數(shù)值離散，而在其余兩方向采用連續(xù)函數(shù)的半數(shù)值半解析方法，該方法能有效地將三維問(wèn)題簡(jiǎn)化為一維問(wèn)題求解，從根本上解決了常用數(shù)值計(jì)算方法在模擬三維地下水運(yùn)動(dòng)時(shí)存在的計(jì)算工作量大，占用內(nèi)存多、耗時(shí)大等缺點(diǎn)。本文基于有限層法的優(yōu)點(diǎn)，推導(dǎo)了以伽遼金法結(jié)合貝塞耳函數(shù)為基礎(chǔ)的層狀非均質(zhì)各向異性承壓含水層的穩(wěn)定流有限層方程，并編制了相應(yīng)的計(jì)算程序。通過(guò)對(duì)兩個(gè)經(jīng)典算例的數(shù)值解與解析解對(duì)比分析，驗(yàn)證了本文方法的正確性。關(guān)鍵詞：有限層法；各向異性；承壓含水層；穩(wěn)定流；降深中圖分類(lèi)

2、號(hào)：P641 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：宏博（1981），男，碩士研究生，主要研究方向?yàn)榄h(huán)境巖土工程及數(shù)值分析。Analysis of Steady Groundwater Flow in the Confined Aquifer by Finite Layer Method ZHU Hong-bo , WANG Xu-dong , ZAI Jin-min (College of Civil Engineering , Nanjing University of Technology, Nanjing , China)Abstract: The finite layer method (FLM

3、) is a quasi-numerical and quasi-analytic method, which is to discretize one dimension of the spatial domain using finite elements, approximating variations in the other two dimensions using continuous function. The method therefore reduces three-dimensional problems to one-dimensional problems effe

4、ctively. The FLM can resolve the disadvantages radically, such as the large amount of works, the much memory and more time, et al, by using common numerical methods to simulate three-dimensional groundwater flows. Taking advantages of the FLM, finite layer formulations are derived for calculating th

5、e groundwater drawdown in the heterogeneous and anisotropic confined aquifer, which based on Galerkin method and Bessel function, and a computer program using FORTRAN language is developed. Two numerical examples are presented and compared with analytical solutions to demonstrate the validity of fin

6、ite layer method. Key words: finite layer method; anisotropic; confined aquifer; steady groundwater flow; drawdown地下水運(yùn)動(dòng)由于受含水層的非均質(zhì)、各向異性和邊界條件等因素的影響，其運(yùn)動(dòng)狀態(tài)十分復(fù)雜，對(duì)于非均質(zhì)、各向異性問(wèn)題很難用解析方法進(jìn)行求解，必須借助于數(shù)值模擬方法求其近似解1-2。國(guó)內(nèi)外學(xué)者不斷嘗試尋找有效的數(shù)值方法用于地下水模擬計(jì)算，薛禹群等將多尺度有限元法應(yīng)用于地下水模擬，節(jié)省了計(jì)算量，提高了計(jì)算精度3。朱學(xué)愚等采用邊界元法對(duì)二維穩(wěn)定流進(jìn)行了分析計(jì)算4。但是，數(shù)值法在進(jìn)行三

7、維問(wèn)題的計(jì)算分析時(shí)，常常會(huì)遇到計(jì)算分析占用內(nèi)存多、計(jì)算工作量大等問(wèn)題。Cheung于1968年首先提出了有限條法，并成功地將有限條法應(yīng)用于結(jié)構(gòu)工程問(wèn)題的分析5。該方法能有效地將三維問(wèn)題簡(jiǎn)化為一維問(wèn)題求解，從而大大減少單元數(shù)、節(jié)省計(jì)算工作量和內(nèi)存的需要量。有限層法正是利用了有限條法的特點(diǎn)，進(jìn)一步發(fā)展和完善起來(lái)的。它是一種對(duì)空間某一方向進(jìn)行數(shù)值離散，而在其余兩方向采用連續(xù)函數(shù)的半數(shù)值半解析方法。Small和Booker等運(yùn)用了有限層法對(duì)土的固結(jié)與流變耦合問(wèn)題做了較為系統(tǒng)的研究6；Stanley、Myron等將有限層理論應(yīng)用于地下水分析7；Southcott也采用該方法對(duì)豎向荷載下的樁土相互作用問(wèn)題

8、進(jìn)行了探討8；支喜蘭、王秉綱利用有限層方法對(duì)彈性地基上受對(duì)稱(chēng)荷載作用的圓形厚板進(jìn)行了力學(xué)分析9；宰金珉、王旭東等將有限層法應(yīng)用于土與結(jié)構(gòu)物的共同作用分析10-11。有限層法已在巖土工程數(shù)值計(jì)算中得到了廣泛的應(yīng)用。本文利用有限層法的優(yōu)點(diǎn)，以伽遼金法結(jié)合貝塞耳函數(shù)為基礎(chǔ)，推導(dǎo)了承壓含水層穩(wěn)定流有限層方程，編制了相應(yīng)的計(jì)算程序，通過(guò)算例驗(yàn)證了該方法的正確性。1 承壓含水層穩(wěn)定流模型的數(shù)學(xué)描述假設(shè)層狀非均質(zhì)、各向異性、等厚承壓含水層，初始水位水平，定水頭邊界條件，如圖1所示。根據(jù)水量均衡原理和達(dá)西定律，柱坐標(biāo)下承壓含水層穩(wěn)定流的數(shù)學(xué)模型可表示為 (1)式中：s(r, z)為降深；Kr、Kz分別為r徑向

9、、z方向的滲透系數(shù)；W(r, z)為源匯項(xiàng)，補(bǔ)給源取正值，排泄匯取負(fù)值；R為影響半徑。圖1 承壓含水層完整井物理模型 Fig.1 Configuration of fully penetrating well physical model in confined aquifer2 有限層格式推導(dǎo)有限層計(jì)算簡(jiǎn)圖如圖2所示。圓形計(jì)算區(qū)域徑向r的計(jì)算尺寸取為R，豎向z的計(jì)算尺寸取為c。同時(shí)沿z軸方向?qū)⒊袎汉畬与x散成L個(gè)層元。圖2 有限層計(jì)算簡(jiǎn)圖Fig.2 Finite layer models2.1 試探函數(shù)設(shè)降深在徑向用解析函數(shù)表示，而z方向采用有限層元離散，層元內(nèi)降深線(xiàn)性分布。則試探函數(shù)可表示

10、為下列分離變量的函數(shù)項(xiàng)級(jí)數(shù) (2)式中：為降深試探函數(shù)；mj為待求系數(shù)；為滿(mǎn)足定水頭邊界的第一類(lèi)零階貝塞耳函數(shù)，其中，為零階貝塞耳函數(shù)的第m個(gè)正零點(diǎn)；M為級(jí)數(shù)項(xiàng)數(shù)；L為有限層層元數(shù)。 Nj(z)為線(xiàn)性插值基函數(shù)，其形式如下： (3)2.2 伽遼金方程根據(jù)加權(quán)余量法原理，可得有限層法伽遼金方程 (4)式中：D為含水層計(jì)算區(qū)域；W為源匯項(xiàng)；為微分算子， ；為權(quán)函數(shù)，；為第一類(lèi)零階貝塞耳函數(shù)；Ni(z)同樣為線(xiàn)性插值基函數(shù)，其形式同式(3)。將含水層離散為L(zhǎng)個(gè)層元，則式(4)可改寫(xiě)為(5)式中：D e為層元計(jì)算區(qū)域；W e為層元源匯項(xiàng)。將式(2)代入式(5)得 (6)式中：，為層元線(xiàn)性插值基函數(shù)。2

11、.3 層元有限層方程取一典型有限層層元l，其節(jié)面分別為l，l+1，層厚cl，如圖3所示。圖3 典型有限層層元Fig.3 Finite layer element models利用貝塞耳函數(shù)Am (r)和Am (r) 的正交性，可得 (7)式中：為一階貝塞耳函數(shù)。故，當(dāng)時(shí)，由式(6)得： (8)結(jié)合基函數(shù)的性質(zhì)，由式(8)積分可得到層元l的有限層方程： (9)式中：為層元滲透矩陣系數(shù)；為水量矩陣系數(shù)。 (10)對(duì)于完整井抽水： (11)其中，為抽水率，Q為井流量。對(duì)于半球狀點(diǎn)源： (12)式(9)用矩陣可簡(jiǎn)化表示為： (13)式中：為層元滲透矩陣；為水量矩陣；為待求系數(shù)矩陣。2.4 整體有限層方

12、程由層元有限層方程組裝可得整體有限層方程： (14)式中：為整體滲透矩陣，；為整體水量矩陣， 。整體矩陣由層元矩陣疊加而成，疊加方法如圖4所示。圖4 Gm矩陣的疊加Fig.4 Matrix assembly for Gm 3 有限層方程求解式(14)整體有限層方程為L(zhǎng)+1階線(xiàn)性代數(shù)方程組，可采用追趕法12求解方程，求得未知系數(shù)mj后，由式(2)計(jì)算任意位置的降深。上述計(jì)算過(guò)程通過(guò)編制相應(yīng)的計(jì)算程序?qū)崿F(xiàn)。4 數(shù)值算例及分析4.1 承壓含水層中完整井假設(shè)承壓含水層厚c=20 m，為了便于對(duì)比分析，滲透系數(shù)取Kr = Kz = 5 m/d，計(jì)算區(qū)域取為R =2 500 m的圓形區(qū)域。在計(jì)算區(qū)域中心布

13、置一完整井，抽水量Q = 400 m3/d，如圖1所示。采用有限層法模擬計(jì)算完整井抽水引起的降深。為了驗(yàn)證承壓含水層穩(wěn)定流有限層分析方法的正確性，引入無(wú)限承壓含水層穩(wěn)定流完整井解析解2進(jìn)行對(duì)比分析。圖5 單井標(biāo)準(zhǔn)降深分布曲線(xiàn)Fig.5 Normalized drawdown curve of single well圖5給出了單井有限層解和解析解的標(biāo)準(zhǔn)降深曲線(xiàn)。其中，將距抽水井的距離和降深進(jìn)行無(wú)量綱化處理。由圖可知，無(wú)論是在抽水井附近，還是在遠(yuǎn)離抽水井的邊界處，有限層解與解析解都有較好的吻合程度，表明了本文方法的正確性。4.2 承壓含水層中點(diǎn)匯為了與半無(wú)限承壓含水層點(diǎn)源解析解13進(jìn)行對(duì)比分析，假

14、設(shè)承壓含水層厚c=800 m，滲透系數(shù)Kr = Kz = 5 m/d，計(jì)算區(qū)域取為半徑R=1 000 m的圓形區(qū)域。在含水層頂面中心布置一點(diǎn)匯，點(diǎn)匯的抽水量Q =200 m3/d，如圖6所示。采用有限層法模擬點(diǎn)匯抽水引起的降深。圖6 點(diǎn)匯物理模型Fig.6 The point sink physical model 圖7、圖8分別給出了以點(diǎn)匯為中心沿水平向和垂直向有限層數(shù)值解和解析解的標(biāo)準(zhǔn)降深曲線(xiàn)。其中，將距點(diǎn)匯的距離和降深進(jìn)行無(wú)量綱化處理。由于解析解中點(diǎn)匯趨于無(wú)窮小，所以在點(diǎn)匯附近（距點(diǎn)匯00.5 m）有限層解與解析解有一定的差異。但隨著距點(diǎn)匯距離的增加，有限層解與解析解則具有較好的一致性，

15、再一次驗(yàn)證了本文方法的正確性。圖7 點(diǎn)匯抽水標(biāo)準(zhǔn)降深分布曲線(xiàn)（水平向）Fig.7 Normalized drawdown curve of point sink（horizontal）圖8 點(diǎn)匯抽水標(biāo)準(zhǔn)降深分布曲線(xiàn)（垂直向）Fig.8 Normalized drawdown curve of point sink（vertical）5 結(jié) 論(1) 以伽遼金法結(jié)合貝塞耳函數(shù)為基礎(chǔ)，推導(dǎo)了層狀非均質(zhì)各向異性承壓含水層穩(wěn)定流有限層方程，并編制了相應(yīng)的計(jì)算程序。 (2) 算例的對(duì)比分析表明，有限層解和解析解具有良好的一致性，驗(yàn)證了地下水穩(wěn)定流分析方法的正確性。(3) 可將該方法進(jìn)一步推廣到非穩(wěn)定流問(wèn)

16、題的求解，并可應(yīng)用于熱傳導(dǎo)的計(jì)算分析。參考文獻(xiàn)1 朱學(xué)愚, 錢(qián)孝星, 劉新仁. 地下水資源評(píng)價(jià)M. 南京:南京大學(xué)出版社, 1987.2 孫訥正. 地下水流的數(shù)學(xué)模型和數(shù)值方法M. 北京:地質(zhì)出版社, 1981.3 薛禹群, 葉淑君, 謝春紅,等. 多尺度有限元法在地下水模擬中的應(yīng)用J. 水利學(xué)報(bào), 2004, (7): 7-13. 4 朱學(xué)愚, 謝春紅. 邊界元法在計(jì)算地下水穩(wěn)定水位和流量中的應(yīng)用J. 南京大學(xué)學(xué)報(bào):自然科學(xué)版, 1987,23(1):64-78.5 Cheung Y K. Finite strip method in structural mechanics M. New

17、York: Pergamon Press, 1976.6 Small J C, Booker J R. Finite layer analysis of primary and secondary consolidationC/ Proceedings of the 4th International Conference of Numerical Methods in Geomechanics. Edmonton:1982(1):365-371.7 Stanley S Smith, Myron B Allen, Jay Puckett. The finite layer method for groundwater flow models J.Water Resources Research, 1992, 28(6): 1715-1722.8 Southcott P H, Small J C. Finite layer analysis of vertically loaded piles and pile groupsJ.Computers and