一種模擬介質(zhì)交界面處達(dá)西流速的分區(qū)多尺度有限單元法

謝一凡，吳吉春，葉 逾，謝春紅，魯春輝，5

（1.河海大學(xué) 水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210098；2.河海大學(xué) 水利水電學(xué)院，江蘇 南京 210098；3.南京大學(xué) 地球科學(xué)與工程學(xué)院，江蘇 南京 210023；4.南京大學(xué) 數(shù)學(xué)系，江蘇 南京 210093；5.河海大學(xué) 長江保護(hù)與綠色發(fā)展研究院，江蘇 南京 210098）

1 研究背景

天然地下水系統(tǒng)常由多種不同的含水介質(zhì)組成，伴隨著裂隙、透鏡體等地質(zhì)特征。在含水介質(zhì)交界面處，根據(jù)折射定律，交界面法向達(dá)西流速和水頭具有連續(xù)性，切向達(dá)西流速則具有不連續(xù)性且需按折射定律呈比例［1-3］。然而，有限元等傳統(tǒng)方法難以保證達(dá)西流速的連續(xù)性，更不能保證交界面處的達(dá)西流速符合折射定律。另一方面，地下水問題具有時(shí)空大尺度特性、非均質(zhì)性等特征。有限元法在模擬地下水問題時(shí)常需要精細(xì)剖分的網(wǎng)格以保證精度，而精細(xì)剖分的網(wǎng)格中含有大量的節(jié)點(diǎn)，需要大量的計(jì)算消耗和計(jì)算時(shí)間進(jìn)行模擬，導(dǎo)致有限元等傳統(tǒng)方法在模擬地下水問題時(shí)的計(jì)算效率較低［4-5］。本文提出的分區(qū)多尺度有限元法（MSFEM-D）能夠有效解決上述問題，不僅能夠令交界面處的達(dá)西流速滿足折射定律，并且具有極高的計(jì)算效率。

連續(xù)達(dá)西流速算法能夠保證節(jié)點(diǎn)達(dá)西流速值唯一，即流入流出某截面的流量守恒，具有重要意義。然而，有限元等傳統(tǒng)算法無法保證節(jié)點(diǎn)水力梯度連續(xù)性，也無法獲得連續(xù)的節(jié)點(diǎn)達(dá)西流速?，F(xiàn)有連續(xù)達(dá)西流速算法具有比傳統(tǒng)有限元等方法更高的精度，主要有兩種模擬方式。第一種是將水流方程和達(dá)西方程進(jìn)行耦合，將水頭和達(dá)西流速同時(shí)作為未知項(xiàng)，迭代求解達(dá)西流速。這種方式具有較高的計(jì)算精度，但需要較大的計(jì)算消耗來進(jìn)行迭代，代表方法為混合有限元法［6-7］。第二種則是將水頭和達(dá)西流速先后分別通過水流方程和達(dá)西方程進(jìn)行模擬。這種方式的算法簡單直接，也具有較高的計(jì)算精度，應(yīng)用比較廣泛，代表方法有Yeh的伽遼金有限元模型［8］、張志輝等的三次樣條法［9］、Batu的雙重網(wǎng)格法［10］等。其中，Yeh的伽遼金有限元模型應(yīng)用最廣泛，還可保證全局質(zhì)量守恒［11］。

雖然Yeh在傳統(tǒng)有限元基礎(chǔ)上改進(jìn)了模型［8］，能保證節(jié)點(diǎn)達(dá)西流速的連續(xù)性，卻無法保證折射定律對(duì)切向達(dá)西流速部分的要求，即無法保證交界面兩側(cè)的切向達(dá)西流速按照折射定律呈比例［12］?？茖W(xué)工作者為此提出了間斷有限元法［13］、改進(jìn)的Yeh的S1有限元模型［3］等交界面處理技術(shù)。其中，Zhou等改進(jìn)Yeh的模型提出的S1有限元模型具有較高的精度，并且原理十分簡單，適用于多種形態(tài)的交界面［3］。該模型通過折射定律構(gòu)造了一種Jump函數(shù)，并運(yùn)用其改進(jìn)了Yeh的模型中交界面節(jié)點(diǎn)部分的方程，從而通過迭代修正達(dá)西流速至交界面某一側(cè)的真實(shí)達(dá)西流速，再通過折射定律得到交界面另一側(cè)的達(dá)西流速；但因其運(yùn)用Jump函數(shù)時(shí)需要迭代，需要比Yeh的模型更大的計(jì)算消耗［3］。

