顯式有限解析法模式於超亞臨界混合流之研究

全文

(1)國立交通大學土木工程學系碩士論文. 顯式有限解析法模式於超亞臨界混合流之研究. 研究生:楊璨仲指導教授:葉克家博士. 中華民國九十五年一月.

(2) 顯式有限解析法模式於超亞臨界混合流之研究 Study on Mixed Supercritical and Subcritical Flows Using Explicit Finite Analytic Model. 研究生：楊璨仲. Student：Tsan-Chung Yang. 指導教授：葉克家. Advisor：Keh-Chia Yeh. 國立交通大學土木工程學系碩士論文. 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 Jan 2006 Hsinchu, Taiwan, Republic of China. 中華民國九十五年一月.

(3) 顯式有限解析法模式於超亞臨界混合流之研究. 學生：楊璨仲. 指導教授：葉克家國立交通大學土木工程研究所. 摘要本研究延續 Hsu and Yeh(1996)之一維顯式有限解析法(explicit finite analytic method，簡稱 EFA)模式，發展出適於超臨界流與超亞臨界混合流流況之數值模式。EFA 求解之特點，乃在於求解水流動量方程式時，以特性法觀念解得其中變量(流量與通水斷面積)之局部解析解，並且遵守可蘭穩定性條件；邊界處理方面則透過水流之連續方程式與動量方程式，利用特性法觀念求解邊界處之變量；而在超亞臨界混合流況之內部邊界，根據福祿數的大小來判斷水躍發生的位置，且利用內部相鄰計算點的水位高程，透過外插的方式來得到超臨界流區域的下游邊界水深。為驗證本模式之正確性，本文針對單一坡降矩形渠道、不同坡降矩形渠道問題進行模擬以及分析，且透過與美國陸軍工兵團所研發之 HEC-RAS 模式的比較，說明本模式模擬之成果與適用性，最後並應用於清水溪草嶺潭堰塞湖上游河段，以探討天然陡坡河道之複雜流況。. I.

(4) Study on Mixed Supercritical and Subcritical Flows Using Explicit Finite Analytic Model. Student : Tsan-Chung Yang. Advisor : Keh-Chia Yeh. Institute of Civil Engineering National Chiao Tung University. Abstract This study extends Hsu and Yeh’s (1996) one-dimensional explicit finite analytic model (EFA) for simulating supercritical and mixed supercritical and subcritical flows. The essence of the EFA is the adoption of the concept of method of characteristics to the momentum equation for solving the local analytic solution of the dependent variables (i.e., discharge and cross-section area of flow). To ensure stability of the scheme, Courant condition should be obeyed. The dependent variables at the upstream and downstream boundaries are obtained through the method of characteristics. For the interior boundary condition at mixed supercritical and subcritical flows, the locations of the occurrences of hydraulic jumps are determined according to the values of Froude Numbers. And water depths for supercritical regime at downstream boundaries were calculated. This was done through the method of external interpolation, by utilizing the water surface elevations of the interior neighboring computational points. To test the accuracy of the model, this study simulates and analyses the cases on both a single and different slopes of rectangular channels. To elaborate the functionality and suitability of this model, these simulated results were then compared with the ones done by the HEC-RAS model, developed by the US Army Corps Engineers. Finally, the study applies the results from the model to the upstream of Chin-Shui River at Tsaoling landslide dam area, and discusses the flow condition of natural river.. II.

(5) 謝誌承蒙吾師葉克家教授孜孜不倦、不厭其煩的悉心指導，使得本論文得以順利完成，在此致上最誠摯的謝意。感謝口試委員國立成功大學蔡長泰教授、國家高速電腦中心副主任蔡惠峰博士之細心匡正與建議，使得本論文得以更趨完整。此外感謝楊錦釧教授、張良正教授、許勝田教授、蔡武廷教授在課業上的教導，讓學生受益良多。感謝璁哥、勇哥、東哥、夢祺、恩添、哲宏、昇學、大槍、仕哲學長在學業與研究上細心指導與建言，感謝胤隆、士偉、小宇和曉萍、秀容學姊及同窗好友峰志、俊毅、于軒、益家、弘恩的互相砥礪與扶持，使彼此共度愉悅的研究時光，亦感謝宣汝、力瑋、仲達、欣瑜、浩榮、偉國、雅婷、宥達學弟妹豐富我的研究生涯，因為有了大家的參與，使我這兩年的研究生活過得充實又多采多姿，永難忘懷。最後謹以本論文獻給我親愛的父母、兄長景翔及摯愛的女友宜靜，感謝你們對我的關心與支持，我所有的一切都來自於你們無止境的付出與關愛，謝謝你們。. III.

(6) 目錄中文摘要 ..................................................................................................Ⅰ 英文摘要 ..................................................................................................Ⅱ 謝誌 ..........................................................................................................Ⅲ 目錄 ..........................................................................................................Ⅳ 圖目錄 ......................................................................................................Ⅵ 符號說明 ..................................................................................................Ⅹ 第一章緒論 ..............................................................................................1 1.1 前言..............................................................................................1 1.2 研究方法......................................................................................1 1.3 文獻回顧......................................................................................2 1.3.1 常用水理定床模式 ..........................................................2 1.3.2 超臨界流...........................................................................4 1.3.3 超亞臨界混合流...............................................................4 1.4 本文架構......................................................................................6 第二章理論基礎 ......................................................................................7 2.1 基本假設(de saint Venant 之假設) .............................................7 2.2 控制方程式..................................................................................8 2.3 顯式有限解析法簡介..................................................................9 第三章數值方法 .................................................................................. 10 3.1 齊次雙曲線型方程式................................................................10 3.2 混合型方程式............................................................................11 3.3 水理控制方程式之離散化........................................................11 3.4 外部邊界條件............................................................................14 3.5 內部邊界條件............................................................................17 3.6 交錯格網....................................................................................18. IV.

(7) 3.7 穩定性分析…………….............................................................. 18 3.8 模式演算之流程………...............................................................19 第四章模式驗證 ....................................................................................21 4.1 超臨界定量流流場分析………...................................................21 4.1.1 單一坡度矩形渠道.....................................................21 4.1.2 單一坡度矩形漸縮及漸擴渠道...................................22 4.2 超亞臨界定量混合流流場分析.................................................. 23 4.2.1 單一坡度矩形渠道......................................................... 23 4.2.2 單一坡度矩形漸擴渠道................................................. 23 4.2.3 陡坡接緩坡之矩形渠道 ..................................................24 4.2.4 緩坡接陡坡再接緩坡之矩形渠道 ..................................25 4.3 超臨界變量流流場分析............................................................ 25 4.4 超亞臨界變量混合流流場分析................................................ 26 4.5 結語............................................................................................ 27 第五章草嶺潭河段模擬 ....................................................................... 28 5.1 草嶺潭概述 ................................................................................ 28 5.2 流域概況及模擬範圍 ................................................................ 28 5.3 模擬所需資料 .............................................................................29 5.3.1 地文資料…………......................................................... 29 5.3.2 河床糙度(曼寧 n 值)…………...................................... 29 5.4 模擬結果與討論 .........................................................................30 第六章結論與建議 ................................................................................33 6.1 結論...............................................................................................33 6.2 建議...............................................................................................34 參考文獻...................................................................................................36. V.

(8) 圖目錄圖 3-1 亞臨界流況下之特性速度................................................... 39 圖 3-2 超臨界流況下之特性速度................................................... 39 圖 3-3 超臨界流況下之外插特性速度............................................ 40 圖 3-4 一維特性曲線交錯格點示意圖............................................ 41 圖 3-5 模式流程圖........................................................................ 42 圖 4-1 EFA&HEC-RAS 比較圖(Q=0.281cms, Dini=0.05m)............... 43 圖 4-2 EFA&HEC-RAS 比較圖(Q=0.843cms, Dini=0.1m)................. 43 圖 4-3 EFA&HEC-RAS 比較圖(Q=1.571cms, Dini=0.15m)............... 44 圖 4-4 EFA&HEC-RAS 比較圖(Q=0.1259cms, Dini=0.05m)............. 44 圖 4-5 EFA&HEC-RAS 比較圖(Q=0.3771cms, Dini=0.1m)............... 45 圖 4-6 EFA&HEC-RAS 比較圖(Q=0.7026cms, Dini=0.15m)............. 45 圖 4-7 EFA&HEC-RAS 比較圖(Q=0.0398cms, Dini=0.05m)............. 46 圖 4-8 EFA&HEC-RAS 比較圖(Q=0.1192cms, Dini=0.1m)............... 46 圖 4-9 EFA&HEC-RAS 比較圖(Q=0.2221cms, Dini=0.15m)............. 47 圖 4-10 EFA&HEC-RAS 比較圖(Q=0.281cms, Dini=0.05m).............. 47 圖 4-11 EFA&HEC-RAS 比較圖(Q=0.843cms, Dini=0.1m)............... 48 圖 4-12 EFA&HEC-RAS 比較圖(Q=1.571cms, Dini=0.15m)............. 48 圖 4-13 EFA&HEC-RAS 比較圖(Q=0.1259cms, Dini=0.05m)........... 49 圖 4-14 EFA&HEC-RAS 比較圖(Q=0.3771cms, Dini=0.1m)............. 49 圖 4-15 EFA&HEC-RAS 比較圖(Q=0.7026cms, Dini=0.15m)........... 50 圖 4-16 EFA&HEC-RAS 比較圖(Q=0.0398cms, Dini=0.05m)........... 50 圖 4-17 EFA&HEC-RAS 比較圖(Q=0.1192cms, Dini=0.1m)............. 51 圖 4-18 EFA&HEC-RAS 比較圖(Q=0.2221cms, Dini=0.15m)........... 51 圖 4-19 EFA&HEC-RAS 比較圖(Q=0.1328cms, Dini=0.05m)........... 52 VI.

