張鵬偉 費宇紅* 郝奇琛 李亞松 朱玉晨 孟素花 郭春艷

（1、中國地質(zhì)科學(xué)院水文地質(zhì)環(huán)境地質(zhì)研究所，河北 石家莊 050061 2、福建省水循環(huán)與生態(tài)地質(zhì)過程重點實驗室，福建 廈門 361001）

地下水是重要的自然資源，作為許多地方工業(yè)、農(nóng)業(yè)和居民生活主要水源或唯一供水水源，其水量與水質(zhì)變化備受社會各界關(guān)注。地下水?dāng)?shù)值模擬是研究地下水流場變化、溶質(zhì)運移和熱量運移的重要方法，目前求解地下水模型的數(shù)值方法包括有限差分法(FDM)、有限元法 (FEM)、邊界元法 (BEM) 和有限分析法(FAM)。MODFLOW 作為有限差分地下水流數(shù)值模型的代表，是目前水資源利用、環(huán)境保護等科研和生產(chǎn)領(lǐng)域最為普及的計算程序。空間離散化是構(gòu)建地下水?dāng)?shù)值模型的重要步驟，MODFLOW 為其提供了結(jié)構(gòu)化網(wǎng)格與非結(jié)構(gòu)化網(wǎng)格兩種剖分方式，為深入理解非結(jié)構(gòu)化版本（MODFLOW-USG） 與 傳 統(tǒng) 結(jié) 構(gòu) 化 版 本 ( 包 括MODFLOW-2000，-2005，-NWT，-LGR) 的區(qū)別和特點，本文針對以上兩種網(wǎng)格剖分方式從求解方法、連接方式以及單元格幾何形狀等方面進行分析對比。

1 求解方法

MODFLOW-USG 基于控制體積有限差分法（Control Volume Finite-difference, CVFD）進行非結(jié)構(gòu)化網(wǎng)格剖分，該方法綜合有限差分法和有限元法的優(yōu)點，采用差分法進行離散，又采用不規(guī)則的網(wǎng)格形狀處理模型邊界和網(wǎng)格疏密的問題，從而保證了流場維度的靈活性[1]。對于單元n 與每一個與之相鄰的單元m，控制體積有限差分方程的一般形式為[2]：

式中：Cnm為單元n 與相鄰單元m 之間的系數(shù)；hm和hn分別為m 單元和n 單元的水頭；HCOFn是所有與單元n 水頭hn相關(guān)的系數(shù)項的總和；RHSn為所有常數(shù)項和已知項，也稱平衡方程右側(cè)項。

式中：anm是單元n 與單元m 公共面在垂直于流向方向上的投影面積；knm為單元n 與單元m 之間的滲透系數(shù)；Lnm、Lmn為單元n 和單元m 中心到公共面的垂直距離；Vn為單元n 的體積；SSn為單元的儲水率；?t 為時間步長；t-1 為上一時間步長。（圖1）

圖1 非結(jié)構(gòu)化網(wǎng)格空間離散示意圖

傳統(tǒng)結(jié)構(gòu)化網(wǎng)格剖分所采用的有限差分方程同樣是基于質(zhì)量守恒和水量平衡方程構(gòu)建，對于單元(i,j,k)與相鄰單元(i,j-1,k)，二者間的流量計算方程如下[3]：

式中：hi,j,k是單元(i,j,k)的水頭；qi,j-1/2,k是兩相鄰單元間的流量；Δci和Δvk分別為單元的寬和高；Δrj-1/2為兩單元中心間距。對于單元(i,j,k)有限差分的一般計算方程為：

式中：KR、KC、KV 分別為沿j、i、k 方向的滲透系數(shù)；QSi,j,k為各類源匯項；SSi,j,k為單元儲水率。（圖2）

圖2 結(jié)構(gòu)化網(wǎng)格空間離散示意圖

2 連接方式

結(jié)構(gòu)化網(wǎng)格剖分是指每個內(nèi)部單元格具有相同數(shù)量的毗鄰單元，除了在邊界處的單元以外，每個單元通常固定連接到6 個主方向上的相鄰單元（圖2）。這種連接方式使其具有網(wǎng)格生成速度快、數(shù)據(jù)結(jié)構(gòu)簡單的優(yōu)點，同時方便數(shù)據(jù)和相關(guān)參數(shù)的準備，便于賦值，最重要的是計算穩(wěn)定性好并且收斂速度快。然而在實際模擬過程中，結(jié)構(gòu)化網(wǎng)格剖分存在明顯的局限性：在刻畫不規(guī)則邊界時，需要將研究區(qū)以外的單元設(shè)置為不活動單元，且研究區(qū)邊界往往容易形成鋸齒形邊界（圖3 (b)）；若采用網(wǎng)格加密的方式細化邊界或重點區(qū)域，加密區(qū)域行和列會擴展至模型域的邊緣，從而產(chǎn)生許多非必要的計算單元，大大影響了模型的計算效率并且容易造成模型不收斂的問題。

