軟岩河床沖刷動床數值模式之研發與應用

國 立 交 通 大 學

土木工程學系

博士論文

軟岩河床沖刷動床數值模式之研發與應用

Development and application of soft bedrock mobile-bed model

研究生：廖仲達

指導教授：葉克家

軟岩河床沖刷動床數值模式之研發與應用

學生：廖仲達 指導教授：葉克家 博士 國立交通大學土木工程研究所

摘要

台灣河川具有坡陡流急、沖淤劇烈之特性，當河床表層具有卵礫 石或砂質粒徑之沖積層受到破壞時，若河川又位於軟弱岩盤帶，常造 成底部之軟岩裸露，進而出現嚴重之河床沖刷與下切現象，加速河床 沖刷、影響河防安全。國內外目前對於河道沖淤之數模研究，多著重 在沖積層與傳統輸砂行為之模擬，而軟岩沖刷問題有別於傳統輸砂， 須有額外機制之建置以描述其沖刷行為。 本研究開發具有模擬軟岩河床沖刷功能之水平二維動床數值模 式，引用水力沖蝕機制之流功門檻沖蝕率公式，將該機制整合於顯式 有限解析法動床模式上，考慮河床分層、分區設定，整合懸浮載源交 換機制，讓模式可自動考慮岩床之沖、淤行為。以濁水溪集集堰下游 與大安溪峽谷河段為例，進行軟岩河道沖刷過程之檢定、驗證與比較 分析。 在現場具有明顯深槽、高灘地地形案例之模擬中，流況劇烈複 雜，模擬過程可維持數值穩定，說明本模式之相關理論架構及乾溼點 處理技巧等可適用於現場案例之模擬；此外，模擬底床縱剖面、橫斷 面之沖刷量與實測趨勢相同，檢定驗證成果良好。水流對河床之沖蝕 能力，與流速、剪應力及流功呈某種正比關係，配合遷急點沖刷機制， 驗證流功門檻沖蝕率公式可適用於模擬軟岩河川之沖刷過程。

Development and application of soft bedrock

mobile-bed model

Student: Chung-Ta Liao Advisor: Keh-Chia Yeh Department of Civil Engineering

National Chiao-Tung University

Abstract

Steep slope and severe bed change are the general characteristics of Taiwan’s rivers. These characteristics often cause the bed armour layer flushed away, make the bedrock exposed, and then increase the channel incision rate. Due to the bedrock exposed and the channel characteristics, it makes the incision problem more seriously. Most mobile-bed models of past few decades aimed at the sediment transport of alluvial channel. However, the soft bedrock channel incision processes are different with the alluvial channel. Some empirical formula and mechanisms are needed to build in the model to simulate the erosion processes.

In this study, a bedrock river evolution mechanism is included in a 2-D mobile-bed model, called the explicit finite analytic model. The model can consider both incision and deposition over the bedrock, by combining a stream power type of bedrock erosion rate formula with the flow and sediment transport modules. The multiple bed layers distributed in bedrock river are also considered. In the theorem of non-equilibrium suspended sediment concentration, the velocity scale is set to adjust the overloaded or under-loaded situation over the bedrock riverbed. The downstream reach of Chi-Chi Weir of Choshui River and the Taan River are chosen as the study sites for the model calibration, validation and

comparison.

In the field cases, the proposed model is stable in the simulation of complex river morphology with the main channel and floodplains. For the stream power threshold bedrock erosion rate formula and knickpoint evolution mechanism, the comparison of the calculated bed changes with field data show the model is capable of predicting bedrock incision, with reasonable accuracy and reliability. The simulated bed changes, thalwegs and cross-sectional profiles agree well with the measured data. It also indicates that the velocity, shear stress and critical stream power of flow could form a proper relation of erosion threshold in bedrock rivers. This model could be served as the prediction tool for the bedrock rivers.

Keywords: soft bedrock channel incision, numerical model, stream

誌謝

承蒙指導老師 葉教授克家碩士班、博士班期間之指導與教誨， 讓學生於理論及實務皆可並重學習，以順利完成學業。感謝博士論文 口試委員 蔡教授長泰、蔡博士惠峰、盧教授昭堯、賴博士進松、廖 教授志中，承蒙各位委員對學生論文之指教與建議，讓學生受益良 多。感謝廖教授志中、潘教授以文與黃博士明萬之研究團隊，讓本研 究得以順利進行。感謝美國 NCCHE 之王教授、賈博士、張博士與丁 博士等人，讓學生對研究領域有不同視野。 感謝葉門與楊門團隊歷屆學長姐與學弟妹所有成員，感謝璁哥與 添哥研發之 EFA 模式雛形，讓本研究可基於此架構延伸；感謝棒哥 時常適時提點研究方向或生活瑣事；感謝貓博的美食日記；感謝冠曄 經濟拮据時之援助；感謝凱哥，有了你 418 總是有很多趣事；感謝老 唐，有了你協助讓我可更專心於論文；感謝邱邱、洲哥，祝你們博班 生涯順利；感謝 T 哥、阿富、岱玲、豪哥、健賓、邱哥，祝你們都可 以一帆風順；感謝奇綺，讓我知道甚麼是感恩師父；感謝 Marco、高 思、冠翰、大原、鈞凱，有了你們 418 充滿了血汗及活力；感謝念樺， 感謝這麼多年來在工作、生活等各方面的協助與關照；感謝葉師母， 研究室所有成員生活的大家長。 最後感謝我的家人、父母與兄長，感謝所有親朋好友於過去經濟 困難時的支持；感謝父母與兄長對我求學生涯期間各方面的關愛與照 顧，也感謝你們能容忍我求學與工作之課業繁忙、聚少離多，希望大 家都能永遠身體健康、平安幸福。 仲達 敬上 2013/7

1.

目錄

摘要 ... I Abstract ... II 目錄 ... V 表目錄 ... VIII 圖目錄 ... IX 符號說明 ... XIII 第一章 緒論... 1 1.1 研究動機與目的 ... 1 1.2 文獻回顧 ... 2 1.2.1 泥砂磨蝕機制 ... 2 1.2.2 水力沖蝕機制 ... 4 1.2.3 穴蝕機制 ... 5 1.2.4 遷急點沖刷機制 ... 5 1.2.5 岩床沖刷數模回顧 ... 6 1.3 研究方法 ... 8 1.4 章節大綱 ... 9 第二章 水理與輸砂模組理論 ... 13 2.1 水理模組 ... 13 2.2 輸砂模組 ... 16 2.3 有限解析法簡介與特點 ... 21 2.4 水理方程式數值方法... 22 2.5 輸砂方程式數值方法... 23 2.5.1 懸浮載質量守恆方程式之離散化 ... 24 2.5.2 河床載質量守恆方程式之離散化 ... 25 2.6 乾濕點處理技巧 ... 26

第三章 岩床模組理論 ... 30 3.1 岩床模組概念與限制... 30 3.1.1 岩床模組概念 ... 30 3.1.2 岩床模組限制 ... 31 3.2 沖蝕指數決定法 ... 31 3.2.1 材料強度參數（Ms） ... 32 3.2.2 尺寸參數（Kb） ... 32 3.2.3 抗剪強度參數（Kd） ... 33 3.2.4 地盤構造條件參數（Js） ... 33 3.2.5 顆粒材料之沖蝕指數 ... 34 3.3 岩床沖蝕率模式 ... 34 3.3.1 流功門檻機制 ... 34 3.3.2 遷急點沖刷機制 ... 35 3.4 岩床模組架構與模式整合 ... 37 第四章 沖刷機制探討與敏感度分析 ... 46 4.1 沖刷機制探討 ... 46 4.1.1 模擬條件 ... 46 4.1.2 模擬成果 ... 47 4.1.3 輸砂模組敏感度分析 ... 48 4.1.4 岩床模組敏感度分析 ... 49 4.1.5 格網密度敏感度分析 ... 50 4.2 遷急點沖刷機制探討... 52 4.2.1 模擬條件 ... 52 4.2.2 模擬成果 ... 52 第五章 現場案例模擬與分析 ... 70 5.1 濁水溪案例 ... 70

5.1.1 濁水溪背景概述 ... 70 5.1.2 濁水溪模擬條件 ... 71 5.1.3 濁水溪模擬成果 ... 73 5.2 大安溪案例 ... 75 5.2.1 大安溪背景概述 ... 75 5.2.2 大安溪模擬條件 ... 76 5.2.3 大安溪模擬成果 ... 78 第六章 結論與建議 ... 112 6.1 結論 ... 112 6.2 建議 ... 114 參考文獻 ... 116

