五、 垂向二維水流與熱流偶合數值模式模擬與驗證
6.2 建議
1. 本研究目前著重於地下水流與熱流問題之模擬，基於細胞自動機 之高度可擴充性，未來亦可以將問題延伸至地下水污染傳輸問題 或地層下陷問題等。
2. 本模式目前僅針對純水之熱流問題，因此水流熱傳導係數及水流 比熱均設定為定值，未來問題若延伸擴充至污染傳輸問題時，兩 參數之給定應依據污染物之特性與濃度進一步計算改變。
3. 本研究目前僅以地下水流模式與熱流模式之偶合模擬，因此對於 土壤壓縮性之分析，僅以線性方程式模擬之，未來若透過細胞自 動機理論建立地層下陷模式，可對土壤壓縮性進行更嚴謹之探討。
4. 本研究目前為基礎模擬工具之開發，未來可以進一步使用各項應 用型問題中，例如溫泉可利用水量管理、冷凍工法模擬分析等。
參考文獻
1. McDonald, M.G., and Harbaugh, A.W., (1988) A modular three-dimensional finite-difference ground-water flow model: U.S.
Geological Survey Techniques of Water-Resources Investigations, book 6, chap. A1, 586 p.
2. Trescott, P.C., Pinder, G.F., and Larson, S.P., (1976) Finite-difference model for aquifer simulation in two dimensions with results of numerical experiments: Techniques of Water-Resources Investigations of the United States Geological Survey, book 7, chap. C1, 116 p.
3. Larson, S.P., (1978) Direct solution algorithm for the two-dimensional ground-water flow model: U.S. Geological Survey Open-File Report 79-202, 22 p.
4. Manteuffel, T.A., Grove, D.B., and Konikow, L.F., (1983) Application of the conjugate-gradient method to ground-water models: U.S. Geological Survey Water-Resources Investigations Report 83-4009, 24 p.
5. Ozbilgin, M.M., and Dickerman, D.C., (1984) A modification of the finite-difference model for simulation of two-dimensional ground-water flow to include surface-ground water relationships: U.S.
Geological Survey Water-Resources Investigations Report 83-4251, 98 p.
6. Weeks, J.B., (1978) Digital model of ground-water flow in the Piceance Basin, Rio Blanco and Garfield Counties, Colorado: U.S.
Geological Survey Water-Resources Investigations 78-46, 108 p.
simulation of three-dimensional ground-water flow: U.S. Geological Survey Open-File Report 75-438, 32 p.
8. Trescott, P. C., and Larson, S. P., (1976) supplement to Open-File Report 75-438--Documentation of finite-difference model for simulation of three-dimensional ground-water flow: U.S. Geological Survey Open-File Report 76-59l, 21 p.
9. Kuiper, L.K., (1987) Computer program for solving ground-water flow equations by the preconditioned conjugate gradient method: U.S.
Geological Survey Water-Resources Investigations Report 87-4091, 24 p.
10. Leake, S.A., and Prudic, D.E., (1991) Documentation of a computer program to simulate aquifer-system compaction using the modular finite-difference ground-water flow model: Techniques of Water-Resources Investigations of the United States Geological Survey, book 6, chap. A2, 68 p.
11. Grove, D.B., and Stollenwerk, K.G., (1984) Computer model of one-dimensional equilibrium-controlled sorption processes: U.S.
Geological Survey Water-Resources Investigations Report 84-4059, 58 p.
12. Konikow, L.F., and Bredehoeft, J.D., (1978) Computer model of two-dimensional solute transport and dispersion in ground water:
Techniques of Water-Resources Investigations of the United States Geological Survey, book 7, chap. C2, 90 p.
13. Goode, D.J., and Konikow, L.F., (1989) Modification of a method-of-characteristics solute-transport model to incorporate decay and equilibrium-controlled sorption or ion exchange: U.S. Geological
Survey Water-Resources Investigations Report 89-4030, 65 p.
14. Dunlap, L.E., Lindgren, R.J., and Carr, J.E., (1984) Projected effects of ground-water withdrawals in the Arkansas River Valley, 1980-99, Hamilton and Kearny Counties, southwestern Kansas: U.S.
Geological Survey Water-Resources Investigations Report 84-4082, 168 p.
15. Lewis, F.M., Voss, C.I., and Rubin, Jacob, (1986) Numerical simulation of advective-dispersive multisolute transport with sorption, ion exchange and equilibrium chemistry: U.S. Geological Survey Water-Resources Investigations Report 86-4022, 165 p.
16. Glover, K. C., (1987) A dual-porosity model for simulating solute transport in oil shale: U.S. Geological Survey Water Resources Investigations Report 86-4047, 88 p.
17. Cooley, R.L., (1992) Derivation of finite-element equations and comparisons with analytical solutions; Part 2 of A MODular Finite-Element model (MODFE) for areal and axisymmetric ground-water-flow problems: U.S. Geological Survey Techniques of Water-Resources Investigations, book 6, chap. A4, 108 p.
18. Torak, L.J., (1993) Design philosophy and programming details; Part 3 of A MODular Finite-Element model (MODFE) for areal and axisymmetric ground-water-flow problems: U.S. Geological Survey Techniques of Water-Resources Investigations, book 6, chap. A5, 243 p.
19. Lin, H. C., D. R. Richards, G. T. Yeh, J. P. Cheng and N. L. Jones..
"FEMWATER: a three-dimensional finite element computer model for simulating density dependent flow and transport.", US Army
Engineer Waterways Experiment Station, Pennsylvania State University, Brighan Young University, 1996.
20. Portela, J. N. and Alencar, M. S. “Outage contour using. a Voronoi diagram”, IEEE, Piscataway NJ, ETATS-UNIS, 2004.
21. Thiessen, A.H. ”Precipitation averages for large areas”. Monthly Weather Review, 39, pp.1082-1084, 1911.
22. Sukumar, N. ”Voronoi cell finite difference method for the diffusion operator on arbitrary unstructured grids”. International Journal for Numerical Methods in Engineering, 57, pp.1-34, 2003
23. Tucker, G.E., Lancaster, S.T., Gasparini, N.M., Bras, R.L., and Rybarczyk, S.M. ”An Object-Oriented Framework for Hydrologic and Geomorphic Modeling Using Triangulated Irregular Networks ” Computers and Geosciences, 27, pp. 959-973, 2001.
24. Iske, A. and Käser, M. ”Conservative semi-Lagrangian advection on adaptive unstructured meshes”. Numerical Methods for Partial Differential Equations, 20:388–411, 2004.
25. Fuhrmann, J and Langmach, H. ”Stability and existence of solutions of time-implicit finite volume schemes for viscous nonlinear conservation laws”. APPLIED NUMERICAL MATHEMATICS , 37, pp. 201-230, 2001.
26. Georgy Voronoi, (1907) Nouvelles applications des paramètres continus à la théorie des formes quadratiques. Journal für die Reine und Angewandte Mathematik, 133:97-178.
27. Andreas Fabri, (2007) Voronoi Diagrams in CGAL, the Computational Geometry Algorithms Library. Voronoi Diagrams in Science and Engineering, 2007. ISVD '07. 4th International
Symposium on 9-11 July 2007 Page(s):8 - 14
28. Fortune, S. ”A sweepline algorithm for Voronoi diagrams, Algorithmica”. ALGORITHMICA, 2, pp.153-174, 198.
29. VANGENUCHTEN, M.T. ”A closed-form equation for predicting the hydraulic conductivity of unsaturated soils”. SOIL SCIENCE SOCIETY OF AMERICA JOURNAL, 44, pp.892-898, 1980
30. Fine, R.A. and Millero, F.J. ”The compressibility of water as a function of temperature and pressure”. Journal Chemical Physics, 59, pp. 5529-5536, 1973.
31. Mendicino, G., Senatore A., Spezzano, G. and Straface , S. ”Three-dimensional unsaturated flow modeling using cellular automata”. Water Resource Research, 42, W11419, 2006.
32. Ruhaak, W., Rath V. and Wolf A. ”3D finite volume groundwater and heat transport modeling with non-orthogonal grids, using a coordinate transformation method”. Advances in Water Resource, 31, pp.513-524, 2008.
33. Parsons, J.A. and Fonstad, M.A. ”A cellular automata model of surface water flow”. HYDROLOGICAL PROCESSES, 21, pp.2189-2195, 2007.
34. Folino, G., Mendicino, G.. and Senatore A. ”A model based on cellular automata for the parallel simulation of 3D unsaturated flow”.
PARALLEL COMPUTING, 32, pp. 357-376, 2006.
35. Coppola, E., Tomassetti, B., Mariotti, L. ”Cellular automata algorithms for drainage network extraction and rainfall data assimilation”.HYDROLOGICAL SCIENCES OURNAL-JOURNAL
36. Jacob B.,”Dynamics of Fluids in Porous Media”, Dover Publications, 1988.
37. 李士哲(2001)
金崙地區溫泉資源調查分析之研究，
國立成功大學 碩士論文。38. 許俊賢(2002)
深層岩體熱力
-水力
-力學偶合行為之研究，
國立成 功大學碩士論文。39. 張育弘(2004)
台灣核廢料地下處置場之熱傳及熱應變分析
，國立 中央大學。40. 宋鴻鈞(2003)
三微細胞自動機之研究
-以海洋油污染擴散為例
， 國立中山大學碩士論文。41. 黃宏明(2005) Voronoi Diagram
互動教學程式之設計與研究
，私 立中華大學碩士論文。附錄 A 總水頭計算公式推導
式 A-9 為積分上之恆等式，意即 至 間的定積分值，可以拆成
P
uP
lW P
c∆ x
∆ y
∆ z
圖附錄 A.1-1 控制體積垂向受力示意圖
附錄 B 土壤壓密係數與傳統方法之關聯推導
式 B-5 為本模式描述控制體積內蓄水質量之方程式，其中因此若
控制體積內水體積之時間變化率，其中 為比出水率，代表拘限含水 層因應壓力水頭變化，土體與水體因壓縮或膨脹所釋放出的水量，式 B-11 為比出水率、土體壓縮性與水流壓縮性三者之關係。此外，式 B-12 定義總水頭為壓力水頭與位置水頭之和。由於傳統方式是以水 體積守恆建立概念模式，因此若欲轉換為水質量守恆，則需再乘上密 度。此外傳統常見之控制方程式均為偏微分形式之方程式，在推導過 程中已使其代表單位體積之變化率，因此若欲與本研究對應，則應乘 上控制體積之大小，因此透過傳統方法所描述之控制體積內蓄水質量 之時間變化率可以描寫為式 B-13。觀察式 B-10 與 B-13，兩者公式極 為類似。因此顯示本模式所採用之方程式組合，在經過推導與簡化 後，可與傳統拘限地下水流模式所採用之方法一致。
ss
β α n
ss = + ...(B-11)
z p
h= + ...(B-12)
( ) ( ) ( )
t n p t Vol
z n p
t Vol St
∂ + ∂
×
∂ = + + ∂
⋅
∂ =
∂ ρ0 α β ρ0 α β ...(B-13)
附錄 C Voronoi Diagram 空間切割法
空間切割在數值模擬問題中，為必要的一個步驟，不同數值方 法則有不同之切割方法，以最常見的有限差分法(FDM)為例，其必須 侷限於矩形網格，因此對於不規則形狀之掌握較差。本研究在此採用 Voronoi Diagram 作為切割演算法。
附錄 C.1 徐昇氏網格(Voronoi Diagram)簡介
Voronoi Diagram 是透過散佈於平面上一群平面點位資料(sites or points)為基礎所切割出來的圖形，其圖形如圖附錄 C.1-1 所示。Voronoi Diagram 具備下列特性：
1. 各點(sites or points)所必對應一個控制區塊(cell)，區塊內之任何一 點與該節點之距離，必定小於此點與其他節點之距離。
2. 控制區塊形狀必為凸包(convex hulls)。
3. 兩節點之連線必定正交於其控制區塊(cells)之邊，此邊稱為控制區 塊邊(edge)
4. 控制區塊邊(edge)之交接點稱為角點(vertex)，令任意角點為圓心，
必可找到一圓通過相鄰之節點(如圖附錄 C.1-1 所示)，該圓稱為空 圓(empty circle)。
5. 將任意空圓上之所有節點連線，可以繪製成三角形，其稱為 Delaunay Triangles，其與 Voronoi diagram 為對偶圖形。圖附錄 C.1-1 上，節點 0、1 與 4 之連線所形成之三角形即為 Delaunay Triangles。
圖附錄 C.1-1 Voronoi Diagram 與其相關圖形
數學定義：若給定一組點位資料(N ={N1,N2,N3LLNn})，其中 為第 m 個節點，n 為空間中節點數量，其數值應為
Nm
∞
≤
≤ n
2 。其控制
區塊(Voronoi Diagram)可定義如下：
( )
N{
D D N D N for j i}
V i = : − i ≤ − j ∀ ≠ ...(C.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
附錄 C.2 CGAL 函式庫簡介
CGAL 全名為 Computational Geometry Algorithm Library，是計 算空間幾何上極具可信度的函式庫，其中包含了許多圖形的資料結 構、建構圖形的演算法、格網產生與處理與圖形分析等功能。為 Geometry Factory 所開發的開放原始碼專案，其目的為 C++環境，提 供有效率並且可信賴的空間計算演算法。目前已經被廣泛應用於電腦 繪圖、科學視算、設計模擬、格網生成、數值方法與微生物模擬等各 項領域上。
然而 CGAL 函式庫提供之 Voronoi Diagram 相關功能，僅針對無 限空間之 Voronoi Diagram 空間切割，亦即外圍節點之控制區塊 (Voronoi Cells)並不為封閉空間。然而本研究是應用於數值模擬上，
網格必須是有限個封閉空間，因此不完全符合需求。
因此本研究是以 CGAL 函式庫作為基礎函式庫，透過物件封裝 技巧，進一步開發適合本研究之封閉空間 Voronoi Diagram，使外圍 之控制區塊與模擬邊界進行切割計算，重新輸出適合本研究數值模擬 之空間切割結果。
附錄 C.3 封閉空間之 Voronoi Diagram 函式庫開發
由於 CGAL 之基本功能僅提供非封閉空間之切割，對於外圍控
封閉空間之切割結果，並與模擬區塊邊界框比較套疊，重新製作為封 閉空間之切割結果。
圖附錄 C.3-1 不同邊界框形成之 Voronoi Diagram
圖附錄 C.3-1 為相同之點位資料，套疊不同之模擬區塊邊界框，
會套疊出不同之控制區塊形狀。例如邊界 1 與邊界 3 對控制區塊 1 而 言，便會形成不同之形狀，對於邊界 3 而言，必須刪除凸出邊界的部 分；對於邊界 1 而言，則由於控制區塊本身即為封閉區塊，所以無須 另行處理。另外，不同邊界對於非封閉型之控制區塊而言，則需以邊 界本身與非封閉型之控制區塊，以聯集方式形成封閉之控制區塊。本 研究之處理邏輯如下所示：
1. 非封閉形式之控制區塊：如圖附錄 C.3-1 之控制區塊 0、2、3、4、
5、6。此類狀況則以邊界框與該控制區塊做聯集判斷，可形成一 個封閉多邊形。
2. 封閉形式之控制區塊：如圖附錄 C.3-1 之控制區塊 1、7。此類狀 況可依據邊界框是否穿越該控制區塊，再做分類。若無穿越，則 該控制區塊為正確解；若有穿越，則該控制區塊仍應與邊界框做 聯集判斷。