【摘 要】介紹了基于有限元原理的FEFLOW軟件,它的應(yīng)用領(lǐng)域覆蓋了水量、水質(zhì)與溫度的模擬,軟件提供地理信息系統(tǒng)數(shù)據(jù)接口,能自動(dòng)產(chǎn)生空間有限單元網(wǎng)格。FEFLOW程序包中配備了快速精確的數(shù)值算法來(lái)控制和優(yōu)化求解過程,結(jié)果輸出具有先進(jìn)的圖形視覺化效果,并能很好地處理相關(guān)的水文地質(zhì)條件。舉例說明了FEFLOW軟件在建立黃河下游影響帶三維地下水流模擬模型中的應(yīng)用。

　　【關(guān)鍵詞】有限元原 理水文地質(zhì)條件 影響帶 模擬模型　

　　利用數(shù)值模擬軟件對(duì)地下水流和溶質(zhì)運(yùn)移等問題進(jìn)行模擬的方法以其有效性、靈活性和相對(duì)廉價(jià)性逐漸成為地下水研究領(lǐng)域的一種不可缺少的重要方法,得到了廣泛的應(yīng)用。模型計(jì)算指的是進(jìn)行地下水流動(dòng)或溶質(zhì)運(yùn)移正反演計(jì)算,常用的方法主要是有限差分法、有限元法、邊界元法等[1]。有限差分法是一種古典的數(shù)值計(jì)算方法,其基本思想是:用滲流區(qū)內(nèi)有限個(gè)離散本論文由

　　點(diǎn)的集合代替連續(xù)的滲流區(qū),在這些離散點(diǎn)上用差商近似地代替微商,將微分方程及其定解條件化為以未知函數(shù)在離散點(diǎn)上的近似值為未知量的代數(shù)方程(稱之為差分方程),然后求解差分方程,從而得到微分方程的解在離散點(diǎn)上的近似值[2]。美國(guó)地質(zhì)調(diào)查局在20世紀(jì)80年代開發(fā)的MODFLOW (Modular three-dimensional finite-difference ground-water flow model)為有限差分的典型代表。該軟件包括水井、補(bǔ)給、河流、溝渠、蒸發(fā)蒸騰和通用水頭邊界六個(gè)子程序包,分別用來(lái)處理相關(guān)的水文地質(zhì)條件[3]。自MODFLOW問世以來(lái),不斷有新的子程序包被開發(fā)出來(lái)。例如用來(lái)模擬水位下降引起地面沉降的子程序包[4],用來(lái)模擬水平流動(dòng)障礙(Horizontal flow-barrier)的子程序包等[5]。新子程序包的加入,使MODFLOW的應(yīng)用范圍不斷擴(kuò)大。實(shí)踐證明,只要合理使用,它也可以用來(lái)解決許多地下水在裂隙介質(zhì)中的流動(dòng)問題。經(jīng)過合理的線性化處理,MODFLOW還可以用于解決空氣在土壤中的運(yùn)動(dòng)問題[6]。如果將它與其他用于溶質(zhì)運(yùn)移模擬的程序結(jié)合起來(lái),可以用來(lái)模擬諸如海水入侵等地下水密度為變量的問題[7]。有限單元法是求解偏微分方程定解問題的另一種有效的數(shù)值方法,它的基本思想早在20世紀(jì)40年代就已經(jīng)提出。有限單元法成功地解決了眾多領(lǐng)域的許多計(jì)算問題,例如結(jié)構(gòu)力學(xué)、結(jié)構(gòu)工程學(xué)中的應(yīng)力分析,結(jié)構(gòu)穩(wěn)定性問題;土力學(xué)、巖石力學(xué)中的應(yīng)力變形與穩(wěn)定性分析等;水利工程的液體流動(dòng)、水工結(jié)構(gòu)和壩體分析以及地下水流動(dòng)、水動(dòng)力彌散問題等。20世紀(jì)70年代末,德國(guó)WASY水資源規(guī)劃和系統(tǒng)研究所開發(fā)了基于有限單元法的FEFLOW ( Finite element subsurface FLOWsystem)軟件,它是迄今為止功能最為齊全的地下水模擬軟件包之一。該軟件包具有圖形人機(jī)對(duì)話、地理信息系統(tǒng)數(shù)據(jù)接口、自動(dòng)產(chǎn)生空間各種有限單元網(wǎng)格、空間參數(shù)區(qū)域化及快速精確的數(shù)值算法和先進(jìn)的圖形視覺化技術(shù)等特點(diǎn)。在FEFLOW系統(tǒng)中,用戶可以很方便迅速地產(chǎn)生空間有限單元網(wǎng)格,設(shè)置模型的參數(shù)和定義邊界條件,運(yùn)行數(shù)值模擬以及實(shí)時(shí)圖形顯示結(jié)果與成圖。

