智慧型可適性計算平台之發展及其於地下水流、熱流與污染傳輸耦合模擬之應用

國 立 交 通 大 學

土 木 工 程 學 系

博 士 論 文

智慧型可適性計算平台之發展及其於地下

水流、熱流與污染傳輸耦合模擬之應用

Development of Intelligent and Adaptive

Computation Framework and its

Application on the Simulation of

Groundwater Flow Coupled with Heat and

Pollutant Transport

研 究 生 ： 蔡瑞彬

指導教授 ： 張良正 博士

智慧型可適性計算平台之發展及其於地下水流、熱流與污染

傳輸耦合模擬之應用

Development of Intelligent and Adaptive Computation Framework and its Application on the Simulation of Groundwater Flow Coupled with Heat

and Pollutant Transport

研 究 生：蔡瑞彬 Student：Jui-Pin Tsai 指導教授：張良正 Advisor：Liang-Cheng Chang 國 立 交 通 大 學 土 木 工 程 學 系 博 士 論 文 A Thesis

Submitted to Department of Civil Engineering College of Engineering

National Chiao Tung University in Partial Fulfillment of the Requirements

for the Degree of PhD

in

Civil Engineering October 2010

Hsinchu, Taiwan, Republic of China

I

智慧型可適性計算平台之發展及其於地下水流、熱流與污染

傳輸耦合模擬之應用

研究生：蔡瑞彬 指導教授：張良正 教授

摘要

為了解決傳統數值模式不易擴充及人工智慧方法難與數值模式 深入整合等兩大問題，本研究提出一新的計算架構，智慧型可適性計 算平台，突破傳統數值模式之開發方法，建立一系統化、自動化之數 值建模平台，應用此系統所開發之模式，可隨著需求逐步增加模式模 擬能力，並整合人工智慧方法於模擬模式核心，而使模式具有類「人 工智慧」能力，如可以因應環境變化調整自身模擬範疇等。模式開發 上分為兩階段進行，第一階段為開發一具可擴充性之建模系統，稱為 可適性計算平台(ACF)，以 ACF 所開發之模式可容易的新增或修改模 式所計算之方程式，第二階段再將 ACF 與人工智慧進行整合，而成 為智慧型可適性計算平台(IACF)，其中在人工智慧的部分本研究整合 了專家系統(expert system)及類神經網路(ANN)。 本研究將發展完成之計算平台應用於發展地下水流、熱流及污 染物傳輸耦合模擬模式，為檢驗模式之正確性與適用性，本研究設計 數個簡例，再應用IACF 所發展模式與 Tough2 及 HST3D 兩模式進行 模擬驗證。其中 Tough2 與 HST3D 以不同之方程式描述地下水流、 熱流及污染傳輸問題之運動機制，由於IACF 平台之擴充彈性，本研 究可在短時間內完成與此兩模式之比較驗證。驗證結果顯示，IACF 所發展模式與前述既有模式之模擬結果，具有高度一致性，此亦可證 明了應用IACF 建模相較於傳統建模方式更有效率。本研究亦將專家

II 系統與 ANN 整合入以 IACF 所發展之模式中，由模擬結果顯示，考 量專家系統之判斷規則，可因應環境變化而自動調整計算模擬，有效 提高模擬效率。系統並成功的將類神經網路實作於地下水流模組中， 證明了IACF 具有與不同人工智慧方法整合之彈性。前述模擬驗證結 果顯示了以IACF 發展模式，可容易的修改模式所計算之方程式或新 增原模式未考量之方程式，因此可提升模式之開發效率，協助研究者 致力於問題本身機制之探討，減少程式開發或修改所費之時間與精 力。

III

Development of Intelligent and Adaptive Computation

Framework and its Application to the Simulation of

Groundwater Flow Coupled with Heat and Pollutant

Transport

Student：Jui-Pin Tsai Advisor：Dr. Liang-Cheng Chang Department of Civil Engineering

National Chiao Tung University

Abstract

Difficult to modify or extend the computational functions of a numerical model and incapable of integrating artificial intelligent methods into the computation kernel are two important issues in the conventional numerical model development and maintenance. Therefore, this study proposed a novel numerical model framework to improve the numerical model developing efficiency. One can easily modify or

increase the computational functions of a numerical model developed by the proposed framework. Furthermore, the model can integrate artificial intelligent methods into the computational kernel. The whole framework is denominated as intelligent adaptive computation framework (IACF) and its development was divided into two steps. The first step was to develop an adaptive computation framework (ACF). The ACF can facilitate the development of a model with the flexibility of adding new functions to simulate new modeling mechanisms easily. The second step is to integrate artificial intelligent methods, expert system and artificial neural network, into ACF and completed the IACF.

After developing the IACF, this study applied the IACF to develop a simulation model for solving groundwater flow coupled with heat and pollutant transport problem. To examine the accuracy of the simulation model, the model simulation results of several hypothesis cases were compared with those of Tough2 and HST3D. The Tough2 and HST3D

IV

have applied different equations to describe the flow, heat and pollutant transport problem. The IACF made the model comparison became possible within a limited time. The simulation results showed good consistence among the three models. The situation not only verified the model accuracy but also, on the other hand, demonstrated the capability of IACF to improve model development efficiency. A rule based expert system and ANN model were integrated into the simulation model. Simulation results show that the model can adjust the computations according to the system states and effectively improve the computation efficiency by activating the expert system. The embedding of ANN model was also successfully verified and that demonstrate the flexibility of IACF to integrate various artificial intelligent methods. The flexibility of modifying an exited function or adding a new function in a developed model by using IACF can increase the model development efficiency and facilitate researchers focus more on the problem mechanism instead of program coding and modification.

V 謝 誌 感謝吾師張良正教授對於本論文的細心指導及求學階段中對學 生待人接物與學問研究的啟發，使得學生在各方面皆能有所成長。承 蒙口試委員海洋大學李光敦教授、台灣大學李天浩教授、中央大學蔡 武廷教授、交通大學單信瑜教授、洪士林教授及吳宗信教授、中州技 術學院蕭金財教授、成功大學詹錢登教授細心認真指正拙文，並於口 試期間提出相當寶貴的意見，使本文更加完備，在此謹致衷心謝意。。 在研究室方面，由衷的感謝陳宇文學長在論文及程式上的協助與 指導，讓我在博士班就學期間獲益良多，葉明生學長值得效法的實事 求是態度，何智超學長帶有顏色的幽默、陳鴻輝學長圓融的處世態度 以及NAPL 小組的革命情感，接著感謝弼舜在我最緊急的時候伸出援 手，讓我論文得以及時完成，還有各位可愛的學、弟妹：浚瑋、昀柏、 婉茹、怡釧、柏成、昱維、偽善、阿牛、瀚聖、佑誠、阿海、韋圻、 雲直、阜峻、冠宇、小瑜、阿布、深惠等，有你們的陪伴，讓我增添 許多研究室的美好回憶。 最後非常感謝我的家人對我的全力支持，在這漫長的六年時光內， 替我加油打氣，以及感謝我的女友賢柔對我始終如一，一路與我患難 與共。僅將此研究獻給我最愛的家人及週遭支持、關心我的人。

VI

目錄

摘要 ... I Abstract ... III 謝 誌 ... V 表目錄 ... IX 圖目錄 ... XII 符號使用對照清單 ... XVI 第一章 前言... 1 1.1 研究動機 ... 1 1.2 研究目的 ... 3 1.3 文獻回顧 ... 4 第二章 系統架構與研究流程 ... 12 2.1 傳統建模流程與可適性計算架構建模流程說明 ... 12 2.2 系統概念架構 ... 18 2.3 研究流程 ... 20 第三章 可適性計算平台建置 ... 22 3.1 可適性計算平台程式架構 ... 22

VII 3.2 可適性計算平台之基本求解原理說明 ... 25 3.3 主控協調層建置 ... 27 3.4 計算核心子系統建置 ... 33 3.5 方程式集合分析模組建置 ... 47 3.6 應用模組儲存容器 ... 62 3.7 變數與參數資料結構 ... 63 第四章 地下水流、熱流及污染物傳輸模式建置 ... 74 4.1 數學模式說明 ... 75 4.2 數值離散說明 ... 113 第五章 可適性計算平台模式驗證 ... 122 5.1 暫態垂向二維地下水流、熱流及污染傳輸耦合模擬驗證(一) ... 123 5.2 暫態垂向二維地下水流、熱流及污染傳輸耦合模擬驗證(二) ... 132 第六章 智慧型可適性計算平台建置 ... 149 6.1 智慧型可適性計算說明 ... 151 6.2 智慧型可適性計算簡例 ... 152 6.3 智慧型可適性計算平台系統架構及開發流程說明 ... 154