(9) 圖 4-20 EFA&HEC-RAS 比較圖(Q=0.3804cms, Dini=0.1m).............. 52 圖 4-21 EFA&HEC-RAS 比較圖(Q=0.684cms, Dini=0.15m).............. 53 圖 4-22 EFA&HEC-RAS 比較圖(Q=0.0594cms, Dini=0.05m)........... .53 圖 4-23 EFA&HEC-RAS 比較圖(Q=0.1701cms, Dini=0.1m).............. 54 圖 4-24 EFA&HEC-RAS 比較圖(Q=0.3059cms, Dini=0.15m)............ 54 圖 4-25 EFA&HEC-RAS 比較圖(Q=0.0398cms, Dn_depth=0.1m)…... 55 圖 4-26 EFA&HEC-RAS 比較圖(Q=0.0398cms, Dn_depth=0.2m)…... 55 圖 4-27 EFA&HEC-RAS 比較圖(Q=0.1192cms, Dn_depth=0.2m)…... 56 圖 4-28 EFA&HEC-RAS 比較圖(Q=0.1192cms, Dn_depth=0.4m)…... 56 圖 4-29 EFA&HEC-RAS 比較圖(Q=0.2221cms, Dn_depth=0.3m)…... 57 圖 4-30 EFA&HEC-RAS 比較圖(Q=0.2221cms, Dn_depth=0.6m)…... 57 圖 4-31 EFA&HEC-RAS 比較圖(Q=0.0188cms, Dn_depth=0.1m)…... 58 圖 4-32 EFA&HEC-RAS 比較圖(Q=0.0188cms, Dn_depth=0.2m)…... 58 圖 4-33 EFA&HEC-RAS 比較圖(Q=0.0538cms, Dn_depth=0.2m)…... 59 圖 4-34 EFA&HEC-RAS 比較圖(Q=0.0538cms, Dn_depth=0.4m)…... 59 圖 4-35 EFA&HEC-RAS 比較圖(Q=0.0967cms, Dn_depth=0.3m)…... 60 圖 4-36 EFA&HEC-RAS 比較圖(Q=0.0967cms, Dn_depth=0.6m)…... 60 圖 4-37 EFA&HEC-RAS 比較圖(Q=0.0398cms, Dn_depth=0.1m)…... 61 圖 4-38 EFA&HEC-RAS 比較圖(Q=0.0398cms, Dn_depth=0.2m)…... 61 圖 4-39 EFA&HEC-RAS 比較圖(Q=0.1192cms, Dn_depth=0.2m)…... 62 圖 4-40 EFA&HEC-RAS 比較圖(Q=0.1192cms, Dn_depth=0.4m)…... 62 圖 4-41 EFA&HEC-RAS 比較圖(Q=0.2221cms, Dn_depth=0.3m)…... 63 圖 4-42 EFA&HEC-RAS 比較圖(Q=0.2221cms, Dn_depth=0.6m)…... 63 圖 4-43 Meselhe(1994)數值模擬矩型渠道示意圖 ........................... 64 圖 4-44 EFA&Meselhe(1994)渠道沿程水面線比較圖....................... 64. VII.

(10) 圖 4-45 超臨界變量流流場分析之上游入流歷線 ............................ 65 圖 4-46 超臨界變量流流場分析之上游水深歷線 ............................ 65 圖 4-47 超臨界變量流流場分析 EFA 之渠道沿程水面線….............. 66 圖 4-48 超臨界變量流流場分析 EFA 之渠道沿程水面線................. 66 圖 4-49 超臨界變量流流場分析 EFA 之渠道沿程水面線................. 67 圖 4-50 超臨界變量流流場分析 EFA 之渠道沿程水面線................. 67 圖 4-51 超臨界變量流流場分析 EFA 之渠道沿程水面線….............. 68 圖 4-52 超亞臨界變量流流場分析之上游入流歷線......................... 68 圖 4-53 超亞臨界變量流流場分析之上游水深歷線......................... 69 圖 4-54 超亞臨界變量流流場分析之下游水深歷線......................... 69 圖 4-55 超亞臨界變量流流場分析之 EFA&HEC-RAS 比較圖 …....... 70 圖 4-56 超亞臨界變量流流場分析之 EFA&HEC-RAS 比較圖 .......... 70 圖 4-57 超亞臨界變量流流場分析之 EFA&HEC-RAS 比較圖.......... 71 圖 4-58 超亞臨界變量流流場分析之 EFA&HEC-RAS 比較圖.......... 71 圖 4-59 超亞臨界變量流流場分析之 EFA&HEC-RAS 比較圖 .......... 72 圖 4-60 超亞臨界變量流流場分析之 EFA&HEC-RAS 比較圖.......... 72 圖 4-61 超亞臨界變量流流場分析之 EFA&HEC-RAS 比較圖.......... 73 圖 5-1 民國 88 年草嶺堰塞湖平面示意圖....................................... 74 圖 5-2 草嶺潭地區清水溪之上游河段距堰塞湖溢流口 0.9km~4.5km 處之 3-D 示意圖................................................................. 75 圖 5-3 200 年洪水潰壩下(方案 C1、C2、C3)不同潰壩延時潰口流量歷線比較圖........................................................................ 76 圖 5-4 距草嶺潭溢流口 0.9 公里處之上游入流歷線........................ 76 圖 5-5 距草嶺潭溢流口 0.9 公里處之上游水位歷線........................ 77 圖 5-6 距草嶺潭溢流口 4.5 公里處之下游水位歷線........................ 77. VIII.

(11) 圖 5-7 清水溪上游段之渠道沿程水面線......................................... 78 圖 5-8 清水溪上游段之渠道沿程水面線......................................... 78 圖 5-9 清水溪上游段之渠道沿程水面線......................................... 79 圖 5-10 清水溪上游段之渠道沿程水面線....................................... 79 圖 5-11 清水溪上游段之渠道沿程水面線....................................... 80 圖 5-12 清水溪上游段之渠道沿程水面線....................................... 80 圖 5-13 清水溪上游段之渠道沿程水面線....................................... 81 圖 5-14 清水溪上游段之渠道沿程水面線....................................... 81 圖 5-15 清水溪上游段之渠道沿程水面線....................................... 82 圖 5-16 清水溪上游段之渠道沿程水面線....................................... 82 圖 5-17 清水溪上游段之渠道沿程水面線....................................... 83 附錄一特性線之推導................................................................... 84 附錄二開放邊界條件處理原則..................................................... 85. IX.

(12) 符號說明 A﹕通水斷面積 B﹕渠道寬度 C﹕單位換算常數 C r ﹕可蘭數 c﹕水面波速 Fr ﹕福祿數 g﹕重力加速度 h﹕水深 n﹕曼寧糙度係數 Q﹕流量 R﹕水力半徑 S 0 ﹕底床坡降 S f ﹕摩擦坡降 t﹕時間 u﹕斷面平均流速 x﹕沿渠道中心線的距離 Y﹕渠道橫斷面距離 Z﹕水面高程 β﹕動量校正係數 Δt﹕時間間距 Δx﹕空間間距 μ﹕邊界條件線性組合係數 Φ﹕時間權重係數. X.

(13) 第一章緒論 1.1 前言台灣地區由於經濟發展迅速，在可利用的土地資源日益匱乏之下，往河川中、上游流域發展已是不可避免的趨勢。但是由於台灣地區河川坡陡流急，且中、上游河段或水工結構物下游常為超臨界流或超亞混合流之流況，故充分瞭解其水理現象對河川整治及水工結構物的設計將有極大的幫助。然而在天然河川或人工渠道坡度陡峻處、斷面束縮段、渠道底床局部明顯突起或河道中設置水工結構物，如溢流堰或動態調節閘門下之高速水流，常容易有臨界流況的發生，此現象在物理上屬於不穩定 (instability) 現象，而在數學上屬於奇異點 (singularity)的問題，所以造成數值模擬上相當大的困難。傳統上對於探討超臨界流或超亞混合流的現象，常使用水工模型試驗，藉動力相似原理之水工設計而求得原體之各物理量；然水工模型試驗之結果雖然精度高，但耗時且昂貴，故本文提供一數值模式期能以經濟有效之方式，迅速求解超臨界流或超亞混合流流況時之水理現象。. 1.2研究方法本文根據連續方程式與動量方程式，發展一維顯式有限解析法模式(explicit finite analytic model，簡稱EFA)，利用交錯格點與特性線的觀念，求解內、外部格點的流量與通水斷面積。此法具有推導容易且在程式撰寫不似隱式法般複雜。在邊界條件部分，亞臨界流之邊界條件給定方式為：上游給定流量、下游給定水深或是利用連續方程式透過特性法求解水深；超臨界流之邊界條件給定方式為：上游給定流量與水深、下游則無須給定。而超亞混合流之邊界條件給定方式為：上 1.

