徐昇氏網格應用於二維地下水流數值模式之建立

國 立 交 通 大 學

土 木 工 程 研 究 所

碩 士 論 文

徐昇氏網格應用於二維地下水流數值模式

之建立

Using Voronoi Diagram on Development of 2D

Groundwater Numerical Model

研 究 生 ： 易正偉

指導教授 ： 張良正 博士

徐昇氏網格應用於二維地下水流數值模式之建立

Using Voronoi Diagram on Development of 2D

Groundwater Numerical Model

研

究

生 ：易正偉

Student : Cheng-Wei Yih

指 導 教 授：張良正 博士

Advisor : Dr. Liang C. Chang

國 立 交 通 大 學

土 木 工 程 學 系 碩 士 班

碩 士 論 文

A Thesis

Submitted to Institute of Civil Engineering

National Chiao Tung University

in Partial Fulfillment of the Requirements

for the Degree of

Master of Science

in

Civil Engineering

July 2008

Hsinchu, Taiwan, Republic of China

徐昇氏網格應用於二維地下水流數值模式之建立

學生：易正偉 指導教授：張良正 博士 國立交通大學土木工程研究所

摘要

數值模式已成為分析地下水相關問題之主要工具，而發展數值 模式常用的數值方法有有限差分法、有限元素法與有限體積法等，有 限差分法的優點在於簡單直接且計算效率高，然其缺點在於格網形狀 僅可是矩形而缺乏彈性，有限元素法雖有網格形狀上的彈性，但其理 論較複雜且計算效率較低，有限體積法之理論複雜度及且計算效率介 於兩者之間，其優點為所得數值解符合質量守恆。因此若能延續有限 差分法的簡單及高效率，克服矩形網格限制，整合有限體積法符合守 恆條件的精神，則可發展出更簡單、精確、高效率及高應用彈性的數 值模式。 本研究的目的在以徐昇氏網格為基礎，配合守恆定律，針對地 下水流問題發展一新的地下水流數值模式。本研究藉由 CGAL 函式 庫建立徐昇氏網格，並以各徐昇氏網格為控制體積，配合相關之方程 式進行守恆運算，再以迭代方式進行數值求解。本研究以十一個案例 驗證模式之正確性，及展現模式應用的彈性如局部加密等。結果顯示 本模式在各案例中均可符合質量守恆的要求。此外，配合徐昇氏網格 之局部格網加密，可在增加少量計算時間的情形下，大幅提昇計算精 度，有效提升模式應用彈性及計算效率。

Using Voronoi Diagram on Development of 2D

Groundwater Numerical Model

Student：Cheng-Wei Yih Advisor：Dr. Liang-Cheng Chang

Department of Civil Engineering

National Chiao Tung University

Abstract

Numerical model simulation technique has become a powerful tool for analyzing groundwater problems. Although several numerical methods, such as the finite difference method (FDM), the finite element method (FEM) and the finite volume method (FVM), are widely used for numerical model development the FDM has the advantage of easy implementation and high computation efficiency. However, because the FDM requires a rectangular grid, solving problems with irregular boundary is difficult, thus limiting its practical application. Unlike the FDM, grid shape of the FEM is flexible and easily applies to solve problems with irregular boundary. However, the FEM is more complicated and computationally less efficient than the FDM. Easy implementation and computational efficiency for the FVM lies between the FDM and FEM. The FVM advantage is its guaranteed solution to satisfy conservation conditions. Regarding the previous discussion, this research proposes a novel numerical method, easy to implement, with high computational efficiency and a flexible grid shape for high practical application.

The current study adopts the grid for the proposed groundwater numerical mode that applies the Voronoi diagram and the Computation

geometry algorithm library (CGAL) to generate the Voronoi grid in the model. For each Voronoi cell (grid), the mass conservation of water and other relevant physical laws form the relationship between a cell and its neighboring cells. The proposed model applies iteration methods for computing solutions cell by cell. Using the proposed model to verify model capability solves several hypothesis problems. Solving problems with point source demonstrates that applying the irregular Voronoi grid uses only one tenth of that using the rectangular grid for similar solution accuracy.

謝 誌

感謝張良正老師這幾年來對於本研究的悉心指導，並且老師交予 學生的求學態度更使學生獲益良多。也承蒙口試委員蔡惠峰副主任、 李光敦教授及李天浩副教授仔細地審閱此文，在口試中提出審閱意 見，使本文更趨於完備，於此致上萬分感謝。 在交大的這幾年，感謝研究室內的學長對我的照顧，特別是陳宇 文學長，從研究之初的系統架構規劃至研究中的程式設計到研究末的 案例驗証部份，各個階段都給予我非常多的指導提點與協助，使本研 究能在限定的時間內順利完成。也感謝蔡瑞彬學長常會與我討論地下 水的相關問題並提供意見，使研究能持續在正確的方向，感謝葉明生 學長與陳鴻輝學長常常鼓勵我，何智超學長總是會回答我關於研究相 關的疑問，朱宏杰學長常提醒我報告中不完善的部分，希望各個厲害 的學長都能早日拿到博士學位。 感謝碩士班的學長浚瑋，在學期間總是鼓勵我與教導我許多這個 領域的相關知識，感謝學弟為善，對於本研究做了相當多的協助，使 得本研究可以持續推進。也感謝碩士班的明穎學長、文鴻學長、昀柏 學長、任馥學長、君儀學長、怡釧、婉茹，昱維、柏成、韋圻、翰聖、 敏威學弟以及程翔助理兼同學，你們總是讓研究室的氣氛很熱鬧很有 活力。 感謝海斌、韻玲還有黃碩傑、陳雅婷、謝鳳書、讓我在想要休息 的時候，有一起從事休憩活動的夥伴，讓求學生涯過得更有色彩。 最後要感謝最關心我的家人，爸媽與哥哥對我的關愛，在我求學 的期間總是鼓勵著我，即使在碩士班的求學時間過長也是一直支持著 我，使我沒有後顧之憂，可以專心於學業。感謝韻嫻妳長久以來的永 不間斷的關愛與支持，求學期間的每一天都有妳在我身邊陪伴我，讓

我的心中總是充滿著溫暖，使我對於未來有著更多的憧憬以及更大的 動力，謝謝你們。

章 節 目 錄

中文摘要... Ⅰ 英文摘要... Ⅱ 謝誌... Ⅳ 章節目錄... Ⅵ 表目錄... Ⅷ 圖目錄... Ⅸ 一、 緒論 ... 1 1.1 前言 ... 1 1.2 研究目的 ... 2 1.3 文獻回顧 ... 2 1.2.1 地下水模式之演進 ... 2 1.2.2 Voronoi Diagram 之演進 ... 4 1.2.3 Voronoi Diagram 結合於數值模式 ... 5 二、 研究步驟 ... 6 三、 研究理論與方法 ... 8 3.1 Voronoi Diagram 空間切割法... 8 3.1.1 徐昇氏網格(Voronoi Diagram)簡介 ... 8 3.1.2 CGAL 函式庫簡介 ... 10 3.1.3 封閉空間之 Voronoi Diagram 函式庫開發 ... 10 3.2 數值演算法 ... 12 3.2.1 計算空間架構... 14 3.2.2 方程式離散化... 16 3.2.3 內迭代處理方法... 17 3.2.4 外迭代處理方法... 19

3.2.5 整體數值模擬流程... 20 四、 垂向二維地下水流數值模式開發 ... 22 五、 垂向二維地下水流數值模式模擬與驗證 ... 27 5.1 穩態規則網格案例驗證與模擬 ... 27 5.2 暫態規則網格案例驗證與模擬 ... 43 5.3 暫態不規則網格案例驗證與模擬 ... 58 六、 結論與建議 ... 64 6.1 結論 ... 64 6.2 建議 ... 64 參考文獻... 參 1 附錄 土壤壓密係數與傳統方法之關聯推導 ... 附 1

表 目 錄

表 5.1-1 解析解與數值解於左右兩側邊界穿越流量比較表... 29 表 5.1-2 外迭代收斂標準、計算效率與守恆誤差結果表... 30 表 5.2-1 完全隱示法相對系統守恆誤差表 ... 45 表 5.2-2 相對系統守恆誤差表(案例 5.2-1 與 5.2-4)... 56 表 5.3-1 案例差異說明表 ... 59 表 5.3-2 案例模擬結果比較表 ... 59

圖 目 錄

圖 2-1 研究步驟流程圖 ... 7 圖 3.1.1-1 Voronoi Diagram 與其相關圖形... 9 圖 3.1.3-1 不同邊界框形成之 Voronoi Diagram... 11 圖 3.2.1-1 空間架構說明示意圖 ... 15 圖 3.2.1-2 節點與連結關係示意圖 ... 15 圖 3.2.2-1 網格示意圖 ... 17 圖 3.2.3-1 內迭代流程圖 ... 19 圖 3.2.4-1 外迭代流程圖 ... 20 圖 3.2.5-1 整體數值模擬流程圖 ... 21 圖 5.1-1 案例 5.1-1 格網與模式配置圖... 30 圖 5.1-2 案例 5.1-1 總水頭等值圖... 31 圖 5.1-3 案例 5.1-1 壓力水頭等值圖... 32 圖 5.1-4 外迭代收斂標準與對數守恆誤差趨勢圖... 33 圖 5.1-5 案例 5.1-2 格網與模式配置圖... 35 圖 5.1-6 案例 5.1-2 總水頭等值圖... 36 圖 5.1-7 案例 5.1-3 格網與模式配置圖... 38 圖 5.1-8 案例 5.1-3 總水頭等值圖... 39 圖 5.1-9 案例 5.1-4 材質設定圖... 41 圖 5.1-10 案例 5.1-3 與案例 5.1-4 之總水頭等值圖 ... 42 圖 5.2-1 抽水點各時刻總水頭變化圖 ... 45 圖 5.2-2 抽水井所在網格之瞬時流率變化圖... 46 圖 5.2-3 抽水井所在網格之流量變化 ... 46 圖 5.2-4 水位剖面上視圖 ... 47