　　1 FEFLOW簡(jiǎn)介從FEFLOW的問世起一直到現(xiàn)在,在理論研究和對(duì)實(shí)際問題的處理上,它經(jīng)過了不斷的發(fā)展、修改、提高,在它的發(fā)展過程中,FEFLOW經(jīng)過了大量的測(cè)試和檢驗(yàn),它成功地解決了一系列與地下水有關(guān)的實(shí)質(zhì)性問題,如判斷污染物遷移途徑、追溯污染物的來(lái)源,地?zé)岬哪M,海水入侵預(yù)測(cè)等。

　　1.1 FEFLOW的應(yīng)用領(lǐng)域水量模擬:模擬水源地開采或者油田注水對(duì)區(qū)域地下水流場(chǎng)的影響、模擬水庫(kù)放水或者河流斷流時(shí),河道沿線地下水流場(chǎng)的變化等;水質(zhì)模擬:模擬污染物在地下水中的遷移過程及其時(shí)間空間分布模式[8]、模擬沿海地區(qū)抽取地下水引起的海水入侵等[9];溫度模擬:模擬非飽和帶以及飽和帶溫度場(chǎng)的分布[10]。

　　1.2 系統(tǒng)輸入特點(diǎn)通過標(biāo)準(zhǔn)數(shù)據(jù)輸入接口,用戶既能直接利用已有的GIS空間多邊形數(shù)據(jù)生成有限單元網(wǎng)格,也可以基于地圖用鼠標(biāo)設(shè)計(jì),網(wǎng)格的數(shù)目可以自己指定,可以方便地調(diào)整網(wǎng)格的幾何形狀,增加和放疏網(wǎng)格大小等。用戶可以視具體情況定義第一、第二和第三類邊界,而且還可以對(duì)邊界條件增加特定的限制條件;能方便地定義復(fù)合含水層中的分層開采井和混合開采井以及注水井。所有邊界條件及其限制條件、滲透系數(shù)、補(bǔ)排量既可設(shè)置為常數(shù),也能定義為隨時(shí)間變化的函數(shù)。FEFLOW提供了克里格(Kriging)、阿基瑪(Akima)和距離反比加權(quán)(IDW)三種方法,對(duì)離散的空間抽樣數(shù)據(jù)進(jìn)行內(nèi)插或外推。輸入數(shù)據(jù)格式既可以是ASCⅡ碼文件,也可以是GIS地理信息系統(tǒng)文件,如ArcView的shp格式。