(14) 游給定流量與水深、下游給定水深或利用連續方程式透過特性法求解水深。EFA模式進行演算時採取可蘭數與時間權重兩參數，此兩參數值必須小於1。根據陳(1998)之模擬試驗成果，當可蘭數大時，其數值震盪可被消除；而時間權重因子越小時，數值震盪越顯著，其建議採用範圍介於0.7~1.0之間，本文將依據其試驗結果採用0.7。. 1.3 文獻回顧 1.3.1 常用水理定床模式以下就較常用之定床水理模式加以介紹： 1.HEC-RAS模式 HEC-RAS (Hydrologic Engineering Center’s River Analysis System) 模式係美國陸軍工兵團水文工程中心研發之一維水理演算模式， HEC-RAS模式除了納入原HEC-2模式之定量流模擬演算外，更加入其他演算功能。 HEC-RAS為一維水面線演算模式，適用於河床坡度小於10％之定量緩變流，其主要演算方法，為利用能量守恆(即標準步推法)來計算包括摩擦損失、突擴或突縮損失、使用者給定之附加損失，而利用動量平衡來計算包括河道連接點(合流)、橋樑、轉換流區域，可以處理亞臨界流、超臨界流及超亞混合流之水面線演算，亦具有模擬變量流的功能，而在處理由亞臨界流況到超臨界流況或超臨界流況到亞臨界流況之快速流路變化(水躍)時，亦利用動量方程式求解之。 2.FLDWAV模式 FLDWAV(flood wave routing model)模式為美國國家氣象局(U.S. National Weather Service, NWS)發展之一維渠道洪水演算模式。 FLDWAV模式發展最主要的目的即結合DWOPER模式與DAMBRK. 2.

(15) 模式的特性，其中DWOPER模式可處理多分支渠道模擬，DAMBRK 模式可處理潰壩水流、堤防溢流以及超臨界流模擬。FLDWAV模式即結合了上述兩個模式的優點，所以可處理一維變量流、超臨界流流況、多分支渠道系統以及堤防溢流所形成的洪氾問題，並且能夠提供潰壩水流、堤防溢流、即時洪水變化以及水流流經橋墩、閘門、溢洪道等流況之模擬。 3.MIKE 11 MIKE 11為丹麥水工試驗所(Danish Hydraulic Institute)發展之商業模式，為一完全視窗化之一維水理演算模式，其中流體動力 (hydrodynamics, HD)模組為MIKE 11模式之水理演算模組，可提供其他模組所需的水理資料。該模組可以模擬一維變量流流況，並且可以處理超臨界流流況以及水流流經堰、涵洞以及不規則形狀之水工結構物，並可以處理多分支或是網路系統之渠道與河川。 4.FLO-2D模式 FLO-2D模式為美國科羅拉多大學(Colorado State University)於 1998年針對該州進行洪災保險研究而發展之模式。FLO-2D模式係同步計算一維變量流及二維漫地流之流況，並且於堤防溢流時，計算堤內及堤外水流互動的流況，可以適用於矩型、梯型或是天然河道斷面，並可模擬水流流經橋墩及涵洞等水工結構物。 5.SOBEK模式 SOBEK模式為荷蘭WL|Delft Hydraulics公司與其他荷蘭顧問公司所共同研發，為一套整合河川、都市排水系統與流域管理之商用模式，並分為Rural、Urban、River三套模式，包含降雨逕流、渠流、水質模式、及時控制(real time control)及漫地流(overland flow)等模組、應用上包含河川、都市下水道系統之水理計算及區域淹水模擬，供防. 3.

(16) 洪及水資源管理者進行管理、決策與分析之用。SOBEK模式中水流 (water flow)模組係模擬一維變量流，可適用於規則斷面或是天然河道斷面，可模擬超臨界流與亞臨界流之混合流況，並可模擬水流流經橋墩、孔口、堰涵、洞及抽水站等水工結構物，亦可以應用於網狀系統渠道 1.3.2 超臨界流在超臨界流數值模擬方面，Kupiers & Vreugdenhill(1973)建立二維水深平均數學模式，並使用有限差分法進行超臨界流之數值模擬； Jimenez & Chaudhry(1988)利用MacCormack法計算且加入人工遲滯項處理超臨界明渠流；Bhallamudi & Chaudry (1992)、Rahman & Chaudry (1997)針對超臨界流束縮段之水理情況，亦利用MacCormack有限差分法進行數值模擬計算。國內研究方面，許銘熙等人(1993、1994)先後完成渠道中超臨界流波之數值模擬及自由水面二維垂直水理模式，並針對渠道中急變流況進行各項數值模擬與模型試驗；在實驗方面，Coles & Shintaku在1943年的超臨界流實驗中，做了超臨界流通過束縮段之水槽試驗；而Rouse在1951年亦做了超臨界流通過擴張段之水槽試驗。Ippen在1951年時，也進行渠道束縮段超臨界流水工模型試驗，並以四象限圖解法決定斜震波之水深與震波角。 1.3.3 超亞臨界混合流流場型態因為渠道坡度及寬度的改變、流量突變或水工結構物操作等因素而發生劇烈變化，造成超亞臨界混合流場並存在震波(shock) 及其傳遞，尤其在非定量流中，超臨界流與亞臨界流之界面會隨時移動，更增加了模擬的困難度，因此模式常需以此類流場為模擬對象，以驗證其實用性。依據Cunge et al. (1980)的分類，文獻中求解具震波水流的計算方式可分為三類：(1)震波擬合(shock fitting)法：此法需要. 4.

(17) 內部邊界條件用以決定震波的位置及震波的流場特性，然而如此的作法往往造成較複雜的計算過程；(2)僞滯性(pseudo-viscosity)法：常應用於非延散(non-dissipative)有限差分法中，利用加入人工黏滯項 (artificial viscosity)來抑制震波不連續水面附近可能產生的數值震盪；(3)震波捕捉(shock-capturing)法：透過求解水流控制方程式來獲得整個流場流況，並直接反應流場中震波的形成與傳遞，此法並不需要額外特殊的數值處理技巧，為文獻中最為常見者，本文亦屬於此方法。關於水躍計算方面之研究，Chow (1959)由上游超臨界流和下游亞臨界流兩方向分別依據緩變量流理論計算水位剖面，再以兩邊比力 (specific force)相等之處，決定水躍的位置。MaCorquodale & Khalifa (1983)利用strip-integral法來計算水躍的長度、水面線和底床壓力； Abbott et al. (1969)利用有限差分法，另外Katopodes (1984)利用有限元素法求解迪聖凡納氏(de Saint Venant)方程式以計算水躍位置。 Rahman & Chaudhry (1995)曾經使用MacCormack及two-four兩種數值模式，模擬不同福祿數流況下所形成之水躍現象。Meselhe et al. (1994) 和顏(1995)等，模擬在恆定流況下，因為渠道坡度改變而產生超臨界流與亞臨界流的相互轉換；而黃(1995)、陳(1999)則模擬在非定量流況下，利用渠道底床的突升或突降來探討超亞混合流流況。 Garcia-Navarro et al. (1992)探討在定量流況下，渠道斷面窄縮造成臨界流況的發生；而鄭(1997)、鄧(1998)、陳(1999)等亦探討非定量流況下，渠道斷面窄縮造成臨界流況的情形。而在實驗方面，Gharangik & Chaudhry在1991年之水躍試驗中，利用下射式閘門產生超臨界流後，調整尾水堰以控制尾水高度，使其產生亞臨界流之流況，並根據不同的流場，模擬不同的福祿數。. 5.

(18) 1.4 本文架構本研究係延續Hsu & Yeh (1996)一維矩形渠道之顯式有限解析法數模成果，發展擴充一維Q-A模式應用於超臨界流及超亞混合流之數值模式。由於台灣地區河川坡度變化劇烈，每逢洪水來臨將造成河床劇烈的沖淤變化，河道斷面不規則，且局部超臨界流及超亞混合流發生的情況甚多，因此研發Q-A模式以模擬河道中洪水波與水位高程變化有其必要性。透過一維顯式有限解析法，簡化渠道中之複雜水理現象，依據質量守恆及動量守恆的法則進行處理，本模式透過和商業模式的比較，針對不同渠道坡降的定量流以及變量流況進行模擬試驗以及說明，並應用於清水溪上游河段之流況分析。本文之主要內容包括：利用連續方程式與動量方程式，在滿足 de Saint Venant 之基本假設下，應用 EFA 方法將動量方程式轉換成全微分式，再利用沿特性線積分的觀念，求解動量方程式，配合適當的斷面處理與差分式，求解連續方程式，且為避免數值震盪的影響，採用交錯格網求解，在處理超亞臨界混合流況時，利用內部相鄰計算點的水位，透過外插的方式來得到超臨界流區域的下游邊界水深，並根據福祿數的大小來判斷水躍發生的位置。為驗證本數值模式的適用性與準確性，利用各種不同的渠道坡降來做定量流以及變量流演算，並和美國陸軍工兵團水文工程中心研發之 HEC-RAS 模式比較，說明本模式模擬之成果與適用性；在模式應用方面，實際應用於清水溪上游段河道，以探討本模式是否能夠合理的運用於複雜自然河川水系之洪流演算。. 6.

(19) 第二章理論基礎渠道中之水流因遵守質量守恆及動量守恆定律，因此由數學式描述之控制方程式在滿足de Saint Venant之基本假設下，隨所選擇的兩個因變數(dependent variable)而有不同型式的一維明渠流數學模式。採用之變數為水深與流速時，多用於特徵方程式之求解中；而本文則以流量及通水斷面積為因變數。 2.1 基本假設(de Saint Venant 之假設). 1.流速均勻分佈: 流速均勻分佈在通水面積上，即每一個通水斷面積僅存在一個流速，且假設動量校正係數為一，此即一維水流。換言之，通水斷面的水位線為一水平線，在此水平線上的各點水位相等。 2.靜水壓分佈: 假設渠道中水流之垂向流線曲率很小而可忽略其垂直加速度，則靜水壓分佈的假設成立。 3.渠道定量流摩擦損失估計: 渠底摩擦與紊流效應對水流所造成的損失，可以定量流摩擦律估算。 4.平均底床坡度甚小：當平均渠底坡度甚小，河床坡降 S 0 ≤ 0.1時(底坡線與水平線夾角 θ ≤ 6 ο )，重力沿渠道所造成的分力將很小，甚至可忽略不計，亦即水. 深可以垂向水面與渠底高程差表示。 5.忽略柯氏力及風力的影響: 通常在大範圍之海洋或是湖泊才會考慮到柯氏力及風力的影 7.