圖 5.2-5 不同數值技巧之水位剖面比較圖 ... 48 圖 5.2-6 完全隱示法之邊界流量變化圖 ... 48 圖 5.2-7 抽水井抽水量變化圖 ... 50 圖 5.2-8 抽水點各時刻總水頭變化圖 ... 50 圖 5.2-9 脈衝函數示意圖 ... 51 圖 5.2-10 不同數值技巧之水位剖面比較圖... 51 圖 5.2-11 案例 5.2-3 水位剖面上視圖... 53 圖 5.2-12 案例 5.2-1 與 5.2-3 之水位剖面比較圖 ... 54 圖 5.2-13 案例 5.2-3 邊界流量變化圖... 54 圖 5.2-14 案例 5.2-4 格網與模式配置圖... 56 圖 5.2-15 水位剖面比較圖(案例 5.2-1 與 5.2-4)... 57 圖 5.2-16 邊界流量變化圖(案例 5.2-1 與 5.2-4)... 57 圖 5.3-1 案例 5.3-1 格網與模式配置圖... 60 圖 5.3-2 案例 5.3-2 格網與模式配置圖... 61 圖 5.3-3 案例 5.3-3 格網與模式配置圖... 62 圖 5.3-4 規則網格與不規則網格水位剖面比較圖... 63

1

一、緒論

1.1 前言

在數值模擬的領域內，基本的建構流程，均是先確立問題的概 念模式，再建立其對應的數學模式，接著再應用如有限差分、有限元 素或其它數值方法得出該數學模式的離散方程式，最後求解此離散方 程式，得出該問題的數值解。由於數值模式大多是以有限的節點代表 問題的整個解空間，因此許多的空間切割法因而產生，從基本的方形 網格、矩形網格演進至三角格網及其他格網形狀。而求解以節點上的 變數來表示的離散方程式，傳統上常見的方式有兩種，即是矩陣求解 或迭代求解。其中矩陣求解為將各運算節點之方程式組合成為矩陣方 程式，再以矩陣求解法，如高斯消去法進行求解，而矩陣求解因應節 點方程式的編排順序與特性，有各種增進時間效益的求解方法，如帶 狀矩陣求解法等；迭代求解則不將各運算節點之方程式組合成為矩陣 方程式，而每次僅以各節點之方程式更新該節點變數之值，更新時周 圍節點上之變數則採用其當下之值，如此重複更新各節點變數之值直 至達到設定之收斂標準為止，此類迭代求解也有各種的改良方式如以 內插為基礎的連續鬆弛法(Successive Over/Under Relaxation)。

解地下水流問題常用的數值方法有有限差分法(Finite Difference Method, FDM)、有限元素法(Finite Element Method, FEM)、有限體積 法(Finite Volume Method, FVM)及邊界元素法(Boundary Element Method, BEM)等，其中過去最常使用的兩種方法為有限差分法與有限 元素法，有限差分法的優點在於簡單直接且計算效率高，因此易於實 作，目前最被普遍應用的 Modular three-dimensional groundwater flow model., (MODFLOW) (1988) 模式即是應用有限差分法，以往有限差 分法最大的缺點，在於其格網形狀僅可是矩型而缺乏彈性，造成許多

實際應用上的不便或精度上的降低。相較於有限差分法，有限元素法 的最大優點則在於其網格形狀上的彈性，理論上可以是任意的多邊 形，可適應不規則邊界與需局部加密之情況，惟相對而言，其理論較 複雜且計算效率較低。有限體積法是針對控制方程式進行空間與時間 上之積分，再以此積分形式之控制方程式配合差分與空間內插方法， 建立有限體積法之離散方程式，其優點是所得數值解必定符合守恆條 件，是近年來較為流行的數值方法，惟目前多數以有限體積法發展的 數值模式多仍採用矩形網格，仍有如有限差分模式般在實際應用上的 不便。因此，若能延續有限差分法的簡單及高效率，克服矩形網格的 限制，並整合有限體積法符合守恆條件的精神，應可發展出更簡單、 精確、高效率及高應用彈性的數值模式。

1.2 研究目的

本研究的目的在以徐昇氏網格為基礎，配合守恆定律，針對地 下水流問題發展一新的地下水流數值模式。相對於有限體積法，本研 究發展的模式，概念上不需先進行控制方程式的積分，而是直接以徐 昇氏網格表達之控制體積，配合相關方程式進行(質量)守恆運算，因 此基本概念相對簡單，且可解決矩形網格實際建模上的不便，有相當 高的推廣應用價值。

1.3 文獻回顧

本節將介紹地下水模式的演進、Voronoi Diagram 的發展以及結 合上述兩者之研究。 1.3.1 地下水模式之演進

在地下水領域內，許多地下水數值模式曾先後被建立，美國地 質測量局(USGS)發展的地下水相關數值模式多是以有限差分法與有 限元素法為主。在有限差分法的發展上，由 Trescott, Pinder, and Larson (1976) 發展出以迭代法為求解方式的二維飽和地下水有限差分模 式，而後 Larson (1978) 針對此數值解法部分建立出四點矩陣直接求 解法以取代先前的迭代求解法減少模擬時間，Manteuffel, Grove, and Konikow (1983) 建立出共軛梯度法，增加了對求解稀疏矩陣的效率， 此處可看出在電腦處理速度並未如現今快速的年代，許多研究均投入 在增進求解效率。Ozbilgin and Dickerman (1984) 則使用了曼寧公 式，將地表水入滲加入了二維地下水有限差分模式的計算，Weeks (1 978) 使用以有限個二維侷限層(Confining Layers)切割含水層的概 念，建立擬三維的有限差分模式。

三維飽和地下水有限差分模式由 Trescott (1975) 發表，發表的 時間早於該作者二維之發表時間，此模式為 1970 至 1980 間的主要模 式之一。而後直到 McDonald and Harbaugh (1988) 發表了模組化的有 限差分模式即 MODFLOW，此模式為最被廣泛應用之地下水模式之 一，並且功能完備，可進行三維之地下水模擬，甚至拘限與非拘限含 水層混合的多層地下含水層系統也可模擬，且除地下水系統本身外， 區域補注、蒸發散、抽水井抽水、區域排水與河床水位變動都可模擬， 其後許多研究均以 MODFLOW 為基礎，如 Kuiper (1987) 於 MODFLOW 三維模式中加入前調式共軛梯度法 (Preconditioned Conjugate Gradient Method, PCG)，加速了 MODFLOW 的求解速度， Leake and Prudic (1991) 則進行因抽水而造成含水層之下陷的研究與 模擬。

輸模式，可模擬吸附與不可逆反應等問題，Konikow and Bredehoeft (1978)、 Goode and Konikow (1989) 推展出的二維污染傳輸模式，除 吸附外,尚可模擬傳流、延散與離子交換等現象。

有限元素法的發展則由 Dunlap, Lindgren and Carr (1984) 使用葛 勒金法 (Galerkin Method) 為求解方法建立出地下水之二維有限元素 模式，其後 Lewis, Voss and Rubin (1986) 發表了二維污染傳輸有限元 素模式，可解決傳流、延散與吸附等物理現象與離子交換與化學反應 平衡等化學現象。Glover (1988) 建立了模擬河流的動態方程式， Cooley (1992)，Torak (1993)則模擬了暫態滲流的情形。三維的有限元 素模式則是由 Lin, Richards, Yeh, Cheng & Jones (1996) 發展出，其模 式稱為 FEMWATER，可模擬飽和、不飽和土壤內之地下水流、不定 邊界條件、海水入侵以及污染傳輸等問題，為目前最主要的有限元素 模式。 上述之模式均必須將各方程式整合出一控制方程式，才能給進 模式求解，而本研究於數值演算法上，以方程組運算並透過最佳化迭 代的方式，搜尋出最符合各點當時刻質量守恆的基本變數。此方法的 優點在於，僅需對守恆方程式演譯出各個組成方程式，各個組成方程 式與守恆方程式即可稱為一控制方程組。此方程組不需人工進行前置 之整合方程式的動作，更替內部其中任一方程式都不需重新進行整合 方程式的處理，使得在方程式的擴充上有顯著的便利性。 1.3.2 Voronoi Diagram 之演進 Voronoi Diagram 為一種自然的切割方式，代表了各節點控制的 區域內，區域中的其餘點必最靠近該控制節點，而區域間的交介面， 必定與控制節點間的連線是正交的，在數學的領域上，這是已被研究