另一方面，上述多種達(dá)西流速算法均受到有限元算法框架的限制，在模擬大尺度、非均質(zhì)地下水問題時(shí)效率較低。對(duì)此，侯一釗等提出多尺度有限元法（MSFEM）［14］是一種有效對(duì)策，它能夠通過構(gòu)造滿足局部水流方程的多尺度基函數(shù)來描述細(xì)尺度信息，從而將細(xì)尺度問題轉(zhuǎn)化為粗尺度問題來大幅降低計(jì)算消耗，具有比有限元等傳統(tǒng)方法更高的計(jì)算效率［4，14-15］。近年來，該方法不斷被國內(nèi)外科學(xué)工作者改進(jìn)，提出了多種新型地下水多尺度算法［16-19］，并成功應(yīng)用于地面沉降等實(shí)際問題［20］，顯示該方法在地下水?dāng)?shù)值模擬領(lǐng)域的應(yīng)用前景。作者也將多尺度有限元法與上述第二種方式算法中的三次樣條法、雙重網(wǎng)格法結(jié)合，提出了多種高效達(dá)西流速算法［21-23］，取得了較好的模擬結(jié)果。

本文提出的MSFEM-D屬于第二種方式的達(dá)西流速算法，可以有效解決上述問題。MSFEM-D不僅能夠令交界面處的達(dá)西流速滿足物理規(guī)律，且具有下面兩方面優(yōu)勢(shì)：（1）該方法不僅能夠通過Yeh的伽遼金有限元模型［8］保證達(dá)西流速的連續(xù)性，而且能夠通過區(qū)域分解技術(shù)將研究區(qū)用交界面分區(qū)，避免不同介質(zhì)在模擬時(shí)的相互影響并提升計(jì)算效率，同時(shí)結(jié)合Jump函數(shù)［3］保證交界面達(dá)西流速符合折射定律，且無需任何迭代過程。（2）該方法通過多尺度有限元框架［14］替代有限元框架，能夠提升水頭和達(dá)西流速的問題尺度，將大量細(xì)尺度節(jié)點(diǎn)替換為少量的粗尺度節(jié)點(diǎn)，從而降低解水頭和達(dá)西流速計(jì)算消耗，具有極高的計(jì)算效率。

2 原理和算法

2.1 MSFEM-D的算法流程 MSFEM-D算法適用于地下水二維問題，算法流程圖如圖1。MSFEMD算法分為兩個(gè)部分，即水流問題模擬部分（圖1左上虛線框部分）和達(dá)西流速模擬部分（圖1右側(cè)）。

圖1 MSFEM-D算法流程

MSFEM-D 的地下水流問題模擬部分采用現(xiàn)有的MSFEM 算法［14］模擬水頭，通過構(gòu)造多尺度網(wǎng)格和多尺度基函數(shù)可以顯著提升計(jì)算效率并保證精度。MSFEM-D的達(dá)西流速模擬部分為本文的主要?jiǎng)?chuàng)新部分，該部分綜合了區(qū)域分解算法、MSFEM［14］、Yeh的伽遼金有限元模型［8］和Zhou根據(jù)折射定律構(gòu)造的Jump函數(shù)［3］的優(yōu)點(diǎn)。首先，MSFEM-D運(yùn)用交界面將研究區(qū)剖分為子區(qū)域，并運(yùn)用區(qū)域分解算法將研究區(qū)上的模擬問題劃分為若干子區(qū)域上的子問題。其次，MSFEM-D 將MSFEM 和Yeh 的伽遼金有限元模型進(jìn)行了組合，即，運(yùn)用MSFEM的多尺度基函數(shù)替換Yeh的伽遼金有限元模型中的有限元基函數(shù)。再次，運(yùn)用MSFEM和Yeh的伽遼金有限元模型的組合方法求解子區(qū)域上的子問題；在獲得首個(gè)子問題上交界面處的達(dá)西流速值后，通過折射定律構(gòu)造的Jump函數(shù)獲得其相鄰子區(qū)域交界面上的達(dá)西流速，并作為下一個(gè)子問題的邊界條件。最后，循環(huán)這個(gè)過程，直到求解完畢。需要指出的是，在實(shí)際工作中，需要根據(jù)實(shí)際地下水問題的條件進(jìn)行建立概念模型、選擇數(shù)學(xué)模型，再應(yīng)用本算法，并進(jìn)行模型校正等后續(xù)步驟以完成地下水問題的模擬工作。