(20) 響，或是緯度相差很大時才會考慮，在小範圍的區域以重力、靜水壓、摩擦力為主。. 2.2 控制方程式對於不可壓縮水流之控制方程式，包括水流連續及動量方程式，為如下形式： 1.連續方程式 ∂A ∂Q + = 0 ∂t ∂x. (2-1). 式中，A為通水斷面積；Q為流量；t為時間；x為沿渠道中心線的距離。 2.動量方程式 ∂Q ∂ Q ∂Z + [ Q ( β )] + gA + gAS ∂t ∂x A ∂x. f. =0. (2-2). 式中，β為動量校正係數；Q/A=u為斷面平均流速；g為重力加速度； Z為水面高程；而摩擦坡降可表為: Sf =. Q Q n2 2. 2. C A R. 4 3. =. u u n2 4. C 2R 3. (2-3). 其中，R為水力半徑；C為單位換算常數；n為曼寧糙度係數； S f 之方向與水流方向相同。在式(2-2)中之非線性項是構成跨臨界流數值不穩定之主因，而且若將此非線性項完全忽略，則原亞臨界流數值模式有時竟能克服臨界流之不穩定特性。因此常見的做法是在可能發生臨界流處，令非線性項中之 β =0，雖可讓數值模擬在遇到超臨界流時得以繼續計算，但所模擬出來之結果則未必可用，故在模式中讓 β 值為使用者可調整之參數，並依數值模擬試驗隨之調整。 8.

(21) 2.3 顯式有限解析法簡介本研究採用之有限解析法為美國佛羅里達大學陳景仁(C. J. Chen) 教授所創，可分為隱式及顯式有限解析法兩種。隱式有限解析法已廣泛運用於計算流體力學而得到極佳的成果(Chen & Chen，1984，Chen & Choi，1990)，此方法的特色為: (1)對於不規則的幾何邊界，仍可在卡式座標使用結構性(structured)格網求解。 (2)由於個別元素體的數值離散化(discretization)係採用局部解析(local analytic solution)來近似，因此可將數值演算的捨入誤差降到最低。 (3)數值解之穩定性極佳，屬於無條件穩定之數值模式。由於隱式有限解析法之解法係針對橢圓或拋物線型(elliptic or parabolic)之偏微分方程式求解，如用於雙曲線型(hyperbolic)偏微分方程(明渠水流即為一例)並不太適合，因此乃有顯式有限解析法的發展及提出(Dai，1994)，但僅針對無自由表面之 Navier-Stokes方程式。在求解對流傳輸方程式中，對流項以特性法觀念解得式中變量之局部解析解，且使依時變量再透過適當給定之初始條件而求得，此即為顯式法之求解精神。顯式有限解析法之穩定性，和一般常見的顯式數值方法一樣，受可蘭(Courant)數小於或等於1的限制，雖不如隱式法般屬於無條件穩定，但解法較簡單，故在應用上仍有其優點。. 9.

(22) 第三章數值方法本研究採用顯式有限解析法求解水理控制方程式，可分為隱式及顯式有限解析法兩種。隱式有限解析法已廣泛運用於計算流體力學 (Chen & Chen, 1984a, 1984b；Chen & Choi, 1990)。由於隱式有限解析法之解法係針對橢圓或拋物線型(elliptic or parabolic type)之偏微分方程式求解，如用於雙曲線型(hyperbolic)偏微分方程(明渠水流即為一例)則不太適合，因此才有顯式有限解析法的發展及提出(Dai, 1994)，以下各節中將陸續介紹顯式有限解析法之求解方式與數值特性。. 3.1 齊次雙曲線型方程式 EFA水理模式採用顯式有限解析法來求解雙曲線型淺水波水流動量方程式。在各計算元素(element)中，該數值方法可用以求得動量方程式中移流項部分之局部解析解。其基本觀念可利用以下之一維度一階線性齊次雙曲線型微分方程式來加以說明： φt + uφ x = 0. (3-1). 其中，u為x-方向之速度，在各計算元素中被視為常數。當啟始條件 φ ( x,0) = ϕ ( x) 被適當給定時，可求得式(3-1)之解析解如下： φ ( xi , Δt ) = ϕ ( xi − uΔt ). (3-2). 其中， xi 為待求點之座標； xη1 = xi − uΔt 則定義一條特性線的運動軌跡。該特性線從與啟始平面之交點 φ ( xi ,0) (圖3-4中所標示之η1 點)出發，在經過 Δt 時間後，正好移流通過點 xi 。從物理觀點說明式(3-2) 的意義為，對於式(3-1)之純移流問題，計算點 xi 於某一時刻( t n +1 )之物理量，等於通過該點之特性線與啟始平面( t n )的交點( η1 點)所具有之物理量。至於 η1 點的座標( xη1 )，則可根據移流速度u藉由特性軌跡來. 10.

(23) 加以推求( xη1 = xi − uΔt )。此外，對於不同正負號的組合，根據上風 (upwind)的特性，其相依區域(domain of dependence)可自動調整。對於非線性移流方程式，如明渠流之動量方程式的移流項部分： u t + uu x = 0. (3-3). 因具有複雜的非線性特性，使得解析解並不存在。為解決此一問題，並獲得如式(3-2)的解析解形式，可透過局部線性化的方法，將定義特性線軌跡之移流速度 u 以特徵速度代替之。儘管移流速度本身隨著時間和空間改變，仍可假設此特徵速度在一計算時間間距內為常數，用以代表此時間間距內平均的移流速度。至此，所剩之問題即在於如何選擇適當的特徵速度，此將留待至方程式離散化時一併說明。. 3.2 混合型方程式對於非齊次的混合型方程式，如明渠流之動量方程式： φt + uφ x = Fg. (3-4). 等式右邊可利用已知的物理量以顯式法直接計算之，所得結果視為源項(source term)，可直接加入式(3-2)中得： φ ( x, Δt ) = ϕ ( x0 − uΔt ) + Fg Δt. (3-5). 以上的計算方法即為Dai所提出之顯式有限解析法，並以此模擬穴洞(cavity)流場。本文則將此數值方法應用於具自由液面之流場中，並在模式中引入疊代計算的流程，藉以修正特徵速度與源項。. 3.3 水理控制方程式之離散化在EFA水理模式中，動量方程式的移流項部分係採用顯式有限解析法來加以離散化，其餘的項次以及連續方程式，則採用有限差分法處理之。首先，考慮水流連續方程式，為保留連續方程式的守恆特性，. 11.

(24) 因此以控制體積法(control volume)將其離散化，空間上採用中央差分式與時間上採用前向差分式 (陳, 1998)。連續方程式在此交錯格網上之差分式如下表示： # # Ain +1 − Ain ⎡⎢ ⎛⎜ Qi(+n1+1) − Qi(−n1+1) + Φ ⎢ ⎜ Δt 2Δx ⎣ ⎝. n n ⎞ ⎟ + (1 − Φ)⎛⎜ Qi +1 − Qi −1 ⎟ ⎜ 2Δx ⎝ ⎠. ⎞⎤ ⎟⎥ = 0 ⎟⎥ ⎠⎦. (3-6). 上式中，下標i表示網格點的位置，上標n與n+1表示兩相鄰之計算時間，n為每個計算時段的啟始時間，n+1為目前程式所欲計算流場變數之時間；上標為(n+1)#者，表示每次疊代運算後所得最新之流場變數，且於每個計算時段的第一次運算中，上標為(n+1)#與n之流場變數設定相同；由於流場變數隨時間而改變，因此引入一時間權重因子(time weighting factor) Φ 用以加權兩時段的流場變數，以考慮時間積分的效應。至於動量方程式的部分，在具自由液面的天然河川中，對流項 (convective terms)為主要影響項，故使用顯式有限解析法將其離散化，以全微分的形式表示之(葉等, 1996): DQ ∂ Q ∂Z − gAS = −βQ ( ) − gA ∂x Dt ∂x A. f. (3-7). 其中， ∂Q DQ Q ∂Q = +β ( ) ∂t Dt A ∂x. (3-8). 上式中， D / Dt 表示對時間之全微分計算；而全微分的特性線速度軌跡為: dx Q =β =u A dt. ( β = 1). (3-9). 由式(3-8)可看出移流項部分為非線性，其移流速度本身即為一待求之流場變數，因此式(3-2)中之特性線運動軌跡無法直接計算求得。為解決此種問題，必須先將式(3-8)局部線性化，為各計算點取其特徵. 12.

