顯式有限解析法模式結合預測­修正數值法於 超亞臨界混合流之研究

顯式有限解析法模式結合預測­修正數值法於 超亞臨界混合流之研究

學生:鍾仁凱 指導教授:葉克家 國立交通大學土木工程學系

摘要

本研究延續 Hsu and Yeh (1996)之一維顯式有限解析法(explicit finite analytic method，簡稱 EFA)模式，發展出適於超臨界流與超亞臨 界混合流流況之數值模式。EFA 求解之特點，乃在於求解水流動量方 程式時，以特性法觀念解得其中變量(流量與通水斷面積)之局部解析 解，並且遵守可蘭穩定性條件；邊界處理方面則透過水流之連續方程 式與動量方程式，利用特性法觀念求解邊界處之變量；而在超亞臨界 混合流況之內部邊界，根據福祿數來判斷水躍發生之位置；內部相鄰 計算點的水位高程，先透過預測-修正數值法(MacCormack scheme)來 得到超臨界流區域的流量及通水面積，接著再將預測-修正法所得流量 及通水面積以特性法觀念求解超臨界流區域之下游邊界水深。此方式 可有效減少地形驟變對數值計算之影響，並可解決因流況改變所引起 數學上屬於奇異點(singularity)之問題。本文針對實驗室水槽及實際河 川陡緩坡交替之情況進行數模與評估。

關鍵詞：顯式有限解析法、預測­修正數值法、地形驟變、超亞臨界混合 流

 

 

 

 

 

 

 

 

 

 

Study on Mixed Supercritical and Subcritical Flows Using Explicit Finite Analytic Method and

MacCormack Hybrid Scheme

Student : Ren-Kai Jhong 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 MacCormack scheme and method of characteristics, by utilizing the water surface elevations of the interior neighboring computational points. The mixed supercritical and subcritical flow fields in laboratory flumes and natural rivers will be simulated and evaluated by the proposed model.

Keywords: explicit finite analytic method, MacCormack Scheme, sudden elevation transition, mixed supercritical and subcritical flow.

 

目錄

摘要 ... I Abstract ... II 目錄 ... III 表目錄 ... V 圖目錄 ... VI 符號說明 ... X

第一章 緒 論 ... 1

1.1 前言 ... 1

1.2 文獻回顧 ... 1

1.2.1 常用定床數模之回顧 ... 1

1.2.2 超臨界流 ... 4

1.2.3 超亞臨界混合流 ... 4

1.3 論文架構 ... 5

第二章 理論基礎 ... 7

2.1 基本假設(de Saint Venant 之假設) ... 7

2.2 控制方程式 ... 7

2.3 顯式有限解析法簡介 ... 8

2.4 內部邊界（超亞臨界混合流）數值方法簡介 ... 9

2.4.1 預測-修正法 ... 9

2.4.2 Lax Method ... 10

2.4.3 Lax-Wendroff Method ... 11

2.4.4 TVD Scheme ... 11

第三章 數值方法 ... 13

3.1 齊次雙曲線型方程式 ... 13

3.2 混合型方程式 ... 14

3.3 水理控制方程式之離散化 ... 14

3.4 邊界條件(參考楊 1996 顯示有限解析法於超亞臨界混合流之研究 國立交通大學碩士論文 ) ... 17

3.5 交匯區邊界條件(內部邊界條件)... 20

3.5.1 處理方式 ... 21

3.5.2 交匯區連續方程式、動量方程式處理 ... 22

3.6 超亞臨界混合流處理方法(內部邊界條件) ... 22

3.7 穩定性分析 ... 23

3.8 交錯格網 ... 23

3.9 模式演算之流程 ... 24

第四章 模式檢定 ... 26

4.1 非線性問題（Burger Equation）分析 ... 26

4.2 緩坡接陡坡再接緩坡之矩形渠道(Meselhe, 1994) ... 27

4.3 淡水河系颱洪檢定案例模擬 ... 27

4.3.1 淡水河系概況 ... 27

4.3.2 模擬範圍概述 ... 28

4.3.3 模式演算所需資料 ... 28

4.3.4 馬瑪莎颱洪檢定結果討論 ... 29

4.4 濁水溪颱洪檢定案例模擬 ... 31

4.4.1 濁水溪流域概況 ... 31

4.4.2 模擬範圍概述 ... 32

4.4.3 模式演算所需資料 ... 32

4.4.4 敏督利颱洪檢定結果討論 ... 33

4.5 小結 ... 34

第五章 模式驗證與結果討論 ... 35

5.1 淡水河系颱洪驗證結果討論 ... 35

5.1.1 艾利颱洪模擬結果討論 ... 35

5.1.2 海棠颱洪模擬結果討論 ... 37

5.1.3 泰利颱洪模擬結果討論 ... 39

5.2 濁水溪主流颱洪模擬結果討論 ... 41

5.2.1 龍王颱洪模擬結果討論 ... 41

5.3 結果討論 ... 42

第六章 結論與建議 ... 43

6.1 結論 ... 43

6.2 建議 ... 44

參考文獻 ... 45

附錄一 特性線之推導 ... 99

附錄二 開放邊界條件處理原則 ... 100

表目錄

表 4-1 淡水河流域曼寧 N 值(檢定前) ... 91

表 4-2 淡水河流域率定後曼寧 n 值(1/6) ... 92

表 4-3 淡水河流域率定後曼寧 n 值(2/6) ... 93

表 4-4 淡水河流域率定後曼寧 n 值(3/6) ... 94

表 4-5 淡水河流域率定後曼寧 n 值(4/6) ... 95

表 4-6 淡水河流域率定後曼寧 n 值(5/6) ... 96

表 4-7 淡水河流域率定後曼寧 n 值(6/6) ... 97

表 4-8 濁水溪本流曼寧糙度 n 值表 ... 98

圖目錄

圖 3-1 亞臨界流況下之特性速度 ... 48

圖 3-2 超臨界流況下之特性速度 ... 48

圖 3-3 合流示意圖... 49

圖 3-4 分流示意圖... 49

圖 3-5 一維特性曲線交錯格點示意圖 ... 50

圖 3-6 模式流程圖... 51

圖 4-1 初始速度的不連續面分佈(Burgers 方程式案例) ... 52

圖 4-2 初始速度的不連續面分佈(Burgers 方程式案例) ... 52

圖 4-3 初始速度的不連續面分佈(Burgers 方程式案例) ... 52

圖 4-4 初始速度的不連續面分佈(Burgers 方程式案例) ... 53

圖 4-5 初始速度的不連續面分佈(Burgers 方程式案例) ... 53

圖 4-6 Meselhe (1994)數值模擬矩型渠道示意圖 ... 53

圖 4-7 EFA & Meselhe (1994)渠道沿程水面線比較圖 ... 54

圖 4-8 淡水河流域圖 ... 54

圖 4-9 瑪莎颱洪上游邊界入流歷線 ... 55

圖 4-10 瑪莎颱洪下游邊界河口潮位歷線 ... 55

圖 4-11 艾利颱洪上游邊界入流歷線 ... 55

圖 4-12 艾利颱洪下游邊界河口潮位歷線 ... 56

圖 4-13 海棠颱洪上流邊界入流歷線 ... 56