　　1.3 FEFLOW配備了若干先進(jìn)的數(shù)值求解法來(lái)控制和優(yōu)化求解過程 快速直接求解法,如PCG,BICGSTAB,CGS,GMRES以及帶預(yù)處理的再啟動(dòng)ORTHOMIN法;靈活多變的up-wind技術(shù),如流線up-wind,奇值捕捉法Shock capturing,以減少數(shù)值彌散;皮卡和牛頓迭代法求解非線性流場(chǎng)問題;自動(dòng)調(diào)節(jié)模擬時(shí)間步長(zhǎng);模擬污染物遷移過程包括對(duì)流、水動(dòng)力彌散、線性及非線性吸附、一階化學(xué)非平衡反應(yīng);為非飽和帶模擬提供了多種參數(shù)模型,如指數(shù)式、Van Genuchten式和多種形式的Richard方程;垂向滑動(dòng)網(wǎng)格(BASD)技術(shù)處理自由表面含水系統(tǒng)以及非飽和帶的模擬[11];開放性外部程序接口,以便用戶在FEFLOW系統(tǒng)中連接和使用自己的程序模塊。1.4 系統(tǒng)結(jié)果輸出及顯示FEFLOW提供了其他任何地下水模擬軟件都無(wú)法比擬的、豐富實(shí)用的圖形顯示和數(shù)據(jù)結(jié)果分析工具。其先進(jìn)的圖形視覺化及數(shù)據(jù)分析技術(shù)表現(xiàn)在:有限單元網(wǎng)、邊界條件和模型參數(shù)的三維可視化及顯示;標(biāo)量數(shù)據(jù)的三維彩色(透明或灰度)等勢(shì)面顯示以及其二維平面彩色或等值線顯示;三維地下水流徑追蹤,流動(dòng)時(shí)間及流速動(dòng)畫顯示(包括二維平面、剖面投影或二維平面追蹤);三維體截段的空間顯示和三維交叉剖面組的空間顯示;三維圖形的任意旋轉(zhuǎn),二維、三維圖形的放大或縮小;總體和局部水量平衡分析(包括任意幾何多邊形內(nèi)的水流通量分析);計(jì)算和圖形顯示通過各種邊界條件的水通量、物質(zhì)通量以及其在特定時(shí)間區(qū)間內(nèi)的積分量。

　　2 FEFLOW中特殊水文地質(zhì)問題的處理

　　2.1 地層的不連續(xù)在FEFLOW中,三維模型的空間結(jié)構(gòu)是由許多層(layer)組成,這些層由片(slice)隔開,每?jī)善芯褪且粚?圖1)。網(wǎng)格、邊界條件、初始條件是對(duì)應(yīng)于片,而參數(shù)(給水度、滲透系數(shù)等)是對(duì)應(yīng)于層。FEFLOW只能處理延伸到整個(gè)研究區(qū)域的片和層,所以,理論上FEFLOW不能處理地層不連續(xù)等特殊情況。

　　2.1.1 地層的缺失地層由于地殼的抬升而受強(qiáng)烈的剝蝕,含水巖層出露地面(圖1),或者含水層缺失。FEFLOW在處理這種情況時(shí),仍然認(rèn)為地層是連續(xù)到整個(gè)研究區(qū)域;只是在地層缺失的地方,認(rèn)為厚度是無(wú)限小,默認(rèn)是0.01 m,其參數(shù)參考相鄰層的參數(shù)值。

　　2.1.2 斷層斷層是地層發(fā)生錯(cuò)位,是一種特殊的水文地質(zhì)條件。到目前為止,FEFLOW還沒有很好的方法來(lái)處理斷層、跌水現(xiàn)象,WASY公司也一直設(shè)法解決這一難題。常用的處理方法是把斷層部分作為一個(gè)單獨(dú)的參數(shù)分區(qū),根據(jù)斷層兩邊的水力聯(lián)系程度來(lái)確定含水層的參數(shù)(主要是滲透系數(shù))。

　　2.1.3 小面積強(qiáng)透水帶的處理當(dāng)遇到天窗或者滲透系數(shù)非常大的區(qū)域時(shí)(這些區(qū)域的面積往往比較小),可以把它隔離出來(lái)當(dāng)成一單獨(dú)的層(圖2)。在這一特殊的層上,強(qiáng)透水帶位置的厚度、滲透系數(shù)等參數(shù)按照實(shí)際賦值;除此之外的地方,參數(shù)考慮相鄰層相應(yīng)的參數(shù)。