(25) 速度 u * ，用以代替原有之移流速度 u 。當然，對於任一計算時段，各計算元素中所取的特徵速度要能適當表示移流速度之大小。，如此將使得(3-7)式變成一線性非齊次雙曲線型偏微分方程式，透過線性化的過程，並藉由特性速度計算其特性線之運動軌跡，如此，式(3-2)便能用以求解式(3-8)中移流項的部分。由於在一時段內的變化量跟總變化量比起來要小很多，因此可以前一時段的值當作已知值代替。為進一步考慮式(3-2)速度項 u 和式 (3-7)源項之影響在時間上的變化，分別將(3-7)、(3-9)式對時間積分，得 η2 ⎡ ⎤ ∂Z ∂ ⎡Q ⎤ + gAS f ⎥dt Qη 2 − Qη1 = − ∫ ⎢ βQ ⎢ ⎥ + gA η1 ∂x ∂x ⎣ A ⎦ ⎣ ⎦. η2. X η 2 − X η1 = − ∫ udt. (3-10) (3-11). η1. 其中，η1和η2分別代表特性曲線與前一時段之空間軸的交點和目前時段之空間軸的交點位置。動量方程式在交錯格網上之差分式可表示如下： # # ⎡ ⎛ ⎛ Q ⎞ (n +1) ⎛ Q ⎞ (n+1) ⎜⎜ ⎟ ⎢ −⎜ ⎟ ⎜ ⎝ A ⎠r Qin+1 − Qηn1 ⎢ # ⎝ A ⎠s (n +1) + ⎢ΦQi β⎜ Δt nΔx ⎜ ⎢ ⎜ ⎢ ⎝ ⎣. # # ⎡ ⎛ (n +1) − Z s(n +1) ( n +1)# ⎜ Z r ⎢ΦgAi ⎜ nΔx ⎢⎣ ⎝. [ ΦgAi(n+1) S (fin+1) #. #. ⎞ ⎛ ⎛ Q ⎞ n ⎛ Q ⎞ n ⎞⎤ ⎟ ⎜ ⎜ ⎟ − ⎜ ⎟ ⎟⎥ ⎟ ⎜ ⎝ A ⎠ r ⎝ A ⎠ s ⎟⎥ ⎟ + (1 − Φ)Qηn1 β ⎜ ⎟⎥ + nΔx ⎟ ⎜ ⎟⎥ ⎟ ⎜ ⎟⎥ ⎝ ⎠⎦ ⎠. n n ⎤ ⎞ ⎟ + (1 − Φ ) gAηn1 ⎛⎜ Z r − Z s ⎞⎟⎥ + ⎜ nΔx ⎟ ⎟ ⎠⎥⎦ ⎝ ⎠. + (1 − Φ ) gAηn1 S nfi. ]= 0. (3-12). 式中，考慮流場可能存在亞臨界流與超臨界流同時發生之混合流況，抑或兩種流況於同一位置交替發生，因此式(3-11)等式右邊非齊次項部分之數值離算方式，可根據各計算點之流場變數判斷，當水流為亞臨界流時，採用中央差分法將其離散化，而當超臨界流時，採用上風. 13.

(26) 法(前項或後項差分法) ，以考慮明渠流中擾動波傳遞的方向。即若為亞臨界流(-1<Fr<1) 則r=i+1，s=i-1，n=2；若為向下游傳遞之超臨界流(Fr＞1) 則r=i，s=i-1，n=1，若為向上游傳遞之超臨界流(Fr＜-1) 則r=i+1，s=i，n=1，其中，上標為 (n + 1)# （包含參數 Q,. Q , A, Z , S f ） A. 代表在目前時段下疊代過程中新計算得之變數值;Φ為權重係數； Δt 為時間間距； Δx 為空間間距；下標為 η1 (包含參數Q、A)以前一時段格點上的值線性內差而得，如 Qηn1 以如下形式表示： Qηn1 = Qin −1. xη1 − xi −1 Δx. + Qin−−11. xi − xη1. (3-13). Δx. 其中，欲計算 Qηn1 的量，首先要決定特性線速度，以確定 η1 點的位置，根據實際模擬經驗，每次疊代運算後所得之最新流速 u i(−n1+1) 適合用來 #. 做為特徵速度。而 xη1 以如下形式近似： xη1 = xi − u i(−n1+1) Δt #. (3-14). 3.4 外部邊界條件邊界條件為數值模式中極為重要的一環，除了所設定的條件與個數必須符合物理意義之外，設定的方式亦必須兼顧模式的功能性與適用性。在一維計算流場中，上、下游水理邊界條件的型態通常有以下三種: (1)以流量歷線為邊界條件，即給定流量與時間之函數關係式，其形式為 Q=f(t)。 (2)以率定曲線為邊界條件，即給定流量與水深之函數關係式，其形式為 Q=f(h)。 (3)以水位歷線為邊界條件，即給定水深與時間之函數關係式，其形式為 h=f(t)。. 14.

(27) 除了給定的邊界條件外，求解邊界上其他未知物理量之處理技巧係透過特性法 (characteristics method) 來估算之 (Garcia-Navarro& Saviron,1992)；本研究即引用特性法求解邊界未知物理量。沿流線方向擾動波傳遞的兩個特性方程式，可以用來近似或推估邊界上未知物理量；此兩種特性方程式，分別存在速度向量(u+c)及(u-c)，其中u為水流流速，c為水面擾動波之波速(wave celerity)。當流場型態發生變化時，擾動波傳遞的方向(向上游或下游)可能會隨著改變。舉例來說，當流速方向指向下游為正值時，特性速度在亞臨界流況下(u+c) 將會指向下游，而(u-c)會指向上游，如圖3-1所示；而在超臨界流況下，特性速度(u+c)及(u-c)皆指向下游，如圖3-2所示。因此，針對不同的流況，需要不同的方式來處理邊界條件的問題。根據Garcia-Navarro & Saviron, (1992)之邊界條件概述，在亞臨界流況下，先就上游邊界而言，由於特性速度(u+c)在時間-空間座標內所形成之特性線落在模擬流場區域的外面，因而無法據此獲得沿著特性線所傳遞的流場訊息，為彌補不足，可以根據實際流況給定上游邊界條件，例如根據流量歷線Q=f(t)做為邊界條件。至於未知的上游邊界水深，則可以根據以下的相對應之方程式(compabillity equation)估算之： gA ∂A DQ ⎫ DA ⎧ = ⎨ gA( S 0 − S f + ⎬÷μ h)− Dt ⎩ B ∂x Dt ⎭. μ = − β u − ( β 2 − β )u 2 + g. A B. (3-15) (3-16). 當上式中之 β 等於 1 時，特性軌跡為 dx A =u− g =u−c dt B. (3-17). 就下游邊界而言，由於特性線速度(u-c)所形成之運動軌跡落在模擬流場區域的外面，因此也必須根據實際流況給定一個下游邊界條 15.

(28) 件，例如水深歷線h=f(t)，至於未知物理量Q，則可利用以下特性方程式估算之： DQ gA ∂A DA = gA( S 0 − S f + h) − μ Dt B ∂x Dt. μ = − β u + ( β 2 − β )u 2 + g. (3-18). A B. (3-19). 當上式中之 β 等於 1 時，特性軌跡為 dx =u+ dt. g. A =u+c B. (3-20). 特性線的詳細推導請見附錄一。若水流流況為超臨界流時，就上游邊界而言，由於特性速度(u+c) 及(u-c)在時間-空間座標內所形成之特性線皆落在模擬流場區域的外面，因而無法據此獲得沿著特性線所傳遞的流場訊息，可以根據實際流況給定上游兩個邊界條件，例如根據流量歷線 Q=f(t)及水深歷線 h=f(t) 做為邊界條件；因兩條特性線皆往下游方向傳遞，故於下游邊界處，以外插特性線之作法求值，如圖 3-3 所示，再仿效亞臨界流時之交錯網格作法，以 C-之特性線利用連續方程式取水深，而以 C+之特性線利用動量方程式求取流量。利用此兩條特性線及其相對應之方程式(3-15)及(3-18)，求解兩個變數。而水流流況為超亞臨界混合流時，將超亞臨界混合流流場的超臨界流與亞臨界流區域分開來討論之。首先，對於超臨界流區域，此兩個邊界條件必須同上述的給定在上游邊界，以符合流場特性波傳遞的方向，然而水躍發生的地方，依據內部邊界條件，可得到一個超臨界流區域的下游邊界水位。再者，對於亞臨界流區域而言，此一內部邊界條件，可當作其上游邊界條件，所以僅須在下游邊界給定一水位條件即可滿足。. 16.

(29) 綜合以上所述，邊界點之計算乃透過特性法之觀念，由兩條特性線交匯而得，沿特性線之擾動波波速可表示成. dx = u ± c ，當亞臨界流 dt. 時， u < c 即兩擾動波波速具相反方向，同交於邊界點，此時位於計算區外之一條特性線須由人為給定此一邊界條件；當臨界流時， u = c 即其中一擾動波波速為零，另一速率由水流決定為正向（向下游）或負向（向上游），若水流方向為正，則在上游需給定一邊界條件，而下游則不須給定，直接由前一時刻之下游邊界變量給定即可。當超臨界流時， u > c 即兩擾動波波速具相同方向，若水流方向為正，則上游處須同時給定兩邊界條件，而下游則不須給定。而當超亞臨界混合流時，則須採用超臨界流及亞臨界流之邊界條件，即上游處同時給定兩個邊界條件，下游處給定一個邊界條件即可，詳細之外部開放邊界條件處理原則如附錄二所示。. 3.5 內部邊界條件當超亞臨界混合流即所謂之跨臨界流發生時，將超亞臨界混合流流場分為超臨界流與亞臨界流區域分開討論之。首先，對於超臨界流區域，此兩個邊界條件必須給定在上游，以符合流場特性波傳遞的方向。而在超臨界流的下游端，即水躍發生的地方，會因為上、下游水深的共軛特性自然形成一個限制條件，對整個流場而言，此限制條件為一內部邊界條件，它的存在使得超臨界流區域多一個邊界條件，導致數值計算發生 over-determined 的現象。為解決內部邊界條件所造成的問題，可利用內部相鄰計算點的水位高程，透過外插的方式來得到超臨界流區域的下游邊界水深。至於超臨界流下游邊界的位置，則定義在福祿數大於 1 的計算格點，並且其下一個計算格點的福祿數小於 1，如此即可在反覆的疊代計算過程中，根據福祿數的大小來判斷水 17.