圖 4-14 海棠颱洪下流邊界河口潮位歷線 ... 56

圖 4-15 泰利颱洪上游邊界入流歷線 ... 57

圖 4-16 泰利颱洪下游邊界河口潮位歷線 ... 57

圖 4-17 馬莎颱洪洪峰時刻鳶山堰下游範圍水面線 ... 57

圖 4-18 馬莎颱洪新海大橋水位歷線比較 ... 58

圖 4-19 馬莎颱洪新海大橋實測值與模擬值比較圖 ... 58

圖 4-20 馬莎颱洪台北橋水位歷線比較 ... 58

圖 4-21 馬莎颱洪台北橋實測值與模擬值比較圖 ... 59

圖 4-22 馬莎颱洪獅子頭水位歷線比較 ... 59

圖 4-23 馬莎颱洪獅子頭實測值與模擬值比較圖 ... 59

圖 4-24 馬莎颱洪土地公鼻水位歷線比較 ... 60

圖 4-25 馬莎颱洪土地公鼻實測值與模擬值比較圖 ... 60

圖 4-26 馬莎颱洪中正橋水位歷線比較 ... 60

圖 4-27 馬莎颱洪中正橋實測值與模擬值比較圖 ... 61

圖 4-28 馬莎颱洪寶橋水位歷線比較 ... 61

圖 4-30 馬莎颱洪五堵水位歷線比較 ... 62

圖 4-31 馬莎颱洪五堵實測值與模擬值比較圖 ... 62

圖 4-32 濁水溪流域圖... 63

圖 4-33 敏督利颱洪上游邊界入流歷線 ... 64

圖 4-34 敏督利颱洪下游邊界水位歷線 ... 64

圖 4-35 龍王颱洪上游邊界入流歷線 ... 64

圖 4-36 龍王颱洪下游邊界水位歷線 ... 65

圖 4-37 敏督利颱洪彰雲大橋水位歷線比較 ... 65

圖 4-38 敏督利颱洪彰雲大橋實測值與模擬值比較圖 ... 65

圖 4-39 敏督利颱洪溪州大橋水位歷線比較 ... 66

圖 4-40 敏督利颱洪溪州大橋實測值與模擬值比較圖 ... 66

圖 5-1 大漢溪、淡水河縱剖面示意圖 ... 66

圖 5-2 艾利颱洪洪峰時刻鳶山堰下游範圍水面線 ... 67

圖 5-3 艾利颱洪新海大橋水位歷線比較 ... 67

圖 5-4 艾利颱洪新海大橋實測值與模擬值比較圖 ... 67

圖 5-5 艾利颱洪台北橋水位歷線比較 ... 68

圖 5-6 艾利颱洪台北橋實測值與模擬值比較圖 ... 68

圖 5-7 艾利颱洪獅子頭水位歷線比較 ... 68

圖 5-8 艾利颱洪獅子頭實測值與模擬值比較圖 ... 69

圖 5-9 艾利颱洪土地公鼻水位歷線比較 ... 69

圖 5-10 艾利颱洪土地公鼻實測值與模擬值比較圖 ... 69

圖 5-11 艾利颱洪中正橋水位歷線比較 ... 70

圖 5-12 艾利颱洪中正橋實測值與模擬值比較圖 ... 70

圖 5-13 艾利颱洪寶橋水位歷線比較 ... 70

圖 5-14 艾利颱洪寶橋實測值與模擬值比較圖 ... 71

圖 5-15 艾利颱洪五堵水位歷線比較 ... 71

圖 5-16 艾利颱洪五堵實測值與模擬值比較圖 ... 71

圖 5-17 海棠颱洪洪峰時刻鳶山堰下游範圍水面線 ... 72

圖 5-18 海棠颱洪新海大橋水位歷線比較 ... 72

圖 5-19 海棠颱洪新海大橋實測值與模擬值比較圖 ... 72

圖 5-20 海棠颱洪台北橋水位歷線比較 ... 73

圖 5-21 海棠颱洪台北橋實測值與模擬值比較圖 ... 73

圖 5-22 海棠颱洪獅子頭水位歷線比較 ... 73

圖 5-23 海棠颱洪獅子頭實測值與模擬值比較圖 ... 74

圖 5-24 海棠颱洪土地公鼻水位歷線比較 ... 74

圖 5-25 海棠颱洪土地公鼻實測值與模擬值比較圖 ... 74

圖 5-26 海棠颱洪中正橋水位歷線比較 ... 75

圖 5-27 海棠颱洪中正橋實測值與模擬值比較圖 ... 75

圖 5-28 海棠颱洪寶橋水位歷線比較 ... 75

圖 5-29 海棠颱洪寶橋實測值與模擬值比較圖 ... 76

圖 5-30 海棠颱洪五堵水位歷線比較 ... 76

圖 5-31 海棠颱洪五堵實測值與模擬值比較圖 ... 76

圖 5-32 泰利颱洪洪峰時刻鳶山堰下游範圍水面線 ... 77

圖 5-33 泰利颱洪新海大橋水位歷線比較 ... 77

圖 5-34 泰利颱洪新海大橋實測值與模擬值比較圖 ... 77

圖 5-35 泰利颱洪台北橋水位歷線比較 ... 78

圖 5-36 泰利颱洪台北橋實測值與模擬值比較圖 ... 78

圖 5-37 泰利颱洪獅子頭水位歷線比較 ... 78

圖 5-38 泰利颱洪獅子頭實測值與模擬值比較圖 ... 79

圖 5-39 泰利颱洪土地公鼻水位歷線比較 ... 79

圖 5-40 泰利颱洪土地公鼻實測值與模擬值比較圖 ... 79

圖 5-41 泰利颱洪中正橋水位歷線比較 ... 80

圖 5-42 泰利颱洪中正橋實測值與模擬值比較圖 ... 80

圖 5-43 泰利颱洪寶橋水位歷線比較 ... 80

圖 5-44 泰利颱洪寶橋實測值與模擬值比較圖 ... 81

圖 5-45 泰利颱洪五堵水位歷線比較 ... 81

圖 5-46 泰利颱洪五堵實測值與模擬值比較圖 ... 81

圖 5-47 濁水溪縱剖面示意圖 ... 82

圖 5-48 龍王颱洪彰雲大橋水位歷線比較 ... 82

圖 5-49 龍王颱洪彰雲大橋實測值與模擬值比較圖 ... 82

圖 5-50 龍王颱洪溪州大橋水位歷線比較 ... 83

圖 5-51 龍王颱洪溪州大橋實測值與模擬值比較圖 ... 83

圖 5-52 海棠颱洪新海大橋絕對誤差分析圖 ... 83

圖 5-53 海棠颱洪台北橋絕對誤差分析圖 ... 84

圖 5-54 海棠颱洪獅子頭絕對誤差分析圖 ... 84

圖 5-55 海棠颱洪土地公鼻絕對誤差分析圖 ... 84

圖 5-56 海棠颱洪中正橋絕對誤差分析圖 ... 85

圖 5-57 海棠颱洪寶橋絕對誤差分析圖 ... 85

圖 5-58 海棠颱洪五堵絕對誤差分析圖 ... 85

圖 5-59 艾利颱洪新海大橋絕對誤差分析圖 ... 86

圖 5-60 艾利颱洪台北橋絕對誤差分析圖 ... 86