2.2 MSFEM-D的網(wǎng)格構(gòu)造 設(shè)水平方向?yàn)閤方向，垂直方向?yàn)閥方向。圖2為MSFEM-D對(duì)含一條水平交界面的示例研究區(qū)Ω的剖分，具體分為水頭模擬網(wǎng)格和達(dá)西流速模擬網(wǎng)格兩部分。

2.2.1 MSFEM-D 的水頭模擬網(wǎng)格 在模擬水頭時(shí)，MSFEM-D 運(yùn)用多尺度網(wǎng)格進(jìn)行剖分［14］，先將研究區(qū)剖分為粗網(wǎng)格，再將粗網(wǎng)格剖分為細(xì)網(wǎng)格。本方法的粗、細(xì)網(wǎng)格的數(shù)量可以根據(jù)數(shù)值模擬具體需要進(jìn)行設(shè)置。圖2為本文的數(shù)值模擬實(shí)驗(yàn)所用的剖分，MSFEM-D先將Ω直接剖分為1800（30×30×2）個(gè)三角形粗網(wǎng)格單元，此時(shí)交界面AB不作為分區(qū)界限。然后，MSFEM-D 再將每個(gè)粗網(wǎng)格單元剖分為9（3×3）個(gè)細(xì)網(wǎng)格單元，圖2右側(cè)小三角形Δijk為一示例粗網(wǎng)格單元及其細(xì)網(wǎng)格單元。

2.2.2 MSFEM-D的達(dá)西流速模擬網(wǎng)格 在模擬達(dá)西流速時(shí)，MSFEM-D先運(yùn)用交界面將研究區(qū)剖分為若干個(gè)僅含單一介質(zhì)子區(qū)域，再進(jìn)行粗、細(xì)網(wǎng)格剖分，具體粗、細(xì)網(wǎng)格的數(shù)量可根據(jù)數(shù)值模擬的需要進(jìn)行設(shè)置。如圖2中，MSFEM-D運(yùn)用AB將研究區(qū)剖分為Zone 1和Zone 2兩個(gè)子區(qū)域，然后延用水頭的模擬網(wǎng)格分別將兩個(gè)子區(qū)域剖分為900（30×15×2）個(gè)粗網(wǎng)格單元，再將每個(gè)粗網(wǎng)格單元剖分為9個(gè)細(xì)網(wǎng)格單元。

圖2 MSFEM-D對(duì)示例研究區(qū)的剖分

2.3 折射定律與Jump函數(shù) 在含水介質(zhì)交界面處，地下水水頭和達(dá)西流速需符合折射定律［1-2］：

式中：V為達(dá)西流速向量；n為交界面的法向單位向量；H為水頭；+、-號(hào)分別代表交界面的兩側(cè)。

2.4 構(gòu)造多尺度基函數(shù) 設(shè)示例粗網(wǎng)格單元為Δijk，其頂點(diǎn)i、j、k的多尺度基函數(shù)分別為Ψi、Ψj、Ψk，Δijk具有p個(gè)內(nèi)點(diǎn)，Mτ，τ=1，2，…，p。以多尺度基函數(shù)Ψi的構(gòu)造過程為例，考慮在粗網(wǎng)格單元Δijk上的簡化的橢圓方程：

在設(shè)定多尺度基函數(shù)的邊界條件后［14］，式（4）適定。

在粗網(wǎng)格單元Δijk中的任意細(xì)網(wǎng)格單元Δabc內(nèi)，Ψi可被線性表示為［14］：

式中：Na、Nb、Nc為有限元線性基函數(shù)；Ψi(a)、Ψi(b)、Ψi(c)為Δabc頂點(diǎn)處的Ψi的值。對(duì)式（4）進(jìn)行伽遼金變換，再將式（5）代入，可以得到關(guān)于Ψi的方程組，求解可得Ψi在Δijk中所有節(jié)點(diǎn)的值。其他多尺度基函數(shù)的構(gòu)造過程與Ψi類似，這里不再贅述。

2.5 運(yùn)用MSFEM-D模擬二維地下水問題