許久的一個圖形，Voronoi Diagram 此圖形最早在 1644 年就已由法國 哲學家笛卡兒 (René Descartes) 所建構，而後由 Dirichlet (1850)在他 自己對於二次多項式研究中，使用了 Voronoi Diagram，再由 Voronoi (1907) 針對此圖形進行更深入的研究，將此圖形延伸往更高的維 度。而後此圖被用於各個領域上，如 Portela (2004) 以 Voronoi Diagram 界定出無線網路基地台的選擇，Overholt (2007) 使用 Voronoi

Diagram 作為機器人路徑規劃的依據，而在氣象與地理學之上， Thiessen (1911)將 Voronoi Diagram 用於區域降雨量的推估，因此在水 文領域也稱此圖形為徐昇式多邊形 (Thiessen＇s Polygon)。Sukumar (2003)將 Voronoi Diagram 應用於延散方程式內之有限差分法，使 Voronoi Diagram 可以做為空間差分之網格。

1.3.3 Voronoi Diagram 結合於數值模式

近年來，已有研究將 Voronoi Diagram 做為數值模式內空間切割 的方法。Gregory E. Tucker (2001)將 Voronoi Diagram 應用於水文地質 模式，針對 Voronoi Diagram 進行探討，將 Voronoi Diagram 內網格間 的連線可視為可能之流線，並在數值解法上使用 Voronoi Polygons 為 有限體積的區塊(cell)進行質量守恆求解。而 Sukumar (2003) 則將 Voronoi Diagram 結合有限差分法發展數值延散模式，並完整推導對 應之有限差分式。A. Iske and M. Käser(2004)以應用 Voronoi Diagram 建立可適性非結構網格，其計算節點會依據傳流項(advection)而移 動，該研究亦針對移動前後的質量守恆問題進行改善。

本研究將以徐昇氏網格 (Voronoi Diagram) 作為垂向二維地下 水流數值模式之空間切割法，以此空間切割法進行局部加密，展示徐 昇氏網格的優點與適用性。

2

二、研究步驟

本研究之研究步驟如圖 2-1 所示，本研究首先開發新型態地下水 模擬模式，並藉由地下水流模擬案例驗証本模式之正確性。 一、 模擬模式開發： 在模擬模式開發方面，可以再進一步細分為「空間切割」與「節 點計算」兩個步驟。在空間切割方面，本研究採用徐昇氏網格(即 Voronoi Diagram)進行模擬格網切割，透過徐昇氏網格，節點與周邊 節點間所形成的線段必定與所在網格邊界垂直，此一部份的功能實做 則藉由 CGAL(Computational Geometry Algorithm Library)函式庫。

在「節點計算」方面，本研究將各個徐昇氏網格視為控制體積， 進行地下水流進出質量平衡，穿越各控制體積邊界的通量則以各節點 上的變量(水頭)所得之差分式近似計算之，藉由此質量守恆可建立各 徐昇氏網格內節點與鄰近節點的關係式，基於此關係式進一步以迭代 方式逐步求解整個模擬區域之解。 二、 模式驗證： 接著本研究設計多個地下水模擬案例，藉由模擬結果驗證模式正 確性。

3

三、研究理論與方法

本研究運用「Voronoi Diagram」空間切割法建立地下水流數值 模式，因此 3.1 節將先針對 Voronoi Diagram 進行說明介紹，之後在 於 3.2 節介紹本研究之數值模擬方法。

3.1 Voronoi Diagram 空間切割法

空間切割在數值模擬問題中，為必要的一個步驟，不同數值方 法則有不同之切割方法，以最常見的有限差分法(FDM)為例，其必須 侷限於矩形網格，因此對於不規則形狀之掌握較差。本研究在此採用 Voronoi Diagram 作為切割演算法。 3.1.1 徐昇氏網格(Voronoi Diagram)簡介

Voronoi Diagram 是透過散佈於平面上一群平面點位資料(sites or points)為基礎所切割出來的圖形，其圖形如圖 3.1.1-1 所示。Voronoi Diagram 具備下列特性： 1. 各點(sites or points)所必對應一個控制區塊(cell)，區塊內之任何一 點與該節點之距離，必定小於此點與其他節點之距離。 2. 控制區塊形狀必為凸包(convex hulls)。 3. 兩節點之連線必定正交於其控制區塊(cells)之邊，此邊稱為控制區 塊邊(edge) 4. 控制區塊邊(edge)之交接點稱為角點(vertex)，令任意角點為圓心， 必可找到一圓通過相鄰之節點(如圖 3.1.1-1 所示)，該圓稱為空圓 (empty circle)。 5. 將任意空圓上之所有節點連線，可以繪製成三角形，其稱為 Delaunay Triangles，其與 Voronoi diagram 為對偶圖形。圖 3.1.1-1

上，節點 0、1 與 4 之連線所形成之三角形即為 Delaunay Triangles。 圖 3.1.1-1 Voronoi Diagram 與其相關圖形 數學定義：若給定一組點位資料(N ={N₁,N₂,N_3LLN_n})，其中 為第 m 個節點，n 為空間中節點數量，其數值應為 m N ∞ ≤ ≤ n 2 。其控制 區塊(Voronoi Diagram)可定義如下：

( )

N

{

D D N D N for j i

}

V i = : − i ≤ − j ∀ ≠ ... (3.1.1-1) 公式描述了 Voronoi Diagram 最重要的特性，就是控制區塊內的 任何點必定與其對應之節點距離最近。代表該節點在空間上，比任何 其他節點更能代表該區塊。因此在水利或氣象領域中，Voronoi Diagram 被應用於將點位形式之雨量資料推估集水區內之總降雨量。

在數值運算上，由於兩節點連線必正交於其控制區塊邊，因此對於穿 越格網邊界之流量估算上，無須進行角度轉換之向量計算。

