行政院國家科學委員會專題研究計畫完整成果報告
台灣環島海岸水位預報系統之建立-總計畫暨子計畫：
台灣環島暴潮預報模式及數值網格產生法之研究(III)
A forecasting system of water elevations around Taiwan：
Studies on Numerical Modeling of Storm Surges and GridGeneration in Coastal Waters around Taiwan
(III)
計畫類別：整合型計畫 計畫編號：NSC92－2625－Z－002－014 執行期間：92 年 08 月 01 日至 93 年 07 月 31 日 總計畫主持人：蔡 丁 貴 計畫主持人：蔡 丁 貴 計畫參與人員：王 鄭 翰、鄭 正 德 處理方式：
可立即對外提供參考 （請打ˇ）□一年後可對外提供參考 □兩年後可對外提供參考 （必要時，本會得展延發表時限） 執行單位：國立台灣大學土木工程學研究所 中華民國 九十四 年 四 月 一 日摘 要
一、總計畫：
「台灣環島海岸水位預報系統之建立」為三年期整合型計畫，整 合之構想是以台灣環島（包括澎湖、金門、馬祖等各離島）海岸水位 之預報系統建立為主題，分為四個主要的部分，包含台灣環島天文潮 汐預報模式、台灣環島暴潮預報模式、數值網格產生法研究以及台灣 環島海岸水位預報資訊系統建立。其主要目的在結合海象及海岸工程 之專精研究人員，促成學術界專業研究人員之積極參與，希能儘早完 成國內海岸溢淹預警作業系統，以早日發揮海岸溢淹災害防治功能。 本計畫為「海岸環島海岸水位預報系統之建立」整合型計畫之總 計畫，主要負責海岸水位模擬示範區之選定(在環島水位計算完成 後)，共用資料之蒐集與格式制定建檔，每年度除擔任各計畫協調工 作與資料相互引用事宜，以及現場調查支援問題，及現場看察工作 外，並召開 3 至 4 次期中檢討會以廣集交換計畫執行心得與需求，了 解各計畫之進度，以期確切達成計畫目標與進度。 本計畫的研究成果包括：(1)潮汐半日型主要分潮的調和常數及 潮流橢圓計算；(2)邊界元素數值計算方法中關於奇異性問題之解 析；(3)運用對數學奇異性與幾何奇異性的研究成果，探討其對邊界 符合網格系統產生的影響。同時，根據複變函數的特性，研究出其進 行複變轉換時的角度限制式；(4)小區域(單連通)以及大區域(複連通) 網格之建立；(5)完成整體地理資訊系統規劃以及初步系統建置；(6) 建立於地理資訊系統可供查詢之現場觀測資料分析結果。一、子計畫：
本計畫為三年期整合計畫「台灣環島海岸水位預報系統之建立」 之子計畫之一，其研究目的在發展台灣本土的暴潮及溢淹預測模式， 期能降低每年颱風侵犯台灣時淹水所造成生命財產上的損失。同時， 配合總計畫，提供暴潮數值計算以及建立邊界符合數值網格予另兩個 子計畫計算所需。本研究除了引進美國聯邦緊急事故處理局發展出的 暴潮預測模式---FEMA 暴潮模式之外，並藉由探討邊界元素法在奇 異性問題上的理論基礎，建立並改善原網格系統為邊界符合座標系 統；另外，也進一步改善 FEMA 模式數值計算方法以增加計算效率 及穩定性。 模式發展初期以數值正交網格之建立、應用改良為主，中期將以 整合正交網格與新模式以及驗證為主，而在末期時並考慮天文潮位之 變化，驗證曲線座標系統之新模式。至於風場模式之修正則將於研究 後期暫以氣象局預報的大區域風場代替移動的圓形風場模式，同時考 量背景氣壓來模擬，待他日氣象界發展出較佳的風場模式時再加以改 良應用。 本模式之發展預期將可提供沿岸暴潮，沿岸河川水位，暴潮溢淹 範圍等預測之所需，乃至進一步建立台灣溢淹區洪災預警及疏散管理 系統。其並可做為本國沿海海岸後退線劃定、洪水保險等之依據。 關鍵詞：台灣、邊界符合座標系統、水位預報、潮汐、暴潮、地理資 訊系統ABSTRACT
1.Main program:
The object of this proposal is to form an integrated research group to make efforts in establishing a forecasting system of water elevations around Taiwan. It includes four components, namely, study of tidal prediction model on seas surrounding Taiwan, numerical modeling of studies on storm surge in coastal waters around Taiwan, numerical modeling of studies on grid generation technology in coastal waters around Taiwan, and studies on an information system of coastal water elevation around Taiwan. This integrated program is to be completed in period of three years. The team of several researchers will devote themselves to the works of numerical modeling, database establishment, and in-site investigations of coastal water elevations around Taiwan. The final goal is the establishment of a predictive and warning system for coastal floods around Taiwan.
This integrated project will make a selection of main regions suffered from coastal flooding after water elevations around Taiwan are obtained. Common data format, file system, and in-site field investigations is arranged throughout this integrated research. Three to four progress meetings was conducted during the passed research period to exchange research results and to share experiences in order to achieve the goal of each individual sub-project.
The study results of this year include: (1) the harmonic constants associate with the current ellipse of the principle constituents of semidiurnal tide are solved on coast and seas around Taiwan, (2) solving singularity problem of BEM and the establishment of the formulation for extending its application to multi-connected regions, (3) the computing of astronomic tide and storm surge around Taiwan, (4) the establishing of small-area (simple connected) boundary-fitted orthogonal grid system and large-area (multiply-connected) boundary-fitted orthogonal grid system, (5) completing the first-stage arrangement of Geographic Information System, and (6) in-site data associated with the data analysis can be found in GIS.
2.Sub-program:
system of water elevations around Taiwan”.The main purpose of this study is to develop a new storm surge model by considering the local boundary effects. And the success of this model is expected to reduce losses of life and property caused by storm surges. The other purpose of this study is to establish and offer the boundary-fitted orthogonal grid systems for another two sub-studies of this program.
The study first introduces the FEMA storm surge model developed by the Federal Emergency Management Agency, U. S. A., and further improves the model theoretically and numerically by replacing the original grid system with the boundary-fitted grid system. Meanwhile, newer viewpoints about the singularity problems in BIEM are discussed.
Finally, a new numerical scheme is established in this storm surge model to improve the computational efficiency and stability.
In the first stage of the development of this new storm surge model, the principle work is to improve the boundary-fitted coordinate system. In the second stage, the main work is to combine the FEMA model with the orthogonal coordinate grid system. In the last stage, variation of the astronomic tide will be included and the verifications and validations of the new model with curvilinear coordinate system will be made.
The development of this model will meet the needs of the coastal storm surge prediction, coastal river water level prediction, and coastal flooding prediction. Moreover, this study can be taken as the base of the inundation prediction and evacuation management system for coastal flooding areas of Taiwan. Also, the development of this model can provide as guidelines of the setback line and flooding insurance.
目 錄
摘要 ...I Abstract ...III 目錄 ... V 表目錄 ...VI 圖目錄 ... VII 第一章 前言 ...1 1.1 研究目的...1 1.2 研究方法...2 1.3 研究步驟及工作重點 ...6 第二章 理論分析 ...10 2.1 文獻回顧...11 2.2 天文潮理論 ...14 2.3 邊界符合座標系統理論 ...16 2.4 建立邊界符合座標系統所遭遇的問題與對策 ...22 2.5 暴潮理論...34 第三章 數值結果與分析 ...45 3.1 保角網格系統範例演算 ...45 3.2 天文潮結果與分析 ...58 3.3 邊界符合座標之建立 ...72 3.4 暴潮偏差數值結果與分析 ...78 第四章 結論與建議 ...97 6.1 結論 ...97 6.2 建議 ...98 參考文獻 ...100表目錄
表 3-1-1-1 Cosine 凹面範例之保角網格品質評估 ...33 表 3-1-2-1 半圓形甜甜圈範例之保角網格品質評估...33 表 3-1-3-1 近梯形範例之保角網格品質評估...33 表 3-1-4-1 四個半圓區域範例之保角網格品質評估...33 表 3-2-1 1996~2002 淡水河口歷年來暴潮偏差表 ...33 表 3-2-2 天文潮各分潮參數...34 表 3-2-3 1999~2002 歷年天文潮水位預報誤差表...35 表 3-2-4 台灣環島潮位站天文潮計算誤差表 ...47表 3-3-1 替代弧形路徑逼近(Contour Approach Method)法計算結果 ...49
圖目錄
圖 1-3-1 台灣環島水位預報系統分工架構圖...5 圖 2-2-1 天文潮模式建立流程圖 ...15 圖 2-3-1 座標轉換示意圖 ...16 圖 2-3-2 複變映射示意圖 ...17 圖 2-4-1 保角網格產生流程圖 ...17 圖 2-4-2 自 z 到 w 平面之複變轉換示意圖 ...17 圖 2-4-3 李曼平面示意圖 ...17 圖 2-4-4 內角平滑化示意圖 ...17 圖 2-4-5 蜂窩狀區域範例，A、D 再各細分成兩個內角...17 圖 2-4-6 封閉多邊形角度限制式示意圖 ...17 圖 2-4-7 以蘭嶼導範例說明角度限制式之使用 ...17 圖 2-4-8 根據角度限制式選擇複變轉換順序。經過 255 次的複變轉換 之後，原區域轉換為一超矩形區域...17 圖 3-1-1-1 a.沒有用平滑係數所做出的網格 b.平滑係數為 0.01.c.平滑 係數為 0.11 以及 d.平滑係數為 0.11.但 c.中最靠近邊界的網 格線與邊界網格距離設定為定值。(Akcelik, 2001) ...37 圖 3-1-1-2 四邊節點數相同之均勻保角網格系統(41x41) ...37 圖 3-1-1-3 四邊節點數相同但網格數加密之保角網格系統 (101x101)，可發現左上與右上角隅之網格線分佈較平滑 ...37 圖 3-1-1-4 不均勻保角網格分佈，靠近左右邊界的網格線較密 ...37 圖 3-1-1-5 邊界點點數降低之保角網格系統 ...37 圖 3-1-2-1 半圓形甜甜圈範例，邊界節點等距 ...37 圖 3-1-2-2 半圓形甜甜圈範例，邊界節點等距時，產生的網格交錯(Luis, 1996)...37圖 3-1-2-3 本研究產生之 41x41 的保角網格系統 ...37
圖 3-1-3-1 360 個邊界點及 41x41 網格線之近梯形保角網格...37
圖 3-1-3-2 a.沒有使用平滑係數，b,c 使用平滑係數 0.01，a,c 使用滑動 邊界法(Sliding boundary)(Akcelik et al., 2001) ...37
圖 3-1-4-1 四個半圓區域，邊界節點等距 ...37
圖 3-1-4-2 四邊皆為 Dirichlet conditions. (B1)RL (Rinskin & Leal, 1983) 系統所產生的網格系統.(B2)RL 系統加上使用平滑函數控 制.(B3)RL 系統加上使用微調係數0.01. (B4)兼用 RL 系 統,微調係數0.002以及控制係數c 0.5。 ...37 圖 3-1-4-3 超矩形平面上的邊界節點分佈 ...37 圖 3-1-4-4 本研究產生之四個半圓區域保角網格系統 ...37 圖 3-1-4-5 四個半圓區域 A 點之放大圖。(1)四個角附近點數加密。 (2)四個物理奇異點附近沒有加密...37 圖 3-2-1 天文潮模式驗證 (一) ...37 圖 3-2-2 天文潮模式驗證 (二) ...37 圖 3-2-3 天文潮模式驗證 (三) ...38 圖 3-2-4 天文潮模式驗證 (四) ...38 圖 3-2-5 天文潮模式驗證 (五) ...39 圖 3-2-6 天文潮模式驗證 (六) ...39 圖 3-2-7 潮位站設置圖 (資料來源：中央氣象局) ...40 圖 3-2-8 竹圍天文潮模式驗證 ...41 圖 3-2-9 新竹天文潮模式驗證 ...41 圖 3-2-10 台中港天文潮模式驗證 ...42 圖 3-2-11 箔子寮天文潮模式驗證 ...42 圖 3-2-12 東石天文潮模式驗證 ...43 圖 3-2-13 高雄天文潮模式驗證 ...43 圖 3-2-14 後壁湖天文潮模式驗證 ...44 圖 3-2-15 成功天文潮模式驗證 ...44
圖 3-2-16 花蓮天文潮模式驗證 ...45 圖 3-2-17 蘇澳天文潮模式驗證 ...45 圖 3-2-18 梗枋天文潮模式驗證 ...46 圖 3-2-19 基隆天文潮模式驗證 ...46 圖 3-3-1 Kisu(1988)問題示意圖 ...48 圖 3-3-2 Kisu(1988)和陳(1994)計算結果與解析解之比較...48 圖 3-3-3 東石、布袋一帶海域之邊界符合正交座標網格 ...51 圖 3-3-4 台灣環島邊界符合正交網格系統之建立(一) ...52 圖 3-3-5 台灣環島邊界符合正交網格系統之建立(二) ...53 圖 3-4-1 網格編號(30,1)所得數值結果 ...55 圖 3-4-2 網格編號(40,2)所得數值結果 ...55 圖 3-4-3 象神颱風淡水暴潮偏差計算結果 ...58 圖 3-4-4 象神颱風基隆暴潮偏差計算結果 ...58 圖 3-4-5 象神颱風高雄暴潮偏差計算結果 ...59 圖 3-4-6 象神颱風花蓮暴潮偏差計算結果 ...59 圖 3-4-7 賀伯颱風淡水暴潮偏差計算結果 ...60 圖 3-4-8 賀伯颱風基隆暴潮偏差計算結果 ...60 圖 3-4-9 賀伯颱風高雄暴潮偏差計算結果 ...61 圖 3-4-10 賀伯颱風東石暴潮偏差計算結果 ...61 圖 3-4-11 象神颱風淡水環島暴潮偏差計算結果 ...62 圖 3-4-12 象神颱風竹圍環島暴潮偏差計算結果 ...62 圖 3-4-13 象神颱風新竹環島暴潮偏差計算結果 ...63 圖 3-4-14 象神颱風箔子寮環島暴潮偏差結果 ...63 圖 3-4-15 象神颱風高雄環島暴潮偏差計算結果 ...64 圖 3-4-16 象神颱風後壁湖環島暴潮偏差計算結果 ...64 圖 3-4-17 象神颱風花蓮環島暴潮偏差計算結果 ...65 圖 3-4-18 象神颱風花蓮環島暴潮偏差計算結果 ...65 圖 3-4-19 象神颱風梗枋環島暴潮偏差計算結果 ...66 圖 3-4-20 象神颱風基隆環島暴潮偏差計算結果 ...66
圖 3-4-21 賀伯颱風淡水環島暴潮偏差計算結果 ...67 圖 3-4-22 賀伯颱風竹圍環島暴潮偏差計算結果 ...67 圖 3-4-23 賀伯颱風新竹環島暴潮偏差計算結果 ...68 圖 3-4-24 賀伯颱風台中港環島暴潮偏差計算結果 ...68 圖 3-4-25 賀伯颱風東石環島暴潮偏差計算結果 ...69 圖 3-4-26 賀伯颱風高雄環島暴潮偏差計算結果 ... 69 圖 3-4-27 賀伯颱風成功環島暴潮偏差計算結果 ... 70 圖 3-4-28 賀伯颱風蘇澳環島暴潮偏差計算結果 ... 70 圖 3-4-29 賀伯颱風梗枋環島暴潮偏差計算結果 ...71 圖 3-4-30 賀伯颱風基隆環島暴潮偏差計算結果 ... 71
第一章
導論
1.1 研究目的 海岸水位之預報準確性不但關係我國已經高強度開發之海岸地 區人民生命安全與財產保障，對於碩果僅存的海岸保護與保育地區之 經營管理策略擬定，具有關鍵成敗的重要性。同時，海岸河口水位更 關係重要河川之洪水演算模式之預報能力，其準確性即為河川洪水演 算下游邊界條件的重要依據，影響河川防洪設計之標準及洪水預報之 操作，其重要性非常明顯。 海岸水位之預報包括天文潮位、暴潮位、波浪抬升水位。其中天 文潮位為地球表面水體受太陽、月亮及其他星體之重力作用在地球自 轉的情形下產生近週期性的水位變化，而暴潮位為特殊氣象條件如颱 風的中心氣壓差及風場剪應力造成水面的變化，不但水位受低氣壓作 用而抬昇，風場在地形影響下更擴大水位抬升。風場產生的波浪在近 岸又因非線性及碎波效應而造成波浪而堆高水位。然波浪之效應相對 較小，在本整合計畫內暫不列入考量。 本計劃之目的在建立台灣環島水位預報模式，以預報天文潮位及 潮流、暴潮，數值網格產生法研究以及環島海岸水位預報資訊系統建 立四部分為工作重點。其中，以天文潮位之即時水深做為預報暴潮水 位及流場之基礎，並建立快速而準確之數值網格，作為提昇水理計算 準確度的努力重點；台灣環島海岸水位預報資訊系統建立之目的則為 現場資料蒐集、整理、分析與呈現，以 GIS 為基礎，呈現即時之環島 海岸水位，並提供數據作為水理數值模式驗證之依據。本整合計畫所 稱台灣環島除台灣本島外，尚包括澎湖、金門及馬祖等離島，以為全 國海岸水位完整考量之設計。 本研究期望藉由台灣海岸水位預報系統之建立，提供海岸防災、 海岸保護之依據，藉以建立台灣安全、永續以及美麗的海岸環境。1.2 研究方法 天文潮位之常用最小週期成分為半日潮，其週期約為 12 小時 30 分，暴潮位發生之時間尺度在觀測地點受颱風影響之時間長度約在 24 小時，而波浪因受風吹而產生，最長風浪週期約在 25 秒以下，颱 風進行中，沿著路徑有時產生更長之湧浪，其週期約在 1～3 分鐘。 海岸水位之變動量當中，我國環島東部天文潮位振幅較西部為小，西 部之潮位振幅分佈由南北向中部逐漸增加，在台中港附近產生最大之 潮差，約有 4.5 公尺。 暴潮潮位一般均由總潮水位扣除天文潮位（以無暴潮條件之預報 天文潮位為基礎）而得，受地形之影響甚大，岬灣凹處，暴潮往往較 為 明 顯 ， 海 岸 結 構 物 如 防 波 堤 之 上 風 處 （ upwind ） 又 較 下 風 處 （downwind）之暴潮位為大。 海岸水位預報系統之建立，其另一關鍵要素在於所使用之數值模 式必須經過驗證（verification）及檢驗（validation）之程序，使預報 之結果，能有信心使用，作為洪災、淹水等之預報與預警之設計、洪 水保險等海岸保護與保育策略措施之依據。海岸水位現場水位站網資 料之取得，分析與整理，是一項基本而不可或缺的工作，我國之海岸 水位站之建立與管理，目前分在各業務單位自行處理，缺乏統一之標 準，參考點之高程可能不一致，亟需進一步加以整理，使之一致而成 為有用之資料，達到掌握海岸水位的正確資訊的目的。而此類資料對 於海岸水位預報系統，更可提供數值模式最重要的檢驗數據來源。 海岸水位預報數值模式之建立，不論是以三維或二維之空間來描 述問題，對於現場實際問題的掌握，模式必須具備的重要功能之一， 就是要能快速而準確產生與現場地形契合（boundary-fitted）的網格 系統，才能進行實際而準確的水位及流場計算，因此，數值模式的前 置處理，也就是計算網格產生模式的搭配也是一個重要而有效的工 具，其效率的好壞也會決定水位數值模式的預報能力的高低。 本整合型計畫包含三個子計畫：(1)台灣環島天文潮汐預報模式
之研究，主持人為莊文傑博士1、(2)台灣環島暴潮預報模式及數值網 格產生法之研究，主持人為蔡丁貴教授2、以及(3)台灣環島海岸水位 預報資訊系統之建立，主持人為蘇青和博士3。 各計畫的研究重點及方法摡述如下： 子計畫一以預報天文潮位及潮流為重點。考量臺灣環島海域的水 深、海岸地形及地球自轉效應，於線性及忽略底床與水體表面摩擦與 擴散效應下，依據淺水長波系統方程式建立含括地轉效應與水深影響 因子的二維潮波水動力計算模式(Tsay, 1991；Juang, 2000)，並進一步 使用有限元素法，配合無反射的局部輻射邊界條件，利用線性等參數 三角形元素，在開放海域邊界上，給定單一分潮潮波的入射角，依據 所建立的潮波水動力模式，計算台灣環島沿岸及近岸海域任意位置上 潮汐半日型各主要分潮的振幅與相位調和常數，同時，並可計算求出 潮流橢圓(current ellipse)。 子計畫二為台灣環島暴潮預報模式及數值網格產生法研究兩部 分，以天文潮位之即時水位（子計畫一）做為預報暴潮水位及流場之 基礎，並建立快速而準確之數值網格，作為提昇水理計算準確度的努 力重點。在環島暴潮預報模式的建立方面，子計畫二引入美國美國聯 邦緊急事故處理局(Federal Emergency Management Agency, FEMA)所 發展出的暴潮及溢淹模式(FEMA , 1983, 1988)，並利用自行發展的邊 界符合曲線正交座標系統(Tsay and Hsu, 1997)改善 FEMA 模式邊界不 符合之缺點，來降低數值計算受地形之影響。同時在矩形計算區域內 修正 FEMA 模式，嘗試用預測－修正法(predictor-corrector method)(林 聰悟、林佳慧，1997)計算以增加準度及效度。 數值網格產生法方面，乃先採用複變映射理論將不規則區域轉換 成標準區域，再配合邊界積分元素法解拉普拉斯方程式 (Laplace equation)，而將標準區域轉換至矩形區域，並在矩形區域中形成網 1 交通部運輸研究所港灣技術研究中心研究員 2 國立台灣大學土木工程研究所
格。其中使用邊界元素法時，角落 (Corner)的多重方向導數(Multiple normal derivatives)，需另外處理，方能使轉換順利進行。最後，再利 用邊界積分法，與柯西里曼條件( Cauchy-Riemann condition)，將矩形 區域中之網格反轉換回標準區域，然後利用複變映射理論反轉換回原 幾何形狀，如此，原幾何形狀之網格便可建立完成。 子計畫三為台灣環島海岸水位預報資訊系統之建立，其目的為現 場資料蒐集、整理、分析與呈現，以 GIS 為基礎，呈現即時之環島海 岸水位，並提供數據作為水理數值模式驗證之依據。在資料庫的資料 品管方面，則配合實測資料繪圖，以統計分析，能譜分析，調和分析 等一連串的分析加以篩選和補遺。 本年度的計畫成果主要在於邊界符合網格系統的改善，此研究基 於 Tsay and Hsu 於 1997 年所完成的邊界符合網格系統加以發展。主 要應用前面所述的兩種奇異性於邊界符合網格系統之產生。在 Tsay and Hsu 的研究中，雖然已經使用複變轉換的技術來克服流線與勢能 線不垂交的問題，可惜並沒有深入討論其中關於物理奇異性的問題。 而這些問題，必須在具退化邊界的區域內產生網格時才會突顯出問 題。同時，在用到邊界元素法產生網格時，亦即，在運用解拉普拉斯 方程式(Laplace Equation)產生保角網格系統(Conformal grids)時，由於 前人並不了解(至少在文獻中沒有發現)物理奇異性問題，因此認定保 角網格的產生極為不可能，進而加入微調項於拉普拉斯方程式中，犧 牲保角性而維持正交性(Thompson et. al., 1977) (保角網格系統可以 使得座標轉換後，拉普拉斯項維持拉普拉斯型態，避免控制方程式複 雜化，增加運算效率)。其實，只要在運用邊界元素法產生網格前， 利用之前所述的複變轉換將物理奇異性去除，產生保角網格並非難 事。然而，即便複變轉換有此優點，但在運用時，卻要注意 branch-cut 問題，其有轉換時計算域壓縮或是拉伸的問題，如果轉換後計算域為 拉伸，超越 branch-cut 的部分便會產生多值(Multi-value)的問題，使 得網格扭曲交錯，這也是文獻中所記載完全保角的網格不易產生的原
因。此研究中找出複變轉換的限制式，作為實際運用時，安排物理奇 異點先後轉換順序之依據，也就是說，將不符合限制式的轉換點放在 較後面的順序來轉換，符合限制式的角度先轉換的話，常常原本不符 合限制式的奇異點也可以成功轉換。在研究論文中，於蘭嶼島網格例 子中，可連續轉換 255 次將 255 個物理奇異點去除，而保證為保角轉 換，便足以說明限制式的可行性，由於轉換所花費的 CPU 時間很少， 因此即使轉緩 255 次速度仍然很快。另外，由於在邊界元素法基礎研 究中成功去除掉數值邊界層，因此利用此產生邊界符合網格系統，將 可產生極為靠近邊界(10-9)的網格系統，對於計算邊界層(Boundary layer)，將是有力的工具(Wang, 2005)。
1.3 研究步驟及工作重點 各研究計劃子題以及其研究綱要與分工合作架構茲說明如下：
台灣環島海岸水位預報系統之建立
蔡丁貴
國立台灣大學土木工程學系教授
台灣環島天文潮 預報模式 莊文傑 港灣技術研究中 心研究員 台灣環島暴潮預 報模式 蔡丁貴 台灣大學土木工 程學系教授 台灣環島海岸水 位預報系統 蘇青和 港灣技術研究中 心研究員 數值網格產生模式 蔡丁貴 台灣大學土木工程 學系教授 GIS資訊系統 Commercial Package 圖 1-3-1 台灣環島水位預報系統分工架構圖 由上圖可以清楚看出各子計劃之間的相關性。以下就各子計劃的 工作內容以及子計劃之間的關係作進一步的說明：一、 台灣環島天文潮數值預報模式：
近一兩年來，台灣海峽及環島天文潮位之物理特性與數值模擬已 具有可觀的成績（Lin et al, 2000a, 2000b; Juang, 2000; Tsay, 2000）。當 某一地點之天文潮型由半日潮主導時，不但潮差之預測可達 20 公分 以內，滿潮時之預測亦可達到約 10 分鐘之誤差內，以此為台灣環島 之海岸水位預報系統之天文潮位預報基礎，應有在較短時間達成建立 整體預報系統的實效。 天文潮位數值模式預報之重點工作包括： 1.半日潮之流場模式建立、驗證及檢證； 2.全日潮之水位模式建立、驗證及檢證； 3.全日潮之流場模式建立、驗證及檢證； 4.合成潮之預報系統之建立、改善、驗證及檢證。 二、 暴潮水位數值模式 暴潮水位與流場之預報一般均利用平面二維之質量守恆關係式 及運動方程式進行數值模式的建立。最為成功的數值模式應屬美國聯 邦緊急事故處理局(FEMA)所發展建立之直角座標網格有限差分，內 置細網格的數值模式，並利用該模式進行沿岸淹水機率評估，做為洪 水保險費率計算的依據。此模式已成功應用於嘉義白水湖海岸之暴潮 計算(郭金棟等，1998)。 暴潮預測技術又分為兩個子計畫，其一為颱風壓力場及風場之預 測；其二為受壓力場及風場而運動之暴潮位水理預測模式。本子計畫 之重點在此一階段將先以水理預測模式之考量為優先，而颱風壓力場 及風場之預測則利用該領域現有之理論為依據。 本計畫將利用 FEMA 模式，引入邊界符合正交座標系統改善其 模式中無法計算任意形狀的不足之處，並以預測之天文潮位水深（子
計畫一）進行颱風暴潮位計算。研究成果需以現場潮汐資料驗證（子 計畫三）。各年度計畫工作之分配情形為： 1.利用 FEMA 暴潮溢淹模式，並以發展之數值正交網格（非卡式 座標）進行大網格計算環島暴潮分佈。 2.以大網格計算之網格點為邊界條件，內置（embed）細網格對 局部海岸地區進行暴潮計算，地點資料選擇順序依其重要性逐 一建檔（子計畫三）。 3.考慮天文潮汐之變化，考慮我國沿岸特有之地貌（如鹽田、魚 塭等），驗證曲線座標系統之新模式，完成建立新的暴潮預報模 式。 三、 數值網格產生技術之改善與應用 數值預報模式之預報能力由其準確性與計算速率決定，數值模式 所依據之數學模式是對所欲描述之水理現象經過假設與簡化而得，有 其數學上的必要性，但對其結果之應用則有其限制，此可由現場現象 觀測之檢驗與計算結果比較而知其適當性。除此之外，數值模式的另 一關鍵問題，則為計算網格前置準備的簡易性，或網格在計算方法上 所產生誤差的降低，均是預報模式能力提昇的重點。 本子計畫重點為改善現有以邊界元素法建立之數值網格產生方 法，提昇此方法之準確性及應用範圍。其主要工作項目分列如下： 1.解決邊界元素法在基礎點（Base point）靠近邊界時之不連續問 題。 2.擴充應用範圍至複連通區域（Multiply-connected domains）。 3.改善座標轉換偏導數值計算方法的準確性。網格點分佈之選擇 策略。 4.邊界符合座標系統之有限差分法數值網格產生。(子計劃二)
5.邊界符合區域之有限元素法計算網格產生。(子計劃一) 四、 海岸水位及流場站網之建立 海岸水位預報系統之成功與否，在於預報數值模式是否能準確並 快速產生答案，因此，不論是在系統發展階段及應用階段，現場資料 之蒐集、整理、分析與呈現均有其重要性，不但必須有歷時之變化， 而且也需有空間上環島不同位置不同的反應。此一工作學理上雖然不 能預期重大突破，但此基礎工作，不容有錯誤，因此資料如何蒐集、 整理分析之方法探討、及彙整後資料之呈現、保管與流通都是日後海 岸管理與研究工作不可或缺的資訊。 本子計畫之工作項目可以列如以下各項： 1.利用 GIS 蒐集水位站網各水位站及流場觀測站之詳細資料。 2.檢驗校正現場水位及流場資料。 3.進行調合分析，建立各水位站各分潮之振幅、相位及流速。 4.預報天文合成潮之水位及流速。 5.由分析之潮位資料，建立環島潮位、潮時及流場之分佈。 6.暴潮發生時段，將預報之合成天文潮位由總水位扣除，求得暴 潮偏差潮位。 若以年度分工而考慮，則 第一年：以現有我國國際港口及臨近重要河川如淡水河等出口之 水位站與流場觀測站為主，建立台灣環島（包括離島） 資訊網站芻形。 第二年：增加以國內商港及漁港現有之水位站及流場觀測站為 主，改善及擴充環島海岸水位資訊網站。 第三年：彙整所有水位站及流場觀測站資料，完成環島海岸水位 資訊網站。
接下來的第二章及第三章將針對子計畫二：「台灣環島暴潮預報 模式及數值網格產生法研究」詳加說明，至於子計畫一以及子計畫三 的 詳 細 研 究 理 論與 成 果 請 參 考 國科 會 NSC90-2625-Z-124-001 及 NSC90-2625-Z-124-002 計畫編號之子計畫成果報告書，於此不贅述。 第四章結論與建議部分則統合三個子計畫目前的研究成果，作一 簡單的結論，並提供後續研究之建議。
第二章
理論分析
2.1 文獻回顧 潮位的計算包括兩個部份：一是由天體運動所引起的天文潮 (Astronomical tide)，一是由異常氣象所引起的潮位，如颱風造成的暴 潮(Storm surge)。當颱風侵襲引起暴潮之際，若再碰上陰曆朔望之天 文大潮，則往往會造成低窪地區海水倒灌。（1994,李賢文） 天文潮的計算濫觴於 1867，William Thomson 於當年推導出調和 分析的方法來分析潮位，其原理乃是基於任何運動皆可以分解成一系 列的調和運動(harmonic motions)，因此可以利用一系列的調和函數 (harmonic functions)來逼近任何運動。雖然到 1867 William Thomson 才運用於潮汐分析之上，但實際上，此原理早在西元前 356 年就被 Eudoxas 發現。當年，他運用一系列的圓周運動來說明星球的不規則 運動。自 William Thomsom 之後，幾百年來陸陸續續又有許多學者專 家針對此技術加以改良，並利用天體運行關係推導出各星球對地球的 引力在地球上所造成的分潮週期公式。 而在暴潮預報模式的發展歷史中，可分成經驗方法、解析方法及 流體動力數值方法三種。 經驗方法是利用統計分析方式，求出海面高程與其他氣象因子的 關係。最簡單的就是直接求出暴潮位與鄰近平均風之關係，以及暴潮 位與中心氣壓之關係，在將兩者所求數值疊加而得暴潮位。1971 年 王博先生對國聖港、鹿港、外傘頂三地暴潮之推算，1975 年黃壽銘 先生研究花蓮港，1976 年魏靖松先生推算高雄港等地暴潮，及 1977 年 王 崇 岳 先 生 研 究 淡 水 河 暴 潮 ， 1992 年 Minoru 利 用 下 式 ： a(1010p)bw2cos 預 測 東 京 (Tokyo) 、 名 古 屋 (Nagoya) 及 大 阪
(Osaka)三個臨海的大城市暴潮。2000 年蔡瀚陞研究淡水河口颱風暴 潮水位；。2002 年 Ou 等人，應用此法，計算台灣新竹、台中、東石、 與花蓮等區域，此皆屬經驗法則。
解析方法則是利用流體動力方程式推算暴潮。流體動力方程式包 括質量守恆及動量守恆式，其為一組非線性方程式。要解這組非線性 方程式則勢必將方程式簡化為線性方程式。而此通常只能適用於簡單 地形，對於自然界複雜的海岸地形就顯得不適宜(1967,Bretschneider)。 流體動力數值模式如同許多其他的數值計算一般，在電子計算機 問世後才發展起來，利用電子計算機快速計算之能力求解非線性方程 式之逼近解，因此毋需簡化條件，在自然界複雜的海域也能應用，是 目前最被廣泛使用的方法。 流體動力數值的先驅是德國海洋學家 Hansen(1956)，他在 1953 年以此法計算北海暴潮，獲得令人鼓舞的結果。1958 年 Platzman 以 此法計算密西根湖的暴潮，1959 年 Fisher 以與 Hansen 不同之格點構 造同樣推算北海之暴潮，同時討論到計算的穩定性問題，以及平滑技 巧和格點間距對演算之影響。1963 年 Platzman 再度以數值法計算北 美 伊 利 湖 之 暴 潮 ， 提 出 共 軛 格 點 模 式 (Conjugated Richardson Lattices)，討論如何處理沿邊界之格點。Harris 和 Jelesnianski(1964) 則發展出暴潮偏微分方程式之近似形式，同時將邊界以定差方式表 示，並討論其誤差。Jelesnianski 在 1965 年更以粗細兩種不同網格系 統計算暴潮，在大海域先用粗網格計算，再代入小海域細網格之邊界 條件中計算暴潮。Sielecki 則於 1968 年以能量不滅觀點求得穩定解。 1979 年李賢文先生以完整的流體動力模式預報台灣海峽沿岸之 暴潮。1984 年李賢文先生更進一步發展環繞台灣周圍海域之颱風暴 潮數值模式。1987 年劉肖孔先生發展一個三維颱風暴潮數值預報模 式，模式中並考慮了溫度和鹽度的變化。 1976 年美國聯邦緊急事故發展署(FEMA)為了洪災保險之需，發 展出一套暴潮預測模式，稱為 FEMA 模式。1983 年美國國家科學院 (National Academy of Sciences, NAS)會議中推薦此模式並建議修改此 模式以因應需要。因此有了 1985 年新版之 FEMA 暴潮預測模式。
在 Jelesnianski 領導下，1992 年美國氣象局發展出 SLOSH(Sea, Lake, and Overland Surges from Hurricane)模式推算海水倒灌之地區，
供民眾避難之用。 1996 年 Watson 等人研究，影響暴潮的因素很多，包括波流、地 形效應、暴風路徑與前進速度等相當重要，因為發生溢淹的位置與跟 當地發生高潮時間有很大的關係，暴風強度是決定暴潮偏差量主要因 素，但在預測上確有困難。而地形效應上顯示，V 型港灣會累積能量， 故暴潮偏差量會增大。 就美國氣象局暴潮模式的發展史來看，在 70 年代，模式計算僅 限於海岸線外的水域，即假定一高牆立於海岸線上，海水只能堆高， 無法穿透。之後，美國氣象局根據 60 年代暴潮觀測資料的努力，用 一組歷史最高水位觀測值來調整風力場係數。80 年代，氣象局推動 海灣與內陸漫灘的模式，利用粗細網格交疊之計算方式求算結果。 目前，SLOSH 模式在颱風路徑的預測精確度仍有限，登陸在 24 小時內，平均誤差仍在正負 100 海哩，而少許路徑的誤差，卻往往有 天壤之別的結果。因此，不適用單一路徑模式做暴潮預報。 SLOSH 和 FEMA 模式相異之處，在於前者用極座標網格，而 FEMA 則是用卡氏座標網格計算。後者並以網格對角線來逼近海岸 線。在風場模式上，SLOSH 分成風場為海岸線以外之”海風”，及海 灣內與淹水區之”湖泊風”。海風取早期外海模式(SPLASH, Special Program to List the Amplitude of Surges from Hurricanes)(Jelesnianski, 1992);湖泊風則根據 1949 年美國陸軍工程部設立的觀測站資料所得 出之摩擦係數修正外海模式而得。 至於在數值網格產生法方面，Thompson 等(1985)利用保角映射 原理 (Comformal mapping) 解橢圓性之偏微分方程式，並配合柯西里 曼條件 (Cauchy-Riemann conditions)將一由四條平滑曲線所圍成，四 個頂點為直角的不規則區域轉換成矩形區域，並在矩形區域中建立網 格，再做反轉換，使得不規則區域內產生網格。 Tsay 等(1989, 1990) 將邊界積分元素法運用於 Thompson 等 (1985) 所提出的區域轉換理論，以做為保角轉換之依據，在特定的 不規則區域中產生網格，而建立一自動網格形成模式。並將此應用於
緩坡方程式 (Mild-slope equation) 來描述波浪運動。另於 1990 年發 表的技術報告中提及，其優點為座標轉換與反轉換之關係在建立之 後，如需調整網格時可以不必重複再解座標轉換與反轉換所需之偏微 分方式。 Tsay(1997)利用複變映射理論，將不規則區域轉換成規則區域， 再利用 Tsay 等(1989,1990)之研究在任意不規則區域中形成正交網 格，由於複變映射理論之應用，將 Tsay 之模式擴展至任意不規則區 域，使 Tsay 之模式應用更為廣泛。 2.2 天文潮理論 天文潮係地球自轉情形下，地表受太陽、月球及其他星體重力作 用所產生近週期性之水位變化。運用前人研究天體運行關係推導出的 各星球對地球之引力在地球上所造成的分潮週期公式，吾人可求出地 球上不同經緯度受到星球引力作用下的潮汐週期。然後再根據這些週 期進行調和參數的推求。調和分析的原理，便是假設潮汐的漲落皆可 以由一系列和星球引力相關的調和項(sin 或 cos 函數)疊加之後逼近達 成。
k r r r r a T t C y t y 1 0 2 cos ) ( (2-1) 其中，ya表天文潮位，y0代表平均水位，為未知數；Cr為調和 項 cos 函數的振幅，為未知數；Tr為各分潮週期，已知數；則為相r 位角，乃未知數。為求得最佳潮位估計，可運用最小平方法原理。假 設t為 t 時間分析資料之誤差，即，t yt ya(t)，yt為實測潮位。 則誤差平方和為 2 1 1 1 0 2 2 cos
t n t n t t k r r r r t t T t C y y E (2-2) 欲使誤差平方和最小，需 0 0 y E ， 0 r C E ， 0 r E 。其中，
2k 1
個方程式可解
2k 1
個未知數，如此便可求得 k 個分潮的調和參數(即振幅與相位角及平均水位)，代回(2-1)即可預測任意時間 t 的天文 潮位。在實際計算上，會遭遇到現場量測水位資料闕漏的問題，因此， 必須應用第一次求出的天文潮位反求引用資料中闕漏的現場水位資 料，然後以此修正過的資料重新再求出新的天文潮分潮參數。利用此 新參數可求出任意時間的天文潮位。而經過非颱風時期的現場潮位資 料驗證後，便能確認天文潮模式的準確性。 天文潮模式的建立流程如圖 2-2-1，為：(1)取得足夠長度且品質 良好的潮位資料；(2)過濾掉明顯錯誤的實測資料，確認資料品質。 例如，由過去紀錄中，河口水位即使在颱洪時期也未曾超過 3 公尺以 上，因此若平日水位有高於或低於 3 公尺者，便可視為不合理值而去 除；(3)進行調和分析，得到天文潮分潮參數，並利用此參數計算天 文潮位；(4)繪圖觀察計算值與觀測值的水位歷線是否相符，如果有 明顯的相位差，表示有重要的分潮參數沒有考慮到，必須回頭檢視天 文潮分潮選擇是否恰當，調整過後，再重新計算一遍。(5)統計計算 值與觀測值的年平均誤差，如果超過設定的門檻值，表示缺漏的資料 已經發生影響，必須利用前一次計算所得之天文潮位將原始實測資料 闕漏部分補遺或修正不合理值，進行第二次的調和分析，而獲得一組 新的天文潮參數及天文潮位。 其中，天文潮分潮的週期可從天文及潮汐相關書籍查表而得。
圖 2-2-1 天文潮模式建立流程圖 2.3 邊界符合座標系統理論 由於現有模式計算偏導數時會出現些許誤差，尚無法準確描述座 標轉換間尺度因子。因而擬對現有的理論基礎重新檢討，針對有誤差 的地方進行修正。以下就現有模式的理論基礎作一說明： 邊界符合座標網格形成，先採用複變映射理論將不規則區域轉換 成超矩形區域(Hyper-rectangular region)(如圖 2-3-1)，再配合邊界積分 元素法解拉普拉斯方程式 (Laplace equation)，而將超矩形區域轉換至 矩形區域而在矩形區域中形成網格。其中使用邊界元素法時，處理角 落 (Corner)的多重方向導數(Multiple normal derivatives)，需另外處 理，而使轉換能順利進行。並須進行反轉換之計算，利用邊界積分法， 否 天文潮分潮參數及天 文潮位計算 蒐集潮位觀測資料，清除不合理之觀測值 調整分潮個數 計算與觀測水位 歷線相位相符? 結束 計算與觀測值年平均 誤差超過門檻值 是 原始觀測資料闕 漏部分或不合理 值以計算值代入 否 是
與柯西里曼條件( Cauchy-Riemann condition)，將矩形區域中之網格反 轉換回超矩形區域，再利用複變映射理論反轉換回原幾何形狀，最 後，原幾何形狀之網格便可建立完成。 圖 2-3-1 座標轉換示意圖 此外，將不規則區域轉換成超矩形區域之過程中，乃利用複變映 射理論推導座標轉換系統偏導數之關係。而超矩形區域轉換至矩形區 域之過程中，則利用邊界積分元素法，求區域內偏導數的值。將上述 所 求 得 之 偏 導 數 組 合 後 ， 可 建 立 原 幾 何 座 標 系 統 轉 換 (Physical coordinate system)至矩形區域座標系統(Transformed coordinate system) 之偏導數。 一、 複變映射理論—建立超矩形區域 由於 Tsay 等(1989, 1990) 之自動網格形成模式，僅適用在如(圖 2-3-1) 的區域：圖中四點 A、B、C、D 夾角均為 2，其餘邊界上任 意點皆為平滑曲線，這樣的區域於現有模式中稱為「超矩形區域」。 為了能銜接 Tsay 等 (1989, 1990)之模式，必須先將各種不規則區域均 能轉換為超矩形區域。而這樣的轉換方法即以「複變映射理論」為基 礎而建立。 如圖 2-3-2，以 W 複數平面夾角之區域轉換至 Z 複數平面夾 角之區域為﹕
圖 2-3-2 複變映射示意圖 1 A W dW dZ (2-3) 式中， A 積分後得 A W W A Z Z 0 1 ( 0) (2-4) 式中，Z0 為 Z 複數平面的支點（Branch point），W0為 W 複數 平面的支點。在式 (2-4) 中，因為 Z 平面與 W 平面皆為複數平面， 且複數函數(WW0) 為多值函數。為保持函數為單值，需建立一條人 為的界線W B0 。參見(圖 2-3-2) ，其中W0點稱為支點，B 點在無窮遠 處，（當然其他以W0為始點的射線都可），這條界線稱為支線（Branch line）或分枝切割（Branch cut）。在區域轉換的過程中，為使得區域 不至被分割，其選取的支線必須在轉換區域之外，如此才能確保轉換 區域之連續性。 在此考慮複變映射理論之反轉換，自 Z 平面夾角之區域轉換 至 W 平面夾角之區域，其反轉換方法如下： A Z Z A W W 1 0 0 [ ( )] (2-5)
因此，利用以上之理論，可將不規則區域，選取邊界上四點轉換 為直角，其餘邊界轉換為平滑曲線；經過複變映射理論多次的轉換， 即可將任意不規則的區域轉換為超矩形區域。而經由複變映射理論之 反轉換，亦可將超矩形區域反轉換回原不規則幾何區域。 複變映射理論正轉換與反轉換之計算方法，利用映射原理轉換過 程中其偏導數之推求方法，在使用映射原理轉換時，座標系統轉換之 偏導數可得如下： -1 1 0)] ( [ Re A Z Z A X U (2-6) -1 1 0)] ( [ Re i A Z Z A Y U (2-7) -1 1 0)] ( [ Im A Z Z A X V (2-8) -1 1 0)] ( [ Im i A Z Z A Y V (2-9) 式中，U 和 V 為 W 平面上的實部及虛部。所求得之偏導數，亦 可驗證其符合柯西里曼條件( Cauchy-Riemann condition ) X V Y V X U Y U , (2-10) 二、 邊界積分元素法--建立矩形區域 利用邊界元素法解拉普拉斯方程式，使得 X-Y 平面之超矩形區 域，經由保角轉換至 平面之矩形區域。如(圖 2-3-1），X-Y 平 面上之封閉區間 ABCD，邊界為四個區段：AB、BC、CD 與 DA， 在平面上選取適當的矩形A B C D，使 A、B、C、D 四點對應 至A、B、C、D四點。且使得A B、C D之值為某一常數，D A、 B C 之值為某一常數，因此，可採用拉普拉斯方程式作保角轉換：
0 2 2 2 2 Y X (2-11) 0 2 2 2 2 Y X (2-12) 且 , 必須滿足柯西—里曼方程式 Y X (2-13) X Y (2-14) 利 用 二 維 拉 普 拉 斯 方 程 式 的 基 本 解 ( 自 由 空 間 格 林 函 數 free-spaceGreen‘sfunction)ln r 以及變分法(variation method)，可將方 程式(2-11)、(2-12)改成邊界積分方程式如下： ds n r n r P) (ln ) ln (
(2-15) ds n r n r P) (ln ) ln (
(2-16) 式中， 為邊界，ds 為沿著邊界之微小線段，r 為區域內或邊界 上一奇異點 P 至邊界之距離。 2 1 2 2 ] ) ( ) [(X Xp Y Yp r (2-17) n r r r n 1 ) (ln (2-18) 當奇異點 P 在區域內時，2；當 P 點在邊界上時，其值 為邊界內部之夾角值。因此可利用(2-15)、(2-16)式，計算出邊界上之 所有 , 值或 n n , 值。 由於在矩形區域邊界上所有 X、Y 值皆已知，但邊界上之所有 n X 、 n Y 值皆未知，在 Tsay 等(1989, 1990)之研究中，於處理矩形區 域邊界之導數時，利用柯西里曼條件(Cauchy-Riemann condition)及採 用線性計算，導致邊界之導數產生誤差。而現有模式採用邊界積分元 素法解拉普拉斯方程式，能精確地求出邊界上之所有 n X 、 n Y 值。值得一提的是，在矩形區域的四個角落(Corner)，重方向導數（Multiple normal derivatives)，可利用柯西里曼條件(Cauchy-Riemann condition) 增加四條方程式處理，即可順利進行。 在利用邊界積分法求解 X2 0過程中，發現邊界若有 N 個節 點，因其邊界條件為德列契列特（Dirichlet）形式，X 值已知，但 n X 未知，且在四個角落（Corner）因有多重方向導數（Multiple normal derivatives ) X n X n , ( )皆未知，如此便只有 N 條方程式，但卻有 N+4 個未知數。同理在求解 2 0 Y 過程中亦相同。因此，本研究 之作法乃採用 2 0 X 與 2 0 Y 同時求解，如此則將有 2*N 條方 程式；另外加上A、B、C、D四點之柯西里曼條件(Cauchy-Riemann condition)，可增加 8 條方程式，則共有 2*N+8 條方程式，而未知數 正好有 2*N+8 個。如此即可將邊界上所有 X n n X n n , Y , ( ), ( Y) 皆可求得。 至目前為止，已經順利完成將不規則幾何區域轉換成矩形區 域，因現有模式之目的為在矩形區域中建立網格，以方便有限差分 （Finite difference method）之計算。所以可視有限差分法之需要，建 立矩形中之網格。 利 用 上 述 已 經 求 得 之 矩 形 區 域 邊 界 上 的 偏 導 數 X n n X n n , Y , ( ), ( Y)與 X、Y 值，做反轉換至超矩形區域，而在 超矩形區域中建立網格。 其反轉換之公式如下： ds n X r n r r X P X( ) ln 2
(2-19) ds n Y r n r r Y P Y( ) ln 2
(2-20) 式中，P 表示為區域內之點 2 1 2 2 ] ) ( ) [( P P r (2-21)在式(2-19)及(2-20)積分方程式中，右邊為邊界積分，且邊界上的 X、Y 值與 X、Y 之偏導數皆已知，利用邊界積分法計算，即可求得 超矩形區域內之 X、Y 值，如此便將矩形區域中之網格點成功轉換至 超矩形區域。為建立超矩形區域轉換至矩形區域的偏導數，利用邊界
積分方程式（Liggett and Liu, 1983）：對 , 分別偏微分得：
ds n r n r r P (ln ) ) 1 ( ) ( 2
(2-22) ds n r n r r P (ln ) ) 1 ( ) ( 2
(2-23) 複變映射原理與邊界積分法將轉換過程中，座標系統之偏導數已 求出。將上述所求出之偏導數利用連鎖律（Chain rule）組合，而建 立原幾何區域轉換至矩形區域之偏導數。 2.4 建立邊界符合座標系統所遭遇的問題與對策 本研究中所產生的邊界符合網格系統，為保角網格系統，也是自 原區域轉換到矩形區域時，座標轉換式中最精簡，限制最多，最不易 產生，但又最完美的轉換 (Thompson et. al., 1985；Ryskin and Leal,1983)。當一x-y座標系統轉換到另一個-系統時，其轉換矩陣可表 為 2 2 2 2 22 12 21 11 y x y y x x y y x x y x g g g g A (2-24) 則在原物理區域內的 Laplace 方程式20,20可改寫為 0 2 2 2 11 2 12 2 2 22 x g x g x g (2-25) 0 2 2 2 11 2 12 2 2 22 y g y g y g (2-26)
當g12 0時，即 0 2 2 22 11 2 x g g x (2-27) 0 2 2 22 11 2 2 y g g y (2-28) 而在計算區域座標系統所產生的矩形網格對應到原區域的 邊界符合座標網格系統即稱為正交網格系統。 若 進 一 步 加 上 g11 g22 條 件 的 話 ， 則 式 (2-28) 可 化 簡 為 0 , 0 2 2 x y 。此時，所產生的邊界符合網格系統稱為「保角網格」 (Conformal grids)。 過去的文獻認為，此保角網格系統的產生有其困難性，穩定性上 也有問題。尤其在邊界形狀過於崎嶇時，常會有網格疊置(Overlapping) 的問題，即網格交錯或者內部網格跑到域外去的現象。稱此為平滑性 (Smoothness)問題，為了製造出平滑的網格系統，常常要犧牲網格的 保角性。 以平滑性考量為主的網格系統以 Thompson 等人於 1976 年提出 的 TTM (Thompson, Thmes, and Mastin)系統為主，而以正交性為考量 的則以 Rinsky 於 1983 年提出的 RL (Rinsky and Leal)系統為代表。前 者以一對相關的橢圓型偏微分方程組控制網格產生 ) , ( 2 2 2 2 P y x 以及 ) , ( 2 2 2 2 Q y x (2-29) 為了平滑網格，不得不犧牲保角性，於方程式等號右側加上控制 項P(,)以及Q(,)。基本上，此二控制項之參數可利用試誤法(Try and error)來獲得最佳的平滑效果。由於網格線多半在靠近邊界時遭遇 到 扭 曲 ， 因 此 可 將 此 二 參 數 設 定 為 自 邊 界 向 內 域 的 指 數 遞 減 (Exponential decay)型態；亦即，越靠近邊界，修正越多，越接近內部， 則越退化成拉普拉斯方程式之保角網格系統。後續相關的研究，有一 大部份是延續 TTM 系統架構，而針對不同的P(,)、Q(,)控制項型 態所產生的網格系統加以討論。
至於 RL 系統，則強調網格的正交性(Orthogonality)以及尺寸因子 (Scale factor)的等稱性。簡單而言，雖不能如保角網格般處處符合柯 西-里曼條件(Cauchy-Riemann conditions)，但卻藉由引入扭曲函數 d f (Distortion function)來維持正交性 y f x d (2-30) x f y d (2-31) 以上兩式分別對、作偏微分的動作後再相加，可得下列二方 程式 0 ) 1 ( ) ( x f x f d d (2-32) 0 ) 1 ( ) ( y f y f d d (2-33) 此處的 fd可定義為 h h fd (2-34) 2 2 11 y x g h (2-35) 2 2 22 y x g h (2-36) (2-32)、(2-33)即被稱為 RL 系統。當 fd 1時，此系統即退化成保 角 網 格 系 統 ， 後 續 許 多 研 究 ， 多 集 中 在 以 不 同 的 方 式 計 算 出 d
f (Luis,1996; Jeng, 1999; Akcelik, 2001; Zhang, 2004)。
由以上的文獻回顧可知，保角網格為最佳的網格系統，但是由於 產生保角網格不易，因此退而追求網格系統的平滑性以及正交性，因 而個別發展出 TTM 以及 RL 網格系統。至於爲何產生保角網格不易，
中並未見到相關的討論。本研究即根據邊界元素法研究中的幾何奇異 性(物理奇異性)加以探討網格扭曲的原因，並提出以複變轉換做為對 策解決此問題，然而此法亦有其先天上的限制，在運用上不可不知， 其限制詳見於本節後文。 在判斷網格系統品質好壞時，可利用 MDO、ADO、MAR、AAR 等四項指數來決定。其定義如下：
|90 |
maxi,j o i,j MDO
1 2 1 2 , | 90 | 2 1 2 1 m i n j j i o n m ADO 其中，下標i以及 j分別代表在x、y方向的網格點位置。m、n則 為此二方向上的網格數。角度定義如下 h h g j i 12 1 , cos
1 2 1 2 , max 2 1 2 1 m i n j j i d f n m AAR (2-37) 為位於
i, j 點的扭曲函數 fd。MDO 所代表的是網格中偏 離正交性最大者，ADO 所代表則為所有網格偏離正交之平均值，MAR 為網格曲線最平滑者，AAR 則為平均的網格平滑度。當網格為完全 正交時，亦即，g12 0，MDO 以及 ADO 的最佳值為 0 度，最差為 90 度；反之，當網格系統為保角網格時，則進一步g11 g22，此時， MAR 以及 AAR 為最佳值 1，所表現出來的網格則最為平滑。 綜合而言，保角網格系統之 MDO、ADO 為 0，MAR、AAR 為 1 ， 不但是最正交者，也是最平滑的網格系統，所以其為網格產生法 中最受歡迎的。而本研究所產生的網格系統，正是此種網格系統。 然而，在保角網格系統產生的過程中，如圖 2-4-1 所示，將會遭 j i d j i d j i f f MAR , , , 1 , max max j di f ,遇到兩個問題：(1)由於邊界元素法中的數學性奇異性，在運用邊界 元素法做映射時，內部網格線不可過於靠近邊界，必須距離至少 6 10單 位長度以上，否則會產生數值發散的現象。關於這點，本研究利用幾 何解析方法，經由正確估計邊界元素法中關於角度的積分項於基礎點 逼近邊界時的角度值，可大大降低數學奇異性所帶來的誤差。以目前 的技術，在沒有幾何奇異性的情況下，可以產生至少距離邊界 9 10以 上的網格線(Wang, 2005)，對於需要計算到極靠近邊界的流體力學問 題，此類網格的精密度將有決定性的影響。(2)由於邊界崎嶇時所帶 來的幾何奇異性(物理奇異性)問題，使得當要計算具物理奇異性的奇 異點時，其勢能值會產生較大的誤差，表現在網格上，就是網格線扭 曲或是網格線重疊交錯的現象。欲解決此問題，可利用複變函數先將 幾何奇異點的奇異性消除掉，最後形成四個角為直角，且四邊平滑的 超矩形區域(Hyper-rectangular region)，之後，再利用邊界元素法進行 平面轉換。然而，複變轉換並非萬能，在轉換之間，需要注意 branch-cut 的問題，以及所衍生的角度限制式。 其原理簡述如下： 如圖 2-4-2 中所示，當使用(2-38)式的複變轉換將 z 上某一點轉換 圖 2-4-2 自 z 到 w 平面之複變轉換示意圖
到 w 上時，由於角度放大的原因，原來在 z 平面上2角度內定義的 點，分佈到 w 平面上的角度範圍變成了4，這超越了複數平面四個 象限所能定義的範圍，因而產生了多值(Multi-value)的問題，亦即， 原來 z 平面上等相位角的兩個點，映射到 w 平面上時成為一點。這 個多值的問題，表現在網格產生上，便是網格邊界的扭曲變形，無法 形成多邊形。因為本研究產生網格的方法中，第一步是將不規則邊界 利用上述複變函數的方法轉換到四個角直角、四邊平滑的超矩形區 域，之後才在此區域內產生網格，只要邊界在轉換中沒有出問題，則 於超矩形區域內切割網格也不會出問題。因此，複變轉換中多值的問 題在此處顯得格外重要。 針對此問題，李曼(Riemann)於 1991 年提出三度空間的連接平面 構想，如圖 2-4-3 所示 圖 2-4-3 李曼平面示意圖 其意義為，利用三維變數來紀錄相位角變化，以達到一對一映射 的目的。 另外一個工程應用上的方法，乃是利用封閉的多邊形內角關係， 做為限制式，來確保 branch-cut 線不穿越到區域內。應用這樣的原理， 又可以下列兩種方法為之： 1. 改變內角
由於複變轉換造成非一對一映射的狀況出於當原來物理區 域上的角度被「放大」映射到超矩形區域上成為 90 度或者 180 度時，有可能使得平面上其他節點相對此點的轉角超過 360 度而 造成多對一映射的情況，因此必須將內角為銳角的情形改善為內 角為鈍角的情況，使其成為「壓縮」而非「放大」計算域的情況。 2 (2-38) 其中，為轉點之外其他的邊界節點內角，為轉換前的角 度，為轉換後的角度。其意義為，轉換倍率不可使其他角度超 越 2。由於本研究所轉換的超矩形區域僅含有 180 度以及 90 度 兩種角度，因此，對於原區域中欲轉換成 180 度的角度而言，即 ，其限制為 2 0 2 (2-39) 對於 2 者 2 0 2 4 (2-40) 假設地形最為崎嶇的情況下，以至於有某個節點相對於轉點的 角度為2，則上面二式可改寫成 ，其限制為 2 (2-41) 對於 2 者 4 (2-42) 因此，為了確保複變轉換不出問題，對於欲轉換到超矩形 區域上成為 90 度的原區域角度，必須限制大於 45 度，而欲轉換成
180 度者，則必須限制大於 90 度。根據以上二式，本研究推展出一 簡便的邊界平滑化技術。如圖 2-4-4 所示，在銳角兩邊極小距離內 增取 B、C 兩點，使得ABAC，則原角度便會被平滑化成鈍角
2 ，此角可滿足以上兩式。由於在計算上，AB與AC可取非常小 的距離，因此，即使邊界改變，在物理上，對於計算應無太大之影響。 圖 2-4-4 內角平滑化示意圖 為了驗證理論的可行性，本研究取一蜂窩型區域產生網格系統 (圖 2-4-5)進行研究。研究發現，若不在 A 以及 D 點採用上述的方法 加以平滑化，將難以產生保角網格系統。也就是將 A,D 的 3 角利用前 面所述之方法於此二點非常靠近的附近各加兩個點以替換此點，造成 兩個 3 2 的角度，節點數也因此增加為八點。此範例顯示上面的理論 是可行的。圖 2-4-5 蜂窩狀區域範例，A、D 再各細分成兩個內角。 2. 利用封閉多邊形的角度限制式進行複變轉換順序的調整
雖然，利用上面的方法可以有效避免複變轉換時角度多值的 問題，然而，當區域的節點距離過大時，複變轉換的保角性誤差 變大(Schinzinger and Laura，1991)，則轉換時的角度控制將更加 困難。因此，本研究發展出第二種角度控制技術，也就是封閉多 邊形的角度限制式，作為複變轉換點順序調整的依據。 如圖 2-4-6 所示，對於一封閉 n 多邊形而言，若各外角為A1、 2 A 、A3…An，則
1
2
3
2 3 2 1 A A An n A (2-43) 其意義為：當每次做一次複變轉換時，便要檢查一次外角和是否 有符合多邊形的角度和2，若不符合，表示中間有邊界交錯的情形出現，則必須更換各角度複變轉換的順序，直到符合角度限制式為 止。若更換順序仍無法達成避免邊界交錯之目的，則需用第一個方法 加以輔助，交互運用，達成目標。則其必須在每一次做複變轉換時加 以判定。 圖 2-4-6 封閉多邊形角度限制式示意圖 以下以蘭嶼島的範例說明此方法。圖 2-4-7 為台灣蘭嶼島，邊界 點數為 255 點。根據式(2-43)作為判斷式，可調整複變轉換之順序進 行轉換，經過 255 次的複變轉換之後，原區域轉換成一四個內角直 角，四邊平滑之超矩形區域。
圖 2-4-7 以蘭嶼導範例說明角度限制式之使用。 圖 2-4-8 根據角度限制式選擇複變轉換順序。經過 255 次的複變 轉換之後，原區域轉換為一超矩形區域。 簡言之，保角網格系統產生過程，在矩形區域轉換至超矩形區域 的過程中，過去若取矩形區域中離邊界小於106距離來模擬邊界，則 會發生數值計算收斂性的問題。今藉由對邊界元素法奇異性問題的研
析，此一困難已得以克服 (Wang & Tsay, 2005)。其原理乃是將產生網 格時遭遇的奇異性問題分成數學奇異性(Mathematical Singularities)以 及幾何奇異性(Geometrical Singularities)兩部分進行。數學奇異性的部 分利用解析幾何的方式求出邊界積分方程式中，當基礎點(Base point) 逼近邊界點，解析解角度項的正確求法。經過解析方法以及範例證 實，若能正確計算角度值，則計算時並不會受到所謂數值邊界層誤差 的困擾。然而，由於推導的過程所使用的邊界元素為線性元素，因此 當邊界較崎嶇時，便需要使用較密的邊界點來模擬，這是它的缺點。 另外，當遇到退化邊界(Degenerated Boundaries)時，也就是內角偏離 180 度過大時，如 cut-off wall 或是 corner problem，在接近尖端(tip) 或轉角(vertices)處會有勢能集中的現象，這種由邊界幾何形狀所引發 的奇異性問題，稱為幾何奇異性問題(Geometrical Singularities)，遇此 問題，則須藉助複變轉換(Complex mapping)的
此成果已經為著名的期刊 Engineering Analysis with Boundary Elements(EABE)所接受，並且預期短時間內將可發表。所以可在矩形 中任意建立規則的網格，使得在應用時，因矩形邊界而能很容易使用 有限差分法做計算，而將問題簡單化，並提高計算之精度。 2.5 暴潮理論 本文所指的暴潮起因，單指為颱風所造成的現象。在計算暴潮 時，可分成流體動力模式與颱風模式兩種模式交互進行，由颱風模式 所計算出的風場與壓力場提供流體動力模式所需。流體動力模式的邊 界條件由颱風模式所計算出來之虹吸水位決定，即流體動力模式的動 力來源為颱風，包括風剪力、壓力梯度、負壓造成海水面上之高程。 下一節所述之颱風模式與流體動力模式係參考 FEMA 模式。 1 控制方程式 (1)颱風模式 颱風規模參數包括颱風中心氣壓深度P(mb)、中心最大風速