2.5.1 運(yùn)用MSFEM-D模擬地下水水頭 以地下水二維非穩(wěn)定流問題為例，考慮如下方程：

式中：W為源匯項(xiàng)；Ω為研究區(qū)；S為貯水系數(shù)。在地下水模擬問題的定解條件確定后，式（6）適定。

在水頭模擬時(shí)，無需運(yùn)用交界面進(jìn)行分區(qū)，可以直接剖分研究區(qū)。在式（6）兩邊乘以多尺度基函數(shù)Ψi，進(jìn)行伽遼金變換，再離散到每一粗網(wǎng)格單元Δijk上，可得：

根據(jù)文獻(xiàn)［14］，在粗網(wǎng)格單元Δijk上，Ψi可被線性表示為：

式中Hi、Hj、Hk分別為水頭H在粗網(wǎng)格單元Δijk頂點(diǎn)i、j、k處的值。

將式（7）離散到細(xì)網(wǎng)格上，并代入式（5）、式（8）可得式（7）的具體表達(dá)式。再運(yùn)用Ψj、Ψk乘以式（6）兩端，并重復(fù)上述過程，可得粗網(wǎng)格Δijk上關(guān)于水頭H的單元?jiǎng)偠染仃?。?lián)立所有粗網(wǎng)格上的單元?jiǎng)偠染仃?，結(jié)合源匯項(xiàng)W、初始條件、邊界條件即可得到右端項(xiàng)，從而得到研究區(qū)上的總方程組，結(jié)合Cholesky分解法等矩陣求解方法即可獲得研究區(qū)上的水頭在研究區(qū)粗尺度節(jié)點(diǎn)上的值。

2.5.2 運(yùn)用MSFEM-D模擬達(dá)西流速 在獲得研究區(qū)的水頭分布后，MSFEM-D可以運(yùn)用這些已知的水頭值，結(jié)合區(qū)域分解技術(shù)、MSFEM［14］、Yeh的伽遼金有限元法［8］以及Zhou的Jump函數(shù)［3］進(jìn)行達(dá)西流速的求解。其中，水頭值可以提供研究區(qū)的全局信息，區(qū)域分解技術(shù)可以分離介質(zhì)避免互相影響、MSFEM基函數(shù)可以提高模擬效率、Yeh的伽遼金法可以保證達(dá)西流速的連續(xù)性。以Vx的求解過程為例，在圖2所示的研究區(qū)Ω上考慮達(dá)西方程：

如2.2.2 節(jié)所述，MSFEM-D 已經(jīng)將研究區(qū)Ω劃分為Zone 1、Zone 2 兩個(gè)子區(qū)域。MSFEM-D 運(yùn)用區(qū)域分解技術(shù)將研究區(qū)上的達(dá)西流速求解問題劃分為Zone 1、Zone 2這兩個(gè)子區(qū)域的子問題。設(shè)先求解Zone 1上的子問題，則需在Zone 1上考慮式（9）。設(shè)子區(qū)域Zone 1上的粗尺度網(wǎng)格數(shù)目為γ1，粗尺度節(jié)點(diǎn)數(shù)目為Nu1，運(yùn)用ΨI，I=1，2，…，Nu1乘以達(dá)西方程的兩端，并在Zone1 上積分，并將Zone 1離散為γ1個(gè)粗網(wǎng)格Δijk，可得：

在任意粗網(wǎng)格Δijk上，根據(jù)現(xiàn)有文獻(xiàn)中的理論［8，14］，Vx可以被多尺度基函數(shù)線性表示為：

式中Vx(xi，yi)、Vx(xj，yj)、Vx(xk，yk)為達(dá)西流速Vx在粗網(wǎng)格Δijk頂點(diǎn)i、j、k上的值。

將式（11）代入式（10）左端，將式（8）代入式（10）右端，可得：

將式（12）、式（13）結(jié)合式（5），可得Zone 1上關(guān)于達(dá)西流速Vx的方程組的具體形式，求解可以獲得Zone 1 區(qū)域上的Vx值，其中包括了在Zone 1 的邊界上，即在Zone 1 一側(cè)的介質(zhì)交界面AB上的Vx值。結(jié)合Jump函數(shù)式（3），即可獲Zone 2一側(cè)的介質(zhì)交界面AB上的Vx值，這將作為Zone 2上的子問題的第一類邊界條件。在Zone 2上，運(yùn)用類似上述達(dá)西流速的模擬過程，即可獲得Zone 2區(qū)域的粗尺度Vx值。在所有子問題求解完畢后，可以得到Ω上的Vx值。在Ω進(jìn)行類似上述Vx的達(dá)西流速的模擬過程，即可以得到Ω上的粗尺度Vy值。