圖 5-61 艾利颱洪獅子頭絕對誤差分析圖 ... 86

圖 5-62 艾利颱洪土地公鼻絕對誤差分析圖 ... 87

圖 5-63 艾利颱洪中正橋絕對誤差分析圖 ... 87

圖 5-64 艾利颱洪寶橋絕對誤差分析圖 ... 87

圖 5-66 泰利颱洪新海大橋絕對誤差分析圖 ... 88

圖 5-67 泰利颱洪台北橋絕對誤差分析圖 ... 88

圖 5-68 泰利颱洪獅子頭絕對誤差分析圖 ... 89

圖 5-69 泰利颱洪土地公鼻絕對誤差分析圖 ... 89

圖 5-70 泰利颱洪中正橋絕對誤差分析圖 ... 89

圖 5-71 泰利颱洪寶橋絕對誤差分析圖 ... 90

圖 5-72 泰利颱洪五堵絕對誤差分析圖 ... 90

圖 5-73 龍王颱洪各水位站模擬水位絕對誤差分析圖 ... 90

符號說明

A﹕通水斷面積 B﹕渠道寬度 C﹕單位換算常數

C

r﹕可蘭數 c﹕水面波速

F

r﹕福祿數 g﹕重力加速度 h﹕水深

n﹕曼寧糙度係數 Q﹕流量

R﹕水力半徑

S

0﹕底床坡降

S

f ﹕摩擦坡降 t﹕時間

U﹕特徵速度 u﹕斷面平均流速

x﹕沿渠道中心線的距離 Y﹕渠道橫斷面距離 Z﹕水面高程

β﹕動量校正係數 Δt﹕時間間距 Δx﹕空間間距

μ﹕邊界條件線性組合係數 Φ﹕時間權重係數

第一章 緒 論

1.1 前言

台灣地區由於經濟發展迅速，在可利用的土地資源日益匱乏下，往 中、上游流域發展已是不可避免的趨勢。但由於台灣地區河川坡陡流急，

且中、上游河段或水工結構物下游常為超臨界流或超亞混合之流況，故充 分瞭解其水理現象對河川整治及水工結構物的設計將有極大的幫助。然而 在天然河川或人工渠道坡度陡峻處、斷面束縮段、渠道底床局部明顯突起 或河道中設置水工結構物，如溢流堰或動態調節閘門下之高速水流，常有 臨界流況的發生，此現象在物理上屬於不穩定（instability）現象，在數學 上是屬於奇異點（singularity）的問題，所以造成數值模擬上相當大的困難。

EFA 水理模式採用顯式有限解析法來求解雙曲線型淺水波水流動量方 程式。在各計算元素(element)中，可用以求得動量方程式中移流項部分之 局部解析解。模式進行規則渠道演算時，可處理陡、緩坡交界之流況，於 水躍流場或其他水位突變，不致有數值震盪現象發生；於網狀河流亦可進 行主支流之模擬。但應用於現場特殊複雜混合流況案例時，仍有需突破的 瓶頸，如計算通水面積小於零等，故本文將尋找能以更有效處理超亞臨界 混合流的方法，解決分析現場渠道超臨界流或超亞混合流流況時之水理現 象。

1.2 文獻回顧

1.2.1 常用定床數模之回顧

以下就較常用之定床水理演算模式介紹：

1. HEC-RAS 模式

HEC-RAS（Hydrologic Engineering Center’s River Analysis System）模 式係美國陸軍工兵團水文工程中心所發展之一維水理演算模式，HEC-RAS 模式除了納入原HEC-2 模式之定量流模擬演算外，更加入其他演算功能。

HEC-RAS 為一維水面線演算模式，適用於河床坡度小於 10%之定量 緩變流，可處理亞臨界流、超臨界流及混合流之水面剖線演算，亦具有模 擬變量流的功能。

2. UNET 模式

UNET（One-Dimensional Unsteady Flow Through a Full Network of Open Channels）主要為應用於網狀系統渠道之一維變量流水理演算模式，

由美國陸軍工兵團水文工程中心支援並推廣本模式。UNET 模式可以應用 於一維變量流之水流模擬，但不支援超臨界流流況的模擬。除了可以模擬 單一渠道，而共最大特色為可模擬多分支渠道，更可模擬完整之網狀系統 渠道。

3. FLDWAV 模式

FLDWAV（Flood Wave Routing Model）模式為美國國家氣象局(U.S.

National Weather Service, NWS) 所 發 展 之 一 維 渠 道 洪 水 演 算 模 式 。 FLDWAV 模式發展最主要的目的即結合 DWOPER 模式與 NWS DAMBRK 模式的特性，其中DWOPER 模式可處理多分支渠道模擬，NWS DAMBRK 模式可處理潰壩水流、堤防溢流以及超臨界流模擬。FLDWAV 模式結合了 上述兩個模式的優點，可處理一維變量流、超臨界流流況、多分支渠道系 統以及堤防溢流所形成的洪氾問題，並且可模擬潰壩水流、堤防溢流、即 時洪水變化以及水流流經橋、閘門、溢洪道等流況。

4. MIKE 11 模式

MIKE 11 為丹麥水工試驗所(Danish Hydraulic Institute)所發展之商業 模式，此模式為一完全視窗化之一維演算模式，其中流體動力（HD，

Hydrodynamics）模組為 MIKE 11 模式之水理演算模組，可提供其他模組 所需的水理演算資料。該模組可以模擬一維變量流流況，並且可以處理超 臨界流流況以及水流流經堰、涵洞以及不規則之結構物，並可處理多分支

或是網路系統之渠道與河川。

5. NewC 法一維河川水理模式

 

Kutija & Hewett於2002年發展NewC法，針對河川在陡緩坡交界的問題 進行處理。NewC法是在交錯網格的配置下以有限差分法進行Saint Venant 方程式的離散，具有無條件穩定的特點。NewC演算式差分離散控制方程 式交錯格點配置，是在整數節點處配置流量變數 Q，在半節點處配置水位 變數 h；通水斷面 A、自由水面寬度 b與輸水容量係數 K 等待定係數均 視為水位h的函數，這些配置於半節點處，並採用時間權重方式使方程式 能穩定收斂。另外，在程式撰寫上透過雙掃法(double sweep)以達到高效率 求解。

6. SOBEK 模式

SOBEK模式為荷蘭WL | Delft Hydraulics公司與其他荷蘭顧問公司所 共同發展之一整合性軟體系統，該模式整合河川、人工渠道與下水道系統 之模擬。SOBEK模式中水流（water flow）模組係模擬一維變量流，可適 用於規則斷面或是天然河道斷面，並可模擬水流流經橋墩、孔口、堰、涵 洞及抽水站等水工結構物，亦可應用於網狀系統渠道。

7. 淡水河模式

淡水河現有三種預報模式，分別是台大顏清連教授所研發淡水河洪水 預報系統、十河局的 REFOR 模式以及台大蔡丁貴教授所發展的洪水預報 與淹水預警系統等三種洪水預報模式。以下就台大顏清連教授所發展之洪 水預報模式進行說明：預報模式具有預測降雨、降雨逕流、河川演算與水 庫操作等功能，可進行即時演算與事前演算之功能；該模式並以水位站觀 測值進行水位回饋演算，以提升模式預報之精確性，唯須於河道多處設立 水位站，方能達此功能；該模式僅適於緩坡亞臨界流況。 