圖3 結(jié)構(gòu)化網(wǎng)格剖分方式

非結(jié)構(gòu)化網(wǎng)格則是指對于所有單元來說，網(wǎng)格單元與相鄰單元的連接數(shù)不是固定值，每一個單元可以與任意數(shù)量的網(wǎng)格相鄰。如圖4 所示，一個單元與周圍單元的連通性取決于公共邊（面）的數(shù)量，而每個單元的不盡相同，這就導(dǎo)致所列的方程組是非結(jié)構(gòu)化的。非結(jié)構(gòu)化網(wǎng)格允許定義不規(guī)則形狀的模型域，在模型區(qū)域之外，不需要設(shè)置不活動單元。

圖4 非結(jié)構(gòu)化網(wǎng)格剖分方式

3 單元格幾何形狀

結(jié)構(gòu)化網(wǎng)格剖分通常對模型進行正交矩形剖分，但以單元幾何形狀區(qū)分結(jié)構(gòu)化與非結(jié)構(gòu)剖分是常見的錯誤認知，需要注意的是規(guī)則三角形或六邊形網(wǎng)格同樣是結(jié)構(gòu)化的，關(guān)鍵在于其每個單元具有相同數(shù)量的毗鄰單元。對于結(jié)構(gòu)化網(wǎng)格剖分，模型層必須在整個模型領(lǐng)域中保持連續(xù)，所有模型層中的單元數(shù)必須相同，這意味著需要插入完整的模型層，且必須在弱透水層尖滅區(qū)域設(shè)置最小厚度。結(jié)構(gòu)化網(wǎng)格單元編號根據(jù)網(wǎng)格的行、列、層屬性編號。

如圖4 所示，非結(jié)構(gòu)化網(wǎng)格剖分允許單獨或組合使用嵌套網(wǎng)格、三角形、矩形或其他多邊形形狀的幾何圖形以適當(dāng)?shù)仉x散區(qū)域，從而有效提高河流和井等對象的模型分辨率。非結(jié)構(gòu)化網(wǎng)格單元體必須是凸柱狀體，且柱體側(cè)面垂直于水平底面，但各層無需在整個模型區(qū)域保持連續(xù)，模型層厚度可設(shè)置為零。各層可獨立離散化，比如對淺層河流處細化剖分，深層保持粗化網(wǎng)格。為方便預(yù)處理和后處理，MODFLOW-USG 保留了分層的概念，但單元格編號不顯示分層信息，單元從上層開始編號，然后依次按網(wǎng)格層級順序編號至最底層單元。

對于計算精度而言，標準的CVFD 方程要求連接相鄰單元中心點的連線垂直于公共面，并且連線與公共面的交點應(yīng)處于公共面的中心[4]。但嵌套網(wǎng)格、四叉樹網(wǎng)格以及具有非正多邊形網(wǎng)格可能會違背此要求，且偏離CVFD 法的要求越大，地下水流數(shù)值解的精度損失也就越大；然而增加細化網(wǎng)格的數(shù)量，提高模型分辨率又會提高模型精度，降低誤差，二者之間的關(guān)系很難量化。有不少學(xué)者已通過理論推導(dǎo)以及模擬實驗的方式已經(jīng)證明，分辨率越高，模型對水位模擬的精度越高，并且相比較于細結(jié)構(gòu)化的網(wǎng)格剖分，其計算時長也會大大降低[5-6]。

4 結(jié)論

總結(jié)來說，結(jié)構(gòu)化網(wǎng)格剖分具有操作簡便、生成速度快、數(shù)據(jù)結(jié)構(gòu)簡單等優(yōu)點，但在研究重點區(qū)域加密、模擬不規(guī)則研究區(qū)邊界以及模擬不連續(xù)弱透水層或地層尖滅等方面具有一定局限性。

非結(jié)構(gòu)化網(wǎng)格剖分允許刻畫不規(guī)則邊界，支持局部加密，計算效率及計算精度更高，穩(wěn)定性更強。