上述過程中，MSFEM-D運(yùn)用了多尺度基函數(shù)模擬地下水流和達(dá)西流速，得到了水頭和達(dá)西流速的粗尺度網(wǎng)格上的解。因而，MSFEM-D的水頭和達(dá)西流速的總方程組的階數(shù)均遠(yuǎn)低于精細(xì)剖分的傳統(tǒng)方法，具其有極高的計(jì)算效率。

3 數(shù)值模擬實(shí)驗(yàn)

采用如下縮寫方式：Vx表示x方向上的達(dá)西流速；Vy表示y方向上的達(dá)西流速；FEM-F表示精細(xì)剖分的有限元法；Model-Yeh-F 表示精細(xì)剖分的Yeh的伽遼金有限元模型；Model-Zhou-F 表示精細(xì)剖分的Zhou等的S1有限元模型；MSFEM-D表示多尺度有限元與區(qū)域分解法的組合方法。

本節(jié)模擬了含水平介質(zhì)交界面的滲透系數(shù)非均質(zhì)各向異性的地下水二維穩(wěn)定流問題以及含水平與傾斜兩個(gè)介質(zhì)交界面的地下水穩(wěn)定流問題的水頭和達(dá)西流速，并將MSFEM-D 所模擬的水頭與FEM-F比較，MSFEM-D所模擬的達(dá)西流速與Model-Yeh-F［8］以及Model-Zhou-F［3］兩種現(xiàn)有達(dá)西流速算法進(jìn)行比較。為了比較計(jì)算消耗與精度，F(xiàn)EM-F、Model-Yeh-F、Model-Zhou-F 的網(wǎng)格數(shù)目均等于MSFEM-D的細(xì)網(wǎng)格數(shù)目，且所有數(shù)值方法均采用了線性基函數(shù)的邊界條件。同時(shí)，所有數(shù)值方法均采用C++編寫，均采用Cholesky 分解法求解總方程的矩陣，沒有應(yīng)用并行計(jì)算技術(shù)。在測量CPU時(shí)間時(shí)，各方法均在同一計(jì)算機(jī)上獨(dú)立運(yùn)行。

3.1 含水平交界面的滲透系數(shù)非均質(zhì)各向異性的二維地下水穩(wěn)定流問題 本例基于Cainelli等在2012年工作［4］中的數(shù)值算例所設(shè)置，具有解析解。本例的地下水水流和達(dá)西流速的控制方程分別為，二維地下水穩(wěn)定流方程和達(dá)西方程（式（9））。研究區(qū)域Ω為［0，1］×［0，1］。研究區(qū)內(nèi)含有兩種不同的各向同性含水介質(zhì)（圖2），分別位于Zone 1和Zone 2兩個(gè)區(qū)域，它們的交界面為AB（y=0.5）。滲透系數(shù)呈非均質(zhì)各向異性變化，Zone 1、Zone 2的滲透系數(shù)分別為：K1x=（1+x）（1+y），K1y=10（1+x）（1+y）和K2x=10（1+x）（1+y），K2y=100（1+x）（1+y）。本例中，地下水水流方程中的源匯項(xiàng)W和研究區(qū)的邊界的第一類邊界條件由如下水頭H的解析解給出：

此外，x方向上的達(dá)西流速Vx的解析解為：

y方向上的達(dá)西流速Vy的解析解為：

本例中，MSFEM-D先將研究區(qū)剖分為1800（30×30×2）個(gè)三角形粗網(wǎng)格單元，再將每個(gè)粗網(wǎng)格剖分為9（3×3）個(gè)三角形細(xì)網(wǎng)格單元，即總計(jì)16 200個(gè)細(xì)網(wǎng)格單元；Model-Zhou-F、Model-Yeh-F 直接將研究區(qū)剖分16 200（90×90×2）個(gè)網(wǎng)格單元。