1.2.2 超臨界流

在超臨界流數值模擬方面， Jimenez & Chaudhry (1988)利用MacCor- mack（預測-修正數值法）法計算且加入人工遲滯項處理超臨界明渠流；

Bhallamudi & Chaudry (1992)、Rahman & Chaudry (1997)針對超臨界流束縮 段之水理情況亦利用MacCormack有限差分法進行數值模擬計算。國內研究 方面，溫至剛(1986)對於預測-修正法在明渠流況適用性有詳細之探討；許 銘熙等人(1993、1994)先後完成渠道中超臨界流之數值模擬，並針對渠道 中急變流況進行各項數值模擬；而古孟晃(2003)也利用總變量削減法(TVD) 法完成渠道中超臨界流況之數值模擬。

在實驗方面，Coles & Shintaku在1943年的超臨界流實驗中，進行超臨 界流通過束縮段之水槽試驗；而Rouse在1951年亦做了超臨界流通過擴張 段之水槽試驗。Ippen在1951年時，也進行渠道束縮段超臨界流水工模型試 驗，並以四象限圖解法決定斜震波之水深與震波角。

1.2.3 超亞臨界混合流

流場型態因為渠道坡度及寬度的改變、流量突變或存在水工結構物等 因素而發生劇烈變化，造成超亞臨界混合流場且導致震波(shock)的形成與 傳遞，尤其在非定量流中，超臨界流與亞臨界流之界面會隨時移動，更增 加了模擬的困難度，因此模式常需以此類流場為模擬對象，以驗證其實用 性。依據Cunge et al. (1980)的分類，文獻中求解具震波水流的計算方式可 分為三類：(1)震波擬合(shock fitting)法：此法需要內部邊界條件用以決定 震波的位置及震波的流場特性，然而如此的作法往往造成較複雜的計算過 程；(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年之水躍試驗中，利用下 射式閘門產生超臨界流後，調整尾水堰以控制尾水高度，使其產生亞臨界 流之流況。

1.3 論文架構

由上述文獻回顧可知，超臨界流與超亞臨界混合流之處理已有許多研 究 ， 故 本 研 究 引 用 前 人 預 測- 修 正 數 值 方 法 ( M a c C o r m a c k method)之觀念，作為一維顯式有限解析法處理超臨界流與超亞臨界混合流 之方法，以避免流況劇烈變化所以引起的數值震盪。

本研究係延續Hsu & Yeh (1996)一維矩形渠道之顯式有限解析法數模 成果，更新一維Q-A模式演算超臨界流及超亞混合流之數值方法。由於台

灣地區河川坡度變化劇烈，每逢洪水來臨將造成河床劇烈的沖淤變化，河 道斷面不規則，且局部超臨界流及超亞混合流發生的情況甚多，因此研發 Q-A模式以模擬河道中洪水波與水位高程變化有其必要性。

本文之主要內容包括：利用連續方程式與動量方程式，在滿足de Saint Venant 之基本假設下，應用EFA方法，將動量方程式轉換成全微分式，再 利用沿特性線積分的觀念，求解動量方程式，配合適當的斷面處理與差分 式，求解連續方程式，且為避免受數值震盪的影響，採用交錯格網求解，

在處理超亞臨界混合流況時，第一步利用預測-修正數值法求取超亞臨界混 合流區域的各斷面通水面積及流量，第二步再利用第一步所得通水面積、

流量代入顯式有限解析法，得到真正所需超亞臨界混合流區域之流量及通 水面積，並根據福祿數的大小來判斷水躍發生的位置。

為驗證本數值模式的適用性與準確性，利用求解Burger Equation及 Meselhe (1994)實驗室案例，說明本模式在模擬數值上的不連續及超亞臨界 混合流之成果；在模式應用方面以淡水河及濁水溪為例，並與商用模式 SOBEK比較，以探討本模式在加入處理超亞臨界混合流數值方法後，是否 能準確模擬各水位站水位歷線的變化，另外，選取流況交換頻繁區域的水 面線，判別是否合理及水面線有無明顯數值震盪所產生的突升或突降。

第二章 理論基礎

渠道中之水流遵守質量守恆及動量守恆定律，因此由數學式描述之控 制方程式在滿足 de Saint Venant 之基本假設下，隨所選擇的應變數而有不 同型式的一維明渠流數學模式。採用之變數為水深與流速時，多用於特徵 方程式之求解中；而本文採用流量及通水斷面積為應變數。

2.1 基本假設(de Saint Venant 之假設) 1.流速均勻分佈:

流速均勻分佈在通水面積上，即每一個通水斷面積僅存在一個流速，

此即一維水流。換言之，通水斷面的水位線為一水平線，在此水平線上的 各點水位相等。

2.靜水壓分佈:

假設渠道中水流之垂向流線曲率很小而且忽略其垂直加速度，因此水 深方向速度梯度為零，以及水面坡度甚小，可忽略垂向加速度，則靜水壓 分佈的假設成立。

3.渠道定量流摩擦損失估計:

渠底摩擦與紊流效應對水流所造成的損失，可以定量流摩擦律估算。

4.平均底床坡度甚小：

當平均渠底坡度甚小時，重力沿渠道所造成的分力將會很小，甚至可 忽略不計，亦即水深可以垂向水面與渠底高程差表示。

5.忽略柯氏力及風力的影響:

通常在大範圍之海洋或是湖泊才會考慮到柯氏力及風力的影響，或是 緯度相差很大的時候才會考慮，在小範圍的區域以重力、靜水壓、摩擦力 為主。

2.2 控制方程式

對於不可壓縮水流之控制方程式，包括水流連續方程式與水流動量方

程式，為如下形式。

1.連續方程式

= 0

∂ ± + ∂

∂

∂

q

l

x Q t

A

(2-1) 式中，A 為通水斷面積；Q 為流量；t 為時間；x 為沿渠道中心線的距離；

q

l為單位渠長之支流側流量，

+ q

_l為分流之處理，

− q

_l為合流之處理。

2.動量方程式

= 0

±

∂ + + ∂

⎥ ⎦

⎢ ⎤

⎣

⎡ ⎟

⎠

⎜ ⎞

⎝

⎛

∂ + ∂

∂

∂

l l

f

q u

x gAS gA Z A

Q Q x t

Q

β (2-2) 式中，β 為動量校正係數；Q/A=u 為斷面平均流速；g 為重力加速度；Z 為水面高程；

q

_l為單位渠寬之支流側流量，

+ q

_l為分流之處理，

− q

_l為合 流之處理；

u

_l為支流在主流方向的速度分量；而摩擦坡降可表為:

3 4 2

2

3 4 2 2

2

R C

n u u R A C

n Q

S

_f

= Q =

(2-3) 其中，R 為水力半徑；C 為單位換算常數；n 為曼寧糙度係數。

2.3 顯式有限解析法簡介

本研究採用之有限解析法為美國佛羅里達大學陳景仁 (Chen C. J.) 教 授所創，可分為隱式及顯式有限解析法兩種。隱式有限解析法已廣泛運用 於計算流體力學而得到頗佳的成果 (Chen & Chen，1984，Chen & Choi，

1990)，此方法的特色為:

(1)對於不規則的幾何邊界，仍可在卡式座標使用結構性(structured)格網求 解。

(2)由於在個別元素體的數值離散化(discretization)係採用局部解析(local analytic solution)來近似，因此可以把數值演算的捨入誤差降到最低。

(3)數值解之穩定性極佳，屬於無條件穩定之數值模式。

由 於 隱 式 有 限 解 析 法 之 解 法 係 針 對 橢 圓 或 拋 物 線 型(elliptic or parabolic)之偏微分方程式求解，如用於雙曲線型(hyperbolic)偏微分方程(明

渠水流即為一例)並不太適合，因此乃有顯式有限解析法的發展及提出 (Dai，1994)，但僅求解無自由表面之 Navier-Stokes 方程式。在求解對流 傳輸方程式中，對流項以特性法觀念解得式中變量之局部解析解，且使依 時變量再透過適當給定之初始條件而求得，此即為顯式法之特色。顯式有 限解析法之穩定性，和一般常見的顯式數值方法一樣，受可蘭(Courant) 數 小於或等於1 之限制，雖不如隱式法般屬於無條件穩定，但解法卻較為簡 單，故在應用上仍有其優點。

2.4 內部邊界（超亞臨界混合流）數值方法簡介

當超亞臨界混合流發生時，在超臨界流的下游端，即水躍發生的地方，

會因為上、下游水深的共軛特性自然形成一個限制條件，對整個流場而 言，此限制條件為一內部邊界條件，它的存在使得超臨界流區域多一個邊 界條件，導致數值計算發散的現象，若要避免發散的情況發生，勢必要在 發生跨臨界流區域採用不同的數值方法來做處理，其所需要的特性應要能 處理非線性問題及具備有數值震盪極小及數值消散程度小。

以下列舉四種處理非線性問題的數值方法：

2.4.1 預測-修正法

在第一章文獻回顧中，可發現許多學者利用此數值方法處理超亞臨界 混合流，主要特色為可處理陡變的大梯度問題而造成的振盪性誤差。預測 -修正法也是本文內部邊界處理所選用的數值方法。

此方法包括預測-修改兩步驟，可以在空間和時間達到二階精度；主要 觀念是將一個時間間距(Δ t)內的計算步驟分為兩個步驟，即預測(predictor step)與修正(corrector step)。

假設存在一系統為

( ) H ( U x t )

x U F t

U = , ,

∂ + ∂

∂

∂

U為任意物理性質

預測

( )

[ ( ) ( ) ]

iⁿ

n i n i n

i n

i l

i

F F F tH

x U t

U − − − − + Δ

Δ

− Δ

= 1

ε ₊₁

1 2

ε ε ₋₁ (2-4) 修正

(

( )

) [

^{( )}

( )

^{( )}

( )

^{( )}

]

i^{( )l} l

i l

i l

i l

i n i n

i

t H

F F

x F U t

U

U 1 2 1 2

2 2

1

1 1

1

+ − + − + Δ

Δ

− Δ +

=

_{+ −}

+ ε ε ε (2-5)

穩定條件

_max

≤ 1 Δ

Δ u x

t

(2-6) 註：

ε = 1

或

ε = 0

從上面式子中可觀察到，在計算下一個時間步驟 (time step) 上的第i 點時會用到前一個時間的 i - 1、i 、i + 1 三點的值，所以只能計算網格中 間域的點，而計算邊界上的點則是需配合其他數值方法求解。在求預測值 與修正值時，前項差分與後項差分的順序是可變換的(Chaudhry,1993)，然 因修正步驟較預測步驟重要，因此修正步驟之方向應與擾動波傳遞的方向 一致，所以本文在處理超亞臨界混合流時，以後項差分求預測值，以前項 差分求修正值。

2.4.2 Lax Method

此方法為顯式法，在空間可達到二階精度，時間達到一階精度；主要 觀念是將非線性系統直接使用泰勒級數展開，化簡為差分式。

假設存在一系統為

( ) ^H ( ^{U x t} )

x U F t

U = , ,

∂ + ∂

∂

∂

U為任意物理性質 差分式

n i n i n i n i n

i

H

x F F t

u

u =

Δ + − Δ

−

_{+ −}

+

2

1 1 1

(2-7) 穩定條件

≤ 1 Δ t

(2-8)

此方法最大特點就是差分式簡易，且不需要額外考慮擾動波傳遞的方 向；但在求解現場複雜問題時，容易產生數值發散。

2.4.3 Lax-Wendroff Method

此方法為顯式法，可以在空間與時間皆可達到二階精度；主要觀念是 將非線性系統先使用泰勒級數展開，再利用Jacobian轉換的概念求取差分 式。

假設存在一系統為

( ) ^H ( ^{U x t} )

x U F t

U = , ,

∂ + ∂

∂

∂

U為任意物理性質 差分式

( )

^{( )}_{( )}

[ ( )( ) ( )( ) ]

iⁿ

n i n i n i n i n i n i n i n i n

i n i n

i n

i

u u F F u u F F H

x F t

x F u t

u + − − + − +

Δ + Δ Δ −

− Δ

=

_{+ − + + − −}

+

1 1

1 2 1

2 1

1 1

2 4

(2-9)

穩定條件

max

≤ 1 Δ

Δ u x

t

(2-10) 此方法最大特點就是對Lax method 作一修正，以求高穩定性；但在求 解的過程中仍會有數值震盪的狀況產生，但模擬結果較Lax method 為佳

^。

2.4.4 TVD Scheme

此方法可以捕捉複雜震波問題，求解不連續解相當精確，且可抑制非 物理性的振盪。但在程式的撰寫上有一定的困難度。

假設存在一系統為

( ) H ( U x t )

x U F t

U = , ,

∂ + ∂

∂

∂

U為任意物理性質 差分式

n i n

i n i n

i n

i

h h H

x u t

u ⎥ +

⎦

⎢ ⎤

⎣

⎡ − Δ

− Δ

=

+ −

+

2 1 2 1

1 (2-11)

h 為數值通量

( )

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡ + − Δ

= + + +

+

n i i n i n i n

i

u F

F h

2 1 2 1 1

2 1 2

1

α

(2-12)

( )

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡ + − Δ

= − − −

−

n i i n i n i n

i

u F

F h

2 1 2 1 1

2 1 2

1

α

(2-13)

其中α 為特徵速度

⎪ ⎪

⎭

⎪⎪ ⎬

⎫

⎪ ⎪

⎩

⎪⎪ ⎨

⎧

≠

− Δ

−

= Δ

=

+ + +

+

+

0

0

2 1 1

1

2 1

2 1

n n i

i n i

n i n i

n i n

i

i

if u

u u

F F

u if

u

α

；

⎪ ⎪

⎭

⎪⎪ ⎬

⎫

⎪ ⎪

⎩

⎪⎪ ⎨

⎧

≠

− Δ

−

= Δ

=

− −

−

−

−

0

0

2 1 1

1

2 1

2 1

n n i

i n i

n i n i

n i n

i

i

if u

u u

F F

u if

u

α

(2-14)

n i n i n

i

u u

u = −

Δ

₊

+ 1

2

1 ； ⁿ _iⁿ _iⁿ

i

u u

u

₁

2

1 −

−

= −

Δ

(2-15) 穩定條件

0 1

2 1

≤ Δ

≤ Δ

+

x

i

t

α (2-16)

第三章 數值方法

本研究採用顯式有限解析法求解水理控制方程式，該法為美國佛羅里 達大學陳景仁教授所創，此方法的特色為對於個別計算元素的數值離散，

係採用局部解析解來近似，因此可以把數值演算的捨入誤差低，且數值穩 定性極佳，屬於無條件穩定之數值模式。然而隱式有限解析法較適用於橢 圓或拋物線型偏微分方程式，後來顯式有限解析法的發展及提出(Dai 1994)，乃針對於雙曲線型偏微分方程式(河川水流即為一例)。以下各節中 將陸續介紹顯式有限解析法之求解方法與數值特性。

3.1 齊次雙曲線型方程式

EFA水理模式採用顯式有限解析法來求解雙曲線型淺水波水流動量方 程式。在各計算元素(element)中，該數值方法可用以求得動量方程式中移 流項部份之局部解析解。其基本觀念可利用以下之一維度一階線性齊次雙 曲線型微分方程式來加以說明：

= 0 +

_x

t

u

φ

φ (3-1) 其中，u為x-方向之速度，在各計算元素中被視為常數。當啟始條件

( ) x , 0 ϕ ( x )

φ =

被適當的給定時，可求得式(3-1)之解析解如下：

( x

₀

, Δ t ) =

ϕ

( x

₀

− u Δ t )

φ (3-2) 其中，

( x

₀

)

為待求點之座標；

x = x

₀

− u Δ t

則定義一條特性線的運動軌跡。

該特性線從與起始平面之交點

φ ( ) x , 0

(圖3-5中所標示之η 點) 出發，在經

1

過

Δ t

時間後，正好移流通過點

( x

₀

)

。從物理觀點說明式(3-2) 的意義為，對 於式(3-1)之純移流問題，計算點

( x

₀

)

於某一時刻

( t

ⁿ⁺¹

)

之物理量，等於通過 該點之特性線與啟始平面

( t

ⁿ⁺¹

)

的交點(η 點)所具有之物理量。至於

1

η 點的

1

座標

( x

_η1

)

，則可根據移流速度

u

藉由特性軌跡來加以推求(

x

_η1

= x

₀

− u Δ t

)。

此外，對於不同正負號的組合，根據上風(upwind)的特性，其相依區域 (domain of dependence)將會自動的調整。

對於非線性的移流方程式，如明渠流之動量方程式的移流項部分：

= 0 +

_x

t

uu

u

(3-3) 因具有複雜的非線性特性，使得解析解形式並不存在。為解決此一問題，

並獲得如式(3-2) 的解析解形式，可透過局部線性化的方法，將定義特性 線軌跡之移流速度

u

以特徵速度代替之。儘管移流速度本身隨著時間和空 間而改變，仍可假設此特徵速度在一計算時間間距內為常數，用以代表此 時間間距內平均的移流速度。至此，所剩之問題即在於如何選擇適當的特 徵速度，此將留待至方程式離散化時一併說明。

3.2 混合型方程式

對於非齊次的混合型方程式，如明渠流之動量方程式

g x

t

+ u

φ

= F

φ (3-4) 等式右邊可利用已知的物理量以顯式的方法直接計算之，所得結果視為源 項(source term)，可直接加入式(3-2)中得：

( x , Δ t ) =

ϕ

( x

₀

− u Δ t ) + F

_g

Δ t

φ (3-5) 以上的計算方法即為Dai所提出之顯式有限解析法。本文則將此數值方法應 用於具自由液面之流場中，並在模式中引入疊代計算的流程， 藉以修正 特徵速度與源項。

3.3 水理控制方程式之離散化

在EFA 水理模式中，動量方程式的移流項部分係採用顯式有限解析法 來加以離散化，其餘的項次，以及連續方程式，則採用有限差分法的方式 處理之。首先，考慮水流連續方程式，為保留連續方程式的守恆特性，因 此以控制體積法(control volume)將其離散化，空間上採用中央差分式與時 間上採用前向差分式(陳, 1998)。連續方程式在此交錯格網上之差分式如 下：

( ) ( ) ( )

0 )