2.

表目錄

表 3-1 岩石材料強度表 ... 39 表 3-2 節理組數參數表 ... 39 表 3-3 節理面粗糙參數表 ... 40 表 3-4 節理風化參數表 ... 41 表 3-5 地盤構造條件參數表 ... 42 表 3-6 本研究軟岩河床沖刷模式相關機制彙整表 ... 43 表 4-1 run24 案例格網密度敏感度分析格網尺寸表 ... 54 表 5-1 濁水溪案例模式選用參數表 ... 81 表 5-2 大安溪案例模式選用參數表 ... 81

3.

圖目錄

圖 1-1 台灣各軟岩河床沖蝕河段現場照片 ... 10 圖 1-2 河床載顆粒沖擊岩盤表面示意圖（修改自 Foley, 1980） .... 10 圖 1-3 流功與沖蝕指數門檻關係圖 ... 11 圖 1-4 遷急點地形對應之剪應力 與底床啟動剪應力_c ... 12 圖 1-5 不同遷急點地形變遷形式示意圖 ... 12 圖 2-1 顯式有限解析法特性線示意圖 ... 28 圖 2-2 一維空間之沈滓移流軌跡示意圖 ... 29 圖 3-1 水流沖蝕能力與剪力、流速或流功之關係曲線圖 ... 43 圖 3-2 遷急點後退機制河床坡度判斷示意圖 ... 44 圖 3-3 本研究軟岩河道沖刷模式演算流程圖 ... 44 圖 3-4 本研究軟岩河道沖刷模式各模組整合關係流程圖 ... 45 圖 3-5 本研究軟岩河道沖刷模式河床分層概念圖 ... 45 圖 4-1 輸砂模組模擬底床高程縱剖面比較圖(t=1hr) ... 55 圖 4-2 輸砂模組模擬底床高程縱剖面比較圖(t=4.5hr) ... 55 圖 4-3 輸砂模組模擬底床高程縱剖面比較圖(t=10hr) ... 55 圖 4-4 輸砂模組模擬底床沖淤變化縱剖面比較圖(t=1hr) ... 56 圖 4-5 輸砂模組模擬底床沖淤變化縱剖面比較圖(t=4.5hr) ... 56 圖 4-6 輸砂模組模擬底床沖淤變化縱剖面比較圖(t=10hr) ... 56 圖 4-7 輸砂模組模擬與實測底床變化誤差分析比較圖 ... 57 圖 4-8 輸砂模組不同時刻流速與底床高程縱剖面圖 ... 58 圖 4-9 岩床模組模擬底床高程縱剖面比較圖 ... 59 圖 4-10 底床累積沖刷深度與流速逐時變化圖(CS11 位置) ... 60 圖 4-11 底床累積沖刷深度與流速逐時變化圖(CS21 位置) ... 60

圖 4-13 輸砂模組作用層厚度敏感度分析比較圖(t=1hr) ... 61 圖 4-14 輸砂模組作用層厚度敏感度分析比較圖(t=4.5hr) ... 61 圖 4-15 輸砂模組作用層厚度敏感度分析比較圖(t=10hr) ... 61 圖 4-16 輸砂模組作用層厚度敏感度分析比較圖(t=1hr) ... 62 圖 4-17 輸砂模組作用層厚度敏感度分析比較圖(t=4.5hr) ... 62 圖 4-18 輸砂模組作用層厚度敏感度分析比較圖(t=10hr) ... 62 圖 4-19 岩床模組無因次沖蝕係數(Kp)敏感度分析比較圖 ... 63 圖 4-20 岩床模組沖蝕指數(Kh)敏感度分析比較圖 ... 64 圖 4-21 無因次沖蝕係數 Kp敏感度分析底床高程逐時變化圖 ... 65 圖 4-22 沖蝕指數 Kh敏感度分析底床高程逐時變化圖 ... 65 圖 4-23 不同格網密度模擬底床沖淤變化縱剖面比較圖(t=1hr) ... 66 圖 4-24 不同格網密度模擬底床沖淤變化縱剖面比較圖(t=4.5hr) ... 66 圖 4-25 不同格網密度模擬底床沖淤變化縱剖面比較圖(t=10hr) .... 66 圖 4-26 疏格網模擬與實測底床變化誤差分析比較圖 ... 67 圖 4-27 一般格網模擬與實測底床變化誤差分析比較圖 ... 67 圖 4-28 密格網模擬與實測底床變化誤差分析比較圖 ... 67 圖 4-29 岩床模組以不同格網密度模擬底床高程縱剖面比較圖 1 ... 68 圖 4-30 岩床模組以不同格網密度模擬底床高程縱剖面比較圖 2 ... 68 圖 4-31 遷急點沖刷案例模擬底床縱剖面圖 ... 69 圖 5-1 濁水溪流域圖 ... 82 圖 5-2 濁水溪集集攔河堰河段航照圖與現場照片 ... 83 圖 5-3 濁水溪案例模擬格網圖 ... 84 圖 5-4 濁水溪集集攔河堰河段 87 ~ 96 年實測河床縱剖面變化圖 .. 85 圖 5-5 濁水溪案例集集堰流量、泥砂量邊界條件歷線圖 ... 86 圖 5-6 濁水溪案例集集堰流量、泥砂量邊界條件歷線圖 ... 86 圖 5-7 濁水溪案例名竹大橋水位邊界條件歷線圖 ... 87

圖 5-8 濁水溪案例名竹大橋水位邊界條件歷線圖 ... 87 圖 5-9 濁水溪案例河床沖蝕指數 Kh分布圖 ... 88 圖 5-10 濁水溪案例模擬洪峰流速分布圖 ... 89 圖 5-11 濁水溪案例模擬洪峰流場分布圖 ... 89 圖 5-12 濁水溪案例模擬與實測底床縱剖面圖(檢定) ... 90 圖 5-13 濁水溪案例模擬與實測底床縱剖面圖(驗證) ... 91 圖 5-14 濁水溪案例模擬與實測斷面圖(檢定) ... 92 圖 5-15 濁水溪案例模擬與實測斷面圖(檢定) ... 92 圖 5-16 濁水溪案例模擬與實測斷面圖(檢定) ... 93 圖 5-17 濁水溪案例模擬與實測斷面圖(檢定) ... 93 圖 5-18 濁水溪案例模擬與實測斷面圖(驗證) ... 94 圖 5-19 濁水溪案例模擬與實測斷面圖(驗證) ... 94 圖 5-20 濁水溪案例模擬與實測斷面圖(驗證) ... 95 圖 5-21 濁水溪案例模擬與實測斷面圖(驗證) ... 95 圖 5-22 大安溪流域圖 ... 96 圖 5-23 大安溪研究範圍航照圖與峽谷段現場照片 ... 97 圖 5-24 大安溪 87 ~ 98 年實測河床縱剖面變化圖 ... 98 圖 5-25 大安溪案例模擬格網圖 ... 99 圖 5-26 大安溪案例上游流量、泥砂量邊界條件歷線圖(89~93 年) ... 100 圖 5-27 大安溪案例上游流量、泥砂量邊界條件歷線圖(93~98 年) ... 100 圖 5-28 大安溪案例下游水位邊界條件歷線圖(89~93 年) ... 101 圖 5-29 大安溪案例下游水位邊界條件歷線圖(93~98 年) ... 101 圖 5-30 大安溪案例河床沖蝕指數 Kh分布圖 ... 102 圖 5-31 大安溪案例模擬洪峰流速分布圖 ... 103

圖 5-32 大安溪案例模擬剪應力 τ 分布圖 ... 103 圖 5-33 大安溪案例模擬福祿數分布圖 ... 103 圖 5-34 大安溪案例模擬洪峰流速與實測底床縱剖面圖 ... 104 圖 5-35 大安溪案例模擬底床沖淤與流速縱剖面圖(檢定) ... 105 圖 5-36 大安溪案例模擬底床沖淤與流速縱剖面圖(驗證) ... 106 圖 5-37 大安溪案例模擬與實測底床縱剖面圖 ... 107 圖 5-38 大安溪案例模擬與實測底床縱剖面誤差分析圖 ... 108 圖 5-39 大安溪案例模擬與實測斷面圖(檢定) ... 108 圖 5-40 大安溪案例模擬與實測斷面圖(檢定) ... 109 圖 5-41 大安溪案例模擬與實測斷面圖(檢定) ... 109 圖 5-42 大安溪案例模擬與實測斷面圖(檢定) ... 109 圖 5-43 大安溪案例模擬與實測斷面圖(驗證) ... 110 圖 5-44 大安溪案例模擬與實測斷面圖(驗證) ... 110 圖 5-45 大安溪案例模擬與實測斷面圖(驗證) ... 110 圖 5-46 大安溪案例模擬與實測斷面圖(驗證) ... 111