如前所述，MSFEM-D 應(yīng)用MSFEM 直接在研究區(qū)內(nèi)計(jì)算全局水頭，而Model-Zhou-F、Model-Yeh-F 均使用FEM-F 進(jìn)行全局水頭模擬，F(xiàn)EM-F 的網(wǎng)格剖分與Model-Zhou-F、Model-Yeh-F 一致。MSFEM-D和FEM-F在本例都能獲得很高的水頭計(jì)算精度，全局平均誤差分別為0.079%和0.024%。

圖3展示了解析解Vx以及各方法的數(shù)值解Vx在交界面AB（圖3（a））以及截面x=0.5 m處（圖3（b））的值。從圖3（a）中可以看出，MSFEM-D與Model-Zhou-F所獲的數(shù)值解以及解析解在AB兩側(cè)均具有兩個(gè)不同的值。這是由于水流在交界面產(chǎn)了折射，引起AB的切向達(dá)西流速Vx在交界面上不連續(xù)。MS?FEM-D、Model-Zhou-F 的曲線與解析解重合，顯示這兩個(gè)方法均具有較高的精度并符合折射定律，并且多尺度網(wǎng)格和精細(xì)剖分的網(wǎng)格都可以較有效的處理非均質(zhì)滲透系數(shù)。另一方面，Model-Yeh-F在交界面兩側(cè)的值相等，且位于解析解兩側(cè)值中間的位置，具有較大的誤差。這是由于Model-Yeh-F在求解達(dá)西方程時(shí)，沒有對(duì)交界面進(jìn)行處理，因而其得到的解為交界面兩側(cè)達(dá)西流速的平均值，不符合折射定律。因此，Model-Yeh-F與其他方法在交界面處的誤差值具有較大的差異，以Zone 2側(cè)的達(dá)西流速結(jié)果為例：MSFEM-D、Model-Zhou-F 的相對(duì)誤差值分別為0.52%和0.15%；Model-Yeh-F對(duì)于Zone 2側(cè)的解析解的誤差為45%。此外，由于Zone 1側(cè)的解析解絕對(duì)值較小，Model-Yeh-F的對(duì)于Zone 1的解析解的相對(duì)誤差甚至達(dá)到了100%。從圖3（b）中可以看出，在交界面附近節(jié)點(diǎn)上Model-Yeh-F 的解出現(xiàn)了振蕩，這是由于Model-Yeh-F 在交界面附近節(jié)點(diǎn)上的解需要“補(bǔ)償”其在點(diǎn)（0.5，0.5）的解的誤差。由于MSFEM-D、Model-Zhou-F對(duì)交界面進(jìn)行了處理，它們的解在交界面附近依然保持了高精度，和解析解重合。在遠(yuǎn)離交界面AB的位置，Model-Yeh-F 可以表現(xiàn)的很好，精度與MSFEM-D、Model-Zhou-F相近，故圖中三個(gè)方法均與解析解重合。事實(shí)上，Model-Yeh-F在遠(yuǎn)離交界面處的值與Model-Zhou-F完全一致，這是因?yàn)镸odel-Zhou-F僅對(duì)交界面處的達(dá)西流速的相關(guān)方程進(jìn)行了修正。

圖3 達(dá)西流速Vx的解析解和MSFEM-D、Model-Zhou-F和Model-Yeh-F所獲的數(shù)值解

圖4展示了解析解Vy以及各方法所模擬的數(shù)值解Vy在研究區(qū)Ω上的分布。由于Vy在交界面AB上連續(xù)，MSFEM-D 可以直接運(yùn)用水頭網(wǎng)格在研究區(qū)Ω上求解，Model-Zhou-F 無需運(yùn)用Jump函數(shù)進(jìn)行修正，因而Model-Zhou-F和Model-Yeh-F所獲的Vy場完全相同。如圖4所示，雖然y方向的滲透系數(shù)在Zone 1、2 內(nèi)不同，但解析解和各方法的數(shù)值解Vy在交界面AB處均連續(xù)，這符合折射定律。同時(shí)，MSFEM-D 和Model-Zhou-F 所獲的數(shù)值解與解析解很接近，兩種方法的全局平均相對(duì)誤差均分別為1.1%和0.2%。這是由于Model-Zhou-F是直接在細(xì)尺度上求解方程的，能夠獲得比MSFEM類方法略高的計(jì)算精度，但這也產(chǎn)生了大量的計(jì)算消耗，降低了其的計(jì)算效率。