1

(

^{1 1}

1 1 1 1

1

1 ^{# # #}

=

⎥ ⎤

⎢ ⎡

⎟⎟

⎞

⎜⎜

⎛ − ±

Φ

−

⎟ +

⎜ ⎞

⎛ − ±

Φ

− +

₊⁺ ₋^{+ +} _{+ −}

+

A Q Q Q Q Q Q

A

_iⁿ _iⁿ _iⁿ _iⁿ _lⁿ _iⁿ _iⁿ _lⁿ

(3-6)

式中，下標

i

表示網格點的位置； 上標

n

與

n + 1

表示兩相鄰之計算時間，

n

為每個計算時段的啟始時間，

n + 1

為目前程式所欲計算流場變數之時間；

上標

n + 1

^#表示每次疊代運算後所得最新之流場變數， 且於每個計算時段 的第一次運算中，上標為

n + 1

^#與

n

之流場變數設定為相同；

Δ t = t

ⁿ⁺¹

− t

ⁿ。 由於流場變數隨時間而改變，因此引入一權重係數(time weighting factor)

Φ

用以加權兩時段的流場變數，以考慮時間積分的效應；由於將（2-1）式轉 換成交錯格網上之差分式，（2-1）式中，

q

_l 為單位渠長之支流側流量以

x q

_l

Q

^l

= Δ

2 ^表示，

Q

^l為支流側流量，Δ

x

為主流方向之空間間距，−

Q

_l為合 流，

+ Q

_l為分流。

動量方程式的部分，在具自由液面的天然河川中，對流項為主要影響 項，故使用顯式有限解析法將其離散化，以全微分的形式表示之(葉等, 1996):

l l

f

q u

x gAS gA Z A Q Q x Dt

DQ − ±

∂

− ∂

∂

− ∂

=

β

( )

(3-7) 其中，

x Q A Q t

Q Dt

DQ

∂ + ∂

∂

= ∂

β (3-8)

上式中，

D / Dt

表示對於時間之全微分運算；

q

_l為單位渠寬之支流側流量

l l

l

B

q

=

Q

，

Q

_l為支流側流量，

B

_l為支流之寬度，

+ q

_l為分流之處理，

− q

_l為

合流之處理；

u

_l為支流在主流方向的速度分量。

全微分的特性線速度軌跡為:

) 1 ( =

=

=

β u β A

Q dt

dx

(3-9)

由式(3-8)可看出移流項部份為非線性，其移流速度本身即為一待求之 流場變數，因此式(3-2)中特性線之運動軌跡並無法直接計算求得。為解決

此一問題，必須先將式(3-8)局部線性化，為各計算點取其特徵速度

u

^*，用 以替代原有之移流速度

u

。當然，對於任一計算時段， 各計算元素中所取 之特徵速度要能適當表示實際移流速度之大小。如此，透過線性化的過 程，並藉由特徵速度計算其特性線之運動軌跡，如此將使得(3-8)式變成一 線性非齊次雙曲線型偏微分方程式，故式(3-2)便能用以求解式(3-8)中移流 項的部份。

由於在一時間段內的變化量跟總變化量比起來要小很多，因此可以前 一時段的值當作已知值代替。在研究中為進一步考慮(3-9)的移流速度 u 和 (3-7)式的源項(source term)之影響在時間上的變化，分別將(3-7)、(3-9)式對 時間積分，得

dt u q x gAS

gA Z A

Q Q x

Q

Q

2

⁻

1

^{= −} ∫

1²

_⎢⎣ ^⎡ _∂ ^∂ _⎢⎣ ^⎡ _⎥⎦ ^{⎤ + ∂} _∂ ⁺

_{f l l}

_⎥⎦ ^⎤

η η η

η

β _m

(3-10)

∫

=

− ²

1 1 2

η η η

η

X udt

X

(3-11) 其中，η1和η2分別代表特性曲線與前一時段之空間軸的交點和目前時 段之空間軸的交點位置。

動量方程式在此交錯格網上之差分式如下表示：

( )

( ) ( )

( ) ( ) ( )

( ) ( )

[

^{(1 )}

]

^{( ) (1 ) 0}

) 1 (

) 1 (

#

#

#

#

#

#

#

#

#

#

1 ) 1 ( 1

1 1

1 1

1 1

1 1

1

1 1 1

⎥=

⎥⎦

⎤

⎢⎢

⎣

⎡

⎟⎟

⎠

⎞

⎜⎜

⎝ Φ⎛

−

⎟+

⎟

⎠

⎞

⎜⎜

⎝ Φ⎛ Φ

− + Φ

⎥+

⎥⎦

⎤

⎢⎢

⎣

⎡ ⎟

⎠

⎜ ⎞

⎝

⎛ Δ Φ −

−

⎟⎟+

⎠

⎜⎜ ⎞

⎝

⎛

Δ Φ −

+

⎥⎥

⎥⎥

⎥

⎦

⎤

⎢⎢

⎢⎢

⎢

⎣

⎡

⎟⎟

⎟⎟

⎟

⎠

⎞

⎜⎜

⎜⎜

⎜

⎝

⎛

Δ

⎟⎠

⎜ ⎞

⎝

−⎛

⎟⎠

⎜ ⎞

⎝

⎛ Φ

− +

⎟⎟

⎟⎟

⎟

⎠

⎞

⎜⎜

⎜⎜

⎜

⎝

⎛

Δ

⎟⎠

⎜ ⎞

⎝

−⎛

⎟⎠

⎜ ⎞

⎝

⎛ Φ

Δ +

−

+ + +

+

+ +

+

+ +

+ +

n l n

l n l

l n

l n l

fi n n

fi n i

n s n n r n

s n n r i

n

s n

n r n

s n

n r i n

n i

B u u Q

B S Q

gA S

gA

x n

Z gA Z

x n

Z gA Z

x n

A Q A

Q x Q

n A Q A

Q t Q

Q Q

m

η

η

η

η β β

(3-12)

上式中，考慮流場可能存在亞臨界流與超臨界流同時發生之混合流況，抑 或兩種流況於同一位置交替發生，因此式(3-12)非齊次項部份之數值離算

差分法將其離散化，而當超臨界流時，則採用上風法(前項或後項差分法)，

以考慮明渠流中擾動波傳遞的方向。即若為亞臨界流(-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 )

^#（包 含參數

A Z S

_f

A

Q , Q , , ,

）代表在目前時段下疊代過程中新計算得之變數值；Φ 為權重係數；

Δ t

為時間間距；

Δ x

為空間間距。

為計算

Q

_ηⁿ₁，首先要決定特徵速度，以確定η 點的位置，根據實際模

1

擬的經驗，每次疊代運算後所得之最新流速

u

_i⁽₋ⁿ₁^+1)#適合用來做為特徵速度。

故

x

_η1以如下形式近似：

x

_η1

= x

_i

− u

_i^{( )}₋ⁿ₁^+1#

Δ t

(3-13) 下標為η (包含參數 Q、A)以前一時段斷面上的值線性內差而得，如

1 Q

_ηⁿ₁以 如下形式表示：

x x Q x

x x Q x

Q

ⁿ _i^{n i} _i^{n i} Δ + −

Δ

= ^{−1 1}− ⁻¹ ₋⁻₁^{1 1}

1

η η

η (3-14) 式(3-6) 、(3-12) 中等式右邊的物理量，皆可用前一時刻之初始條件(上 標

n

)，與每次疊代運算所得之最新流場變數(上標

( ^{n + 1} )

^#)來直接計算求 得，繼而獲得所求時刻之最新流況(上標

( n

+¹

)

)。

3.4 邊界條件(參考楊 1996 顯示有限解析法於超亞臨界混合流之研究 國立交通大學碩士論文 )

邊界條件為數值模式中極為重要的一環，除了所設定的條件與個數必 須符合物理意義之外，設定的方式亦必須兼顧模式的功能性與適用性。在 一維計算流場中，上下游水理邊界條件的型態通常有以下三種：

(i) 以流量歷線為邊界條件，給定流量與時間之函數關係式，其形式為

)