(30) 躍發生的位置。此外，對亞臨界流區域而言，上游邊界處已存在一個邊界條件，僅需在下游邊界處給定一水位條件，即可滿足亞臨界流流況條件。. 3.6 交錯格網為避免數值計算發生所謂”checkerboard”之震盪現象，因此採用交錯格網來進行差分，如圖 3-4 所示，Q 代表流量，A 代表通水斷面積。所謂的交錯格點即為在奇數點利用動量守恆方程式計算流量，在偶數點則利用連續方程式計算通水斷面積；也就是每一格點只利用一條特性線軌跡，再利用內差的方法求得偶數點上的流量與奇數點上的通水斷面積。在交錯格網中 Qη1 即為動量控制方程式中的 Qξn ，以格網中的特性線為例， Qη1 是由前一時刻的 Qin−1 與 Qin 內差而得，而 Qη1 的位置是由 Qin−1 格點上的流速 U in−1 乘上一開始所假設的△t 值來決定。再將 Qη1 與 Qin +1 做一比較，相差範圍在誤差容許內，則允許跳出計算。而在. 每個時間間距中，反覆的疊代計算為必要的，透過此種過程，可以不斷的修正移流項部分因為局部線性化，以及源項與離散化連續方程式因為顯式求解所造成的誤差，而讓流場變數同時趨於收斂。. 3.7 穩定性分析顯式有限解析法和隱式法不同處為須滿足穩定條件，並非為無條件穩定。一維水流之偏微分方程式係屬雙曲線型，以 Courant, Friedrichs, and Lewy(CFL)於1967年對雙曲線型偏微分方程式進行數值解之穩定性條件分析為發韌。Dai(1994)證明當條件. u Δt Δx. ≤ 1. 成立. 時，顯式有限解析法為穩定的數值方法。有關顯式有限解析法應用在淺水波方程式之穩定性分析，詳見陳氏(1998)。 18.

(31) 3.8 模式演算之流程模式之演算流程如圖 3-5 所示，茲扼要敘述如下：步驟一：給定初始條件初始條件包括：渠道坡度、可蘭數、時間權重、交匯角度、交匯格點、Δt、總模擬時間、模擬類型、疊代次數、初始水深與流量、上游入流流量、下游水深。根據初始條件，可得到每一斷面之初始通水斷面積、濕周與河面寬。步驟二：判斷Δt 值由初始條件所給定之流量、通水斷面積及計算所得之速度，經由 CFL 穩定條件判斷之副程式判定初始條件所給定之Δt 是否過大，若給定之Δt 值比經由 CFL 穩定條件下所得之Δt 值小，則採用給定之值，若大於經由 CFL 穩定條件下所得之Δt，則以 CFL 穩定條件下所得之Δt 取代之。步驟三：邊界條件本文中，若流況為超臨界流，上游給定條件為流量及水深，下游無須給定條件，若流況為超亞混合流，上游給定條件為流量及水深，下游給定條件為水深，配合連續方程式與動量方程式，以一線性組合係數 μ 將兩方程式合而為一，最後形式如下所示： ⎡ DA ⎤ ⎡ DQ ⎤ ⎢⎣ Dt ⎥⎦ + μ ⎢⎣ Dt ⎥⎦ + F = 0. (3-21). 根據所給定之條件即可求得未知數（Q 或 A）。步驟四：內部斷面計算配合連續方程式、動量方程式與交錯格網，分別依據所給定之初始條件求得未知數(Q 或 A)，此時在交錯格網中的 Qη1 位置是由 U in−−11 乘上Δt 所得，每一斷面上之速度項則是由流量除以通水斷面積而得。 19.

(32) 步驟五：CFL 穩定判斷由每一斷面計算所得之新的流量、通水斷面積及流速，計算 CFL 穩定條件下所得之Δt 值，透過有如步驟二的方法給定新的Δt 值，此一步驟將影響內插值的位置(譬如 Qη1 的位置將會隨著改變)，但是若採用新的Δt 值，則所計算得到的新流量與通水斷面積也將回復到舊的流量與通水斷面積（也就是前一時刻的值），亦即此一Δt 的過程必須重新計算。步驟六：內部穩定條件由計算所得之新的流量與通水斷面積，將與前一內部疊代的流量與通水斷面積作一比較，若兩者的絕對誤差值小於容許的誤差範圍，則內部疊代即可停止；若非如此，則給予適當的內部疊代次數。步驟七：收斂條件將此一時刻與前一時刻的流量和通水斷面積做一比較，若兩者的絕對誤差小於容許的範圍，則判定此一渠道水流已達到穩定狀態。步驟八：判斷模擬類型經由給定的模擬類型：穩定流況、非穩定流況，決定模擬結束條件，若是模擬穩定流況，則程式將模擬至整個河道達到穩定狀態為止，若是模擬非穩定流況，則視所欲模擬之時間而定。. 20.

(33) 第四章模式驗證本章節首先針對定量流方面進行模擬及演算，並分別以單一坡度矩形渠道、單一坡度矩形漸縮及漸擴渠道，進行超臨界流流場的分析，說明模式於單一坡度矩形渠道中的模擬，接著再以單一坡度矩形渠道、單一坡度矩形漸擴渠道、陡坡接緩坡之矩形渠道、緩坡接陡坡再接緩坡之矩形渠道，模擬超亞臨界混合流流場，接著再針對超臨界流流況以及超亞臨界混合流流況作變量流分析，並透過與美國陸軍工兵團HEC-RAS模式的比較，說明本模式模擬之成果與適用性。 HEC-RAS模式為西元1967年由美國陸軍工兵團水文工程中心 (The Hydrologic Engineering Center, U.S. Army Corps of Engineers)首創，為一視窗化整合型軟體系統，此系統包含使用者圖形介面 (Graphical User Interface,簡稱GUI)、水理分析單元(定量流水面剖線、變量流模擬、輸砂/動床演算)、資料儲存與管理、圖表製作與資料彙整等功能，功能強大且下載方便，在工程界具有一定的公信力，故將模擬結果與HEC-RAS模式進行比較。. 4.1 超臨界定量流流場分析 4.1.1 單一坡度矩形渠道今考慮三種不同坡降的矩形渠道，其坡降依序分別為0.5、0.1、 0.01，渠道曼寧粗糙係數皆假設為混凝土材質(n≒0.016)，數值格點皆為51點，數模試驗採用可蘭數(Cr)為0.7，而其邊界條件方面，則皆控制上游水深分別為0.05公尺、0.1公尺、0.15公尺進行模擬。當渠道坡降為0.5時，假設渠道長為10公尺，寬為1公尺，上下游端底床高度差5公尺之矩形渠道，而其上游流量則依照曼寧公式，分別依序推估. 21.

(34) 而得Q = 0.281 cms、0.843cms、1.571cms，而其渠道沿程水面線如圖 4-1至圖4-3所示。當渠道坡降為0.1時，假設渠道長為10公尺，寬為1 公尺，上下游端底床高度差1公尺之矩形渠道，而其上游流量則依照曼寧公式，分別依序推估而得Q = 0.1259 cms、0.3771cms、0.7026cms，而其渠道沿程水面線如圖4-4至圖4-6所示。而當渠道坡降為0.01時，假設渠道長為100公尺，寬為1公尺，上下游端底床高度差1公尺之矩形渠道，而其上游流量依照曼寧公式，分別依序推估而得Q = 0.0398 cms、0.1192cms、0.2221cms，其渠道沿程水面線如圖4-7至圖4-9所示。上述三種坡降之矩形渠道數值演算，皆與HEC-RAS之計算結果相當一致。 4.1.2 單一坡度矩形漸縮及漸擴渠道依上述三種不同坡降的矩形渠道來做漸縮及漸擴的模擬以及檢定，並沿用其數值模擬假設，首先就矩形漸縮渠道而言，依據三種不同坡降的渠道特性，皆分別自上游端寬為1.0公尺處，每一斷面漸縮寬度0.01公尺，即每一斷面寬自左右兩端各縮減0.005公尺，漸縮至最下游斷面寬0.5公尺，假設其上游水深同上皆固定為0.05公尺、0.1公尺、0.15公尺，因其上游端通水斷面積不變，故仍延用上述之不同坡降邊界流量；而其渠道坡降為0.5之渠道沿程水面線如圖4-10至圖4-12 所示，坡降為0.1之渠道沿程水面線如圖4-13至圖4-15所示，坡降為 0.01之渠道沿程水面線如圖4-16至圖4-18所示。接下來就矩形漸擴渠道而言，由於渠道完全超臨界流現象不易存在於底床坡降為0.01之矩形漸擴渠道處，故依據底床坡降為0.5及0.1 之渠道特性來作模擬及分析，其漸擴斷面自上游端寬為0.5公尺處，每一斷面漸擴寬度0.01公尺，即每一斷面寬自左右兩端各加寬0.005 公尺，漸擴至最下游斷面寬1.0公尺處，其上游水深同上皆固定為0.05. 22.