4.

符號說明

Ad：集水區面積 s b ：二次流比例係數 C：懸浮質濃度 Ccoeff：泥砂磨蝕機制檢定參數 Cf：乾溼點函數係數 Cq：泥砂磨蝕機制覆蓋項修正係數 Cs：泥砂磨蝕機制懸浮項修正係數 a c ：參考高程沈滓濃度 e c ：水深平均濃度 n d ：前一時刻計算之水深 E：岩床沖蝕速率 Ek：遷急點後退沖刷率 Em：作用層厚度 Ep：流功門檻沖蝕率 Es：泥砂磨蝕機制岩床沖蝕率 Et：剪力門檻機制岩床沖蝕率 F：體力（body force） Fe：岩床面泥砂覆蓋層厚度之相關比例係數； f ：Darcy Weisbach 係數 g：重力加速度 mn g ：逆變之公制係數 m n g ：協變之公制係數

h：水深 hw：乾溼點參考水深 Js：地盤構造條件參數 Kb：尺寸參數 Kd：抗剪強度參數 Kh：沖蝕指數 Kp：流功門檻無因次沖蝕係數 Kr：沖蝕阻力係數 Kt：剪力門檻無因次沖蝕係數 k：von Karman 係數 v k ：軟岩強度參數 s L ：泥砂顆粒躍動長度 Ms：材料強度參數 n：曼寧糙度係數 P：流功 p：孔隙率 Pcr：臨界流功 Pp：壓力 q：單位寬度流量 qb：河床載通量 qs：單位河寬之供砂量 qt：單位河寬之輸砂能力 R：水力半徑 ra：曲率半徑 S：河床坡度

Scr：遷急點後退臨界坡度 f S ：作用層源 Sr：懸浮載源 s：砂比重 t：時間 U：水深平均流速 u：x 座標方向流速 * u ：剪力速度 u：水深平均 x 座標方向流速 ns u ：自由水面之二次流流速 Vmd：乾溼點最小乾床速度 v：y 座標方向流速 v：水深平均 y 座標方向流速 s w ：泥砂顆粒之沖擊速度 f w ：泥砂顆粒沈降速度 Y：軟岩彈性模數 b Z ：底床高程 s Z ：水面高程 T  ：軟岩張力強度 α：流功門檻機制權重指數  ：粒徑百分比  ：二次流函數中距底床之深度 ρ：流體密度 ρs：泥砂密度

c  ：底床臨界剪應力 ：動力黏滯係數 ：水流方向 ：水流側方向 λ：遷急點後退指數 x  ：空間間距 t  ：時間間距

第一章 緒論

1.1 研究動機與目的

台灣河川由於具有坡陡流急、沖淤劇烈之特性，當河床表層具有 卵礫石或砂質粒徑之沖積層受到破壞時，若河川又位於軟弱岩盤帶， 常造成底部之軟岩裸露，進而出現嚴重之河床沖刷與下切現象，對跨 河構造物與堤防等之安全造成威脅。台灣眾多河川中，如大漢溪後村 堰下游、頭前溪隆恩堰與中正大橋下游、大甲溪石岡壩下游、大安溪 峽谷河段、濁水溪集集攔河堰下游、八掌溪仁義潭攔河堰下游等，河 床底層之軟弱岩盤皆已裸露，且受到水流持續沖刷，形成上下游河床 嚴重高程落差或下切現象，如圖 1-1 所示。 由於岩床沖刷問題有別於傳統沖積層河川之輸砂，根據不同地質 材料與水砂條件常造成不同沖刷行為，僅考慮一般輸砂過程之動床模 式在應用上恐嫌不足，須有額外岩床沖刷機制之探討與建置，以適當 反映其劇烈沖刷過程。然而，國內外在河道沖淤之數模研究，多針對 沖積層與一般輸砂行為進行探討，即僅考慮泥砂顆粒運移與底床之交 換，鮮少針對軟岩沖刷機制與水理輸砂模式同時進行探討。 本研究在學理創新與貢獻方面，為結合水理輸砂模式與岩床沖刷 機制，模式可考慮岩床表面之沖刷，並整合輸砂過程中懸浮載濃度剖 面超載時之淤積效應，目的為研發具有實用價值之軟岩沖刷水理輸砂 模式，探討流功門檻軟岩沖刷機制之特性，利用現場案例之實測資 料，進行模式之檢定、驗證、比較，並評估岩床沖刷機制之合理性與 適用性，以作為岩床沖刷過程與工程應用之參考。

1.2 文獻回顧

根據早期文獻對岩床河道（bedrock channel）之定義（Gilbert, 1877;

Howard, 1994; Montgomery et al., 1996; Whipple, 2004），普遍認為岩質

河床為缺乏連續覆蓋之沖積層泥砂，且長期輸砂條件為輸砂通量 （sediment flux）小於輸砂能力（sediment capacity）之河道。由於以 缺乏覆蓋層泥砂之定義方式較為主觀，常受到輸砂條件之不確定性影 響，且套用在不同國家之河川與實際情形有很大差異，因此，Turowski et al.（2008）回顧定義岩床河道之相關文獻，以台灣 81 條河川為例， 重新對岩床河道進行定義，認為若一條河川之拓寬、底床降低與側向 變遷本質上與岩床沖刷有關，即可定義為岩床河川。 一般學理上認為水流作用為岩床沖刷之主要原因，主要作用為水 流在河道流動過程中對岩床所造成之切割與沖擊。就岩床沖刷之物理 過程機制來說，目前較有系統之沖刷、沖蝕成因大致可分為三種型

態，分別為泥砂磨蝕（sediment abrasion）、水力沖蝕（hydraulic erosion）

及穴蝕（cavitation）機制，茲說明如下。 1.2.1 泥砂磨蝕機制 關於泥砂磨蝕概念，基本假設為岩床表面之磨蝕為水流帶動之泥 砂滾動、跳動所造成。Bitter（1963a、1963b）曾分析泥砂沖擊河床 表面所造成之磨蝕（wear）效應，區分為低角度切削（cutting）磨蝕 及高角度變形（deformation）磨蝕；Foley（1980）則將此概念應用 於小區域河床面之河床載磨蝕模型（bed-load abrasion model），並提 出磨蝕率公式，概念如圖 1-2 所示。

Sklar and Dietrich（2004）修正 Foley（1980）之河床載磨蝕公式， 將磨蝕率以顆粒沖擊影響之平均岩石體積、沖擊率及裸露區塊三種因

子 表 示 之 ， 針 對 其 磨 蝕 經 驗 公 式 進 行 敏 感 度 分 析 ， 分 析 供 砂 量 （sediment supply）、無因次相對剪力、粒徑大小等與磨蝕率之變化趨 勢，其泥砂磨蝕率公式如下： _        t s s T v s s s q q L k Y w q E ₂ 1 2  (1-1) 式中，Es 為磨蝕率(m/s)；qs 為單位寬之輸砂率(kg/m/s)；qt 為單位寬 之輸砂能力(kg/m/s)；w_s為泥砂顆粒之沖擊速度(m/s)；Y 為岩體彈性 模數(Pa)；k_v為岩體強度參數，須進行檢定（根據以往試驗，其值介 於 1012~1013）；T為岩體張力強度(Pa)；Ls為泥砂顆粒躍動長度(m)。 上式中(1-qs/qt)項反映岩體表面之覆蓋效應，qs、qt 為泥砂非平衡

與平衡狀態下之輸砂量，前者為輸砂率（sediment transport rate），後 者為輸砂能力（sediment transport capacity）。以一河道之上游邊界來

說，其輸砂率可視為上游之供砂量（sediment supply），輸砂能力則為

該水流條件下之理論輸砂量，可由一般平衡輸砂公式計算，當輸砂率 大於輸砂能力時，水流為超載，河床將產生淤積，岩體表面覆蓋泥砂 顆粒；反之，河床則磨蝕沖刷；若輸砂率等於輸砂能力，則為不沖不 淤狀態。