Voronoi Diagram 發展至今已有相當多的演算法計算，如 half plane intersection、Divine-and-Conquer、Fortune’s Algorithm 等等，其 中 Steve Fortune 於 1986 年發表的掃線演算法，由於在最差的狀況下 也有 的計算複雜度，為建構 Voronoi Diagram 最有效率的演 算法之一。 ) log (n n O 3.1.2 CGAL 函式庫簡介

CGAL 全名為 Computational Geometry Algorithm Library，是計 算空間幾何上極具可信度的函式庫，其中包含了許多圖形的資料結 構、建構圖形的演算法、格網產生與處理與圖形分析等功能。為 Geometry Factory 所開發的開放原始碼專案，其目的為 C++環境，提 供有效率並且可信賴的空間計算演算法。目前已經被廣泛應用於電腦 繪圖、科學視算、設計模擬、格網生成、數值方法與微生物模擬等各 項領域上。

然而 CGAL 函式庫提供之 Voronoi Diagram 相關功能，僅針對無 限空間之 Voronoi Diagram 空間切割，亦即外圍節點之控制區塊 (Voronoi Cells)並不為封閉空間。然而本研究是應用於數值模擬上， 網格必須是有限個封閉空間，因此不完全符合需求。 因此本研究是以 CGAL 函式庫作為基礎函式庫，透過物件封裝 技巧，進一步開發適合本研究之封閉空間 Voronoi Diagram，使外圍 之控制區塊與模擬邊界進行切割計算，重新輸出適合本研究數值模擬 之空間切割結果。 3.1.3 封閉空間之 Voronoi Diagram 函式庫開發

由於 CGAL 之基本功能僅提供非封閉空間之切割，對於外圍控 制區塊需要另行處理，因此本研究透過物件封裝技巧，將其提供之非 封閉空間之切割結果，並與模擬區塊邊界框比較套疊，重新製作為封 閉空間之切割結果。 圖 3.1.3-1 不同邊界框形成之 Voronoi Diagram 圖 3.1.3-1 為相同之點位資料，套疊不同之模擬區塊邊界框，會 套疊出不同之控制區塊形狀。例如邊界 1 與邊界 3 對控制區塊 1 而 言，便會形成不同之形狀，對於邊界 3 而言，必須刪除凸出邊界的部 分；對於邊界 1 而言，則由於控制區塊本身即為封閉區塊，所以無須 另行處理。另外，不同邊界對於非封閉型之控制區塊而言，則需以邊 界本身與非封閉型之控制區塊，以聯集方式形成封閉之控制區塊。本 研究之處理邏輯如下所示： 1. 非封閉形式之控制區塊：如圖 3.1.3-1 之控制區塊 0、2、3、4、5、

6。此類狀況則以邊界框與該控制區塊做聯集判斷，可形成一個封 閉多邊形。 2. 封閉形式之控制區塊：如圖 3.1.3-1 之控制區塊 1、7。此類狀況可 依據邊界框是否穿越該控制區塊，再做分類。若無穿越，則該控 制區塊為正確解；若有穿越，則該控制區塊仍應與邊界框做聯集 判斷。

3.2 數值演算法

在地下水流的問題中，其控制方程式是基於水的質量守恆概念 所建立。根據 3.1 節所建立的控制區塊(Voronoi cells)可以視為一個有 限體積的控制體積，因此基於 Voronoi cell 所建立的水的質量守恆方 程式(即連續方程式)可寫為式 3.2-1a 與 3.2-1b，分別為穩態模擬與暫 態模擬，如下所示： 0 = +

∑

∈ i Num j ij q M i α & _{... (3.2-1a)} i Num j ij t i _{M q} t S i + = ∂ ∂

∑

∈α & _{... (3.2-1b)} 其中 代表節點 i 所代表之控制體積內的蓄水質量， 代表 節點 i 之相鄰節點數量， 則代表節點 i 與節點 j 之間之穿越流量， 則代表節點 i 所屬控制區塊內的源流項(或沈流項)。後續將進一步 定義穿越流量與蓄水質量之計算。 t i S Num_i ij M i q

(

)

⎪ ⎩ ⎪ ⎨ ⎧ ∂ ∂ − = = s darcy ij darcy f ij e s h K V rea A V M ˆ r r r & ρ ... (3.2-2)

(

)

P V P A P A B V f 0 2 2 1 0 − + + = ρ ... (3.2-3)

z P h= + ... (3.2-4) 式 3.2-2 用以定義式 3.2-1 之穿越流量，亦即等號右方之第一項， 穿越流量為水流密度、控制表面面積與達西流速的乘積，達西流速又 為水力梯度與水力傳導係數之乘積。另外，式 3.2-3 為水流密度之轉 換方程式，其為壓力與溫度之函數，其中 0、 、 與 V A₁ A₂ B均為溫度之 函數，為 Rana A. Fine(1973)提出之轉換公式。式 3.2-4 為壓力水頭與 總水頭之關係式，總水頭為壓力水頭與位置水頭之和。 Vol nS St =(ρ_{f d})t ... (3.2-5) )] ( 1 [ ₀ 0 / 0 P P n n n= +α − ... (3.2-6) 式 3.2-5 為描述控制體積內的蓄水質量，其為水流密度、飽和度、 孔隙率與控制體積大小之乘積。其中水流密度隨壓力之變化，可以透 過式 3.2-3 計算。在飽和度部分，若所探討之問題為未飽和層之地下 水流議題，其數值將隨壓力水頭(稱負壓或張力)變化而改變，因此需 要透過土壤之特性曲線加以定義；若所探討之問題為飽和層之地下水 流議題，由於飽和與未飽和之差異在於水流是否完全充滿整個孔隙， 因此飽和度不隨壓力而改變。由於本研究在此僅以飽和問題進行探 討，因此設定為 1。 另外根據地下水理論，拘限含水層之水量進出與壓力變化關 係，係受到水的壓縮性與土的壓縮性所造成，亦即土壤孔隙率會隨壓 力變化，式 3.2-6 為土壤孔隙率隨壓力變化方程式，其中 為未受水 壓影響下之土壤孔隙率， 為土壤之壓密係數，其隨土壤種類不同而 改變。 0 n / α

3.2.1 計算空間架構 在計算上，空間架構可分為節點、連結與控制體積三個部分， 如圖 3.2.1-1 所示，黑點即代表節點，意即數值模擬上的計算點；涵 蓋節點的多邊形稱為控制體積，意即該計算點的代表範圍；兩節點之 連線稱為連結，連結之存在代表該兩節點互為相鄰關係，每個連結必 穿越一個控制表面。 由於本計算模組是將模擬區域配置有限個數的計算節點，如同 其他數值方法一樣，例如：有限差分法將空間切割成有限個網格區 塊。因此，僅在計算節點上擁有空間之場變數，因此若需要非節點位 置之場變數數值，則需要搭配空間推估方法進一步推估。 前述已經列出所需之方程式，其中部分方程式僅探究節點本身 之變化，例如：控制體積內蓄水質量變化相關方程式，此類方程式稱 為節點類型之方程式。其他方程式則與控制表面有關，例如：達西公 式即是用以描述穿越控制表面流速之相關方程式，其亦作用於連結 上，在此稱為連結類型之方程式。 圖 3.2.1-2 為節點與連結關係示意圖，兩相鄰節點可以形成一連 結，其與控制區塊的邊界交點，由於 Voronoi Diagram 的特性，連結 必定與邊界垂直，且交點必定位於連結中點。 連結類型方程式可能需要連結中點之參數資訊，以達西公式為 例，若欲計算穿越控制表面之流速，其為連結中點位置(即控制表面 上)之水力梯度與水力傳導係數的乘積，水力梯度透過中央差分方式 可以藉由兩端節點數值來計算，而水力傳導係數則需要進一步透過空 間推估方法來估算，以水力傳導係數為例，其數值可以藉由兩節點水 力傳導係數之調和平均進行推估。

圖 3.2.1-1 空間架構說明示意圖

3.2.2 方程式離散化 圖 3.2.2-1 為網格示意圖，C、E、W、N 與 S 代表計算節點，因 此符號 CE 代表 C 與 E 兩點之連結中點，以下將以此網格示意圖來進 行描述。在 3.2 節中提及之方程式是以連續方程式為起始(式 3.2-1a 與 3.2-1b)。在暫態模擬上，對該式進行時間上之積分，使得等號左 式代表在某段時間之蓄水改變量，等號左式代表某段時間之穿越總量 (如式 3.2.2-1 所示)。若進一步將積分方程式進行離散化，可將式 3.2.2-1 改寫為式 3.2.2-2，等號右方可以拆成時刻τ 與時刻τ +∆τ兩項， 其中透過係數α ，來切換兩者之比重。若α等於 1，則代表完全受到 時刻τ 之控制，亦為傳統之顯示法(explicit method)；若α等於 0，則 代表完全受到時刻τ +∆τ之控制，亦為傳統之完全隱示法(fully implicit method)。 dt q Q dt t S C Num j C Cj t C i

∫

∑

∫

+∆ +∆ _∈ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + = ∂ ∂ τ τ τ _α τ τ τ / ... (3.2.2-1)

(

α

)

τ α τ τ α τ α τ τ τ ∆ ⎪⎭ ⎪ ⎬ ⎫ ⎪⎩ ⎪ ⎨ ⎧ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + − + ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + = − ∆ + ∈ ∈ ∆ +

∑

∑

C Num j C Cj C Num j C Cj C C q Q q Q S S i i / / 1 ... (3.2.2-2) 式 3.2.2-3 為達西公式(式 3.2-2)之離散型式，水力梯度的微分項 部分可藉由中央差分進行離散，其數值即以節點 C 與節點 E 之總水 頭差值除上兩點距離，其含意代表連結中點。水力傳導係數部分則必 須透過空間推估方法進一步推估，最簡單的推估方式包含算數平均 (式 3.2.2-4a)、幾何平均(式 3.2.2-4b)與調和平均(3.2.2-4c)等，其中變 數x可代表任意待處理之變數，以達西公式中的水力傳導係數(K)而 言，變數x即為水力傳導係數(K)，傳統上水力傳導係數多建議以調 和平均進行計算。另外，上風法亦可作為連結中點數值推估之基礎，

意即透過水流方向之判斷，取其上風方位數值作為連結中點數值。 ⎪⎩ ⎪ ⎨ ⎧ − − = = CE E C CE CE darcy C CE CE CE darcy CE f C CE L h h K V N Area V Q , / , , / ρ ... (3.2.2-3) 2 E C CE x x x = + ... (3.2.2-4a) E C CE x x x = ... (3.2.2-4b) E C CE x x x 1 1 2 + = ... (3.2.2-4c) 控制表面 CE 之穿越流量則為水流密度(ρ_{f ,CE})、達西流速 ( )、穿越面積( )與方向係數( )的乘積。上述達西流速 的定義，其流向是以節點 E 往節點 C 時為正，其中 代表方向係 數的乘積，其用意用以將穿越量數值改為流入控制體積 C 為正，其 數值為 1 或-1，在此方向係數為-1。反之，若欲計算控制體積 E 之穿 越量，其方向係數( )，其數值為 1，其穿越量恰與控制體積 C 之 穿越量數值正負相反。 CE darcy V _, Area_CE N_CE/C C CE N _/ E CE N _/ 圖 3.2.2-1 網格示意圖 3.2.3 內迭代處理方法

本研究在此採用節點運算概念，亦即各運算節點均可獨立運 算，運算時則假設周遭運算節點代表之壓力水頭值為正確值，而本身 之壓力水頭則透過最陡坡降法等微分形式之最佳化方法進行求解。 圖 3.2.3-1 為內迭代執行流程圖，在地下水流議題中，壓力水頭 為待解變數，由於微分形式最佳化方法需要給予初始解，並以初始解 作為初始搜尋點，在此必須給予初始壓力水頭作為初始猜值，代入方 程式集合中，節點之壓力水頭與相鄰節點之壓力水頭可以進一步計算 出穿越流量，代入連續方程式中，令控制表面的總穿越量與控制體積 內的變化量之和為守恆誤差，透過差分方式可以求得壓力水頭值對守 恆誤差之微分值，透過微分資訊之指引，可以逐步往最佳解方向收 斂。當守恆誤差趨近於零時，代表控制表面的總穿越量與控制體積內 的變化量相等，意即節點本身之壓力水頭求得。 在收斂條件方面，其收斂條件共有三種，滿足其一即停止內迭 代演算。首先是迭代次數之最大限制量，其次是守恆誤差ε之一階微 分值趨近於零，最後則是守恆誤差ε本身小於設定之內迭代收斂標 準。

圖 3.2.3-1 內迭代流程圖

3.2.4 外迭代處理方法

內迭代流程是透過最佳化方法求得各節點數值，其概念是假設 周遭節點之數值為正確，然實際僅有透過邊界條件之設定（Dirichlet

B.C.），方可確定其數值為正確。因此，外迭代流程如圖 3.2.4-1 所 示，反覆執行各節點之內迭代流程，當每次外迭代流程可將整體數值 往結果逼近，當節點數值每次之迭代改變量漸漸趨緩，若其改變量小 於外迭代收斂標準，則可認定為外迭代收斂。 圖 3.2.4-1 外迭代流程圖 3.2.5 整體數值模擬流程 整體數值模擬流程如下圖 3.2.5-1，首先依據模式設定檔，讀入 空間切割相關資訊、水文地質參數、邊界條件與方程式集合等資訊， 接著則依據模擬之模擬型態、起始時刻、結束時刻與模擬間距，開始 進行模擬，流程中則當外迭代流程收斂後，則進行是否進行下一時刻 之模擬，若否則完成模擬。若穩態形式之模擬，則僅執行外迭代流程 一次；若非穩態形式之模擬，則依據起始時刻、結束時刻與模擬間距 等資訊進行判斷。

4

四、垂向二維地下水流數值模式開發

本章節將利用前述之說明建立垂向二維地下水數值模式，前述 章節已經建立控制方程式對應離散方程式集合，本章節將進一步說明 依據前述離散方程式集合，實做前述的方程式集合。最後再說明本模 式呼叫程式函式之順序。前述章節曾提及本空間數值架構分成節點、 連結與控制體積三個部分，在方程式中亦可分類為節點類型之方程式 與連結類型之方程式兩種。本節將詳細描述因應需要所實做之各類方 程式，並說明各方程式所屬之類別與特性： 1. Func_Character_PS (斜體字代表方程式實做函式名稱) 式 4-1 是以 van Genuchten(1980)所提出之壓力水頭與飽和度的特 性曲線關係式，修改為可計算橫跨飽和與未飽和之地下水流問題。未 飽和問題中，可透過壓力水頭與土壤參數計算對應之土壤飽和度；飽 和問題中，其數值等於土體孔隙率。飽和問題中該數值則等於土體孔 隙率。本函式為節點類型之方程式。 ⎩ ⎨ ⎧ ≥ = < + = − + = − 0 1 0 ] ) ( 1 [ ) )( ( h for h for h n n p S e e r e r d θ α θ θ θ θ γ β ... (4-1) 其中S_d為飽和度、θ_r為殘餘含水量、θ_e為有效含水量、 為孔隙 率、 為壓力水頭、 n p α 、β、γ 則為土壤係數。 2. Func_Calc_Water_Density 式 4-2 可依據水體所受壓力水頭與當時溫度，計算因應之水密 度。本函式為節點類型之方程式。

(

)

p V p A p A B V f 0 2 2 1 0 ₋ + + = ρ ... (4-2)

其中

(

)

(

)

9 6 15 5 12 4 9 3 6 2 3 1 0 5 5 4 4 3 3 2 2 1 0 6 0 10 15 . 18 C , 10 2952 . 393 C , 10 7562 . 149 C , 10 44846 . 55 C , 10 92221 . 7 C , 10 224944 . 18 C , 99998396 . 0 C 1 − − − − − − × = × − = × = × − = × − = × = = + + + + + + = T C T C T C T C T C C T C V 4 9 3 6 2 4 3 1 3.2891 2.391 10 T 2.8446 10 T 2.82 10 T 8.477 10 T A = − × − + × − − × − + × − 4 12 3 10 2 8 6 5 2 10 299 . 3 10 942 . 7 10 449 . 3 10 913 . 3 10 24 . 6 T T T T A − − − − − × − × + × − × − × = 4 5 3 2 2 10 2789 . 2 10 0478 . 1 21554 . 2 037 . 147 32 . 19654 T T T T B = + − + × − − × − 3. Func_Set_Temperature_Const (

節點方程式

) 本函式用以設定地下水流溫度，由於本問題僅探討地下水流問 題，在此可將溫度設定為 。未來若搭配熱流模式，溫度之設定則 可隨熱流模式結果而變化。本函式為節點類型之方程式。 C ° 20 4. Func_WaterCapacity 式 4-3 則用以將土壤飽和度轉換為控制體積內的蓄水質量。本函 式為節點類型之方程式。本函式為節點類型之方程式。 Vol nS S f d t t ) (ρ = ... (4-3) 5. Func_DarcyLaw_Flux 式 4-4 首先透過達西公式，以水力梯度與實際水力傳導係數乘積 計算斷面穿越流速，其次再乘上控制表面面積與水流密度，轉換為斷 面之穿越水質量。本函式為連結類型之方程式。

(

)

(

)

⎪ ⎩ ⎪ ⎨ ⎧ ≠ ∆ + − + − = = j i where x z p z p K V rea A V M j j i i darcy darcy r r r & ρ ... (4-4)

6. Func_Calc_Porosity_T20 式 4-5 用以定義壓力水頭與實際孔隙率關係，本函式用以呈現土 體之壓縮性。未飽和問題中，其數值固定不變；飽和問題中，其數值 隨壓力增加而線性增加。本函式為節點類型之方程式。 0 0 )] ( 1 [ 0 0 0 ' 0 ≤ = > − + = P for n n P for P P n n n α ... (4-5) 7. Func_Interpolate_Density 本函式透過連結之兩端點水流密度值，以算數平均推估連結中 央之水流密度，其數值可代表控制表面上的水流密度。本函式為連結 類型之方程式。 8. Func_InterporlateDensity_UpWind 本函式透過連結之兩端點水流密度值，以上風法之概念推估連 結中央之水流密度，其數值可代表控制表面上的水流密度。本函式為 連結類型之方程式。 9. Func_InterporlateK 本函式透過連結之兩端點實際水力傳導係數，以調和平均推估連結中 央之實際水力傳導係數，其數值可代表控制表面上的水力傳導係數。 本函式為連結類型之方程式。 10. Func_InterporlateMoisture 本函式透過連結之兩端點含水量，以算數平均推估連結中央之 含水量，其數值可代表控制表面上的含水量。本函式為連結類型之方

程式。 11. Func_InterporlatePorosity 本函式透過連結之兩端點實際孔隙率，以算數平均推估連結中 央之實際孔隙率，其數值可代表控制表面上的實際孔隙率。本函式為 連結類型之方程式。 12. Func_Continuity_Universal 本函式為地下水模式中之守恆方程式，依據模擬問題可分為穩 態與暫態兩種模擬方式，其中暫態模擬又可依據顯示法(explicit method)、完全隱示法(fully implicit method)與混合隱示法

(Crank-Nicholson method)，建立不同時間項之處理方式。 在穩態模擬上，連續方程式如式 4-6 所示，總穿越量與源項之 總和應為零，意即進出平衡。 0 / + =

∑

∈ C Num j C Cj q M i α ... (4-6) 在暫態模擬上，連續方程式如式 4-7 所示，前述顯示法、完全 隱示法與混合隱示法之差異僅在係數α 的選定上，若係數α設定為 1，表示所使用之數值處理方法為顯示法；而係數α設定為 0，則表示 所使用之數值處理方法為完全隱示法；而係數α 設定為 0.5，則表示 所使用之數值處理方法為混合隱示法。

(

α

)

τ α τ τ α τ α τ τ τ ∆ ⎪⎭ ⎪ ⎬ ⎫ ⎪⎩ ⎪ ⎨ ⎧ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + − + ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + = − ∆ + ∈ ∈ ∆ +

∑

∑

C Num j C Cj C Num j C Cj C C q M q M S S i i / / 1 ... (4-7) 上述函式為本研究為地下水流議題所實做之函式，其中方程式

Func_Interpolate_Density 與 Func_InterporlateDensity_UpWind 都屬於

水流密度之空間推估函式，為相同功能但內部處理方法不同之函式， 因此可以依據需求擇一使用。此外，方程式 Func_Character_PS 為 van Genuchten 所提出之特性曲線模型，因此未來亦可因應不同之需求， 選用其他特性曲線模型，例如 Brook 等人所提出之方法。

5

五、垂向二維地下水流數值模式模擬與驗證

前述章節已經建立完成地下水流議題之運算模組，以下將透過 十一種不同類型之地下水案例進行驗證，其中將針對均質與非均質、 穩態與非穩態，以及與規則網格與不規則網格之差異，建立不同的模 擬案例。並針對各個案例的模擬結果進行誤差分析與圖形化輸出，以 探討模式的正確性。以下介紹各個案例設定、模擬結果與相關誤差分 析。

5.1 穩態規則網格案例驗證與模擬

本節為穩態案例之展示，將分為四個不同的案例進行數值驗證 與討論。案例 5.1-1 為無源流項且邊界設定對稱之均質案例；案例 5.1-2 則在邊界條件形狀上改為不對稱；案例 5.1-3 則另外增加一個抽水 井，用以驗證源流項(source)或沈流項(sink)；案例 5.1-4 除了增加抽 水井外，另外在含水層材質方面，則改變為非均質案例。 1. 案例 5.1-1 案例說明： 整體模擬區域如圖 5.1-1 所示，為 11 公尺乘上 13 公尺之垂向二 維方形薄板。網格切割上，由原點規則地每 1 公尺便配置一個計算節 點，透過 Voronoi Diagram 可以產生規則形狀之 11×13 矩形網格，因 此每格長與寬均為 1 公尺。在材質方面，整體區塊設定為相同材質， 其中土壤孔隙率設定為 0.38，水力傳導係數設定為 0.01( )，在 土壤壓密係數( )方面，本案例與後續案例均設定為 0。 day m / ' α 在邊界條件的設定上，左方邊界之總水頭設定為 120(m)，右方 邊界之總水頭設定為 80(m)，上下方邊界則設定為無流量邊界(No

Flow Boundary)，因此本案例水位將由左往右逐步遞減，亦即水流方 向應為由左往右流動。本案例為穩態模擬，全體網格之初始猜值雖可 任意設定，但過於懸殊之設定值易影響模擬收斂性與效率，因此設定 為 100(m)，較為鄰近模擬結果。透過模式求解，各網格數值必逐步 依據邊界條件之設定值，而收斂於前述之水流型態。 本案例將探討不同之外迭代收斂標準對於模式模擬精度之影 響，以下將分別針對當外迭代收斂標準為 、 、 、 與 時進行探討。 -4 10 1× 1×₁₀-5 1×₁₀-6 1×₁₀-7 -8 10 1× 在內迭代收斂條件方面，最大內迭代次數設定為 30。此外，由 於內迭代收斂標準極為重要，其數值同時會影響精度與計算效率，因 此本研究並不直接設定一數值。第一次外迭代時，所有節點之內迭代 均以極小值作為內迭代收斂標準，因此其停止條件是以最大迭代次 數，其後的外迭代則採首次外迭代所有節點之平均守恆誤差，作為後 續內迭代收斂標準。 數值結果： 圖 5.1-2 為總水頭之模擬等值圖，從圖面看來總水頭由左向右逐 步遞減，且其等值線分佈均勻，顯示遞減趨勢趨於直線。圖 5.1-3 為 壓力水頭之模擬等值圖，由於壓力水頭為總水頭扣除位置水頭，因此 其數值由上而下遞增，因此等值線略呈傾斜。由於案例邊界條件設計 極為簡單，因此套用達西公式即可簡單計算左右兩方邊界之穿越流 量，其為質量流率，因此達西流速應再乘上穿越面積與水流密度。 表 5.1-1 為解析解與數值解於左右兩側邊界之穿越流量比較表， 數值解方面是將左右兩側網格所計算之穿越流量累加所得。由表上所 示，隨著外迭代收斂標準逐步降低，數值解與解析解之差異亦逐步降

低，當外迭代收斂標準為 時，其數值誤差已可滿足 0.01%之門 檻。 -7 10 1× 另外，表 5.1-2 用以展示不同外迭代收斂標準下，其計算效率與 守恆誤差結果，此表的計算時間是以中央處理器(cpu)為 AMD

Opteron™ Processor 246，記憶體(ram)為 2GB 的個人電腦執行模式所 得出的時間。表中顯示隨著收斂標準之降低，數值模式需要以更多的 計算時間，換得更加的計算精度。在此均勻網格中，每個網格均與另 外四個網格相鄰，而兩相鄰網格會因為水位差而形成穿越流量，在此 穩態模擬中，四面之穿越流量總和理論上應該要為零，然在數值上則 無法完全為零，其數值稱為該網格之守恆誤差，因此平均守恆誤差是 將所有網格守恆誤差取絕對值再加總平均，而對數守恆誤差則是將前 述之平均守恆誤差取對數值。 圖 5.1-4 是以計算時間與對數守恆誤差所繪出之圖形，由圖面顯 示隨著計算資源之增加，其計算精度的確可逐步降低，然在圖面上顯 示，當對數守恆誤差約為-4 時，其斜率變為和緩。若以工程經濟學觀 點來看，成本增加量(cost)與利益增加量(benefit)可以用以決定最佳之 投資規模，因此亦可決定該點為最佳之外迭代收斂標準，其對應之外 迭代收斂標準為 -7，因此後續案例將直接以此為收斂標準。 10 1× 表 5.1-1 解析解與數值解於左右兩側邊界穿越流量比較表 外迭代收斂標準 左側邊界流量(kg/day) 右側邊界流量(kg/day) -4 10 1× 452.4700(95.847%) -491.3430(104.081%) -5 10 1× 470.1291(99.587%) -474.0252(100.413%) -6 10 1× 471.8834(99.959%) -472.2711(100.041%) -7 10 1× 472.0580(99.996%) -472.0970(100.004%) 數值解 -8 10 1× 472.0753(99.999%) -472.0790(100.000%) 解析解 無 472.077(100%) -472.077(100%)

表 5.1-2 外迭代收斂標準、計算效率與守恆誤差結果表 外迭代收斂標準 計算時間(s) 外迭代次數 平均守恆誤差 對數守恆誤差 -4 10 1× 12 79 2.71841×10-1 -0.57 -5 10 1× 16 134 2.725×10-2 -1.56 -6 10 1× 19 188 2.711×10-3 -2.57 -7 10 1× 23 243 2.716×10-4 -3.57 -8 10 1× 47 297 2.9×10-5 -4.54 圖 5.1-1 案例 5.1-1 格網與模式配置圖

0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10 11 12 圖 5.1-2 案例 5.1-1 總水頭等值圖

0 2 4 6 8 10 0 2 4 6 8 10 12 圖 5.1-3 案例 5.1-1 壓力水頭等值圖

-5 -4 -3 -2 -1 0 0 10 20 30 40 50 迭代時間(秒) 對 數 守恆 誤差 對數守恆誤差與迭代時間比較圖 圖 5.1-4 外迭代收斂標準與對數守恆誤差趨勢圖 2. 案例 5.1-2 案例說明： 案例 5.1-1 為上下對稱的案例，水流流動極為簡單，案例 5.1-2 則以非對稱案例進行探討。整體模擬區域如圖 5.1-5 所示，模擬範圍、 模擬網格、材質設定等與案例 5.1-1 一樣。在邊界條件的設定上，左 下方邊界長度為 3 公尺的區塊，其總水頭設定為 10(m)，右上方邊界 長度為 3 公尺的區塊，其總水頭設定為 30(m)，其餘邊界則設定為無 流量邊界(No Flow Boundary)，因此本案例水流將由右上方往左下方 流動。本案例為穩態模擬，全體網格之初始猜值雖可任意設定，但過 於懸殊之設定值易影響模擬收斂性與效率，因此設定為 20(m)，較為 鄰近模擬結果。透過模式求解，各網格數值必逐步依據邊界條件之設 定值，而收斂於前述之水流型態。 由於案例 5.1-1 已經檢驗不同之外迭代收斂標準對於模式模擬精 度之影響，以下案例之外迭代收斂標準設定為 -7。 10 1×

數值結果： 圖 5.1-6 為總水頭之模擬等值圖，從圖面看來總水頭分佈由右上 向左下逐步遞減，且以對角線略顯對稱。若觀察等值線鄰近邊界部 分，於無流量邊界條件區域之等值線，則往往顯得與邊界近似垂直； 而在定水頭邊界條件區域之等值線，則鄰近設定之水頭數值，顯示模 擬結果的確依據邊界條件之設定，而有正確的模擬數值。 本案例因邊界條件之設定較為複雜，解析解較難求得，因此本 案例並無解析解配合驗證。本案例可以應用穩態模擬之進出水量需達 平衡之特點進一步檢驗，本案例右上方邊界之流入總量為 80.7653 (kg day)，左下方邊界之流出總量為 80.6871(kg day)，兩者之相對誤 差為 0.096871%，由此可知本案例符合進出平衡之條件，符合水質量 之守恆特性。

0 2 4 6 8 10 0 2 4 6 8 10 12 圖 5.1-6 案例 5.1-2 總水頭等值圖 3. 案例 5.1-3 案例說明： 案例 5.1-1 為無源流項案例，亦即無內部抽水的案例，案例 5.1-3 則進行源流項之探討。整體模擬區域如圖 5.1-7 所示，模擬範圍、模 擬網格、材質設定等與案例 5.1-1 一樣。在邊界條件的設定上，左右 邊界之總水頭均設定為 80(m)，上下邊界則設定為無流量邊界(No Flow Boundary)。於點位(5, 5)處配置抽水井，並以 500(kg/day)之抽水

量進行抽水。因此本案例本身並無背景水流流動，受到抽水行為，因 此總水頭將在以抽水井為中心，形成一個洩降錐。本案例為穩態模 擬，全體網格之初始猜值雖可任意設定，但過於懸殊之設定值易影響 模擬收斂性與效率，因此設定為 80(m)，較為鄰近模擬結果。透過模 式求解，各網格數值必逐步依據邊界條件之設定值，而收斂於前述之 水流型態。 由於案例 5.1-1 已經檢驗不同之外迭代收斂標準對於模式模擬精 度之影響，以下案例之外迭代收斂標準設定為 -7。 10 1× 數值結果： 圖 5.1-8 為總水頭之模擬等值圖，從圖面看來的確如預期一般， 形成一個橢圓形的洩降錐，其洩降錐頂的洩降量為 30(m)。此外，本 研究亦進一步作水平衡分析，由於系統中有個抽水井，且上下邊界設 定為無流量邊界，因此在穩態系統中，必須透過左右邊界之定水頭邊 界條件來補充水量，且其補充量必須等於抽水量。左右邊界之穿越流 量為 499.5716(kg/day)，與內部抽水量相比，其計算誤差約為 0.086%， 因此可以判定通過水平衡分析。

0 2 4 6 8 10 12 圖 5.1-8 案例 5.1-3 總水頭等值圖 4. 案例 5.1-4 案例說明： 相較於案例 5.1-3，本案例將均質含水層改為非均質含水層，兩 者之網格切割、邊界條件與抽水量等均一致。圖 5.1-9 為本案例之材 質設定圖，其中材質 1 與案例 5.1-3 之材質相同，其孔隙率為 0.38、 水力傳導係數為 0.01(m/day)；材質 2 之孔隙率維持為 0.38，其水力 傳導係數則更改為 0.001(m/day)，可以視為一層垂直的阻水層。其位

置位於抽水井之左方，厚度約一公尺。 由於案例 5.1-1 已經檢驗不同之外迭代收斂標準對於模式模擬精 度之影響，以下案例之外迭代收斂標準設定為 -7。 10 1× 數值結果： 圖 5.1-10 為案例 5.1-3 與案例 5.1-4 之總水頭之模擬等值圖，圖 面上實線為案例 5.1-3 之數據，而虛線是案例 5.1-4 之模擬結果。兩 者在抽水井位置均產生洩降錐，然案例 5.1-4 因受到抽水井左方阻水 層之阻礙，在抽水井左方之水力坡降較案例 5.1-3 為陡，意即阻水層 左方區域受到抽水井之影響極小，直至阻水層右方方迅速下降。此 外，有阻水層案例之最大洩降量約為 40(m)，而無阻水層案例之最 大洩降量約為 30(m)，其原因應為受到左方阻水層之影響，水源供 應來源僅剩右方，透過較深的洩降錐與較大的水力坡降，方可供應 相同之水量。 本研究亦進一步針對案例 5.1-4 作水平衡分析，由於系統中有個 抽水井，且上下邊界設定為無流量邊界，因此在穩態系統中，必須透 過左右邊界之定水頭邊界條件來補充水量，且其補充量必須等於抽水 量。左右邊界合計之穿越流量為 499.9744(kg/day)，與內部抽水量相 比，其計算誤差約為 0.00512%，因此可以判定通過水平衡分析，其 中左側(含阻水層側)流量為 137.4873(kg/day)，右側(無阻水層側)流量 為 362.4871(kg/day)，兩側流量因阻水層而相差了 225(kg/day)。

0 2 4 6 8 10 0 2 4 6 8 10 12 圖 5.1-10 案例 5.1-3 與案例 5.1-4 之總水頭等值圖

5.2 暫態規則網格案例驗證與模擬

本節則展示暫態(transient)案例，以下將以四個不同的案例進行 驗證與討論。案例 5.2-1 具有一個源流項，左右邊界設定對稱，且其 整體材質為相同材質，與前述案例 5.1-3 配置相同，本案例除了驗證 模式的基本正確性外，亦針對完全隱示法(fully implicit method)與混 合隱示法(Crank Nicholson method)兩種對時間處理上的不同數值技 巧(numerical scheme)進行探討與分析。案例 5.2-2 與案例 5.2-1 大致均 相同，不同之處則是採取時變性的源流項。案例 5.2-3 則加入一個垂 直分佈的阻水層，藉此探討材質非均值之影響。案例 5.2-4 與案例 5.2-1 大致相同，但採取更細之規則網格進行模擬，藉此比較不同網格尺寸 對於模擬結果之影響。 5. 案例 5.2-1 案例說明： 本案例之網格配置與邊界條件均與案例 5.1-3 相同，為案例 5.1-3 之暫態案例，其左右邊界設定為。在邊界條件的設定上，左右邊界之 總水頭均設定為 80(m)，上下邊界則設定為無流量邊界(No Flow Boundary)。在初始條件上，所有位置之初始總水頭為 80(m)，意即初 始水位代表未經抽水時狀態，處於靜水壓分佈。於點位(5, 5)處配置 抽水井，並以 500(kg/day)之抽水量進行抽水。另外，模擬間距為 0.01 天(約 14.4 分鐘)，總模擬時刻數為 6 個時刻。因此單一時刻之抽水量 為 5(kg)。在時間離散上的處理，本案例將分別針對完全隱示法與混 合隱示法進行探討。 數值結果：

圖 5.2-1 為抽水井位置總水頭之時間變化圖，圖中同時顯示完全 隱示法與混合隱示法兩者。從圖面看來完全隱示法之數值結果是由初 始水位逐漸往下遞減，至總水頭 50(m)時趨於穩定，洩降錐約為 30(m) 深，相較於穩態案例 5.1-3 而言，穩定時之水位與穩態水位相近。相 較之下，混合隱示法則有震盪之問題，時刻 1 洩降錐頂水位約為 25(m)，時刻 2 洩降錐頂水位則躍升為 70(m)，極為不合理。 此外，圖 5.2-2 為透過達西公式計算抽水井所在網格與周遭網格 之瞬間穿越流率，顯現由混合隱示法所得之瞬間穿越流率亦有震盪情 形。圖 5.2-3 則是各時刻區間之穿越流量變化圖，然而從圖面顯示兩 數值技巧所得之穿越流量卻幾近一致，兩者一致的原因是因為控制方 程式必須符合水量守恆，因此穿越流量必須等於抽水量。然而在完全 隱示法中，穿越流量之計算是完全以後一時刻之瞬間流率乘上時間間 距長度；而在混合隱示法中，則是綜合前一時刻與後一時刻之瞬間流 率，以兩者之平均乘上時間間距長度。因此在混合隱示法中，為了使 時刻 1 的穿越流量等於該時刻之抽水量，必須透過時刻 0 與時刻 1 之 平均瞬間穿越流率，然本案例之初始條件為完全無流動之水流條件， 因此時刻 1 必須以更大之瞬間穿越流率，方可符合質量守恆之條件， 此為數值震盪之由來。 圖 5.2-4 為抽水井附近之水位剖面上視圖，圖 5.2-5 是線段CC'之 水位剖面圖，完全隱示法之水位剖面則逐步往抽水井位置下降，且時 刻 2 之水位低於時刻 1 之水位，其數值結果符合預期。然混合隱示法 之模擬結果則與震盪之情形，且時刻 2 之水位高於時刻 1 之水位。由 本案例結果顯示，在本問題中，完全隱示法是較適合的數值技巧。 圖 5.2-6 為完全隱示法之不同時刻邊界流量變化圖，由圖面顯示 邊界流量由低逐步增大，到時刻 3 之邊界流量趨近於抽水井之 5(kg)

之抽水量，其變化顯示抽水之影響範圍逐步往外延伸，直到時刻 3 方 趨於穩態。 式 5.2-1 為相對系統守恆誤差計算式，BCF 代表邊界流量，q 代 表抽水量，因此系統內部變化量、邊界流量與抽水量應維持守恆，若 不為零則代表系統之絕對守恆誤差。表 5.2-1 為應用完全隱示法的相 對系統守恆誤差表，各時刻之相對系統守恆誤差均極小，因此從水量 平衡的觀點顯示本模式在暫態模擬上亦有不錯的表現。 q q BCF S S_t − _t − − = +1 ε ... (5.2-1) 表 5.2-1 完全隱示法相對系統守恆誤差表 時刻編號 _{1 2 3 4 5 6} 守恆誤差(%) -0.023 0.030 0.022 0.038 0.014 -0.009 20 25 30 35 40 45 50 55 60 65 70 75 80 0 20 40 60 80 100 時間(分) 總水頭 (m) 完全隱示法 混合隱示法 圖 5.2-1 抽水點各時刻總水頭變化圖

0 100 200 300 400 500 600 700 800 900 1000 0 20 40 60 80 100 時間(分) 瞬時 流率 (k g/ da y) 完全隱示法 混合隱示法 圖 5.2-2 抽水井所在網格之瞬時流率變化圖 0 1 2 3 4 5 14.4 28.8 43.2 57.6 72 86.4 時間(分) 流 入 流量 (kg) 完全隱示法 混合隱示法 圖 5.2-3 抽水井所在網格之流量變化圖

20 25 30 35 40 45 50 55 60 65 70 75 80 0 1 2 3 4 5 6 X座標 總水 頭(m ) 完全隱示法_時刻1 完全隱示法_時刻2 混合隱示法_時刻1 混合隱示法_時刻2 圖 5.2-5 不同數值技巧之水位剖面比較圖 3 3.5 4 4.5 5 14.4 28.8 43.2 57.6 72 86.4 時間(分) 邊界 流量 (k g ) 邊界流量 圖 5.2-6 完全隱示法之邊界流量變化圖 6. 案例 5.2-2 案例說明：

本案例與案例 5.2-1 之差異在於抽水量之設定，本案例抽水量為 由無抽水逐步遞增，抽水量設定方式如圖 5.2-7 所示，抽水量隨時間 逐步增加。 數值結果： 圖 5.2-8 為本案例抽水井位置之各時刻水位變化圖，兩數值技巧 所得之結果均逐步下降，兩者水位相近，其下降趨勢亦相近，僅混合 隱示法之下降趨勢較不圓滑。相較於前述案例 5.2-1 混合隱示法之震 盪結果，在此有不同之結果，其原因在於案例 5.2-1 之初始水位為完 全沒有抽水之水流狀態，意即在時刻 0 之前的抽水量為 0，開始模擬 後則瞬間躍升為設定之抽水量，其抽水型態為脈衝函數(如圖 5.2-9 所 示)，因此對於混合隱示法而言形成震盪。在案例 5.2-2 而言，其抽水 量由小逐步遞增，減弱起始抽水所造成的衝擊，因此使用混合隱示法 不至於造成數值震盪。 圖 5.2-10 為本案例不同數值技巧在時刻 6 之水位剖面，完全隱 示法與混合隱示法之水位剖面均極為類似，亦可作為混合隱示法可應 用於本類型問題之依據。

0 0.5 1 1.5 2 2.5 3 3.5 14.4 28.8 43.2 57.6 72 86.4 時間(分) 抽水 量 (kg) 抽水量 圖 5.2-7 抽水井抽水量變化圖 20 25 30 35 40 45 50 55 60 65 70 75 80 0 20 40 60 80 100 時間(分) 總水頭 (m) 完全隱示法 混合隱示法 圖 5.2-8 抽水點各時刻總水頭變化圖

0 0.5 1 1.5 0 1 2 3 4 5 6 7

x

y

圖 5.2-9 脈衝函數示意圖 50 55 60 65 70 75 80 0 1 2 3 4 5 6 X座標 總水頭 (m ) 完全隱示法_時刻6 混合隱示法_時刻6 圖 5.2-10 不同數值技巧之水位剖面比較圖 7. 案例 5.2-3 案例說明： 本案例為穩態案例 5.1-4 之暫態版本，因此與案例 5.2-1 之差異，

在於附加的垂向阻水層。圖 5.1-9 為本案例之材質設定圖，其中材質 1 與案例 5.1-3 之材質相同，其孔隙率為 0.38、水力傳導係數為 0.01(m/day)；材質 2 之孔隙率維持為 0.38，其水力傳導係數則更改為 0.001(m/day)，因此可以視為一層垂直的阻水層。 其餘初始條件、邊界條件與抽水量則與暫態案例 5.2-1 一致。初 始條件設定為所有位置之初始總水頭為 80(m)。抽水井抽水量設定為 500(kg/day)。另外，前述已經比較完全隱示法與混合隱示法之優劣， 顯現完全隱示法較適合本案例型態之模擬，後續案例均以完全隱示法 進行模擬。 數值結果： 由於本案例因為材質分佈並非左右對稱，因此在水位剖面之擷 取上則分別取出抽水井左右兩方，圖 5.2-11 為本案例展示之水位剖 面。 圖 5.2-12 為本案例與案例 5.2-1 之水位剖面比較圖，由於系統於 時刻 3 之後漸趨穩態，因此本圖亦僅展現前三個時刻的水位剖面，由 圖面顯示案例 5.2-1 之水位剖面均左右對稱，本案例數據則顯現左右 不對稱，左方兩個節點之水力梯度極陡，從本案例之材質與網格配置 圖(圖 5.1-9)看來，前述的左方兩個節點均位於阻水層中，因此為了滿 足抽水之消耗，左方節點必須透過更大的水力梯度方能維持一定的穿 越流量。 圖 5.2-13 為左右邊界之穿越流量變化圖，從圖面顯示左方流量 約為右方流量之半。綜合前述之探討，本案例可以顯現本模式的確可 以應用於暫態非均質問題之模擬上。

20 25 30 35 40 45 50 55 60 65 70 75 80 0 2 4 6 8 10 12 X座標 總水頭 (m ) 案例5.2-3_時刻1 案例5.2-3_時刻2 案例5.2-3_時刻3 案例5.2-1_時刻1 案例5.2-1_時刻2 案例5.2-1_時刻3 圖 5.2-12 案例 5.2-1 與 5.2-3 之水位剖面比較圖 0 5 14.4 28.8 43.2 57.6 72 86.4 時間(分) 總水頭 (m ) 總邊界流量 左側邊界流量 右側邊界流量 圖 5.2-13 案例 5.2-3 邊界流量變化圖 8. 案例 5.2-4 案例說明： 本案例所模擬的問題與暫態案例 5.2-1 相同，僅在網格尺寸上改

變，其網格為 0.33(m)見方的方形網格，其計算節點數量將增加 9 倍， 圖 5.2-14 為本案例的網格配置圖。本案例將與案例 5.2-1 進行比較探 討，進一步瞭解網格尺寸對於模擬結果之影響。 數值結果： 圖 5.2-15 為本案例與案例 5.2-1 之水位剖面，兩案例之差異在 於網格尺寸大小，案例 5.2-1 為粗網格案例，本案例為細網格案例， 其水位剖面上視圖如圖 5.2-4 所示，由圖面所示大部分區域之水位均 吻合，然在抽水井鄰近區域，細網格案例之水位低於粗網格案例，其 兩者最大洩降量約差異 10(m)。 此外，圖 5.2-16 為兩案例之邊界流量變化圖，不同網格尺寸均 可得到相同的邊界流量。表 5.2-2 為兩案例之相對系統守恆誤差表， 兩案例之相對系統守恆誤差均極小，表示各時刻之邊界流量、抽水量 與系統蓄水變化量均符合質量守恆定律。 這表示在不同網格尺寸案例中，模擬結果雖均可符合質量守恆 定律，但在水位的呈現上可能會有不同的結果。由於抽水井附近水位 變化較為劇烈，傳統模擬技巧多會建議配置較細之網格方可精確掌握 水位變化，因此可認定本案例之模擬結果較接近真值，而案例 5.2-1 之模擬結果則在洩降錐上有約 10(m)的誤差。 然在計算效率上，本案例的計算節點為案例 5.2-1 之九倍，因此 可以預期其需要耗費更多的計算時間，方可達到收斂。案例 5.2-1(粗 網格案例)需要耗費 82 秒，而本案例(細網格案例)則需要耗費 3945 秒，需要增加約 50 倍之計算時間。因此對於網格尺寸應該依據所需 精度與可接受之計算時間來衡量，藉此決定模式網格尺寸。

表 5.2-2 相對系統守恆誤差表(案例 5.2-1 與 5.2-4) 時刻編號 _{1 2 3 4 5 6} 案例 5.2-1 守恆誤差(%) -0.023 0.030 0.022 0.038 0.014 -0.009 案例 5.2-4 守恆誤差(%) -0.103 -0.146 -0.148 -0.134 -0.124 -0.083 圖 5.2-14 案例 5.2-4 格網與模式配置圖

20 25 30 35 40 45 50 55 60 65 70 75 80 0 1 2 3 4 5 6 X座標 總水 頭(m) 案例5.2-4_時刻1 案例5.2-4_時刻2 案例5.2-4_時刻3 案例5.2-1_時刻1 案例5.2-1_時刻2 案例5.2-1_時刻3 圖 5.2-15 水位剖面比較圖(案例 5.2-1 與 5.2-4) 0 1 2 3 4 5 14.4 28.8 43.2 57.6 72 86.4 時間(分) 流 入 流量 (kg) 案例5.2-4 案例5.2-1 圖 5.2-16 邊界流量變化圖(案例 5.2-1 與 5.2-4)

5.3 不規則網格之暫態飽和侷限含水層模擬

前述 5.2 節之案例 5.2-1 與 5.2-4 已經用規則矩形網格來驗證並 顯示網格尺寸對於模擬結果之影響，雖然計算精度大幅提昇，但是其 計算量亦暴增至將近 10 倍上下。此外，由於抽水井之影響，鄰近抽 水井之水流流速遠高於遠離抽水井位置之水流流速，其水位變化亦 然。因此僅需在抽水井周遭配置較細網格，在外圍區域仍可以較粗網 格節省計算量，如此可在一定的精度需求下，維持低計算量。 表 5.3-1 為本節案例之差異說明表，其中分為極座標佈點型態與 矩形佈點型態兩種，極座標佈點是以抽水井為圓心，環狀往外延伸佈 點；矩形佈點則在鄰近抽水井附近區域，以更密的矩形佈點附加在原 本網格上。圖 5.3-1 至圖 5.3-3 分別為三個案例的配置網格，加密後 之計算節點數量列於表 5.3-1 上。 表 5.3-2 為各案例的模擬結果比較表，若以案例 5.2-1 與 5.2-4 而 言，兩者網格數量將近十倍，所消耗之計算時間則將近五十倍，而兩 者之洩降錐深度差異將近 8 米深，顯現較細網格對於加強計算精度有 明顯效果，但其耗費之計算量亦大幅提昇。而本節所建立之加密格網 在計算節點方面分別為 147、255 與 255 個計算節點，在計算時間方 面則分別為 90、415 與 315 秒，其洩降錐誤差分別為 6.87、0.70 與 0.21 米深。由上述數據顯示，若在關鍵區域加強配置一定比例之計算 節點，並透過 Voronoi Diagram 的空間切割方法，可以在一定的增加 成本下，有效提升計算精度，以案例 5.3-3 而言，其計算誤差僅約 0.21 米深，計算時間則不到案例 5.2-4 的十分之一。 圖 5.2-4 為水位剖面之上視圖，圖 5.3-4 則是依據圖 5.2-4 所述之 剖面，繪出各案例之水位剖面比較圖，從圖面顯示除案例 5.3-1 外， 其餘案例剖面與案例 5.2-4 剖面趨勢一致。