(35) 公尺、0.1公尺、0.15公尺，而上游流量則依照曼寧公式推求，在底床坡降為0.5時，依序推估而得Q = 0.1328 cms、0.3804cms、0.684cms，其渠道沿程水面線如圖4-19至圖4-21所示，而在底床坡降為0.1時，流量依序推估而得Q = 0.0594 cms、0.1701cms、0.3059cms，其渠道沿程水面線如圖4-22至圖4-24所示。根據上述之各種超臨界流水理演算案例，不論是不同坡降的矩形渠道或漸縮及漸擴矩形渠道問題，EFA水理模式皆與HEC-RAS之計算結果符合，有相當不錯的超臨界流流況模擬結果。. 4.2超亞臨界定量混合流流場分析 4.2.1 單一坡度矩形渠道本節將以底床坡降為0.01之矩形渠道來作模擬及分析，其渠道特性及數值模擬假設同上述方式，並沿用上述之超臨界流上游邊界條件，即控制上游水深分別為0.05公尺、0.1公尺、0.15公尺，而其上游流量則亦同上述依照曼寧公式推估而得，並假設其所對照之下游水深皆為2倍上游水深及4倍上游水深進行模擬。當上游端水深控制為0.05 公尺時，其不同下游水深之渠道沿程水面線如圖4-25及圖4-26所示；上游端水深控制為0.1公尺時，其不同下游水深之渠道沿程水面線如圖4-27及圖4-28所示；而上游端水深控制為0.15公尺時，其不同下游水深之渠道沿程水面線如圖4-29及圖4-30所示。上述三種固定上游端水深且給定不同下游水深之矩形渠道數值演算，皆與HEC-RAS之計算結果相當一致。 4.2.2 單一坡度矩形漸擴渠道接下來就單一坡度矩形漸擴渠道而言，同上述依據底床坡降為 0.01之渠道來作模擬及分析，其渠道漸擴斷面特性及數值模擬假設亦. 23.

(36) 同上述方式，上游水位同上皆固定為0.05公尺、0.1公尺、0.15公尺而其上游流量則依照曼寧公式，依序推估而得 Q = 0.0188 cms 、 0.0538cms、0.0967cms，而其所對照之下游水位皆假設為2倍上游水深及4倍上游水深進行模擬。當上游端水深控制為0.05公尺時，其不同下游水深之渠道沿程水面線如圖4-31及圖4-32所示；上游端水深控制為0.1公尺時，其不同下游水深之渠道沿程水面線如圖4-33及圖4-34 所示；而上游端水深控制為0.15公尺時，其不同下游水深之渠道沿程水面線如圖4-35及圖4-36所示。上述三種固定上游端水深且給定不同下游水深之矩形漸擴渠道數值演算，由EFA與HEC-RAS之數值演算結果觀察發現，在模擬混合流流場時，HEC-RAS容易出現數值震盪的現象，而EFA在無須藉助任何人工黏滯係數的情況下，亦不會有數值震盪的現象發生，仍有不錯的模擬解果。 4.2.3 陡坡接緩坡之矩形渠道接下來就陡坡接緩坡之矩形渠道進行模擬，而其數值模擬假設為一矩形渠道總長100公尺，渠寬1公尺，上下游端底床高度差7.5公尺之矩形渠道，渠道曼寧粗糙係數皆假設為混凝土材質(n≒0.016)，數值格點為101點( Δx =1m)，數模試驗採用可蘭數(Cr)為0.7；此渠道分為兩個不同底床坡降，依序分別為slope=0.01、0.005，而其長度範圍各為50公尺；而其上游水位同上皆固定為0.05公尺、0.1公尺、0.15公尺，上游流量則依照曼寧公式，分別依序推估而得Q = 0.0398 cms、 0.1192cms、0.2221cms，並假設其所對照之下游水位皆為2倍上游水深及4倍上游水深進行模擬。當上游端水深控制為0.05公尺時，其不同下游水深之渠道沿程水面線如圖4-37 及圖4-38所示；當上游端水深控制為0.1公尺時，其不同下游水深之渠道沿程水面線如圖4-39及圖 4-40所示；當上游端水深控制為0.15公尺時，其不同下游水深之渠道. 24.

(37) 沿程水面線如圖4-41及圖4-42所示。上述三種固定上游端水深且給定不同下游水深之陡坡接緩坡之矩形渠道數值演算，由 EFA 與 HEC-RAS之數值演算結果觀察發現，在模擬混合流流場時，HEC-RAS 之渠道沿程水面線於靠近下游端附近，容易出現數值震盪的現象，而 EFA亦不會有數值震盪的現象發生。 4.2.4 緩坡接陡坡再接緩坡之矩形渠道(Meselhe, 1994) 本研究模擬Meselhe於1994年所做的數值試驗，在Meselhe的數值試驗中，考慮水躍發生區域附近，水深與流速梯度急遽變大所造成的非靜水壓分佈的特性，需要加入人工黏滯性(artificial viscosity)來消除或抑制水躍附近所產生的高頻率震盪，並且需要一些額外的參數來設定其大小，過大的值將會導致不必要的數值延散(dispersion)，對於不同類型的流況，這些參數必須個別的加以率定，以便能在獲得抑制數值震盪效果的同時，不至於發生太大的數值擴散現象，以求得最佳的模擬結果，然而從應用的角度來看，這卻降低了模式的實用性，並造成使用上的不便。而其數值模擬假設為一矩形渠道總長1000公尺，渠寬8公尺，曼寧粗糙係數(n≒0.015)；此渠道分為三個不同底床坡降，依序分別為slope=0.002、0.009、0.001，而其長度範圍如圖4-43所示，分別為300公尺、300公尺、400公尺，其數值格點為51點( Δx =20m)，數模試驗採用可蘭數(Cr)為0.7；在邊界條件方面假設上游流量為20 cms，而上游起始水位為4.5m，下游水位為2.5m進行模擬。其數值模擬後與Meselhe所做的數值模擬及利用直接步推法推估之定量流水面線值(Henderson, 1966)之渠道沿程水面線，如圖4-44所示，水躍的位置與共軛水深皆與Meselhe所做的數值模擬試驗相吻合。. 4.3 超臨界變量流流場分析. 25.

(38) 在天然河川中，流況會隨著時間而改變，水流並非真正定量流況，故本文亦發展出可模擬超臨界變量流流場以及超亞臨界變量混合流流場之數值模式，以期能適用於模擬天然河川；本文假設於起始模擬時刻河道水流為定量流狀態，以作為變量流之起始演算流況。而其模擬之假設為一直線渠道長10公尺，寬為1公尺，上下游端底床高度差1公尺之矩形渠道(slope=0.1)，渠道曼寧粗糙係數假設為混凝土材質 (n≒0.016)，數值格點為51點，數模試驗採用可蘭數(Cr)為0.7，由於超臨界變量流上游端流況易漸變成亞臨界流流況，故縮短其模擬時距為50分鐘；其數值模擬之邊界條件設置，在上游輸入一三角形流量歷線，如圖 4-45 所示，而上游水深歷線則由流量歷線之曼寧公式 (Manning formula)推估而得，如圖4-46所示，進行模擬超臨界變量流流場。其不同時刻之渠道沿程水面線，如圖4-47~4-51所示；由於 HEC-RAS之數值變量流模擬演算，在下游端需設置邊界條件，且其超臨界變量流模擬後之流量與起始條件差異過大，並不符合質量守恆條件，故在此未採納比較，而EFA在模擬超臨界變量流演算時，亦無此種問題產生。. 4.4超亞臨界變量混合流流場分析而在模擬超亞臨界變量混合流流場中，其假設於起始模擬時刻河道水流狀態亦同上為為定量流狀態，以作為變量流之起始演算流況。其模擬之假設為一直線渠道長100公尺，寬為1公尺，上下游端底床高度差1公尺之矩形渠道(slope=0.01)，渠道曼寧粗糙係數假設為混凝土材質(n≒0.016)，數值格點為51點，數模試驗採用可蘭數(Cr)為0.7，模擬時間為6小時；而其數值模擬時之邊界條件設置，在上游端輸入一三角形流量歷線，如圖4-52所示，而上游端水深歷線則由流量歷線. 26.

(39) 之曼寧公式(Manning formula)概估而得，如圖4-53所示，下游水深歷線假設為一亞臨界流流況之三角形水位歷線，如圖4-54所示，進行模擬超亞臨界變量流流場。其不同時刻之渠道沿程水面線，如圖 4-55~4-61所示；由EFA與HEC-RAS之數值演算結果觀察發現，在模擬超亞臨界變量混合流流場時，EFA之模擬後水位略大於HEC-RAS 之數值演算結果，且其發生水躍現象的位置也較HEC-RAS之數值演算結果明顯。. 4.5 結語根據上述諸多水理演算案例，不論是超臨界流或超臨界流與亞臨界流同時存在或交替發生之複雜混合流流場問題，EFA水理模式在無須藉助任何人工黏滯係數的情況下，在模擬水躍流場，或是其他水位高程突變，水面曲線急遽變化的流場時，亦不會有數值震盪的現象發生，仍有相當不錯的模擬結果，因此初步驗證了模式的模擬能力與廣泛的適用性。. 27.

(40) 第五章草嶺潭河段模擬 5.1 草嶺潭概述草嶺地區位於雲林縣古坑鄉草嶺村阿里山山脈之西支嶺線地區，地形屬於山岳地帶，坡度甚陡，地形發育受地質構造與岩層性質影響甚劇。清水溪為本地區主要河流，由東而西橫切草嶺地區。於九二一集集大地震時，發生歷年來最大規模之崩山事件，崩塌之大量土石阻斷清水溪，河道阻斷長度約5公里，高度約40~50公尺，並於崩塌區上游形成一堰塞湖。其崩塌土體估計約有1億2千萬立方公尺，其中僅約2千5百萬立方公尺之崩坍土堆積於坡趾處之清水溪河谷中；另外大部分之坍方塊體(約1億立方公尺)則堆積於清水溪左岸，到交山高程也因而增高了約六、七十公尺。有關九二一地震造成之草嶺堰塞湖災害之基本資料如下(經濟部水利署，民國91年)： a、崩塌堆積堰塞湖上游高度約50公尺(標高約540公尺) b、集水區面積162平方公里 c、崩塌面積約400公頃 d、崩塌長度(沿河上下游)約5公里 e、崩塌土方約1億2千萬立方公尺 f、總積水容量約4,300萬立方公尺. 5.2 流域概況及模擬範圍草嶺潭崩塌地位於濁水溪支流清水溪之上游，距清水溪與濁水溪合流口二十七公里處，亦即雲林縣古坑鄉草嶺村西方二公里處，如圖 5-1所示，草嶺潭天然壩壩趾至濁水溪匯流口全長約30公里，流域面積約有四百二十一平方公里，流經嘉義、雲林、南投等縣，支流有田仔溪、加走寮溪、阿里山溪等，於林內觸口山下匯於濁水溪。全流域 28.