Sklar and Dietrich（2004）之泥砂磨蝕理論原先僅就河床材質本身 與河床載過程進行探討，並無考慮河床之覆蓋效應，而後因相關磨蝕 試驗結果呈現供砂量到一定程度後，磨蝕率有下降之趨勢，因此才導 入覆蓋效應項修正其磨蝕率。惟岩床淤積與磨蝕具有不同之物理機 制，建議磨蝕率與覆蓋效應兩者應各別分析判斷後再整合，在參數檢 定可更單純，降低其不確定性。

1.2.2 水力沖蝕機制

水力沖蝕機制之基本概念認為軟岩沖蝕率與水流作用力、底床剪 應力、能量消散等因子相關，總體來看，可包括泥砂磨蝕或穴蝕產生 之沖蝕量。Howard and Kerby（1983）認為軟岩河道之下切速率隨剪 應力增加而增加，並提出河床沖蝕速率為河道坡降與集水面積之關係 式，如式(1-2)所示： m1 m2 d rA S K dt dz _ (1-2) 式中，z 為底床高程；t 為時間；Kr為沖蝕阻力係數；Ad為集水區面 積；S 為河道坡降；m1 與 m2 為常數。 Annandale（1995）曾分析可定性描述岩塊抽離與河床材質沖蝕成 因之沖蝕關係，而後以明渠水力學常見之四種流況，包含溯源沖刷

（headcuts）、水躍（hydraulic jumps）、底床坡度變化（changes in bed

slope）、明渠流（open channel flow），計算不同流況下之能量消散率

（流功）以特徵化水流之沖蝕能量，再利用 150 組河床材質範圍自凝 聚性粒狀材料到巨大之岩塊不等之現地觀測資料，計算河床材質之沖 蝕指數 Kh（erodibility index），提出水流流功與沖蝕指數具有某種沖 蝕門檻函數關係。 Annandale（2006）更有系統地將水流流功與沖蝕指數迴歸，進而 獲得河床材料發生沖蝕時所需之水流臨界能量門檻值，如圖 1-3、圖 1-4 所示。圖中之 Kh值可分為粒狀材質與岩盤材質，其臨界流功 Pcrit 與沖蝕指數 Kh之門檻關係如下： 1 . 0 , 1 . 0 , 48 . 0 75 . 0 44 . 0       h h h h crit K K K K P (1-3)

岩質河床之岩塊因水流作用而產生較大尺度之崩落或抽離河床 現象，稱為岩塊抽離（rock-block plucking），其為相當複雜之物理與 化學風化崩解過程。Whipple et al.（2000）藉由現地調查，認為岩性、 弱面間距、節理、層面等為決定岩質河床沖刷機制之因素，岩塊抽離 過程中小裂縫經水力作用擴大為破裂面、隨著泥砂顆粒逐漸地透過磨 蝕作用磨蝕弱面、及物理或化學風化崩解作用，讓弱面完全擴展連 通，最後導致獨立岩塊之脫離。此外，Whipple et al.（2000）發現局 部區域大尺度之沖刷量來自抽離機制之貢獻較多，大流域（大於 20 km2）之軟岩沖刷機制以懸浮載磨蝕為主，而較小流域之沖刷機制則 為岩塊抽離或河床載磨蝕。 1.2.3 穴蝕機制 穴蝕機制係指岩床表面受到渦流與紊流作用產生之氣泡沖擊效 應，使岩床產生大規模之溝槽或孔洞（Whipple et al., 2000）。從物理 過程觀察，穴蝕機制可視為小型或大型渦流對岩床產生之磨蝕效應。 1.2.4 遷急點沖刷機制 劇烈沖刷河道常伴隨發生遷急點（knickpoint）沖刷，並加速河道 變遷過程，導致河床沖蝕範圍加長與加劇。遷急點根據地質學之定 義，係指河床面地貌突然落差之陡降點。Gardner（1983）曾採用均 質材料（砂、粉砂或黏土、高嶺土拌合固結）進行水槽試驗，探討遷 急點變遷形式，認為底床剪應力 與底床臨界剪應力c為遷急點變遷 之重要因子，圖 1-5 為遷急點地形對應之底床及沖蝕啟動剪應力變化 圖，其剪應力在遷急點區域相對上下游河段較高。Gardner（1983） 綜合其試驗結果及前人研究成果，提出三種遷急點變遷形式，包括傾

retreat），其變遷形式則由 與c之大小、河床載不連續性、河床材質

抵抗水流沖刷之空間變異性所決定，變遷形式示意圖如圖 1-6 所示。

1.2.5 岩床沖刷數模回顧

(一)CCHE2D 模式

岩床河道沖淤數模之研究方面，美國國家計算水科學與工程中心 （National Center for Computational Hydroscience and Engineering, NCCHE）之 CCHE2D 模式，為一套具有變量流、輸砂、沖淤變化、 污染傳輸、岸壁沖刷模擬功能之水平二維動床模式（Jia and Wang,

1997；Wu and Wang, 2004a），近年更考慮了軟岩沖刷機制（水規所，

2008）。該模式參考 Sklar and Dietrich（2004）之泥砂磨蝕機制，針對