(t

f

Q =

；

(ii) 以率定曲線為邊界條件，給定流量與水深之函數關係式，其形式為

)

(h f

Q =

；

(iii) 以水深歷線為邊界條件，給定水深與時間之函數關係式，其形式為

)

(t f h =

。

除了給定的邊界條件外，求解邊界上其它未知物理量之一般化處理技 巧係透過特性法(characteristics method)來估算之(Garcia-Navarro & Saviron , 1992)；本研究即引用特性法求解邊界未知物理量。沿流線方向擾動波傳遞 的二個特性方程式，可用來近似或推估邊界上的未知物理量。這二個特性 方程式，分別存在特性速度向量(

u

+

c

)及(

u

−

c

)，其中

c

為水面擾動波之 波速。當流場型態發生變化時，擾動波傳遞的方向(上游或下游)可能會跟 著改變。舉例來說，當流速方向指向下游為正值時，特性速度(

u

+

c

)將會 指向下游，而(

u

−

c

)在亞臨界流況下會指向上游，如圖3-1所示；而在超臨 界流況下，特徵速度(

u

+

c

)及(

u

−

c

)則會指向下游，如圖3-2所示。因此，

針對不同的流況，需要不同的方式來處理邊界條件的問題。

根據Garcia-Navarro & Saviron (1992)之邊界條件概述，在亞臨界流況 下，先就上游邊界而言，由於特性速度(

u

+

c

)在時間-空間座標內所形成 之特性線落在模擬流場區域的外面，因而無法據此獲得沿著特性線所傳遞 的流場訊息，為彌補不足，可根據實際流況給定上游邊界條件，例如給定 流量歷線Q=f(t)。至於未知的上游邊界水深，則可利用以下的相對應之方 程式(compabillity equation)估算之：

μ

÷

∂ + + ∂

−

= { (

₀

) }

Dt DQ x

A B S gA S Dt gA

DA

h

f (3-15)

B

g A u

u − − +

−

= β ( β

²

β )

²

μ

(3-16) 當上式中之β 等於 1 時，特性線軌跡為

c B u

g A dt u

dx = − = −

(3-17) 再者，就下游邊界而言，由於特性速度(

u

−

c

)所形成之運動軌跡落在 模擬流場區域的外面，因此也必須根據實際流況給定一個下流邊界條件，

例如水深歷線 h=f(t)，至於未知物理量 Q，則可利用以下的特性方程式估 算之：

Dt DA x

A B S gA S Dt gA

DQ

h

f

− μ

∂ + ∂

−

= (

₀

)

(3-18) 下游邊界的處理類似於上游邊界處，唯其

B g A u

u + − +

−

= β ( β

²

β )

²

μ

(3-19) 當上式中之β 等於 1 時，特性軌跡為

u c B

g A dt u

dx = + = +

(3-20) 在進行特性線推導時，因為邊界並無側流量加入，所以在動量方程式 中並不需要考慮，詳細推導請見附錄一。若水流流況為超臨界流時，就上 游邊界而言，由於特性速度(

u

+

c

)及(

u

−

c

)在時間-空間座標內所形成之 特性線皆落在模擬流場區域的外面，因而無法據此獲得沿著特性線所傳遞 的流場訊息，可以根據實際流況給定上游兩個邊界條件，例如根據流量歷 線Q=f(t)及水深歷線 h=f(t) 做為邊界條件；因兩條特性線皆往下游方向傳 遞，故於下游邊界處，可以利用此兩條特性線及其相對應之方程式求解流 量及通水面積兩個變數。

而水流流況為超亞臨界混合流時，將超亞臨界混合流流場的超臨界流 與亞臨界流區域分開來討論之。首先，對於超臨界流區域，此兩個邊界條 件必須同上述的給定在上游邊界，以符合流場特性波傳遞的方向，然而水 躍發生的地方，依據內部邊界條件，可得到一個超臨界流區域的下游邊界 水位。再者，對於亞臨界流區域而言，此一內部邊界條件，可當作其上游 邊界條件，所以僅須在下游邊界給定一水位條件。

綜合以上所述，邊界點之計算乃透過特性法之觀念，由兩條特性線交 匯而得，沿特性線之擾動波波速可表示成

u c

dt

dx = ±

，當亞臨界流時，

u < c

即兩擾動波波速具相反方向，同交於邊界點，此時位於計算區外之一條特 性線須由人為給定此一邊界條件；當臨界流時，

u = c

即其中一擾動波波速 為零，另一速率由水流決定為正向（向下游）或負向（向上游），若水流 方向為正，則在上游需給定一邊界條件，而下游則不須給定，直接由前一 時刻之下游邊界變量給定即可。當超臨界流時，

u > c

即兩擾動波波速具相 同方向，若水流方向為正，則上游處須同時給定兩邊界條件，而下游則不 須給定。而當超亞臨界混合流時，則須採用亞臨界流及超臨界流之邊界條 件，即上游處同時給定兩個邊界條件，下游處給定一個邊界條件即可，詳 細之外部開放邊界條件處理原則如附錄二所示。

3.5 交匯區邊界條件(內部邊界條件)

本模式將主流與支流分為兩獨立之渠道，在交匯點處，主流的動量方 程式必須加入側流的動量，連續方程式也必須考慮支流所匯入或分出的流 量影響，而支流下游或上游的邊界條件則由匯流點所求得之水位給定之

（並利用邊界的處理方法，求得另一未知物理量Q），在滿足疊代精度的 條件下，將可得到主、支流各自的上下游水位，方程式如下表示：

= 0

∂ + ∂

∂

∂

x Q t

A

^ζ (3-21)

其中

Q

^ζ=

Q

⁽_n⁺¹⁾

− Q

⁽_n⁻¹⁾

m Q

_l表示考慮主、支流匯流後之流量變化情形，n 表示匯流點，

^Q

⁽ⁿ⁺¹⁾表示匯流點下一個格點流量，

^Q

⁽ⁿ⁻¹⁾表示匯流點上一個格 點流量，

^Q

^l表示支流流量（−

Q

_l為合流，+

Q

_l為分流）。

在匯流處附近主、支流斷面流速在一維模式中是假設分佈均勻的，且 在主支流各種不同交匯角度皆滿足動量守恆的關係式，其方程式表示如

0 )]