(41) 幾全為山區，夏季濕熱，冬季較冷，夏季西南季風盛行，清水溪流域因處於迎風面，降雨較多，冬季則多屬晴朗天氣之乾季，降雨量少。本研究將以88年清水溪之上游河段草嶺潭地區，坡度變化較陡之河段做為案例模擬對象，模擬範圍從距草嶺潭溢流口0.9公里至4.5公里處，總長約3.6公里，其中，在2.1公里~3.9公里段為921地震後所形成第二水潭之下游出口端，根據經濟部水利署水利規劃試驗所之堰塞湖引致災害防治對策之研究指出，此渠段形成亞臨界流與超臨界流之流況，故以此河段進行現地之超亞臨界混合流水理分析。. 5.3 模擬所需資料 5.3.1 地文資料本模式模擬草嶺潭地區清水溪之上游段河道，其斷面資料係利用民國89年全河道大斷面資料進行模擬，模擬範圍從距草嶺潭溢流口 0.9公里至4.5公里處，共13個河道斷面，斷面間距皆為300公尺，總長約3.6公里，並因本模擬河段兩處之河道底床坡降變化劇烈，故依模擬所需分別在距草嶺潭溢流口1.8公里至2.4公里處及距草嶺潭溢流口 3.0公里至3.6公里處內插4個斷面，故共計17個河道斷面進行模擬，其模擬河道之3-D示意圖如圖5-2所示。 5.3.2 河床糙度(曼寧n值) 一維河川水理模式因不考慮河道沖淤的情況，故假設河槽之底床不會改變，為定床河渠。其曼寧n值採用與經濟部水利署水利規劃試驗所之堰塞湖引致災害防治對策研究相同值，約在0.045~0.06之間。 5.3.3 水文資料 1.上、下游邊界條件本模式所需之邊界條件，均參考「堰塞湖引致災害防治對策之研. 29.

(42) 究」中各項研究成果，作為包括上游河段之入流歷線、水位歷線以及下游水位歷線為邊界條件進行演算。其中上游入流歷線，採用該研究中利用FLAWAV模式模擬方案 C1(200年洪水潰壩延時0.5小時下，潰口流量歷線)，模擬草嶺潭天然壩潰口出流歷線之分析結果，如圖5-3所示，本文取其模擬時間前4小時之流量歷線，時間間距為每3分鐘一筆流量資料，作為本模式之上游入流歷線，如圖5-4所示。而其上下游水位歷線，乃採用該研究中利用Sobek模式模擬草嶺潭堰塞湖後之抬升水位進行模擬，如圖5-5、 5-6所示，模擬水位歷線之時間為4小時，時間間距為每3分鐘一筆水位資料。 2.起始條件在進行洪流演算時，必須輸入起始流量與起始水位資料，方能進行模擬，通常在天然河川中，水流並非真正定量流況，但由於無足夠的資料能提供演算之起始條件，故本文假設於起始模擬時刻河道水流為定量流狀態作為變量流之起始演算流況。首先根據模擬時刻之上、下游邊界條件為起始流量及水位資料，利用線性內差或是迴水演算求得每一斷面的起始水位值。再來將起始演算求得之水位值，利用一維顯式有限解析法的定量流模式進行演算，以求得每一斷面之起始水位。本模式進行演算時採取可蘭數與時間權重兩參數，此兩參數值必須小於1。根據陳(1998)之模擬試驗成果，當可蘭數大時，其數值震盪可被消除；而時間權重因子越小時，數值震盪越顯著，其建議採用範圍介於0.7~1.0之間，故本文將依據其試驗結果採用0.7。. 5.4 模擬結果與討論. 30.

(43) 以200年重現期距，潰壩延時0.5小時之入流洪峰流量13,070cms 為初始流量，漸減至流量1,333cms，每隔9分鐘之渠道沿程水面線，如圖5-7~5-17所示；由上述圖中及圖5-2草嶺潭地區清水溪之上游河段 3-D示意圖可看出上游端之坡度較為陡峭，其高程從草嶺潭溢流口0.9 公里處的534.94公尺驟降至距離草嶺潭溢流口1.2公里處的482.94公尺，其坡降約為17.3 %，且其河道斷面寬度也較溢流口之狹窄斷面寬闊，河道斷面為漸擴狀態，極易形成超臨界流轉變為亞臨界流之流況。而圖中距上游端0~1200公尺處，恰為草嶺潭之第二工區水潭，其蓄水或滯洪之現象甚為明顯，而圖之1200~2700公尺處，為一連續性之陡坡，其平均坡降為9.45 %，且其通水斷面積較小，流速加快為超臨界流之流況，尤其是圖中距上游端2400 ~ 2700公尺處，為一地形起伏變化更劇烈之陡坡，其高程從433.71公尺處驟降達到370.45公尺，其平均坡降更高達21.09 %，與圖中距上游端1200~2100公尺處形成連續陡坡，為本模擬案例河床渠槽坡度最為陡峭之處。而其下游端為2700~3600公尺處，坡度較為平緩，其平均坡度約為3 %，河道流況為亞臨界流之狀態，故本章模擬之清水溪河段，其流況由上游至下游呈現超臨界流交替為亞臨界流，在流經一連續陡坡漸變為超臨界流再交替為亞臨界流之複雜水理河段，由模擬結果可看出由於下游端坡度較緩，流場會受到迴水的影響而水面線逐漸抬升，而其流速也會相對的遞減，其流況為亞臨界流，而上游端從距草嶺潭溢流口0.9公里處至草嶺潭溢流口2.1公里處，即圖5-7~5-17所示的 0~1200公尺處，恰為一陡坡河道銜接一蓄水槽，但因為此蓄水槽之容積甚大，依據「草嶺潭堰塞湖水工模型動床試驗及動床水理分析專題報告」中指出，其蓄水可達3百多萬立方公尺水量，故相較之下其發. 31.

(44) 生水躍現象並不明顯，而本模擬案例的中、下游河段，即圖中5-7~ 5-17 所示的1200~3600公尺處，為一陡坡接更陡坡再接緩坡之渠道坡降特性，其較陡坡銜接緩坡交會處，如圖中所示距上游端2700公尺處，其水位相較於下一個河道斷面之水位淺，而其下一個河道斷面，即圖中距上游端3000公尺處，其水位有明顯的抬升現象，故應存在有局部水躍現象之水理流況。自然河道地形急遽變化，河道斷面的束縮、突擴或是坡度的突升、突降易造成複雜的流況，水面受到超臨界流震波的影響，使模擬結果極易發散不易達到收斂之穩態解，尤其本模擬案例之崩塌區地形變化起伏甚大，與一般天然河道截然不同，更增加了模擬之難度，然草嶺潭崩塌地皆為地形陡峭、變化極大區域，故於水理流況上皆屬於超臨界流、近超臨界流及亞臨界流相互交錯區域，一般水理模式並不具備此種功能，而本文驗證所採用之HEC-RAS模式適用於河床坡度小於10％之緩量變流，而本模擬案例地形變化甚劇，多屬急變量流，所以亦無法模擬，但由本章節之模擬結果顯示，本模式已可模擬於自然河道之複雜河段，並可適用於超流界流與亞臨界流並存之水理情況。. 32.

(45) 第六章結論與建議 6.1結論本研究以一維定床超臨界流與超亞臨界混合流為主要研究對象，模式以亞臨界流演算架構為基礎，尋求各種可能解決跨臨界流問題之方法。在處理超亞臨界混合流流況時，加入一內部邊界條件，利用內部相鄰計算點的水位高程，透過外插的方式來得到超臨界流區域的下游邊界水深，而超臨界流下游邊界的位置，則定義在福祿數大於 1的計算格點，並且使其下一個計算格點的福祿數小於1，如此即可在反覆的疊代計算過程中，根據福祿數的大小來判斷水躍發生的位置。本文針對單一坡降矩形渠道、不同坡降矩形渠道、清水溪上游段河道現地模擬問題進行模擬分析，並透過與美國陸軍工兵團 HEC-RAS的比較，說明本模式模擬之成果與適用性，並獲得結論如下： 1. 本研究發展之EFA水理模式可用以模擬定量流及變量流流況，對於超臨界流問題或是超臨界流與亞臨界流同時存在或交替發生之複雜流場，根據諸多水理演算與美國陸軍工兵團水文工程中心所研發之一維水理演算模式HEC-RAS相比較下有不錯的結果，驗證了EFA法對於具自由水面流場之模擬能力與廣泛的適用性。 2. EFA水理模式在無須藉助任何人工黏滯係數的情況下，在模擬水躍流場，或是其他水位突變，即水面線急遽變化的流場時，亦不會有數值震盪的現象發生，仍有相當不錯的模擬結果。 3. 本研究發現EFA法對於移流特性顯著之流場，尤其在流線曲率大的流況下仍具適用性，本文並實際模擬應用於清水溪草嶺潭地區之上游段河道，對於極其陡峻之河床渠道坡降仍可模擬，顯示出. 33.