沖蝕率公式進行泥砂懸浮與覆蓋效應項之修正，以指數曲線化方式使 沖蝕率與輸砂量之變化關係更加平滑化，反映其非線性過程，修正後 之沖蝕率公式表示如下： ) ( 36 . 0 2 cov ) 1 ]( ) 1 [( 08 . 0 C Csus c s T v s coeff s g e k Y q C E          (1-4) ) ( cov t s q _q q C C   _(1-5) 2 * ) ( 5 . 1 f s sus _w u C C   _(1-6) 式中，Es為指數曲線化之泥砂磨蝕率(m/s)；為流體密度(kg/m 3 )；s 為砂密度(kg/m3 )；g 為重力加速度(m/s2)； 為水流作用於底床之剪應 力 (N/m2 ) ；c 為 底 床 臨 界 剪 應 力 (N/m 2_{) ；} * u 為 剪 力 速 度 (m/s) ，  / *  u ；w_f 為泥砂顆粒沈降速度(m/s)；Ccoeff 為檢定係數；Ccov 為 磨蝕機制覆蓋項；Cq為修正係數(Cq 0.1，需進行檢定)；Csus為磨蝕

機制懸浮項；Cs為修正係數(Cs 0.01，需進行檢定)。 式(1-5)中，岩床表面之覆蓋效應以(qs/qt)之變化來反映，若岩床表 面覆蓋泥砂顆粒，則由顆粒跳動所引起之沖刷將減少；若輸砂量大於 或等於輸砂能力，則沖刷停止。岩床表面之懸浮效應以(u*/wf)項來反 映，當水流強度夠大，所有顆粒皆為懸浮，泥砂跳動將被懸浮行為取 代，則沖刷停止。 但泥砂磨蝕機制之啟動條件僅適用於有輸砂量之條件，無法反映 清水沖刷問題（由式(1-4)知輸砂率或供砂量為零，沖刷量即為零）， 此外，泥砂磨蝕機制之檢定參數頗多，常因輸砂條件與軟岩參數之不 確定性，造成整體沖刷量高估或低估之缺點。CCHE2D 模式亦無法反 映不同岩性強度之岩盤分佈於河道高灘地、深槽等位置抵抗沖蝕能力 之差異。 (二)SRH-2D 模式

美國內政部墾務局（Bureau of Reclamation, US Department of Interior）所發展之 SRH-2D 模式，過去在與水利署水利規劃試驗所合 作下，曾針對台灣軟岩河川之劇烈沖刷問題進行動床模式開發（Lai et al., 2011）。該模式採用之軟岩沖刷機制包括泥砂磨蝕（Sklar and Dietrich, 2004）與水力沖蝕（Annandale, 1995）兩種，並嘗試合併此 二機制進行計算，表示如下： s e t e E F E F E    (1-7)         1 c t t KU E   (1-8)

5 . 1 2 * 5 . 0 ) ( 1 1 1 1 08 . 0 _                     _   f c s v s s w u k gq E     (1-9) 3 2 6 1 66 . 7                R k P_{crit s} c   (1-10) 1 . 0 , 1 . 0 , 48 . 0 75 . 0 44 . 0       h h h h crit K K K K P (1-11) 式中，E 為岩床沖蝕率(m/s)；Fe為岩床面泥砂覆蓋層厚度之相關比例 係數；Et 為剪力門檻機制岩床沖蝕率(m/s)；Es 為泥砂磨蝕率(m/s)； Kt為無因次沖蝕係數，需要配合試驗及現場資料進行檢定；U 為水深 平均之水流流速(m/s)；R 為水力半徑(m)；ks 為粗糙高度（roughness

height）；臨界流功 Pcrit引用 Annandale（2006）之門檻關係，如式(1-3)。

SRH-2D 模式應用於台灣濁水溪之軟岩河道案例中，由於缺乏現 地調查資料，相關參數如底床臨界剪力_c、岩床沖蝕指數 Kh 等之給 定多採用假設值，對於不同機制與物理參數關係之探討較為缺乏。目 前國內外對於二維水理輸砂模式之演算技術已相當成熟，岩床沖刷數 模研發之關鍵在於岩床機制本身，包括機制之選用、模式整合、特性 探討、調查資料與現地驗證等，此為本研究探討之重點。

1.3 研究方法

本研究開發具有模擬軟岩河床沖刷功能之水平二維動床數值模 式，引用流功門檻（廖志中等，2012）沖蝕機制公式，將該機制整合 於顯式有限解析法（explicit finite analytic, EFA）動床模式上，以濁水 溪集集堰下游與大安溪峽谷河段為例，進行軟岩河道沖刷過程之檢 定、驗證與比較分析，探討軟岩質河床之劇烈沖刷特性，該模式可作

為軟岩河川治理規劃應用之參考。

1.4 章節大綱

本研究第一章為「緒論」，對於研究之緣起、目的、方法及文獻 回顧進行說明。第二章為「水理與輸砂模組理論基礎」，針對本研究 模式之水理與輸砂模組進行說明，包含基本假設、控制方程式、數值 方法、輔助關係式等。第三章為「岩床模組理論基礎」，對本研究模 式岩床模組採用之軟岩沖刷機制、理論與模式建置方式進行說明。第 四章為「沖刷機制探討與敏感度分析」，藉由數模針對輸砂模組與岩 床模組之沖刷機制進行探討與參數敏感度分析。第五章為「現場案例 模擬與分析」，採用濁水溪集集堰下游與大安溪峽谷段為案例，進行 模式之檢定、驗證與比較分析。第六章為「結論與建議」，針對研究 成果進行結論說明，並提出後續建議。

(a)頭前溪中正橋河段 (b)大甲溪石岡壩下游

(c)濁水溪集集堰下游 (d)八掌溪仁義潭攔河堰下游

圖 1-1 台灣各軟岩河床沖蝕河段現場照片

1.00E-05 1.00E-04 1.00E-03 1.00E-02 1.00E-01 1.00E+00

1.00E-11 1.00E-09 1.00E-07 1.00E-05 1.00E-03 1.00E-01

S tr e a m P o w e r (K W /m 2) Erodibility Index scour Threshold

圖 1-3 流功與沖蝕指數門檻關係圖 1 （redrawn from Annandale, 2006）

圖 1-5 遷急點地形對應之剪應力 與底床啟動剪應力_c （redrawn from Gardner, 1983）

圖 1-6 不同遷急點地形變遷形式示意圖 （redrawn from Gardner, 1983）

5.

第二章 水理與輸砂模組理論

本研究研發具有模擬軟岩河道沖刷功能之水平二維動床模式，相 關軟岩模組之理論基礎係建構在顯式有限解析法動床模式上，該模式 由交通大學葉克家教授研究團隊研發多年而成，詳細內容可參考 Hsu

et al.（2000）、Lin et al.（2006）、許至璁（2002）、林恩添（2005）、

廖仲達（2006）等文獻，有關模式水理與輸砂模組之理論說明如下。

2.1 水理模組

對於不可壓縮流之控制方程式，其張量表示式如下： 連續方程式 0 ,  i i V (2-1) 動量方程式 i i mn mn m p im i m m i F V g P g V V t V _{   }   , , ,    (2-2) 上式中， i V 為速度向量V在i_{座標方向之分量；t 為時間；} mn g 、gm n為 逆 變 （ contravariant ） 和 協 變 （ covariant ） 之 公 制 係 數 （ metric

coefficient）；ρ 為流體密度；Pp為壓力；為動力黏滯係數；F 為體 力（body force），變數上下標之值為 1~3，代表不同的座標方向。 對一般水深平均河道模式而言，假設其水深方向變化遠小於水平 方向，並忽略風力、科氏力之影響，同時假設靜水壓分佈，將壓力轉 換成水深形式等，可簡化原三維控制方程式之複雜性。對於任意水深 函數 f 而言，水深平均可定義為：



 s b Z Z f z t dz h t f(,, ) 1 (,, , ) (2-3) 式中，Z_b為底床高程；Z_s為水面高程；h為水深。在推導水深平均式 前，式(2-2)中左邊流速與流速梯度之乘積項可利用連續方程式轉換為 保守形式，以方向為例： uv g g g g g g g u g g g g g g g z uw uv g u g z u w u v g u u g                                            ] 2 [ 2 1 ] 2 [ 2 1 ] [ ] [ ] [ 22 22 12 12 11 11 5 . 0 22 2 22 22 12 12 11 11 5 . 0 11 5 . 0 22 2 5 . 0 11 5 . 0 22 5 . 0 11           (2-4) 轉換後不再有流速與流速梯度之乘積項，而流場變數偏微分項可由 Leibnits 法則推導，式(2-2)之水深平均表示式如下： i m m i i mn mn m im i m m i F v u V g P g V V t V _{    }   , ) ( , , ,    (2-5) 在水深平均控制方程式推導過程中還需要 Leibnits 法則、自由水面邊 界條件、底床邊界條件、流速剖面等條件，說明如下： A. Leibniz 法則 s b s b s b s b z z Z Z b b s s Z Z Z Z z t z f h f dz f z t z f z t z f dz f dz t z f ] ) , , , ( [ ) ( ) , , , ( ) , , , ( ) , , , (                                   







(2-6) 式中為函數 f 之任意自變數。 B.自由水面邊界條件

假設流體無法穿過自由水面，可得到自由水面邊界條件： 0 5 . 0 22 5 . 0 11 _              s s s s s s z v g z u g t z w (2-7) 式中 Zs代表水面位置。 C.底床邊界條件 和自由水面邊界條件類似，底床邊界條件為： 0 5 . 0 22 5 . 0 11 _              b b b b b b z v g z u g t z w (2-8) 式中 Zb代表底床位置。 D.流速剖面 水深平均模式中不考慮主流及橫向流速在水深方向之分佈情 形，在水深平均條件下其主流及橫向流速如下（Zimmermann & Kennedy, 1978）： N h u N N u 1 ()1/ (2-9) ) 1 2 (    h u v v _ns  (2-10) 式中，u為主流流速，即水流方向之流速；u為水深平均主流流速；v 為橫向流速，為水深平均流速與彎道二次流流速之和，若在直線道， 二次流流速為零（本研究所指之二次流，為明渠流因彎道離心力所產 生者，一般二次流水深上半部指向凹岸、下半部指向凸岸，並會造成 泥砂由凹岸向凸岸之側向運移）；v為水深平均側向流速； 為距底床 之深度；h為水深；u_ns為自由水面之二次流流速，u_ns b_shu/r_a，b_s為

可由格網點間距與其相對座標所決定；N 為常數，Nk 8/ f ，k為

von Karman 係數， f 為 Darcy Weisbach 係數。

2.2 輸砂模組

對於某一粒徑懸浮載質量守恆方程式 q v C t C _  _{ }    1 ) ( (2-11) 對於某一粒徑作用層質量守恆方程式 0 ) ( ) 1 (        m _{b r f} s q S S t E p    (2-12) 對於整體河床輸砂質量守恆方程式 0 ) ( ) 1 (        b



_{b r} s q S t z p   (2-13) 式中，C 為懸浮載濃度；q為擴散通量及由重力作用下之通量；_s 為泥砂密度；p為孔隙率； 為粒徑百分比；E_m為作用層厚度；q_b 為某種粒徑之河床載通量；S_r為懸浮載源；Sf 為作用層源；zb為底床 高程。 輸砂模組之控制方程式中，懸浮載質量守恆方程式求解水體不同 粒徑大小之懸浮載移流、擴散過程，同時考慮河床載與懸浮載間之交 換機制；作用層質量守恆方程式則求解作用層中，不同粒徑組成與懸 浮載源、作用層源間之交換機制；整體河床輸砂質量守恆方程式則將 上述兩者整合，並反映為底床高程之沖淤變化。 求解式(2-11)至(2-13)時，尚需一些輔助關係式來決定q_b、Em、Sr、 f S 等變數，以下針對非黏性沈滓輔助關係式予以說明。

A.河床載通量（qb） van Rijn（1984a）係以泥砂之中值粒徑為代表粒徑來計算非均勻 質砂之河床載體積通量。其後，Spasojevic（1990）針對每一粒徑（D_k）， 將 van Rijn 之公式稍加修正之。假設河床載運移僅發生在作用層內， 其內某粒徑之百分組成表示為。在一般非均勻之河床質中，較細顆 粒可能被隱藏在較粗顆粒之間，而不易被水流帶動，故 Karim, Holly and Yang（1987）提出一個簡單之經驗因子，稱之為隱藏因子（hiding factor，），對河床載通量予以修正。綜合上述之影響因子而得其河 床載通量如下： 3 . 0 * 1 . 2 ) 1 ( 053 . 0 ) 1 ( ) ( k k k k s k z k t b k z b D T D gD s B D q B q           (2-14) 式中， 13 2 * ] ) 1 ( [  g s D D_k  _k  為無因次顆粒粒徑； ₂ * 2 * 2 * ) ( ) ( k c k c k u u u T   為輸送參 數； c g u u_*  為有效河床剪力速度； ) 3 12 log( 18 90 D d c  為顆粒蔡司係數；  _s s 為砂比重； 為運動滯度；u_*c為臨界剪應力； 0.85 50 ) ( D D_k k   ； 4 . 0 10 4 . 0 0 0 ) ln( 325 . 0 25 . 0 1 * * * *                k k k k w u w u w u w u  ；w_k為粒徑k之沈降速度。 B.懸浮載源（Sr）

對於水深平均模式而言，由於水深方向之維度忽略不考慮，因此

懸浮載源（Sr）可視為聯繫懸浮載與河床載之源項。許多模式採用式

(2-15)計算懸浮載源 Sr，如 Lin & Shen（1984）定義：

) (c c S_r _{s a} (2-15) 式中， 為水體所承載之參考高程沈滓濃度（c_a）與水深平均濃度（c_e ） 之比值，定義為： e a c c   (2-16) 在缺乏濃度剖面之情況下， 之大小並不容易決定，此外，由 Rouse 平衡濃度剖面公式可知，沈滓特性與水理流況對於 之比值會 有直接影響。且水深平均模式利用式(2-15)計算懸浮載源 Sr會有高估 之結果，導致水體承載之濃度瞬間達到平衡濃度之狀態，為改善此問 題，本模式採用平衡濃度剖面積分法。 一般來說，水流在超載情況下，超過平衡濃度的部分會逐漸沈降 至底床；反之，在減載情況下，水流會沖刷河床以補充不足平衡濃度 之部分，超載與減載濃度剖面與底床交換機制之示意圖如圖 2-1 所 示。所謂平衡濃度剖面積分法，係指利用水體承載濃度剖面會趨於平 衡濃度之概念，積分計算兩者剖面之差異，以求得懸浮載與河床載之 交換速率，即為懸浮載源（Sr），反映非平衡之輸砂過程。 關於平衡濃度剖面c_e(z)，本模式採用 van Rijn（1984b）之平衡濃 度剖面公式，表示如下： ) 5 . 0 ( 4 ) ( ) ( _          zh a a a e e h c z c     5 . 0  h z for (2-17)

φ a a a e _δ h z z h δ c z c _         ) ( ) ( ) ( 0.5 h z for (2-18) 式中，c_e(z)為平衡濃度剖面，_a為參考高程， 為懸浮沈滓參數 （ws/cu_*），其中 ws為泥砂沈降速度，c為泥砂擴散與水流擴散 係數之比值， 為 von Karman 係數，u_*為剪力速度。 所謂水體承載濃度剖面，係指參考高程面以上之懸浮沈滓實際分 布情形，本模式假設水體承載濃度剖面與平衡濃度剖面具有相同之無 因次濃度剖面，水體承載濃度剖面可表示如下： ) 5 . 0 ( 4 ) ( ) ( _          zh a a b e h c z c     5 . 0  h z for (2-19)             ) ( ) ( ) ( a a b h z z h c z c 0.5 h z for (2-20) 式中，c (z) 為水體承載濃度剖面，c_b 為參考高程之水體承載濃度 ) ( a c  。當c_b大於c_a時，水流為超載；反之，則為減載狀態。 在得到平衡濃度剖面c_e(z)與水體承載濃度剖面c(z)後，定義一作用 高度如下： t w A_w _* (2-21) 式中，Aw 為作用高度；w*為速度尺度，超載時定義為沈滓沈降速度 ws，減載時定義為河床質躍起速度 wl（若為岩床，躍起速度 wl為零， 作用高度 Aw為零，懸浮載源 Sr積分後為零，即表示岩床表面無河床 減載時產生之懸浮泥砂顆粒）；t為輸砂計算時間間距。河床質躍起

速度 wl 定義為發生跳躍時離開底床之瞬間垂直速度（Hu & Hui,

         3.2 4.5log 1.2 2 . 1 1 . 3 * u w_l (2-22) 式中，為水流強度=_b (_s )gD，D 為河床質粒徑。作用高度定義 了參考高程向上之作用區間，在此作用區間內經由時間t後，水體承 載濃度調整為平衡濃度。根據此假設，懸浮載源可表示為： t dz z c z c S h a a A z z e r _  



     [ ( ) ( )] (2-23) C.作用層厚度（E_m） 作用層厚度代表河床表層與底層泥砂交換之作用範圍，通常與水 深或河床粒徑大小有關，其將影響每單位時間內河床表層與底層之泥 砂交換量。過去文獻中已有不少決定作用層厚度之方法，沖刷現象發

生時，根據 Bennet and Nordin（1977）之研究，E_m可以下式表示：

) ( _bt 1 _bt

m L z z

E    (2-24)

式 中 ， L 為 數 值 參 數 。 當 河 床 表 面 接 近 護 甲 條 件 時 （ armored

condition），作用層厚度接近零，在這種情況下，可用 Borah et al.

（1982）所提出護甲層之厚度（armored-layer thickness），予以修正： p D z z C E m k m k k t b t b m    _



  1 1 ) ( 1 _max  (2-25) 式中，D_{m為不產生移動的最小顆粒粒徑。另外，作用層在淤積期間可} 定義為： ) ( 1 1 t b t b t m t m E z z E      _(2-26)

若不採用上述理論方式計算作用層厚度，一般其範圍約介於 1/10 ~ 1/100 倍水深之間，可針對不同流況與床質條件做調整，或利用實 測資料做率定與校正以決定其值。 D.作用層源（S_f） 作用層源為底部母層（active stratum）與頂面之升降而產生，當 其下降時， )] ( [ ) 1 ( _{s b m} s f z E t p S         (2-27) 其中，s為母層內某一粒徑之百分組成比例。當作用層底部上升時， 式中之s則改為作用層粒徑之百分組成。

2.3 有限解析法簡介與特點

在計算流力領域中常見的數值方法有許多種，諸如有限差分法、 特性法、有限元素法、有限體積法、有限解析法等，各有其優缺點與 適用性。有限差分法在推導及離散化之過程較為容易，但通常需建構 在正交結構型格網（orthogonal structured grid）上模擬，當遭遇不規 則邊界問題時需另加處理，且處理時計算區域（computational domain） 與物理區域（physical domain）座標轉換過程過於繁雜，通常適用於 較簡易規則問題；有限元素法與有限體積法之特點在於可使用非正交 非結構型格網（non-orthogonal unstructured grid）模擬，其格網多為 三角型使模擬邊界形狀容易建構，對於需要更細部之流場描述區塊， 可只加密該區塊之格網數而不必增加總格網維度，儘管離散化與推導 過程相較其他方法複雜，且計算時間亦較多，仍為目前廣泛使用於計 算流力之數值方法之一。

有限解析法可分為顯式法（explicit）與隱式法（implicit）兩種， 該法最早為陳景仁（C. J. Chen）教授所創，創立時為隱式法，其特 色有：(a)可以結構性格網在卡式座標系統處理不規則之邊界；(b)對 於個別計算元素數值離散採局部解析解（local analytic solution）來近 似，可把數值演算之捨入誤差降低；(c)數值穩定性（stability）佳， 為無條件穩定。但隱式有限解析法通常適用於橢圓（elliptic）與拋物 線型（parabolic）偏微分方程式，對於河川水流之雙曲線型（hyperbolic） 偏微分方程式不太適用，因此才有顯式有限解析法之發展（Dai, 1994），初期僅應用在求解無自由表面之 Navier-Stokes 方程式。 顯式有限解析法在求解對流傳輸方程式中，對流項以特性法概念 求得式中變量之局部解析解，並依時變量再透過適當給定之初始條件 求得，為顯式法之特色。此外，該法和一般常見顯式數值方法一樣， 受到可蘭數（Courant number）小於或等於 1 之穩定性限制，模擬時 間間距t無法自行給定，但解法較為簡單，故應用上仍有其優點在。

2.4 水理方程式數值方法

EFA 水理模式採顯式有限解析法來求解雙曲線型淺水波水流動量 方程式，可求得移流項部分之局部解析解，以下用卡式座標中二維一 階線性齊次雙曲線型微分方程式做說明： 0    _{x y} t u v  (2-28) 式中，u, v 分別為 x, y 方向上之速度，當起始條件(x ,y ,0)(x ,y )被 適當給定時，可求得上式(2-28)之解析解如下： ) , ( ) , , (x₀ y₀ t  x₀ut y₀vt  (2-29) 式中(x₀,y₀)為待求點之座標，xx₀ut與yy₀vt定義一條特性線

運動軌跡。該特性線由起始平面 D 點出發，如圖 2-2，經過t時刻後， 交於(x₀,y₀)位置，該特性線上具有相同之物理量，而 D 點之座標可根 據移流速度 u, v 由特性軌跡x_D x₀ut與y_D y₀vt加以推求。但對 於非線性之移流方程式，可透過局部線性化的方法，將移流速度 u, v 以特徵速度代替，儘管移流速度是非線性且隨著時間與空間做改變， 在某一計算時間內仍可假設其為常數，代表某一計算時間內之平均移 流速度。 對於明渠流之動量方程式，其為非齊次混合型形式： g yy xx y x t u v v(  )F  (2-30) 等式右邊利用已知物理量以顯式法直接計算，所得結果視為源項，直 接加入式(2-30)中，可得： t F v t v y t u x t y x , , ) (   ,   )[ ( _xx _yy) _g] ( _{0 0}  _{0 0}    (2-31) 以上計算方法即為 Dai（1994）所提出之顯式有限解析法，本模式將 此方法應用於求解具自由液面之流場，並在模式中引入疊代計算流程 以修正特徵速度與源項。

2.5 輸砂方程式數值方法

求解之輸砂方程式包含了：懸浮載質量守恆方程式(2-11)、作用 層質量守恆方程式(2-12)、整體河床輸砂質量守恆方程式(2-13)三種類 型，懸浮載質量守恆方程式為雙曲線 -橢圓型方程式，具有移流 （advection）與擴散（diffusion）之特性，可表現出懸浮載在水體中 移流與擴散之行為，模式中採用特性法求解；作用層質量守恆方程式 可表現作用層與河床載通量質量守恆之特性，同理整體河床輸砂質量

採用有限差分法離散之。此三種類型之輸砂方程式彼此之物理量有高 度相關，藉由懸浮載源S_r做連結，故以結合演算法同時求解此三種類 型之方程式為佳，再利用 Newton-Raphson 疊代聯立求解懸浮載濃度 C、粒徑百分組成 、底床高程z_b之變動量。 2.5.1 懸浮載質量守恆方程式之離散化 懸浮載質量守恆方程式以特性法解析求解，移流項部分可以全微 分形式表示如下： Dt c D c V t c i i     , i=1,2 (2-32) 其移流軌跡為： i i V t    i=1,2 (2-33) 式中，Vi為懸浮沈滓之移流速度，假設與水流速度相同，而上標 i 則 代表不同座標方向。圖 2-3 為一維空間之沈滓移流軌跡示意圖，假設 移流軌跡介於計算時刻 n1 t 與tn之間，同時定義tn1時刻之端點為到達 點 A（arrival point）， n t 時刻之端點為離開點 D（departure point），離 開點 D 之懸浮質濃度可由起始條件求得，但由於模式採固定格點， 離開點 D 不一定落在格點上，因此需由該點鄰近之點做內插求得濃 度。 一般來說，由於懸浮質變化尺度遠大於河床質，因此兩者在時間 尺度上有所差異，而在要聯立求解之前提下，懸浮載方程式可蘭數將 相對較大，使得移流軌跡穿越若干個計算格點空間，因而需要採取分 段處理之方式以求得較佳之移流軌跡。若移流軌跡由離開點 D 到達 點 A 共跨越 LNS 個格點空間，將軌跡進入及離開各計算格點空間依

序編號為（LNS+1）個節點，各節點相對位置可表示為： ) ( * 2 ) ( 1 1 1 l l l l l l t t u u x x        l=1, 2,…, LNS (2-34) ) ( * 2 ) ( 1 1 1 l l l l l l t t v v y y        l=1, 2,…, LNS (2-35) 式中，下標 l 為節點編號，l=1 代表離開點 D，l=LNS+1 代表到達點 A。利用上式推求節點位置時必須知道各節點上之移流速度，而移流 速度又與節點位置有關，因此可用疊代收斂方式推求正確移流軌跡。 推得正確移流軌跡後，可沿移流軌跡積分式得其離散化式子如 下：



          LNS l l l l a l l a l D A t t h S h S c c 1 1 1 1 _] 2 * ) ( ) ( [   (2-36) 2.5.2 河床載質量守恆方程式之離散化 作用層質量守恆、整體河床輸砂質量守恆方程式以控制體積概念 進行離散化，時間與空間項採後項差分與中央差分進行離散化，並引 入權重係數加權兩時段之變數，其離散化式子如下： 作用層質量守恆離散式： 0 ) ( ) 1 ( )] ) ( ) ( [ ) 1 ( )] ) ( ) ( [ )] ) ( ) ( [ ) 1 ( )] ) ( ) ( [ )] ) ( ) [( ) 1 ( 1 2 2 2 1 2 1 1 2 1 1 2 1 2 1 1 2 1 2 2 1 1 1 2 1 1 2 2 1 1                         p f n rp n rp n s b w n n b n p p n s b s n n b n p p n w b w n e b e p p n w b w n e b e p p n p m n p m s S S S q h q h h h q h q h h h q h q h h h q h q h h h E E t p          (2-37) 整體河床輸砂質量守恆離散式：





[ (1 ) ] 0 )] ) ( ) ( )[ 1 ( )] ) ( ) ( [ )] ) ( ) ( )[ 1 ( )] ) ( ) ( [ )] ) ( ) [( ) 1 ( 1 2 2 2 1 1 2 1 1 2 1 1 2 1 2 1 1 2 1 1 2 2 1 1                 





      n rp n rp n s b w n n b n n s b s n n b n n w b w n e b e n w b w n e b e p p n p b n p b s S S q h q h q h q h q h q h q h q h h h z z t p         (2-38) 式中，為權重係數，模式中採用 0.7。

2.6 乾濕點處理技巧

模式乾濕點處理技巧之詳細理論可參考廖仲達（2006）之研究， 其乾濕點處理係參考蔡智恆（2000）之處理方式進行修正。模式設定 完初始條件後，進入動量方程式演算速度分量前，先判斷模擬範圍格 網點之乾濕狀態，並給予初步修正；動量方程計算完速度分量後，進 入連續方程式計算水深 h，在計算水深同時，再判斷乾濕點狀態，若 為乾點，進入乾濕點修正程式，修正其水深與乾濕點速度分量；動量 及連續方程式內部反覆疊代收斂穩定後進入下一時刻計算。在判斷乾 濕點狀態時，將蔡智恆（2000）之處理方式精簡為八種類型，由於本 研究模式採用顯式法，無解矩陣時之限制，修正方式較隱式法方便。 進行乾濕點修正前，首先要決定參考水深 hw、乾床速度 Vd

（dry-bed transmitted speed）與最小乾床速度 Vmd（minimum dry-bed

transmitted speed）等三個變數。參考水深 hw為模式所允許之最小水 深，其決定依模擬案例尺度不同而異，根據數值經驗採案例水深維度 之 0.01 倍；乾床速度 Vd則為遭遇乾床狀態時所給定之傳遞速度，其 給定時機與該計算點周圍之乾濕狀況有關，其值為：



w



f d C g h V  2* * (2-39) 式中，Cf為係數，其值介於 0 ~ 1 之間，g 為重力加速度，hw為參考 水深；最小乾床速度 Vmd為一個非常小但不為 0 之值，當該計算點判

斷為乾點且不需做速度修正時需給定，根據數值經驗約為模式收斂精 度之 0.1 倍即可。 此修正技巧之理論基礎在於當格網點遭遇乾點時，藉由格網點四 周格點之乾濕狀態，判斷是否給定乾床傳遞速度及參考水深，此速度 及水深為一種假設值，實際物理上並不存在，但在數值計算上藉由此 虛擬值可讓模式不發生奇異點問題，且其值不影響實際計算結果。

z=0 z=a z=h under-loaded concentration profile ca a cb bed level reference level equilibrium concentration profile ca free surface overloaded concentration profile cb t l h A  t s h A  (a)減載 (b)超載 圖 2-1 懸浮載濃度剖面減載與超載示意圖 圖 2-2 顯式有限解析法特性線示意圖

第三章 岩床模組理論

3.1 岩床模組概念與限制

3.1.1 岩床模組概念 本研究根據水深平均條件下所計算之水理因子，包括水深、流 速、剪應力與流功等，導入岩床模組計算，以反映岩床與沖積層之沖 淤分布。在一般明渠流河道中，常認為水力沖蝕為岩床之主要沖刷原 因，泥砂磨蝕所造成之岩床沖刷量則受到現地泥砂條件影響，相對較 複雜且無法掌握。目前文獻中較有系統之泥砂磨蝕機制（Sklar and Dietrich, 2004），其沖刷率公式係根據實驗室磨蝕試驗所獲得，與天 然岩床節理、構造等不可轉移之特性相違，因此泥砂磨蝕機制無法量 化岩床節理與界面風化所造成之影響，現階段應用於明渠流河道尚待 進一步研究及修正。 水力計算中，常用流速大小描述水流流況為湍急或緩慢，剪應力 反映水流作用於單位面積河床之強度特性，而流功可反映水流之能量 特性。過去研究學者多認為，水流對河床之沖蝕能力，應與流速、剪 應力及流功呈現某種正比關係，如圖 3-1 所示，其水流沖蝕能力隨著 相關指標變數增加而增加。本研究之岩床模組，採用流功門檻機制， 由水流條件反映岩床之沖刷量。 根據 Annandale（2006）所提之水力沖蝕概念，提出可定性描述 岩塊抽離與河床材質沖蝕成因之沖蝕門檻模式，岩床適用範圍自粒徑 為無凝聚性之粒狀材料（cohesionless granular）至大塊石（massive hard

rock）等皆適用，沖蝕指數 Kh變化範圍自 10

-2_{至 10}4_{。其中，沖蝕指}

研究關於軟岩之定義，係根據國際岩石力學協會(International Society for Rock Mechanics, ISRM)之分類，其單壓強度介於 0.5MPa ~ 25MPa 間。 3.1.2 岩床模組限制 本研究之軟岩河床沖刷模式適用於一般具有沖積層或軟岩層之 河道，可藉由給定河床質粒徑與軟岩岩性進行沖淤模擬，以反映河床 中不同區塊之沖淤特性。對於跨河構造物如橋墩、攔河堰、丁壩壩頭 與固床工交界處等局部流場與沖刷問題並不適用。雖相關跨河構造物 可藉由格網產生進行幾何形狀與高程之設定，但其局部區域之流場與 沖刷受到水深方向流場與局部紊流影響，屬三維水理現象，使用水平 二維模式計算時有其先天上之限制。

3.2 沖蝕指數決定法

剪力門檻與流功門檻機制具有臨界啟動門檻之概念，其所造成之 岩床沖刷與岩性具有高度相關，其中岩體之岩性可藉由沖蝕指數表 示。沖蝕指數最早由 Kirsten（1982）提出，為四項岩體材料參數之 乘積，各參數於現地調查後依不同條件配合相關經驗式給予不同值， 由四項值之乘積得到沖蝕指數，以作為量化岩床可沖蝕能力之標準， 表示如下： Kh MsKb Kd Js (3-1) 式中，Kh為沖蝕指數；Ms為材料強度參數；Kb為尺寸參數；Kd為抗 剪強度參數；Js為地盤構造條件參數。各參數可由試驗資料或現地觀 測獲得，以下說明四項岩體材料參數之給定方式。

3.2.1 材料強度參數（Ms） 岩體材料強度參數 Ms 值可根據現地調查資料，配合迴歸分析獲 得之經驗公式推求，表示如下： MPa UCS UCS MCR M_s 0.78( )( )1.05, 10 (3-1) MPa UCS UCS MCR M_s ( )( ) 10 (3-2)

式 中 ， UCS 為 岩 體 之 無 圍 壓 縮 強 度 （ unconfined compressive

strength）；MCR 為相對密度係數， 3 10 27  g r MCR



_，



_{r為岩體密度，} g 為重力加速度。彙整不同硬度之岩石材料強度如表 3-1 所示，可看 出岩體材料強度與 UCS、Ms 值存在正相關，強度越高則材料抗沖蝕 能力越強。 3.2.2 尺寸參數（Kb） 尺寸參數 Kb值可量化顆粒或塊體間之尺寸影響，以瞭解不同顆粒 或塊體大小之抗沖蝕能力，由岩石品質指標 RQD （rock quality designation）與節理組數參數 Jn來定義，表示如下： Kb RQD/Jn (3-3) 其中，岩石品質指標之範圍為5 RQD100，節理組數參數範圍為 5 1J_n  ，因此顆粒/塊體尺寸參數範圍為1K_b 100。當無現地鑽孔 取樣資料時，RQD 值可用以下經驗公式獲得：             _0.33 ) ( 10 105 z y x J J J RQD _(3-4) 式中

，

Jx

、

Jy

、

Jz分別為 x、y、z 三個方向上每公尺出現之節理組數。

另外節理組數參數 Jn非直接使用節理組數值，必須參考表 3-2 給定。

因此，RQD 越大、節理組數越少、Kb越大，即岩體越完整（塊體尺

寸越大）抗沖蝕能力越高。

3.2.3 抗剪強度參數（Kd）

抗剪強度參數 Kd值可量化岩體弱面或顆粒間之抗剪強度，由岩體

之節理面粗糙參數（joint roughness number, Jr）與節理改變參數（joint

alteration number, Ja）來定義，表示如下： Kd Jr /Ja (3-5) 節理面粗糙參數 Jr表示不連續面間之粗糙度，受到不連續面形狀 影響，判斷岩石挖掘時節理兩面是否保持緊密狀態或分離，再依不同 接觸面形狀、糙度等獲得此參數。表 3-3 為 Kirsten（1982）根據岩石 節理分離程度、節理面狀態等所分類獲得之節理粗糙數值表。 節理改變參數 Ja與不連續面之空隙填充材料有關，如植生、凝聚 性或非凝聚性材料填充，會影響不連續面之摩擦力，如表 3-4 所示， 其判斷岩石節理分離程度後，再檢視空隙之填充材料性質可獲得此參 數。 3.2.4 地盤構造條件參數（Js） 地盤構造參數 Js和水流流向、岩床較密節理傾向、節理傾角、以

及岩塊形狀有關。岩塊形狀因素採用節理間距比（ratio of joint spacing,

RJS）來代表岩塊形狀比與抗沖蝕程度之影響，從水流流過岩床之縱

剖面觀察岩塊長度比值以計算節理間距比 RJS。

河向或逆著河向）將影響地盤構造條件參數，依照上述資料可以由表

3-5 獲得 Js值。

3.2.5 顆粒材料之沖蝕指數

由於沖蝕指數之決定甚為複雜，且需要大量現地資料與試驗進行 分析以決定各材料參數值，在實務應用上，若缺乏現場調查資料時， Wittler et al.（1998）曾提出在顆粒材質（granular material）河川中，

沖蝕指數 Kh可簡化成以下形式： 3 50 3 50 59 0 . 1 84 . 0 1000 07 . 0 d K J K d K M h s d b s            (3-6) 式中，d50為河床質之中值粒徑，單位為公尺。 沖蝕指數 Kh值為四項岩體參數指標之乘積，係根據現場地表地質 調查後決定，不同區域由於岩性分布不一定均勻，因此實務上 Kh值 通常具有範圍，應用時可取其平均值給定，並根據調查之平面範圍建 立圖層資訊，以供後續應用。

3.3 岩床沖蝕率模式

3.3.1 流功門檻機制 在明渠流條件下，可從水流作用於河床之剪應力



、水深平均流 速 U，獲得平均流功 P（P=



U）。此外，在紊流流況中，瞬間剪應 力與流速變化將隨著平均流功改變，此瞬間作用於岩床之壓力差或剪 應力為造成岩床破壞之主因。 水流流功與其河床材質抵抗能力所換算之臨界流功，為合乎岩床