　　2.2 第三類邊界混合邊界河流在FEFLOW中是作為第三類邊界處理。河流通過底積層和含水層發(fā)生水力聯(lián)系(圖3)。根據(jù)達(dá)西定律有當(dāng)動(dòng)態(tài)地模擬一條河流時(shí),需要知道河流在典型水文年的水位動(dòng)態(tài)曲線。例如已知兩個(gè)水文站的水位歷時(shí)曲線,采用一維線性插值方法,可以得到沿主流道任意結(jié)點(diǎn)的水位歷時(shí)曲線。然而這種插值法用來(lái)模擬短期內(nèi)洪水對(duì)地下水的影響時(shí)就會(huì)出現(xiàn)較大的誤差,因?yàn)樗J(rèn)為上游和下游水位變化是同步的,沒有考慮時(shí)間滯后的影響;也就是說上游和下游的水位同時(shí)到達(dá)歷時(shí)曲線的波峰和波谷(圖4),這種結(jié)果并不是所需要的。為了解決上述問題,FEFLOW提供了一種考慮時(shí)間滯后的一維線性插值方法,通過這種插值(圖5),可以非常準(zhǔn)確地模擬一場(chǎng)暴雨后,沿河附近區(qū)域地下水位的變化。

　　3 FEFLOW軟件的缺點(diǎn)FEFLOW最大的缺點(diǎn)是用于處理源匯項(xiàng)的In Or Out Flow菜單功能過于集中,各個(gè)補(bǔ)排項(xiàng)沒有單獨(dú)的子程序包。這樣雖然節(jié)省了許多源代碼,但是卻給調(diào)參帶來(lái)麻煩。為了節(jié)省調(diào)參所需的時(shí)間,在前處理中必須對(duì)各個(gè)補(bǔ)排項(xiàng)進(jìn)行累計(jì),求出其代數(shù)和(借助自編的程序)。這樣在調(diào)參過程中,即使對(duì)補(bǔ)排項(xiàng)進(jìn)行較小的改動(dòng),只要運(yùn)行自編的程序,補(bǔ)排項(xiàng)的代數(shù)和就能瞬間得出,然后按一定的格式存儲(chǔ)于固定的文件,運(yùn)行FEFLOW時(shí)只要調(diào)用這一文件即可。

　　4 應(yīng)用實(shí)例

　　4.1 研究區(qū)水文地質(zhì)概況研究區(qū)位于黃河的下游,河南省的東北部,地勢(shì)平坦,地貌上屬于黃河沖積平原。地下水主要賦存于第四系松散沉積物之中,總厚度超過300 m。根據(jù)地下水的埋藏條件,可分為4個(gè)含水層,分別對(duì)應(yīng)著Q4,Q3,Q2和Q1層。結(jié)合地下水動(dòng)力條件和開采情況,可將以上4個(gè)含水層歸并為淺層潛水、微承壓水含水層組(Q4和Q3)和深層承壓含水層組(Q2和Q1)。地下水的主要補(bǔ)給方式為大氣降水入滲、渠系滲漏、灌溉入滲和黃河側(cè)滲補(bǔ)給,主要排泄方式為蒸散發(fā)和人工開采。地下水的開采主要集中在黃河兩側(cè)以及黃河古河道,以開采淺層水為主,部分地段開采深部承壓水。目前地下水年總開采量為1.152×109m3,其中淺層水開采量為1.043×109m3,深層水開采量0.109×109m3。

　　4.2 研究區(qū)邊界的確定研究區(qū)的邊界就是黃河影響帶的范圍,即黃河水的側(cè)滲寬度和循環(huán)深度。為了準(zhǔn)確刻劃出黃河影響帶的范圍,于2000年11月和2001年5月進(jìn)行了兩次淺層地下水位的統(tǒng)一調(diào)查,并取得地下水同位素樣品43個(gè),其中黃河水水樣5個(gè)、潛水樣26個(gè)(埋深60 m以內(nèi))、承壓水樣12個(gè)(埋深150~400 m)。采用淺層地下水流場(chǎng)分析法、水動(dòng)力學(xué)解析法、同位素分析法結(jié)合具體的水文地質(zhì)條件,確定了黃河影響帶的范圍,約10~25km[12]。根據(jù)同位素取樣分析和鉆孔資料初步確定了黃河水的循環(huán)深度,小于350 m。黃河概化的正確與否直接影響黃河和地下水之間側(cè)滲量的準(zhǔn)確性。由于黃河為“地上懸河”,且河床底部存在穩(wěn)定的弱透水層;同時(shí)考慮FEFLOW的特點(diǎn),將黃河概化成第三類邊界。黃河側(cè)滲量為黃河水位的函數(shù),黃河主河道上水文站的水位歷時(shí)曲線,通過時(shí)間滯后的一維線性插值方法,可以得到主河道上任意結(jié)點(diǎn)的水位歷時(shí)曲線。因此,只要已知任一年內(nèi)的黃河水位,理論上都可以非常準(zhǔn)確地模擬黃河的側(cè)滲量。