(

[ + =

∂ + ∂

∂ + ∂

∂

∂

l l

f

q u

x gAS gA Z U x Q t

Q

β m (3-22) 式中q 為單位長度側流量（−

q

時為合流，

+ q

時為分流）、

u

_l為支流在主 流方向的速度分量。

3.5.1 處理方式 1.合流:

以 i 表示主流上游段最後斷面位置，i+1 表示主流下游段第一斷面位 置，j 表示支流最後斷面位置，水流流向及示意圖如下所示：

i

i+1

j

主流

支流

圖3-3 一般合流交匯區斷面示意圖

本模式主、合流交匯區採用處理方式為將 i 點水位給予 j 點，因為一 維定床河道的處理方式為上游給定流量，下游給定水位計算，當匯流時，

支流之下游水位未知，故假設 i 點水位和 j 點水位相同，如此合流之邊界 給定，可由特性線方法，算出支流下游邊界（j 點）之流量。

2.分流:

若以i 表示主流上游段最後斷面位置，i+1 表示主流下游段第一斷面位 置，j 表示分流最初斷面位置，水流流向及示意圖如下所示：

i i+1

主流

j 支流

圖3-4 一般分流交匯區隔點示意圖

本模式主、分流交匯區採用處理方式為i+1 點水位與 j 點水位相同，

因為分流之支流上游邊界條件缺少，即利用特性線下游給水位求流量的方 式，在處理上游流量過程中，上游流量經過特性速度之修正，算出分流之 起始流量。

3.5.2 交匯區連續方程式、動量方程式處理

為符合質量守恆方程式與動量守恆方程式，由控制方程式中得知，在 每一斷面所欲求之未知數，在使用交錯斷面的原則下，每一個斷面只需要 一個方程式求解一個未知數；亦即只需要一條特性線。在一般的流況下，

主、支流交匯處水流阻滯現象明顯，所以本模式處理時，利用動量方程式 計算斷面流量，連續方程式計算通水斷面積。本文假設在匯流處附近主、

支流斷面流速在一維模式中是分佈均勻的，而僅考慮支流匯入主流時的交 匯角度，以決定速度在主流方向的分量，在合流將支流在主流方向的動量 分量加入主流的動量守恆的關係式中。在分流則將主流方向的動量扣除支 流的動量分量維持守恆的關係式中。所以在本模式中，於主流匯流處除了 必須滿足質量守恆的條件外，還必須滿足動量守恆方程式，而此動量守恆 方程式也一併滿足支流匯入主流所產生的動量影響，其方程式如(2-2)式所 示。

3.6 超亞臨界混合流處理方法(內部邊界條件)

當超亞臨界混合流即所謂之跨臨界流發生時，在超臨界流的下游端，

即水躍發生的地方，會因為上、下游水深的共軛特性自然形成一個限制條 件，對整個流場而言，此限制條件為一內部邊界條件，它的存在使得超臨 界流區域多一個邊界條件， 導致數值計算發生over-determined 的現象。

為解決內部邊界條件所造成的問題，可利用內部相鄰計算點的水位高程，

先透過預測-修正數值方法求取超亞臨界混合流區域的各斷面通水面積及 流量，以減少地形驟變對模式演算的影響且降低數值震盪幅度、增加之後