圖4 達(dá)西流速Vy的解析解以及MSFEM-D和Model-Zhou-F所獲的數(shù)值解

本例中各方法的計(jì)算消耗對(duì)比如表1 所示。其中，Model-Zhou-F 所消耗的4236 s 的CPU 時(shí)間中，204 s為運(yùn)用FEM-F模擬水頭的時(shí)間；3789 s為模擬Vx場的時(shí)間，因處理交界面而迭代了15次；243 s為模擬Vy場的時(shí)間，因Vy在交界面連續(xù)而無需迭代。Model-Yeh-F則使用了204 s模擬水頭，并分別運(yùn)用兩個(gè)243 s模擬Vx場和Vy場。和其他方法不同，MSFEM-D 所構(gòu)造的模擬水頭、達(dá)西流速的方程組的階數(shù)都很低，故MSFEM-D僅使用不到1 s即可完成本例的模擬，且能獲得與Model-Zhou-F幾乎相同的效果，具有最高的計(jì)算效率。

表1 例3.1中各方法的計(jì)算消耗對(duì)比

3.2 含水平與傾斜兩個(gè)介質(zhì)交界面的二維地下水穩(wěn)定流問題 本例是基于Zhou等在2001年工作中的數(shù)值算例［3］。本例的地下水水流和達(dá)西流速的控制方程與上例相同。研究區(qū)域Ω為［0，120 m］×［0，120 m］。如圖5所示，研究區(qū)所含的兩條交界面為水平交界面AB（y=20 m）和傾斜交界面CD（x+y=160 m）。這兩條交界面將研究區(qū)分為Zone 1、Zone 2、Zone 3 三個(gè)區(qū)域，相應(yīng)的滲透系數(shù)分別為Kx=Ky=100 m/d、K2x=K2y=10 m/d 和K3x=K3y=100 m/d。研究區(qū)上下邊界為隔水邊界，左右邊界為定水頭邊界，左右水頭值分別為0 m和1 m。

本例采用MSFEM-D、Model-Zhou-F 和Model-Yeh-F 進(jìn)行模擬。MSFEM-D 將研究區(qū)剖分為1152（24×24×2）個(gè)粗網(wǎng)格單元，再將每個(gè)粗網(wǎng)格單元剖分為4個(gè)細(xì)網(wǎng)格單元（圖5右側(cè)）；MSFEM-D在模擬達(dá)西流速時(shí)，在每個(gè)子區(qū)域中延用了水頭網(wǎng)格；Model-Zhou-F、Model-Yeh-F 直接將研究區(qū)剖分4608（48×48×2）個(gè)網(wǎng)格單元。本例沒有解析解，故將運(yùn)用Model-Zhou-F 將研究區(qū)剖分為18 432（96×96×2）個(gè)單元獲得的解作為“參照解”。

圖5 例3.2中的研究區(qū)Ω示意圖

各方法所獲的水頭值如圖6所示?？梢钥闯觯芯繀^(qū)的水頭場是按照?qǐng)D6所示的區(qū)域進(jìn)行分布的，能夠較明顯的看到研究區(qū)內(nèi)所包含的交界面。MSFEM-D運(yùn)用MSFEM所獲的水頭場以及Model-Zhou-F和Model-Yeh-F運(yùn)用FEM-F所獲的水頭場的精度都很高，水頭平均相對(duì)誤差分別為0.58%和0.15%。

圖6 水頭參照解以及MSFEM-D和FEM-F所獲的數(shù)值解水頭

圖7 展示了參照解以及各方法的數(shù)值解Vx在交界面AB（圖7（a））以及截面x=80 m 處（圖7（b））的值。從圖7（a）中可以看出，MSFEM-D 和Model-Zhou-F 能夠滿足折射定律，并與參照解重合，它們?cè)诮唤缑鍭B的Zone 1側(cè)的平均相對(duì)誤差分別為0.016%和0.007%。Model-Yeh-F在交界面處依然無法滿足折射定律，具有較大的誤差，其對(duì)Zone 1 側(cè)解析解的誤差為45%。從圖7（b）中可以看出，MS?FEM-D和Model-Zhou-F在交界面附近依然可以取得精確的解，而Model-Yeh-F在AB、CD兩個(gè)交界面附近均出現(xiàn)了“補(bǔ)償”現(xiàn)象，具有較大的誤差。在遠(yuǎn)離交界面處，三種方法均和參照解重合，具有較高精度。