　　4.3 三維地下水流模擬模型及其識(shí)別研究區(qū)內(nèi)水文地質(zhì)概念模型概化為非均質(zhì)、各向同性、空間三維結(jié)構(gòu)、非穩(wěn)定地下水流系統(tǒng)。采用FEFLOW軟件對(duì)地下水流進(jìn)行模擬,采用三角形剖分,沿黃河兩岸和水源地附近加密網(wǎng)格,共剖分17 230個(gè)結(jié)點(diǎn)和26 304個(gè)單元。模擬時(shí)期選擇1999年1月到12月,時(shí)間步長(zhǎng)由模型自動(dòng)控制。源匯項(xiàng)包括降水、灌溉回滲、蒸發(fā)、人工開采等,各項(xiàng)均換算成相應(yīng)分區(qū)的開采強(qiáng)度,然后分配到相應(yīng)的單元格。各水文地質(zhì)單元分區(qū)的參數(shù)初值主要來(lái)自于相關(guān)的水文地質(zhì)報(bào)告。運(yùn)行FEFLOW,可得到在給定水文地質(zhì)參數(shù)和各均衡項(xiàng)條件下的地下水位時(shí)空分布,通過擬合同時(shí)期的流場(chǎng)和長(zhǎng)觀孔的歷時(shí)曲線,識(shí)別水文地質(zhì)參數(shù)、邊界值和其他均衡項(xiàng),使建立的模型能更準(zhǔn)確地定量描述研究區(qū)地下水系統(tǒng)。

　　4.4 研究成果研究表明,1999年研究區(qū)地下水系統(tǒng)總補(bǔ)給量為2.639×109m3,總排泄量2.656×109m3。淺層含水層組總補(bǔ)給量為2.622×109m3,其中降水入滲量占47.11%,河渠入滲量占22.41%,灌溉回滲量占15.02%,黃河側(cè)滲補(bǔ)給量占12.86%,總排泄量為2.630×109m3。多年平均淺層地下水補(bǔ)給資源量約2.835×109m3/a;75%降水保證率下,淺層地下水資源總補(bǔ)給量為2.600×109m3/a。淺層可開采資源量為1.941×109m3/a,深層可開采資源量為0.249×109m3/a。傍河水源地開采條件下,黃河側(cè)滲增量占水源地開采量的65%~75%;8個(gè)擬建水源地同時(shí)開采時(shí),黃河側(cè)滲增量占水源地開采量的67.34%,水源地開采10年后地下水位趨于穩(wěn)定[13]。截滲墻埋深為20 m和45 m,模型運(yùn)行一年后,黃河側(cè)滲量分別較1999年減少了0.175×108m3和0.238×108m3,占黃河側(cè)滲量的5.50%和7.48%[14]。距離黃河大堤10 km以外,流場(chǎng)基本不受截滲墻的影響。最近上市的FEFLOW 5.0,其顏色已經(jīng)不再局限于256色。隨著處理蒸發(fā)和斷層的子程序包的開發(fā),FEFLOW功能更加完善。

　　綜上所述,FEFLOW軟件能很好地處理相關(guān)的水文地質(zhì)條件,具有簡(jiǎn)便的操作、快速精確的數(shù)值算法和先進(jìn)的圖形視覺化技術(shù)等特點(diǎn)。隨著FEFLOW的日益普及,它必將成為水文地質(zhì)學(xué)上一個(gè)強(qiáng)有力的數(shù)值模擬工具。

　　參考文獻(xiàn)：

　　[1]丁繼紅,周德亮,馬生忠.國(guó)外地下水模擬軟件的發(fā)展現(xiàn)狀與趨勢(shì)[J].勘察科學(xué)技術(shù),2002,(1):37-42.

　　[2]孫訥正.地下水流的數(shù)學(xué)模型和數(shù)值方法[M].北京:地質(zhì)出版社,1981.