非結構化與結構化網格剖分在地下水數值模擬中對比分析

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

（1、中國地質科學院水文地質環境地質研究所，河北 石家莊 050061 2、福建省水循環與生態地質過程重點實驗室，福建 廈門 361001）

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

1 求解方法

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

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

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

圖1 非結構化網格空間離散示意圖

傳統結構化網格剖分所采用的有限差分方程同樣是基于質量守恒和水量平衡方程構建，對于單元(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 方向的滲透系數；QSi,j,k為各類源匯項；SSi,j,k為單元儲水率。（圖2）

圖2 結構化網格空間離散示意圖

2 連接方式

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

圖3 結構化網格剖分方式

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

圖4 非結構化網格剖分方式

3 單元格幾何形狀

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

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

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

4 結論

總結來說，結構化網格剖分具有操作簡便、生成速度快、數據結構簡單等優點，但在研究重點區域加密、模擬不規則研究區邊界以及模擬不連續弱透水層或地層尖滅等方面具有一定局限性。

非結構化網格剖分允許刻畫不規則邊界，支持局部加密，計算效率及計算精度更高，穩定性更強。