VIII 6.4 智慧型可適性計算成果說明及驗證 ... 171 第七章 結論與建議 ... 185 7.1 結論 ... 185 7.2 建議 ... 186 參考文獻 ... 188

附錄A Voronoi Diagram 與 CGAL 函式庫 附錄B 總水頭計算公式推導 附錄C 土壤壓密係數與傳統儲水係數之相關推導 附錄D 細胞自動機 附錄E Voronoi Diagram 與可適性計算應用於熱流及污染傳輸 之擴散計算 附錄F 水之熱動力性質公式說明

IX

表目錄

表3.1 圖3.7(B)之節點、相鄰節點、相鄰連結與相鄰邊界虛擬點對應 表... 44 表 3.2 簡例之變數方程式關係矩陣 ... 51 表 3.3 節點與局部參數隸屬關係對應表 ... 71 表 4-1A 地下水流變數與方程式關係矩陣 ... 88 表 4-2 地下水流變數列表及自變數與應變數統計表 ... 89 表 4-3 地下水流函式求解順序表 ... 90 表 4-4A 參考TOUGH2 所建模式熱流傳輸變數與方程式關係矩陣 ... 99 表 4-4B 參考HST3D 所建模式之熱流傳輸變數與方程式關係矩陣 100 表 4-5A 參考 TOUGH2 之熱流傳輸變數列表及自變數與應變數統計表 ... 101 表 4-5B 參考HST3D 之熱流傳輸變數列表及自變數與應變數統計表 ... 102 表 4-6 熱流傳輸函式求解順序表 ... 103 表4-7A 參考 TOUGH2 所建模式之溶質傳輸變數與方程式關係矩陣 ... 109 表 4-7B 參考HST3D 所建模式之溶質傳輸變數與方程式關係矩陣 ... 110 表 4-8A 參考TOUGH2 所建模式之溶質傳輸變數列表及自變數與應變

X 數統計表 ... 111 表4-8B 參考 HST3D 所建模式之溶質傳輸變數列表及自變數與應變 數統計表 ... 112 表 4-6 溶質傳輸函式求解順序表 ... 113 表 5.2 地下水流模擬局部參數設定表 ... 125 表 5.3 地下水流模擬全域參數設定表 ... 125 表 5.4 熱流模擬局部參數設定表 ... 127 表 5.5 污染傳輸模擬局部參數設定表 ... 128 表5.6 案例一模擬初期及模擬結束時，ACF 與 TOUGH2 之最大相對誤 差百分比整理表 ... 132 表 5.7 地下水流模擬局部參數設定表 ... 134 表 5.8 熱流模擬局部參數設定表 ... 135 表 5.9 污染傳輸模擬局部參數設定表 ... 135 表5.20 案例二模擬初期及模擬結束時，ACF 與 HST3D 之最大相對誤 差百分比整理表 ... 140 表5.21 案例三模擬初期及模擬結束時，ACF 與 HST3D 之最大相對誤 差百分比整理表 ... 144 表 5.22 案例二與案例三之PECLET NUMBER計算結果 ... 146 表 6.1 智慧型可適性計算與可適性計算之模擬結果分析表 ... 174

XI

表 6.2 ANN 訓練及驗證資料之範圍說明表 ... 182 表 6.3 ANN 訓練過程表 ... 182 表 6.4 模擬結果最大誤差百分比整理表 ... 183

XII

圖目錄

圖1「傳統數值模擬方法」開發流程圖 ... 3 圖2.1「可適性計算架構建模流程」與傳統建模流程之比較圖 ... 17 圖2.2 可適性計算架構建模流程圖 ... 18 圖2.3 智慧型可適性計算概念架構圖 ... 20 圖2.4 研究流程圖 ... 21 圖3.1 可適性計算平台之系統示意圖 ... 24 圖3.3 可適性計算平台擴散與傳流迭代流程圖 ... 30 圖3.4 可適性計算架構擴散項外迭代模組流程圖 ... 33 圖3.5 可適性計算架構擴散項內迭代流程圖 ... 36 圖3.6 可適性計算架構資訊同步示意圖 ... 38

圖3.7(a)凡諾依圖與 Delaunay Triangle 示意圖 ... 43

圖3.7(b)凡諾依圖之網格架構圖 ... 43 圖3.8 材質與網格分布示意圖 ... 45 圖3.9 材質與網格分布示意圖 ... 46 圖3.10 源匯點與網格配置示意圖 ... 47 圖3.11 方程式相依關係圖 ... 49 圖3.12 應用樹狀搜尋法於方程式解路徑分析示意圖 ... 54

XIII 圖3.13 樹狀搜尋法分析邏輯流程圖 ... 59 圖3.14 方程式樹示意圖 ... 60 圖3.15 多應用模組方程式集合分析流程圖 ... 62 圖3.16 資料儲存結構示意圖 ... 64 圖3.17 節點變數資料結構與函數運作關係圖 ... 65 圖3.18 連結變數資料結構與函數運作關係圖 ... 66 圖3.19 邊界穿越點變數資料結構圖 ... 67 圖3.20 空間變數或參數示意圖 ... 71 圖3.21 局部參數與隸屬節點關係示意圖 ... 72 圖3.22 插入式 (Plug In) 數值模式建構系統 ... 73 圖5.1 案例一之模擬區域及空間切割示意圖 ... 126 圖5.3 熱流傳輸模擬驗證－溫度場分布 ... 130 圖5.5 地下水流模擬驗證－壓力水頭分布 ... 131 圖5.7 污染物傳輸模擬驗證－濃度場分布 ... 132 圖5.9 熱流模擬驗證－溫度場分布 ... 137 圖5.10 污染物傳輸模擬驗證－濃度場分布 ... 138 圖5.11 地下水流模擬驗證－壓力水頭分布 ... 138 圖5.13 污染物傳輸模擬驗證－濃度場分布 ... 139 圖5.14 壓力水頭模擬驗證－壓力水頭分布 ... 141

XIV 圖5.16 污染物傳輸模擬驗證－濃度場分布 ... 142 圖5.18 熱流模擬驗證－溫度場分布 ... 143 圖5.20 案例二有無考量延散機制之污染傳輸模擬結果比較 ... 147 圖5.22 案例三有無考量延散機制之污染傳輸模擬結果比較 ... 148 圖6.1a 初始時刻污染物與熱流分佈圖... 153 圖6.1b 第二個時刻初之污染物與熱流分佈圖 ... 153 圖6.2 智慧型可適性計算架構圖 ... 155 圖6.3 專家系統與計算元溝通示意圖 ... 160 圖6.4 倒傳遞類神經網路架構 ... 163 圖6.5 以神經網路建置方程式示意圖 ... 165 圖6.6(a) 計算元推論法則示意圖 ... 169 圖6.6(b) 計算元狀態分析推論鏈 ... 170 圖6.7 IACF 於時刻 9000 秒時，啟動計算污染傳輸之計算元分布圖 ... 173 圖6.8 IACF 於時刻 2700 秒時，啟動計算污染傳輸之計算元分布圖 ... 176 圖6.9 IACF 於時刻 27900 秒時，啟動計算污染傳輸之計算元分布圖 ... 177 圖6.10 IACF 於時刻 31500 秒時，啟動計算污染傳輸之計算元分布圖 ... 177 圖6.11 計算元狀態分析模組推論過程示意圖 ... 178

XV

圖6.12 ANN 與方程式(4-3A)之計算密度比較圖 ... 183 圖6.13 地下水流模擬驗證－壓力場分布 ... 184 圖6.14 熱流傳輸模擬驗證－溫度場分布 ... 184

XVI

符號使用對照清單

符號使用對照表 符號 對照含意 符號 對照含意 α′ 土壤壓密係數 α 縱向延散性 α 橫向延散性 α 溶質膨脹係數

Van Genuchten 公式參數 Van Genuchten 公式參數

純水黏滯係數 溶液黏滯係數 水流之單位質量熱容量 土體之單位質量熱容量 溶液密度 純水密度 常溫常壓下之純水密度 土體密度 鹽水濃度 參考鹽水濃度 殘餘含水量 飽和含水量 有效含水量 溶液中之溶質濃度(重量 百分濃度) 濃度最大值 濃度最小值 D 機械延散係數張量 D 熱延散張量 D∗ 水動力延散係數張量 重力加速度 位置水頭 ℎ 總水頭 位置高程 熱穿越量 ℎ 單位面積擴散項熱穿越量 ℎ 單位面積對流項熱穿越量 , 等效熱傳導係數 , 水流熱傳導係數 , 土體熱傳導係數 溶液之實際水力傳導係數 相對水力傳導係數 土壤之飽和水力傳導係數 溶質穿越流率 水流質量流率 , 單位面積擴散項溶質穿越 流率 , 單位面積對流項溶質穿越流 率 單位面積水流質量流率 土壤孔隙率 接下頁