圖7 達(dá)西流速Vx的參照解以及MSFEM-D和Model-Zhou-F所獲的數(shù)值解

圖8展示了交界面CD處達(dá)西流速Vx（圖8（a））和Vy（圖8（b））的參照解和MSFEM-D、Model-Zhou-F和Model-Yeh-F所獲的數(shù)值解。如圖8所示，CD為傾斜交界面，不平行于任何坐標(biāo)軸，故達(dá)西流速Vx和Vy均在CD上不連續(xù)，具有兩個(gè)不同的值。本圖中，MSFEM-D、Model-Zhou-F 能夠滿足折射定律，且均能接近參照解，它們?cè)赯one 2側(cè)的Vx的平均相對(duì)誤差分別為7.3%和1.6%（圖8（a））。和上例相比，兩種方法在交界面處的誤差均有增長，這是由于Vx在邊界處變化較快而引起的，可以通過加密網(wǎng)格提高精度。由于Model-Zhou-F是直接在全局細(xì)尺度網(wǎng)格（48×48×2）上進(jìn)行達(dá)西流速求解，而MS?FEM-D則是在（24×24×2）的粗尺度網(wǎng)格上進(jìn)行達(dá)西流速求解；Model-Zhou-F能夠獲得高于MSFEM-D的精度，但所需消耗也遠(yuǎn)大于MSFEM-D，為其的200倍以上。Model-Yeh-F 無法滿足折射定律，在交界面CD處具有較大的誤差，其對(duì)Zone 2側(cè)Vx參照解的平均相對(duì)誤差為28.3%（圖8（a））。

圖8 達(dá)西流速Vx和Vy的參照解以及MSFEM-D、Model-Zhou-F和Model-Yeh-F所獲的數(shù)值解

在本例中，MSFEM-D 僅消耗不到1 s 的CPU 時(shí)間即可完成粗、細(xì)尺度的水頭和達(dá)西流速的模擬，且精度接近Model-Zhou-F；而Model-Zhou-F 則需要迭代14次，共需要207 s才能完成相同的工作；Model-Yeh-F雖然僅需要18 s CPU時(shí)間，但其無法保證達(dá)西流速在交界面上符合折射定律。這一結(jié)果顯示，MSFEM-D具有較高的計(jì)算效率，并能適用于傾斜交界面和多條交界面的情況。

4 結(jié)論與展望

本文提出分區(qū)多尺度有限元法（MSFEM-D）用于模擬交界面處的地下水達(dá)西流速。在構(gòu)造水流和達(dá)西流速方程時(shí)，該方法以MSFEM的框架為基礎(chǔ)，進(jìn)行了尺度提升，提高了計(jì)算效率。在模擬交界面達(dá)西流速時(shí)，該方法運(yùn)用區(qū)域分解技術(shù)分離不同介質(zhì)避免互相影響，并運(yùn)用Yeh 的伽遼金模型保證達(dá)西流速的連續(xù)性，結(jié)合Jump函數(shù)，從而令交界面處的達(dá)西流速符合折射定律。

多種情形下含交界面的地下水?dāng)?shù)值模擬實(shí)驗(yàn)結(jié)果表明：（1）MSFEM-D 不僅能夠在含單條、多條、水平、傾斜交界面的地下水問題中保證交界面處達(dá)西流速滿足折射定律，還可以處理滲透系數(shù)為非均質(zhì)各向異性的情況。（2）MSFEM-D 的全局水頭精度與FEM-F 相近，全局達(dá)西流速的精度與Model-Zhou-F接近，其在交界面處的達(dá)西流速比Model-Yeh-F更加精確。（3）MSFEM-D無需任何迭代過程，其計(jì)算消耗低于Model-Yeh-F，更是遠(yuǎn)低于需要迭代的Model-Zhou-F，能夠節(jié)約大量的計(jì)算成本。

在今后的工作中，將繼續(xù)改進(jìn)MSFEM-D，研究含交界面的三維地下水問題的模擬方法；結(jié)合鉆孔和地質(zhì)圖等實(shí)際信息，提出交界面的描述方法，研究實(shí)際案例的模擬方法；開發(fā)MSFEM-D軟件，并與現(xiàn)有的FEFLOW常用有限元軟件進(jìn)行對(duì)比。