XVII 符號 對照含意 符號 對照含意 流線之正交方向 常壓下之土壤孔隙率 壓力 以 bar 為單位之壓力數值 壓力水頭 以 cm 為單位之壓力水頭數 值 常壓之壓力水頭數值 飽和蒸汽壓力 抽注溶質質量 抽注熱容量 抽注水質量 ̂ 流線方向 飽和度 _, 溶液比熱 , 土體比熱 攝氏溫度 絕對溫度 華氏溫度 達西流速在 方向之分量 達西流速在 ̂方向之分量 達西流速 控制表面 控制體積 ℑ () 代表單位轉換運算子 網格內平均流速

1

第一章

緒論

1.1 研究動機 隨著電腦科技的長足進步，電腦數值模式之計算能力亦愈趨強大， 並已成為科學與工程領域不可或缺的研究與設計工具。數值模式最大 的好處之一為其可重複利用及可操作性，因此可節省大量人力與物力， 加速研究分析的進行。惟隨著數值模式的日趨複雜，其開發與維護的 複雜度亦快數增加，而需大量的資源投入。造成此種現象的原因，除 了問題本身難度增加外，傳統上數值模式開發的缺乏彈性亦為其主要 因素之一。 傳統上數值模式之開發流程如圖1 所示，一般而言可分為 4 個階 段： 1. 概念模式階段(Conceptual model)：對所探討的問題，定義其 涵蓋範圍，進行必要的簡化與假設並對其基本現象進行描 述。 2. 數學定義階段(Mathematical formulation)： 基於前述之概念 模式，進一步以量化的數學式定義問題，通常包含兩個步 驟： a. 對於問題中之各種基本現象或行為以一個或多個數學式 （包括變數及方程式）描述之。 b. 將前述多個變數及方程式以人為方式推導整合為數目較 少之變數及控制方程式，較普遍的為偏微分方程式(PDE）。

2 3. 數值離散階段(Numerical discretization)：定義空間格點，將 原時空上連續之解函數以時空上不連續之離散型變數近似， 再選用適當的數值方法，如有限差分法（FDM），有限元素 法（FEM），有限體積法（FVM），邊界元素法(BEM)及有限 解析法（FAM）等，將前述之控制方程式(組)，進行空間及 時間上之數值積分推導，而將原來之控制方程式(組)表達成 空間格點上之離散型代數方程式。 4. 程式開發階段(Program development)：前述離散型代數方程 式各變數之數值，目前常以計算機(電腦)計算，因此需進行 相關程式的開發，主要有兩大部份；先將所有空間格點上之 離散型代數方程式，組合成矩陣方程式，此矩陣方程式描述 了所有空間格點上離散變數彼此間之時空關係。前述之矩陣 方程式，仍需再以適當之矩陣解法求解各變數之數值。 基於上述建構流程所開發之數值模式本研究在後續文章中稱之 為「傳統數值模式」。由前述開發流程可知，傳統數值模式之開發乃 環環相扣，若欲加入一原模式未考量之變量或現象，則須加入描述此 行為之規則（方程式），如此則須重新推導並修正原問題之控制方程 式、數值解及矩陣方程式，相關電腦程式亦需要重新修改，修改工作 往往相當繁複，亦即前述開發流程雖然定義清楚問題明確，然而卻難 以修改，因此限制了原模式擴充或修正的彈性。。因此，傳統數值模 式，只能模擬其原先控制方程式所描述之現象，而無法擴充本身之模 擬功能，且由於缺乏擴充或修正之彈性，而不易與人工智慧方法結 合。

3 圖1「傳統數值模擬方法」開發流程圖 1.2 研究目的 本研究旨在提出一個新的數值建模方法與實作架構，以解決傳統 數值建模費時費力與所開發之計算模擬軟體，難以更新、擴充及累積 計算功能的問題，並將數值模式進一步結合人工智慧，使其成為具有 調適能力之智慧型計算模擬系統。

4 1.3 文獻回顧 由前述研究目的可看出本研究包含數值建模方式的突破及人工 智慧方法在數值模擬上的應用等兩大層面的議題，以下將今別就此兩 層面進行相關研究的探討。 1.3.1 數值建模方式的改進 在各種建立模式的方法中，以微分方程式建模的方式歷史最為悠 久，已有三百多年的歷史，使用微分方程式建模的好處是在推導微分 方程式的相關技巧已發展成熟，然而在求解解析解的過程皆需靠人力 來進行所有的計算且只有少數的微分方程式有Closed-form solution， 其求解相當不易。近幾十年來數值計算蓬勃發展，研究學者紛紛利用 數值方法求解微分方程式，尤其在偏微分方程式(PDE)之求解方面， 成效卓著，因此數值方法已被認定是數學不可或缺的一部份，但若要 使用數值計算求解一個問題，則至少必須經過三個步驟: (1)將物理行 為轉化為微分方程式(2)使用適當的數值離散方式將這些微分方程式 進行時間與空間的離散，如有限差分法(3)撰寫數值程式 Toffoli(1984)[1]。 Toffoli(1984)[1]認為以上述三個步驟建立數值模式 是一個迂迴而不直接的方法。此外，周成虎等(1999)[2]亦認為雖然以 偏微分方程式(PDE)可靈活的進行複雜數學公式運算，得到精確的定

5 量解，但當使用電腦解決PDE 問題時，由於電腦的運算是架構於離 散基礎上，因此PDE必須對算式本身進行時空離散化，將其展開成冪 次方程式並可能忽略較高階之項次，建立差分方程式，或採用其他轉 換方式，建立離散結構來表示其連續變化量。這些轉換過程不僅是繁 複的，甚至有可能是無法解決的。但最重要的是在這個過程中，PDE 也失去了它最重要的特性－精確性、連續性。Ioannis(1997)[3]亦說明 使用偏微分方程式(PDE)建立模式時，若考慮過多複雜的因子時，則 會大幅增加使用PDE 描述其行為的困難度，無法完整的描述問題中 之所有行為，如海洋中之油污傳輸(宋鴻均，2005) [4]。 不少國內外學者為了試圖解決上述以傳統偏微分建模搭配數值 方法求解之種種問題或限制，因此嘗試以細胞自動機建立數值模型 (Cellular Automata)( 金傳春等人, 2007[5]；Machta, J., 2006[6]：宋鴻 均, 2005[4]：Fox, 2003[7]；Batterman, 2002[8]；Toffoli and Margolus, 1990[9]；Morgan and Morrison, 1999[10]；Papatzacos, 1989[11])。Toffoli and Margolus (1990)[9]指出CA乃直接將巨觀現象背後之基本局部互 動行為(Basic local interaction) 轉換為CA規則，以此建模即可演化產 生原巨觀現象，因此無須推導微分方程式，故不需進行數值離散。在 建立數學模式方面，CA僅需將所有運動機制轉換為CA之規則即可， 因此相較於PDE更易於對原物理系統進行建模 (Toffoli,1984[1]；

6 Omohundro, 1984[12])。宋鴻均(2005)[4]選擇CA建立一海洋油污傳輸 模式，其研究指出若選擇以PDE方式建模，則須先建立足以描述該問 題之偏微分方程式，而其問題囊括之現象十分複雜，若要全部加入考 量則其PDE推導將十分困難，而透過細胞自動機強大的適應性，可較 PDE 更方便的加入邊界、風、流、岸邊附著等影響因子，成功達成 對油污染現象的模擬。金傳春等人( 2007)[5]以CA建立流行病傳染模 式，該研究亦指出與PDE建模方式相比，細胞自動機模型結構簡潔靈 活，計算方式簡單且，模擬過程直觀，易於分析。陳宇文(2009)[13] 發展一數值計算平台，嘗試直接以方程式集合取代傳統偏微分方程式 進行求解，雖其模式成功進行地下水流、熱流與污染傳輸模擬，但其 模擬結果未經驗證，且其在方程式之計算順序分析上亦只能考量單一 物理量問題，當需考量求解多物理量耦合問題時，除需求得各物理量 之方程式求解順序外，仍需物理量間之變數傳遞關係，該研究並未提 及此部份之解決方法，此外，該研究在方程式分析上乃利用人工分析， 未建立自動化分析功能，故該研究仍有相當改善空間。 由前述研究可知，許多研究皆試圖破除以偏微分方程式建模的方 式，使得模式得以同時考量更完整的機制，以較偏微分方程式簡易的 方式建模，惟這些研究皆只致力於如何解決對當下已定義好的問題更 有效率的建模，而未進一步考量如何將原先開發好的模式，更容易的

7 進行擴充修正。例如以前述CA之應用研究為例，當模式開發者要在 原模式加入新的機制時，模式開發者除了增加描述該機制之方程式外， 還須調整該模式各方程式之計算順序及其對應之程式結構；若需增加 考量其他物理量，則模式開發者除了需增加描述該物理量之多條方程 式外，還需人工分析耦合求解時，物理量間變數如何互相傳遞，而這 些皆需以人工方式進行，亦即問題若有任何修正，則皆仍然要如傳統 建模方式般耗費相當人力進行程式修改。總之，前述研究並未考量如 何建立一個一般化的建模流程或平台，使得模式可以漸進、有效率且 系統化的增加模擬機制。 1.3.2 人工智慧方法在數值模擬的應用 此部份又可分為兩大類，一類為人工智慧方法與模式乃是互相結 合的型態，此種方式人工智慧方法與數值模式乃共同存在互相輔助， 另一類為直接以人工智慧方法取代數值模式的角色進行模擬。兩者以 往皆已有相當的研究成果。 在數值模式結合人工智慧方面，其常用的人工智慧方法有專家系 統及啟發式演算法等。專家系統與數值模式的結合方式，通常為兩層 式架構，外層為專家系統之交談式介面(Interactive interface)，其常為 輔助使用者挑選符合問題的數值模式、輔助建立數值模式或進行參數

8

檢定等(Uzel et al., 1988[14]; Abbott, 1991[15]; Chau and Chen, 2001[16]; Chau and Albermani, 2002[17], 2003[18]; Chau, 2004[19]; Kim, 2007[20])，如Uzel[14] 等人(1988)建立之專家系統可以輔助與訓 練使用者建立計算流體力學模式(PHOENICS)，並對於後續模擬結果 進行解釋；Chau and Chen(2001)[16]建立之專家系統可輔助海岸水量 與水質等相關模式建置並可協助所建置模式進行參數檢定；Kim等人 (2007)[20]使用專家系統(Expert system for calibration of HSPF,

HSPEXP)協助檢定水文模式(Hydrologic Simulation Program Fortran, HSPF)；Madsen等人(2002) [21]應用專家系統幫助檢定降雨逕流模式； Chau(2003[21], 2004[23])使用專家系統協助檢定海岸及其相關模式； 陳韋圻(2008)[24]應用專家系統於地下水流模式，同時檢定水力傳導 係數(K值)及儲水係數(Sy值)；王雲直(2010)[25]應用專家系統檢定濁 水溪沖積扇地下水模式。 啟發式演算法與數值模式之結合，通常亦為兩層式架構，其底層 為數值模式，外層為啟發式演算法，而啟發式演算法一般皆作為最佳 化工具(Optimization tools) (Haddadet al., 2006[26])。啟發式演算法如 基因演算法(Genetic Algorithm)、模擬退火演算法(Simulated

Annealing)、禁忌搜尋法(Tabu Search)、蟻群演算法(Ant Algorithm)及 蜜蜂演算法(Honey Bees Algorithm)等，目前已被廣泛的應用於科學、

9

工程及商務等各領域中(McKinney and Lin, 1994[27]; Gen and Cheng, 1997; Wardlaw and Bhaktikul, 1999[28]; Castanedo et al., 2006[29])。在 地下水相關研究領域中，啟發式演算法與地下水模式結合之應用包括 地下水參數檢定(Zheng and Wang, 1996[30]; Abbaspour et al., 2001[31]; Lin et al., 2001[32]; Tsai et al., 2003[33])、地下水管理 (McKinney et al., 1994[27]; Bhallamudi et al., 2003[34]; Chu and Chang, 2009[35])、地下 水復育(Ahlfeld et al., 1998[36]; Guan and Aral, 2004[37]; Chan et al., 2000[38], 2005[39]; Espinoza et al., 2005[40])及地下水污染源推求等 (Datta and Chakrabarty , 2003[41]; Chang, and Yeh, 2005[42]; Singh and Datta, 2006[43] ;Dhar and Datta, 2007[44])。

前述人工智慧與數值模式整合方面的研究，專家系統與啟發式演 算法大都位於數值模式外層，並未將人工智慧整合入數值模式之核心 運算中，因此並未實質增加數值數式本身的模擬或其他能力，無法使 數值模式易具有自我調適等類似「智慧」的能力。 在以人工智慧方法直接進行的模擬方面，則以類神經網路為代表。 神經網路因計算快速及具有建構非線性模型能力之優點，往往作為取 代數值模式之工具。通常作法之一乃藉由數值模式之輸入、輸出及相 關參數等資料進行神經網路訓練，取得原數值模式之特性並取代之， 以利後續加值應用(Malekmohamadi et al., 2008[45];Chu and Chang,

10

2009[46]; Ioannis, et al., 2010[47]) ，在此數值建模將只是一個過程。 Malekmohamadi等人(2008)應用類神經網路取代波傳數值模式

(Numerical wave model)，可藉由風速預測波高。Chu與Chang(2009)[46] 應用類神經網路代替地下水數值模式並結合微分動態規劃模式優選 動態地下水操作策略及成本最小化問題；Ioannis等人(2010)[47]以類 神經網路取代地下水數值模式並結合微分進化演算法(Differential evolution algorithm)求解一區域不違反安全出水量前提下之最大抽水 量問題，黃浚瑋( 2008 )[48]應用類神經網路代替地下水數值模式並結 合遺傳演算法優選地表地下聯合營運之操作策略。另一種情況則直接 以觀測資料，對未釐清的物理、化學或生物行為，透過類神經網路學 習模擬其行為特性Chau(2006)[49]。Coppola (2003)[50] 和Nayak (2006)[51]分別採用坦帕灣、孟加拉灣當地地下水含水層實際水文資 料來訓練類神經網路，並以訓練好之類神經網路來進行水位連續預測， 完成訓練之類神經網路即為具當地代表性之地下水模式。 為了解決前述文獻所指出的模式擴充不易及雖然與人工智慧方 法整合但模式仍不具智慧功能等兩大問題，本研究提出一新的計算架 構，突破傳統數值模式的開發方式，建立一系統化、自動化之數值建 模平台，其中包括發展一分析方程式之方法以及建置一個具可擴充方 程式之建模系統。此方程式分析方法，乃可分析多條方程式之一致性

11 及計算順序，包括多物理量問題變數之傳遞關係，其為使系統具有可 擴充性的核心機制，應用此可擴充模擬功能的系統所開發之模式，可 隨著需求逐步增加模式模擬能力。並在此基礎下將人工智慧方法如專 家系統與類神經網路等，整合入模擬模式核心，而使模式具有「人工 智慧」能力，如可以因應環境變化調整自身模擬範疇等。

12

第二章

系統架構與研究流程

為了突破傳統數值建模之限制，本研究提出一新型態之計算架構， 稱為「可適性計算架構」，而此新計算架構的提出，乃是源自於對傳 統數值建模的重新思考，因此本章在說明新計算架構之前，將仍先以 數值建模發展流程的觀點來比較說明傳統方式與本研究的差異，接著 再說明支持本研究新的模式發展型態，所需系統平台之資訊系統架 構。 2.1 傳統建模流程與可適性計算架構建模流程說明 圖 2.1 為「可適性計算架構」與傳統數值模式建構流程之差異比較， 由圖可知本研究與傳統方式差異甚大，是整個開發流程的改變。圖 2.2 則經整理後本研究提出之數值模式建立流程。圖 2.1 之說明如下： 1. 概念模式(Conceptual model):此階段乃在定義問題及其中之 基本行為描述，因此兩者必須相同。 2. 數學定義(Mathematical formulation)：此階段又可分為兩大步 驟，首先為以各種基本方程式描述第一階段所定義問題中的 各種變量之變化行為，以正確量化描述問題，此本研究與傳 統方式皆相同。而本研究將與傳統方式不同之處為，本研究

13 將保留前述之基本方程式組，並以其做為後續計算之基礎而 不採用如傳統方式般，將前述所得之多條方程式進行人為的 整合推導，以盡量減少方程式及變數之數目。傳統方法的好 處為最後待解的方程式與變數數目較少。傳統上認為方程式 及變數愈少對後續解題將愈有利，惟以本研究的觀點而言， 此種推導與簡化需以人力為之而無法自動化或程序化，因此 傳統方法若欲考慮新的現象或變量而需新增或修改前述之基 本方程式時，則需人為重新推導，而此點將先天上限制了後 續所開發之模式擴充模擬功能的可能。本研究為維持後續增 減方程式的彈性，將不再對前述多條基本方程式進行進一步 整合。如此本研究後續將需處理較多的方程式與變數，惟由 於各方程式並未再經人為整合推導，將可維持其原先較基本 而簡單的型式，而不像傳統方式般，雖然最後面對的方程式 數目較少，惟這些方程式已經人為整合後所得，其型式將較 為複雜。傳統方法常需面對二階以上之微分方程式，惟本研 究的方法，所處理的方程式絕少高於一階微分方程式。方程 式的複雜度亦將直接影響下一步驟數值離散之難易，一階微 分以下之方程式可以很簡單的方式進行離散，若是二階微分 以上之方程式則其數值離散之難度將大為提高。除了不進行

14 方程式整合外，為維持將來增減方程式的彈性，這些基本方 程式將來亦不採同時聯立求解，而是依變數間之相互關係循 序逐條求解，如未來若有方程式的增減，則只需重新定義方 程式計算順序即可。爰此，本研究將在第三章提出方程式分 析的實作方式，以檢驗方程式與變數關係的一致性，同時決 定方程式之計算順序。 3. 數值離散(Numerical discretization) ：在此步驟傳統方法由於 常需面對較高階的微分方程式，因此需要複雜的數值離散方 法，如有限元素法(Finite Element Method, FEM)，有限差分法 (Finite Difference Method, FDM)等皆是。本研究由於絕少面對 一階以上之微分式，因此只需簡單的一階差分，即可對原方 程組進行離散化。另外，為維持增減方程式的彈性，本研究 為了更易於進行離散化，並建議以凡諾依多邊形進行空間格 網之定義。 4. 程式開發 (Program development)：延續步驟 3，傳統方式為 照顧各節點間之互相影響，必需將所有節點上之離散方程式 與變數，再組合成矩陣方程式，接著發展各種不同矩陣解法 之計算程式，對所組成的矩陣方程式，計算其數值解。本研 究將不若傳統數值方法組合所有節點上的變數形成矩陣方程

15 式，反之，本研究將以單一節點上之多條離散方程式為主進 行計算，搭配方程式分析後所得方程式計算順序，再以各節 點互相交換資訊並配合疊代計算的方式，求得整個問題中所 有節點的解。 圖 2.1 中之數學模式建立乃以飽和地下水問題為例，圖中所描述之 方程式為一維穩態的地下水流相關方程式，K 為水力傳導係數，假設 為常數，v 為達西速度，h 為總水頭，若以傳統數值方式建模則需將 其推導成偏微分方程式，其形式為一Laplace 形式之偏微分方程式， 之後再選用適當數值方法(如有限差分法)對此偏微分方程式進行時 空離散，得到所有節點之離散方程式後再建立矩陣方程式，選擇適當 之矩陣解法後即可撰寫數值程式進行求解。若使用者欲建立一未飽和 地下水流問題，則K 為水頭之函數(K=f(h))，其偏微分方程式則需重 新修正，原本為Lapalce 形式之偏微分方程式改為式(2.1)，其形式為 一非線性偏微分方程式，實務上經常會再以式(2.1)為基礎， 進一步假設與推導而得常用之理查方程式(Richard Equation)。因其偏 微分方程式已經修正，故在數值離散與數值模式兩個部分皆需進行大 幅修改，由此可知，即便使用者僅需改變一條方程式，卻需進行整個 模式的重新開發。反觀本研究所提出之開發流程，因直接採用方程式 集合進行求解，修改一條方程式僅需調整方程式之求解順序，以及針

16 對該方程式再進行離散並撰寫其之程式碼即可完成模式修改，此即本 研究所提之「模式擴充性」。此外，就數學模式本身的特性而言，傳 統數值方法則需將相關方程式全部整合入守恆方程式中再推導得控 制程式，使控制方程式中只有待解變數及相關參數，因此若要抽換或 增減任何一條方程式必定要重新推導偏微分控制方程式，並大幅修正 整個數值模式。而本研究無需整合相關方程式，只需將各方程式計算 後之值代入守恆方程式中，因此方程式可依使用者需求進行抽換或增 減，如前述水力傳導係數K 的修正即是，故本研究在模式之修正或 擴充上較傳統數值模式更具彈性。 ( ) ∙ + (ℎ) (2.1)

17

18 概念模式(Conceptual model) 包含問題定義及問題中基本行為之描述。 數學定義(Mathematical formulation) 1.對於問題中之基本行為以一個或多個數學量化的法則 （方程式）描述之。 2. 程式開發(Program development)建立 數值離散(numerical discretization) 以簡單差分方法，配合空間格網定義，建立計算節點上 前述多條方程式之離散方程式(組)。 圖2.2 可適性計算架構建模流程圖 2.2 系統概念架構 傳統上數值計算皆需在待求變量(變數)所在的空間上，畫分格網， 且將各變數的值表達在各格點上，惟最後仍以矩陣的方式整體計算， 換句話說，各格點的意義主要是空間上的座標及其對應的變數。而本 研究主要的計算皆在格點上，因此可將格點代表的意義擴充，各格點 不僅具有空間座標及變數，亦為一個具有獨立運算能力的個體，本研

19 究將其名為「計算元」。本系統為一兩層式架構，包括「主控協調層」 以及「計算層」。為了使數值模式具有智慧，本研究將可適性計算架 構與人工智慧結合，稱之為「智慧型可適性計算架構」 (圖 2.3)。進 一步說明如下： a. 計算層：本研究計算層由計算元組成，各計算元與其鄰近計 算元的相鄰關係，亦如以往由格網定義，計算元的功能為最 底層的數值計算，各計算元皆為獨立計算，其計算方式乃直 接計算多條分散的方程式，本研究將此多條分散的方程式稱 之為「方程式集合」。另外，本研究輔以「多應用模組方程式 集合分析工具」，用以分析方程式集合是否可解以及決定方程 式計算順序，使得各計算元可容易的以人工輔助的方式增加 或改變計算能力，此亦為本研究所提「可適性」之意義，此 部分將於第三章進一步說明。 b. 主控協調層：本研究主控協調層由「協調委員」組成，「協調 委員」可視為另一種計算元，惟其位階乃位於計算層之上， 如圖2.3 所示。各協調委員皆有其職責，如調整系統資源分配 或調整計算元之計算行為等；相較於計算層之主要功能為計 算求解，主控協調層之主要功能為監督及分析所有計算元之 計算狀態，並依分析結果適時調整其計算行為。主控協調層 亦可結合專家系統與機器學習等人工智慧方法增強其功能。

20 圖2.3 智慧型可適性計算概念架構圖 2.3 研究流程 本研究之研究流程圖如圖 2.4 所示，首先為開發智慧型可適性計算 平台，此部分本研究將其分為兩階段進行，第一階段先開發可適性計 算平台，第二階段再開發智慧型可適性計算平台。接著為應用前述步 驟所開發之兩種平台開發地下水流、熱流傳輸與污染質傳輸之耦合求 解數值模式。最後再以伯克利勞倫斯國家實驗室(Lawrence Berkeley National Laboratory) 研發之 Tough2 以及美國地質調查局(USGS)所開

21 發之HST3D 計算模擬軟體進行驗證。 可適性計算平台建置的部分將於第三章說明，智慧型可適性計算平 台的建置將於第六章說明。本研究開發數值模式之方式為將描述問題 所需之方程式儲存於計算平台上即可完成建模，而本研究應用前述計 算平台所開發之地下水流、熱流傳輸與污染質傳輸模式所需之方程式 集合將於第四章進行說明。可適性計算平台所開發之數值模式將於第 五章進行驗證，而智慧型可適性算平台所發展之數值模式將於第六章 說明。 圖 2.4 研究流程圖

22

第三章

可適性計算平台建置

依據前述第二章所提系統架構(如第二章之圖 2.3)可知，系統包含 三大部分，分別為計算層、主控協調層及人工智慧，其中計算層及主 控協調層之實作本研究稱之為可適性計算平台，本章將說明如何實作 此計算平台，另外，人工智慧的實作將於第六章說明。 3.1 可適性計算平台程式架構 可適性計算平台之程式系統圖如圖 3.1 所示，包括計算層與主控協 調層。計算層可再細分為「核心計算子系統」、「應用模組容器」以及 「方程式集合分析模組」。 主控協調層之位階乃在計算層之上，其內包含數個協調委員，協 調委員數量並無限制，可依模式開發者需求擴充協調委員數目。各協 調委員皆有其職責，如調整系統資源分配或調整計算元之計算行為等； 目前可適性計算平台所建置之協調委員主要任務為負責跨計算元間 之相關迭代計算。 在計算層的部分，「核心計算子系統」為計算平台之主要計算工具， 包括數個通用化之模組，如空間離散模組、擴散項內迭代模組、資訊 同步模組，傳流項計算模組，方程式計算順序設定及儲存容器、材質 參數儲存容器、邊界條件設定與儲存容器及源匯項設定與儲存容器。

23 「應用模組儲存容器」為預留之資料夾以儲存應用本平台於特定問題 時，開發者所需建立之方程式集合。本研究共建置了「地下水流」、「熱 流傳輸」及「溶質傳輸」三個應用模組，各自儲存相對應之方程式集 合，將於第四章中說明。「方程式集合分析模組」主要目的為輔助模 式開發者能有效率的建立數學模式。本研究解題乃以描述原問題之多 條方程式直接求解，而不需再推導成一整合型態之偏微分方程式 (PDE)，因此需以「方程式集合分析模組」驗證各應用模組之方程式 集合是否可解，以及定義方程式之求解順序，當多個應用模組須同時 求解時，需再提供不同應用模組間變數之傳遞資訊，此部分將於3.5 節中詳細說明。由圖 3.1 可知，此分析工具連結「核心計算子系統」 與「應用模組容器」兩子系統，為系統關鍵模組。「方程式集合分析 模組」在讀取各應用模組之資訊後，接著產生變數-方程式關係矩陣， 並以此關係矩陣進行後續分析，若各應用模組之方程式集合皆可解， 則可得各應用模組之方程式集合求解順序，之後再將方程式集合求解 順序儲存於核心計算子系統之「方程式計算順序設定及儲存容器」中。 在「主控協調層」與「計算層」外層有一殼層(shell)將此兩子系統包 圍，此殼層之職責為模式輸入與輸出之處理介面。 3.2 節將對可適性計算平台之基本計算原理進行說明，以方便後續 對各子系統實作之進一步之說明。3.3 節將進一步說明主控協調層。

24

3.4 節說明「核心計算子系統」。3.5 節說明「方程式集合分析模組」。 「應用模組容器」將於3.6 節說明。3.7 節則說明可適性計算平台之 「變數與參數資料結構」。

25 3.2 可適性計算平台之基本求解原理說明 在前述第二章之 2.1 節中已說明了本研究之模式開發流程與傳統方 法之差異，由第二章圖 2.1 可知，本研究在建立數值模式時，乃以方 程式集合進行計算而非以傳統偏微分方程式進行建模，其計算流程為 先將各方程式進行時空離散後，再依方程式分析結果決定求解順序， 之後以內外迭代法進行求解。在方程式分析的部分，本研究參考樹狀 搜尋法建立一分析方程式求解順序之演算法，此部分於3.5 節說明。 此外，內外迭代之流程概念說明如圖3.2 所示，首先為步驟(1)，本研 究在求解時乃以單一計算元為基本計算單位，各計算元獨立求解方程 式集合，在求解時先假設由相鄰計算元取得之資訊為真值，目標函數 設定為守恆方程式之誤差最小，之後再利用最佳化方法迭代求解該計 算元之狀態變數，此單一計算元進行最佳化求解之步驟稱為「內迭代」， 其詳細說明參見3.4.1 節。當所有計算元皆完成內迭代計算後，接著 需判斷各計算元在數值上是否已達穩定狀態，在步驟(2)處，本研究 設置「外迭代」步驟分析所有計算元是否皆已達到收斂狀態，分析方 法為比較各計算元前後兩次內迭代狀態變數之差值是否小於一預設 門檻值，若所有計算元皆已收斂，本研究稱之為「外迭代收斂」，此 部分說明詳見3.3.2 節。當系統達到外迭代收斂狀態時即可進入下一 時刻計算，若未達到外迭代收斂，所有計算元需重新啟動進行內迭代

26 計算，在計算之前，各計算元需先接步驟(4)，重新取得相鄰計算元 內迭代後之資訊，而啟動所有計算元進行內迭代計算之步驟本研究稱 之為「外迭代計算」。 此內外迭代方法不同於傳統數值方法之顯示法或隱式法 ，在一計 算元進行內迭代計算時，本研究先取得其相鄰計算元內迭代前之資訊 進行計算，此方式近似於顯示法，但之後若需重新進行外迭代計算時， 本研究則先重新取得相鄰計算元內迭代後之資訊再進行計算，此方式 則近似隱式法，然若以一標準顯示法之定義，乃指在計算該計算元之 狀態變數時，直接取用上一時刻之資訊進行計算不需迭代，而隱式法 乃直接取當下計算中之資訊進行計算；本研究在各時刻之第一次外迭 代計算時雖以上一時刻之資訊作為初始值，但在後續外迭代計算進行 時，則會更新相鄰計算元之資訊，故此內外迭代法並非傳統之顯示法， 而較近似於隱式法。

27 圖3.2 本研究之內外迭代流程概念說明圖 3.3 主控協調層建置 對流體型態之物理量進行守恆方程式計算時，控制體積之邊界流通 量可分為擴散型(diffusion)及傳流型(advection)之流通量，為了數值計 算之穩定與方便，本研究將前述3.2 節之內外迭代求解概念再加以改 進，將其分為擴散項與傳流項之計算，在後續實例應用時，本研究將 考量熱流與污染物等物理量。前述擴散項之計算流程，仍應用原先 3.2 節所定義之外迭代流程，本研究將之稱為擴散項外迭代模組，並 將於3.3.2 節說明，至於傳流項計算之說明請參閱 3.4.3 節。本研究先

28 將擴散項與傳流項分開計算，此兩項再互相迭代求解，而控制此運算 程序之模組，本研究稱之為擴散與傳流迭代模組，此部分將於底下 3.3.1 節說明。 3.3.1 擴散與傳流迭代模組 本節將說明本研究如何計算傳流項與擴散項，計算流程如圖 3.3 所 示。起始時間設定為t =tstr，第一個步驟為進行擴散項計算，包括水 流流場與熱流及污染物擴散項之計算，完成擴散項外迭代計算後，第 二個步驟為以水流流場進行熱流及污染物傳流項之計算，此部分由傳 流項計算模組負責計算，之後即得到考量擴散與傳流兩種機制的熱流 及污染物場分布。但在計算熱流及污染物之傳流項時，所採用之水流 流場仍是先前未考量傳流作用之熱流及污染物影響下所得之流場，故 在傳流項計算模組計算完後，需以新的熱流及污染物場重新計算水流 流場，意即再重新啟動外迭代計算流程，本研究目前此外迭代模組只 啟動水流流場計算，此步驟對應於流程圖之步驟 “水流擴散外迭代 2”， 而在計算水流流場時，水流之相關變數設定皆回到時刻初之值。當得 到新的水流流場時，熱流及污染物在傳流項的部分則需再作修正，因 此需重新啟動傳流項計算模組重算傳流項，當以此新的流場重新計算 熱流及污染物之傳流項時，熱流及污染物之相關變數皆須回復到首次

29 完成擴散項計算後之值。為判斷擴散與傳流迭代是否完成，本研究設 置一「擴散與傳流迭代指標」，其計算方式為將各計算元之狀態變數 在前後兩次迭代之差值除以初始條件中狀態變數之最大值，若所有計 算元上之狀態變數指標皆小於10-5即視為收斂，完成一個時刻之計算， 接著檢查是否抵達最終時刻(t_end)，如未滿足則時刻遞增(t ₌t_+Δt)，執 行下一個時刻之模擬；反之，則完成暫態模擬。

30 模擬開始 是否抵達結束時刻? 模擬結束 是 否 水流、熱流及污染物擴散項外迭代1 起始時刻( t₌ )t_str t t t= +Δ end t t≥ 傳流項計算 水流擴散項外迭代2 傳流項計算 是否達收斂標準? 是 否 反覆運算開始 1 = n 圖 3.3 可適性計算平台擴散與傳流迭代流程圖

31 3.3.2 擴散項外迭代模組 本研究在求解問題時，將擴散項與傳流項分開計算，其中水流、熱 流及污染物之擴散項計算皆為擴散項外迭代模組所掌控。當開始擴散 項外迭代計算時，此模組會先讓所有計算元取得相鄰計算元之資訊， 如相鄰計算元之水頭、溫度及污染物濃度等資訊，接著擴散項外迭代 模組啟動各計算元進行方程式集合之求解。本研究採用最佳化方法做 為各計算元求解方程式集合之工具，此以最佳化方法求解方程式集合 之步驟本研究稱為擴散項內迭代，詳細說明請見3.4.1 節。因各計算 元皆使用最佳化方法求解各自負責的方程式集合，因此在計算開始時， 所有計算元皆須先給定待解變數一初始猜值，再進行最佳化計算，而 在所有計算元皆完成第一次最佳化計算求解後(即完成擴散項內迭代 計算)，所有計算元之待解變數皆已改變，此時擴散項外迭代模組則 會檢驗各計算元之待解變數變化幅度，若所有計算元的變化幅度皆小 於設定的收斂標準則此擴散項外迭代收斂；若有一計算元之待解變數 變化幅度超過收斂標準，則外迭代重新啟動各計算元，在各計算元取 得相鄰計算元最新資訊後，重新進行擴散項內迭代。前述變化幅度之 衡量，為計算元之狀態變數在前後兩次外迭代之差值，除以所有計算 元中初始狀態變數之最大值。擴散項外迭代模組將如此重複執行上述 步驟，直至所有計算元在前後兩次外迭代後之待解變數變化幅度皆小

32 於收斂標準。 上述擴散項外迭代模組之執行流程為處理單一應用模組計算時之 處理方式，若需同時求解多個應用模組之耦合問題時，各計算元則須 計算多個應用模組之方程式集合，各方程式集合求解方式仍為擴散項 內迭代，在求解程序上則為依序進行各應用模組方程式集合之求解， 因此在所有應用模組皆完成擴散項內迭代計算後，計算元內部之各應 用模組待解變數皆已改變，此時計算元需將內部各應用模組之最新資 訊整合，此單一計算元內部應用模組間之資訊更新，本研究稱之為資 訊同步，詳細說明請見3.4.2 節，資訊更新後各應用模組再重新啟動 擴散項內迭代計算，如此往復直至收斂為止。而擴散項外迭代模組在 檢驗耦合問題時，則須檢驗所有計算元各應用模組之待解變數，若有 一計算元之任一應用模組的待解變數變化量超過外迭代收斂標準，則 須重新啟動所有計算元全面重新進行擴散項內迭代，而此時各計算元 無論是相鄰計算元之資訊或是內部各應用模組之資訊皆已更新至最 新狀態。上述步驟之執行流程請見圖3.4。

33 圖3.4 可適性計算架構擴散項外迭代模組流程圖 3.4 計算核心子系統建置 本研究所建立之核心計算模組包括空間離散模組、擴散項內迭代模 組、資訊同步模組、傳流項計算模組以及方程式計算順序設定及儲存 容器、材質參數儲存容器、邊界條件設定與儲存容器及源匯項設定與 儲存容器等四個儲存容器，如圖3.1 所示，以下將分為八個小節進行 說明。 3.4.1 擴散項內迭代模組 相較於傳統數值方法，本研究不採用整合式之控制方程式，而改採

34 用分散式之方程式集合。本研究可直接以分散之各條方程式進行計算。 每次擴散項內迭代運算，均是針對單一節點進行運算，其中待解變數 為xi，在假設周遭節點之變數數值為正確值之前提下，應用牛頓法等 最佳化方法，求解符合守恆定律之xi數值。 在 3.4 節中以「方程式集合分析模組」決定各方程式之計算順序， 守恆形式之雷諾傳輸定律必為最終計算之方程式，如式(4-31) 為守恆 方程式之離散型態(請詳閱第四章 4.2 節數值離散說明)，為符合守恆 之概念，其等號右方必為 0，因此錯誤的xi數值代入式(4-31)其等號右 方則不為 0，意即不會守恆，該數值則以誤差值 (εi)代表之。本研究 在此以最佳化方法搜尋正確之數值(x*i )，式(3-1)為目標函數，是以誤 差值之平方最小為目標。其中，變數 t i f 代表控制體積內某物理量之外 顯性質，而 t i Γ 則代表穿越控制體積之穿越率，兩者均為決策變數xi之 函數。

{

(

)

[

(

)

]

}

i x t i t t i CV t i t t i i i f f V t Z 1 0 min _=ε2₌ +Δ _{− +ωΓ}+Δ _{+ −ω Γ ×Δ} 2 ₌ (3-1) 圖3.5 為擴散項內迭代執行流程圖，圖中之變數( (n) i x )代表問題於節點 i上第 n 次迭代之計算變數數值，其以牛頓法進行修正(如式(3-2)所示)。 ) 0 ( i x 為變數xi 之初始解，首先以xi(0)進行第一次擴散項內迭代計算，意 即將 (0) i x 代入方程式集合中計算，透過守恆方程式，可以算出在節點i

35 上變數數值( (0) i x )對應之目標函數值(Zi)。接著再透過差分方式，進一 步計算該目標函數之一階微分值( ' i Z )。最後，判斷是否已經滿足收斂 標準，若尚未收斂則依據微分值修改變數數值(式(3-2a))，反覆執行此 一動作，直至完全滿足擴散項內迭代收斂標準為止，所得之變數( (n) i x ) 即為正確之數值(x*i )。 ' ) ( ) 1 ( i i n i n i Z Z x x + = − _(3-2a) C Z_i' _≤ c i Z Z ≤ 30 _ Max iter₌ (3-2b) 在擴散項內迭代之收斂條件方面，收斂條件共有三種，滿足其一即 停止擴散項內迭代演算，如式(3-2b)所示。首先是守恆誤差之一階微 分值趨近於零(Z_i' ≤C )，程式處理上則以取一極小數值作為門檻值； 其次則是守恆誤差本身小於設定之擴散項內迭代收斂標準( c i Z Z ≤ )； 最後，為避免在鄰近最佳解時，因一階微分值趨於和緩，修正幅度極 小，導致迭代次數過大之問題，本研究則訂定迭代次數之最大限制量， 本研究在此先設定為30 次，亦即擴散項內迭代次數若大於 30 則視為 收斂。在實際運算上，若計算線性問題，擴散項內迭代流程可於一次 迭代即找得最佳解；若計算非線性問題，則視問題之非線性程度而有 不同之收斂效果，在本研究訂定之三個問題中，均可於數次至十數次

36 間達到收斂。 ( )0 i x i Z ( )n i x ' ∂ ∂ = i i Z Z ( ) ( ) ' 1 i i n i n i Z Z x x + = − ( )0 i x 圖3.5 可適性計算架構擴散項內迭代流程圖 3.4.2 資訊同步模組 如整體模擬問題包含數個物理量問題，不同物理量問題間可能需使 用相同變數，而這些變數只在某些特定物理量進行計算，因此在進行 多物理量問題之耦合計算時，方程式與變數定義不完整的物理量問題， 需先由這些特定物理量問題取得各自計算時所需變數之值，才能進行

37 後續迭代計算，此取得各自計算時所需變數值之步驟，本研究稱為資 訊同步，其相關事宜由資訊同步模組負責處理。在本研究中，定義一 物理量問題所需求解之方程式集合儲存於「應用模組」中，因此底下 文章將以應用模組為單位代表一物理量問題之方程式集合。 圖3.6 為資訊同步模組之示意圖，其中包含三個不同之應用模組，分 別為應用模組A、應用模組 B 與應用模組 C，各自有本身之待解變數 及應由其它模組給定之變數，為簡化計算流程，本研究採各模組先各 自計算再交換更新變數值的方式迭代計算。如在應用模組 A 中， a1 A V 、 2 a A V 與 a3 A V 為應用模組A 之待解變數，並需提供給應用模組 B 與應 用模組C 之使用，為簡化計算流程，A、B 及 C 模組先各自計算，再 藉由資訊同步模組來將數值複製到應用模組B 與 C；相同地，應用模 組A 中需要由其他模組提供之變數如 b1 B V 、 b2 B V 、 b3 B V 與 c4 C V 等四個變數， 亦藉由資訊同步模組複製至模組A，資訊同步模組同步資訊後，接著 再重複各模組之計算，如此往復計算直至收斂為止。如以地下水流與 熱流問題耦合求解為例，地下水流問題需要溫度變數以計算水流密度 與黏滯係數之變化，而熱流問題則需水流速度進行傳流項之計算。在 開始計算時，兩問題先取得對方之初始資訊各自進行內迭代計算，當 兩問題都收斂後再更新彼此資訊，如熱流問題重新取得下水流問題之 水流流速，地下水流問題重新取得熱流問題之溫度，之後再重新各自

38 計算，亦即在各次外迭代計算前，各計算元皆會更新自己計算之所有 物理量問題所需變數資訊，再進行計算。 3 2 1_{, ,} b B b B b B V V V 5 4_, b B b B V V 1 a A V 3 2_, a A a A V V 3 2 1_{, ,} c C c C c C V V V 4 c C V 圖3.6 可適性計算架構資訊同步示意圖 3.4.3 傳流項計算模組 傳統上，應用有限元素法或有限差分法等 Eulerian 方法求解溶質 傳輸或熱流傳輸等問題時，當對流項效應較大時，純以 Eulerian 類 方法求解往往不易收斂，在此參考傳統方法結合質點追蹤 (Particle Tracking Method) 等 Lagrange 類方法，分別求解擴散項與對流項， 本研究稱此對流項之求取方法為傳流項計算模組。擴散項計算之擴散 項內迭代與外迭代，僅計算擴散項部分而將對流項忽略不計，等擴散 項計算收斂後，再以傳流項計算模組計算對流項部份，並需往復迭代

39 計算。式(3-3)用以示意說明對流項之計算，其中vi則為對應之待解變 數，Si代表計算控制體積i內之物理量總量，為待解變數vi之函數， PT i S 則是受到對流項影響後控制體積之物理量總量，而S_adv,j則是代表控制 表面 j之對流項穿越量，其中v_j則代表控制表面上的待解變數數值， 其應以控制表面兩端的數值之平均，仍為待解變數vi之函數，Vj  則代 表水流之流動速度，假設其不為待解變數vi之函數，應先給定。如以 溶質傳輸為例，Si即代表控制體積內的溶質總質量，vi則代表溶液濃 度，而S_adv,j則是溶質之對流穿越量，N 代表環繞計算節點i之控制表 面數。



= + = N j j j j adv i i PT i S v S v V S 1 , ( , ) ) (   ) ( i i f v S ₌ ) , ( ,j j j adv g v V S =  (3-3) 式(3-3)僅可計算受對流項影響後控制體積內之物理量總量，仍須透過 式(3-4)或(3-5)求得問題之待解變數。式(3-4)是直接以逆運算子將物理 量總量轉換為待解變數，因此在實作上不同問題除了需要運算子 f 外， 尚需要實作逆運算子 _f −1_{，且對於非線性之運算子 f 而言，其逆運算} 子較難推導。而式(3-5)則是以最佳化演算法求得最接近之待解變數， 對不同問題而言，可直接沿用已定義之運算子 f 作為問題中的等號限 制式，再以最佳化演算法求解。

40 ) ( 1 PT i PT i f S v ₌ − _(3-4) 2 * v Z ( ) min_* i i PT i S S − = t s, ) ( * * i i f v S = (3-5) 3.4.4 方程式計算順序設定及儲存容器 方程式集合分析模組完成分析後，會將各種應用問題耦合求解時所 需方程式集合及求解順序儲存於本儲存容器中，一種耦合問題即對應 一組方程式集合及其求解順序。在使用者決定模擬問題後，即能以對 應於該模擬問題之方程式集合進行計算，且所有計算元皆計算相同之 方程式集合。 3.4.5 空間離散模組 空間切割在數值模擬問題中，為必要的一個步驟，必須將連續之模 擬區域切割成有限區塊之計算節點或計算格網。不同數值方法均有其 不同之切割方法，以最常見的有限差分法(FDM)為例，其空間離散方 法為矩形網格。本研究採用凡諾依圖 (Voronoi Diagram)作為切割演算 法，並以Delaunay Triangulation 訂定計算節點與相鄰節點之關係 (Delaunay Triangulation 與凡諾依圖之細部說明請詳見附錄 A)，使得 空間網格切割上具備更大之彈性。凡諾依圖的最大特性，在於兩個相

41

鄰計算節點之連線必與對應之穿越邊界正交，如此將便於後續邊界穿 越量之計算。

圖 3.7(a)為凡諾依圖、 Delaunay Triangle 、相鄰節點與連結之示意 圖，其中 Ni、NA、NB與NC為計算節點，取最相鄰的三個節點依序 繪成三角形，由計算節點群所構成之三角形稱為 Delaunay Triangle， 其與凡諾依圖為對偶圖形。圖上之NiNANB三角形即為 Delaunay Triangle，該三角形各邊中垂線交點(V2 )為三角形之外心，將不同三 角形之外心所連區塊即為凡諾依圖，圖上之V1V2V3區塊即為 Ni代表 之凡諾依圖。 「可適性計算平台」將凡諾依圖演算法內建於架構中，模式開發者應 用「可適性計算平台」開發數值模式，在空間離散上僅需在模擬區域 配置計算節點(此計算節點位置即為計算元所在位置)，「可適性計算 平台」會依據輸入之計算節點與整體模擬區域邊界，決定計算節點間 的相鄰關係，以及計算節點對應之控制區塊，該控制區塊即代表該計 算節點之控制體積。 圖 3.7(b)是在模擬區域上規則地等間距配置九個計算節點，透過凡諾 依圖演算法可矩形之規則格網，圖上N0、N1至 N8代表配置之九個計 算節點，包圍計算節點之實線區塊即為該計算節點之控制體積 (Voronoi Diagram)，兩相鄰控制體積共用一個內部控制表面，並以一

42 個長虛線穿越該內部控制表面，並串連兩個計算節點的連結可與控制 表面相交稱為控制表面穿越點，舉例說明：L0代表其中一連結，串連 計算節點N0與 N1，與對應之控制表面相交於一控制表面穿越點，可 以標註其為L0,1，圖上L0至L11合計共有12 個內部相控制表面。此外， 所有控制體積涵蓋整個模擬區域，因此外圍之控制表面即為模擬區域 之邊界，並可各自定義一個邊界穿越點，圖上共有BN0、BN1至 BN11， 合計共有12 個邊界穿越點。如以計算節點 N0為例，其對應之控制體 積為四方形，因此有四個環繞之控制表面，其中兩個為內部控制表面， 分別為 L0與 L2，另外兩個則為邊界控制表面，其上之邊界穿越點分 別為BN0與BN11，邊界條件是設定於邊界穿越點上。圖中符號 V0至 V15為各控制體積(Voronoi Diagram)之角點。

43

圖3.7(a)凡諾依圖與 Delaunay Triangle 示意圖

44 表3.1 圖 3.7(b)之節點、相鄰節點、相鄰連結與相鄰邊界虛擬點對應 表 3.4.6 材質參數儲存容器 空間上各種物理量的變化常受到承載此物理量的介質特性的影響， 而此介質的特性常以各種參數代表之，本研究統稱為材質參數。以水 流、熱流及溶質傳輸為例，其分布會受到多孔介質之水文地質參數影 響，圖3.8 為土壤材質與網格分布示意圖，圖面上虛線描繪出不同材 質區塊，實線描繪出網格，以圖3.8 為例，圖上共分三種不同之土壤 材質，包含Material #1 、#2 與 #3 三區，因此不同網格都應隸屬於 對應的材質區塊。圖上顯示，節點N0與N1位於材質區塊 #2 中，而節 點N2則位於材質區塊 #3 中，藉由空間幾何關係之判定可知節點與材 質區塊之對應關係。

45 圖3.8 材質與網格分布示意圖 3.4.7 邊界條件設定與儲存容器 圖 3.9 為邊界條件與網格配置示意圖，其中粗線為設定之邊界條件， 圖上共有七個線段，包含AB、

BC

、

CD

、DE、EF、

FG

與

GA

，線 段上各自之邊界條件，分別可為 BC #1、#2、#3、#4、#5、#6 與#7。 由3.4.1 節之可適性計算架構網格結構得，從節點往模擬邊界之垂直 投影稱為「邊界穿越點」，各邊界穿越點對應之邊界條件，以BN0為 例，在空間上位於線段AB上，因此其對應之邊界條件為 BC #1。

46 圖3.9 材質與網格分布示意圖 3.4.8 源匯項設定與儲存容器 圖3.10 為源匯點與網格配置示意圖，實線則為網格之控制體積， 網格中央是其計算節點，圖上有三個源匯點，分別為 Source #1、#2 與 #3，分別位於計算節點N0、

N

4與N6之控制體積內部。

47 圖3.10 源匯點與網格配置示意圖 3.5 方程式集合分析模組建置 為了檢驗分析應用模組中，各應用模組的方程式集合所定義之變數 與方程式關係是否完整(well-defined)、及定義方程式集合之求解順序 以及同時求解多個應用模組時，橫跨應用模組間之變數傳遞關係，本 研究建置方程式集合分析模組進行分析。分析結果將儲存於核心計算 模組之「方程式計算順序設定與儲存容器」中。以下將依序說明分析 單一應用模組之方程式集合之分析概念、分析方法以及多應用模組之 方程式集合分析方法。