CCHE2D模式應用於蜿蜒複式河槽變遷之研究

(1)

國立交通大學

土木工程學系

碩士論文

CCHE2D 模式應用於蜿蜒複式河槽變遷之研究

Study on Channel Migration in Meandering Compound

Channel Flow Using CCHE2D Model

研究生：李岱玲

指導教授：葉克家博士

(2)

I CCHE2D模式應用於蜿蜒複式河槽變遷之研究 學生：李岱玲指導教授：葉克家國立交通大學土木工程學系

摘要

台灣的蜿蜒河川於汛期間，洪水漫過高灘地時，在兩岸堤防之間形成蜿蜒複式河槽，有關蜿蜒複式河槽主深槽變遷及河床沖淤方面需要更深入的分析。本研究探討蜿蜒複式河槽溢岸流及未溢岸流對河道穩定及深槽變遷之影響。藉由二維動床模式 CCHE2D 模擬實驗室彎道試驗，了解二次流機制在彎道段之特性，並考慮岸壁沖刷機制，模擬正弦連續彎道河道側向變形過程；溢岸流的研究中，針對三組不同流量及不同上游來砂條件進行模擬，並應用於曾文溪下游之複式斷面河道，探討不同流況下整體河床沖淤趨勢，提供未來河川治理之參考。實驗室及現地案例模擬結果皆與實測資料趨勢吻合。由溢岸流及未溢岸流之模擬結果可知，在未溢岸流情況下，深槽流量對主深槽之側向變形影響最大。當流量增加為滿灘流後，洪水平原流量愈大，深槽蜿蜒產生的 離心力相對減弱，水流受蜿蜒深槽影響愈低。 關鍵詞：CCHE2D 模式、蜿蜒複式河槽、岸壁沖刷、河道變遷

(3)

II

Study on Channel Migration in Meandering Compound Channel Flow Using CCHE2D Model

Student: Tai-Ling Lee Advisor: Keh-Chia Yeh Department of Civil Engineering

National Chiao Tung University

Abstract

Most of the main channels of alluvial rivers are meanders. The main channel becomes the compound channel confined by the embankments during rise of flood. Assessment of channel migration and bed erosion in meandering compound channel is necessary.

The aim of this research is to investigate channel migration in meandering compound channel with inbank and overbank flows using CCHE2D model. Numerical simulations of two experimental bends show the characteristics of the secondary flow. Validations of the model, using sine-generated experimental data with bank erosion, have been performed. The model was also applied to study the morphological changes of experimental meandering compound channels and field case in Tseng-Wen Creek under changing conditions of flow and sediment.

Good agreements between the simulated and available laboratory measurements were obtained. According to simulation results, the discharge in the main channel affects the stability of the main channel. The discharge in the floodplain of compound channel increases as the water depth increases, and the effect of the centrifugal force decreases. The results also point out the importance of the flow interaction between main channel and floodplain under overbank flow condition.

Keywords：CCHE2D Model, meandering compound channel,

(4)

III

誌謝

研究所期間承蒙恩師葉教授克家學業上的諄諄教誨以及生活上的關心，使得兩年的研究生涯過得充實且有意義。感謝成功大學的蔡教授、國網中心蔡惠峰博士，願意將各位寶貴的時間及心力用在我的論文審查上，使得拙作更臻完善。感謝昇學學長、仲達學長、仁凱學長、柏傑學長、彥瑜學姐、唯泰學長、東洲學長平常不吝於回答我的問題，以及生活上的照顧；感謝相處兩年的研究夥伴家榮、韋豪、信富、昀直、舒勤、芳綺、建翔，有了你們，我的研究生活才會如此的精采豐富；感謝還有于軒、健賓的幫忙，令人非常感動，辛苦你們了。感謝每一位曾經幫助過我的人，我都銘記在心。最後要感謝我最愛的父母及家人，給我溫暖與支持。李岱玲 2012.10

(5)

IV 目錄摘要 ... I ABSTRACT ... II 誌謝 ... III 目錄 ... IV 表目錄 ... VII 圖目錄 ... VIII 符號說明 ... III 第一章緒論 ... 1 1.1 前言 ... 1 1.2 研究動機與目的 ... 1 1.3 文獻回顧 ... 4 1.3.1 彎道動床數模 ... 4 1.3.2 岸壁沖刷理論 ... 5 1.3.3 複式河槽 ... 7 第二章 CCHE2D 模式理論及數值方法 ... 13 2.1 水理模式 ... 13 2.2 輸砂模式 ... 17 2.2.1 二次流機制 ... 19 2.2.2 輸砂公式 ... 20

(6)

V 2.2.3 輸砂參數特性說明 ... 25 2.3 岸壁沖刷理論 ... 26 2.3.1 岸壁破壞模組 ... 27 2.3.2 岸壁面沖刷模組 ... 29 2.3.3 格網延展模組 ... 30 2.3.4 岸壁沖刷模組功能限制 ... 32 2.4 模式功能與限制 ... 32 2.5 代表性檢定驗證案例 ... 32 第三章模擬案例 ... 50 3.1 實驗室設計案例 ... 50 3.1.1 單一彎道動床 ... 50 3.1.2 正弦連續彎道 ... 52 3.1.3 蜿蜒複式河槽溢岸流 ... 53 3.2 現場案例 ... 55 第四章結果分析與討論 ... 71 4.1 實驗室案例模擬結果分析 ... 71 4.1.1 單一彎道動床模擬 ... 71 4.1.2 正弦連續彎道側向變形模擬 ... 73 4.1.3 蜿蜒複式河槽溢岸流動床模擬 ... 74

(7)

VI 4.2 現場案例模擬結果分析 ... 80 4.3 小結 ... 82 第五章結論與建議 ... 118 5.1 結論 ... 118 5.2 建議 ... 119 參考文獻 ... 121

(8)

VII

表目錄

表 2.1 Ackers & White’s formula 公式係數(1973, 1990) ... 37

表 2.2 岸壁相關參數建議範圍(Simon et al., 2003) ... 37 表 2.3 CCHE2D 模式功能彙整 ... 38 表 2.4 CCHE2D 模式功能與限制 ... 38 表 2.5 CCHE2D 模式代表性檢定驗證案例彙整表 ... 38 表 2.6 180 °U 型彎道水理模擬參數表 ... 39 表 2.7 正弦彎道水理模擬參數表(θ0=30°) ... 39 表 2.8 正弦彎道水理模擬參數表(θ0=110°) ... 39 表 2.9 180 °U 型彎道動床模擬模擬參數表 ... 39 表 3.1 DHL 及 DLFM 彎道動床模擬各項參數表 ... 57 表 3.2 蜿蜒渠流岸壁侵蝕模擬參數表 ... 57 表 3.3 蜿蜒複式斷面渠道模型幾何條件表 ... 58 表 3.4 複式斷面蜿蜒渠流溢岸流動床模擬參數表 ... 59 表 3.5 代表粒徑與河床質代號 ... 59 表 3.6 模擬現場案例曾文溪中下游主要颱洪事件 ... 60 表 3.7 模式使用泥砂參數設定表 ... 60 表 4.1 CCHE2D 模式單一彎道動床案例模擬與實驗之均方根誤差 ... 84

(9)

VIII

圖目錄

圖 1.1 三種類型蜿蜒複式河槽示意圖(Wormleaton, 2004) ... 11 圖 1.2 蜿蜒複式渠流水流結構圖(Sellin et al., 1993) ... 12 圖 2.1 泥砂在底床的運移(Engelund, 1974) ... 40 圖 2.2 CCHE2D 模式岸壁沖刷計算流程 ... 40

圖 2.3 岸壁土堤破壞示意圖(Osman and Thorne, 1988) ... 41

圖 2.4 岸壁臨界剪應力 ... 41

圖 2.5 CCHE2D 模式未延展前之格網 ... 42

圖 2.6 CCHE2D 模式延展後之格網 ... 42

圖 2.7 懸臂土體破壞形式(模式之限制) ... 43

圖 2.8 180 °U 型彎道水理模擬格網圖 ... 43

圖 2.9 180 °U 型彎道水理模擬與實測流速比較圖(De Vriend, 1979) ... 44

圖 2.10 180 °U 型彎道水理模擬與實測水位比較圖 ... 44 圖 2.11 正弦彎道水理模擬格網圖(θ0=30°) ... 45 圖 2.12 正弦彎道模擬與實測流速比較圖(θ0=30°) ... 45 圖 2.13 正弦彎道水理模擬格網圖(θ0=110°) ... 45 圖 2.14 正弦彎道模擬與實測流速比較圖(θ0=110°) ... 46 圖 2.15 180 °彎道變量流動床模擬(a)實測,(b)計算床形等高線圖(zb h0 ) ... 46

圖 2.16 180°彎道變量流動床模擬底床變化斷面剖面圖(Yen and Lee, 1995)) .... 47

圖 2.17 CCHE2D 模擬蜿蜒渠道底床變化過程圖 ... 48

圖 2.18 內卡河之變量流模擬 ... 48

圖 3.1 Struiksma (1983)DHL 140°彎道實驗幾何形狀圖 ... 61

圖 3.2 DHL 實驗模擬案例使用格網示意圖 ... 61

圖 3.3 Sutmuller and Glerum (1980)DLFM180°彎道實驗幾何形狀圖... 62

圖 3.4 DLFM 實驗模擬案例使用格網示意圖 ... 62 圖 3.5 蜿蜒渠流岸壁侵蝕模擬初始地形圖 ... 63 圖 3.6 蜿蜒渠流岸壁侵蝕模擬初始格網圖 ... 63 圖 3.7 蜿蜒深槽全波長幾何地形圖(cm) ... 64 圖 3.8 複式斷面蜿蜒渠流案例使用格網圖 ... 64 圖 3.9 蜿蜒複式河槽溢岸流模擬結果展示位置圖 ... 65 圖 3.10 複式斷面蜿蜒渠流案例結果斷面位置圖 ... 65 圖 3.11 蜿蜒複式河槽溢岸流案例橫斷面圖 ... 66 圖 3.12 現場案例曾文溪模擬範圍 94 年航拍圖 ... 66 圖 3.13 曾文溪模擬範圍 94 初始 DEM 數值地形分布圖 ... 67 圖 3.14 曾文溪斷面 26 大斷面實測圖 ... 67 圖 3.15 現場案例模擬河床質及曼寧 n 值設定分布圖 ... 68 圖 3.16 麻善大橋含砂量與流量率定曲線 ... 68

(10)

IX 圖 3.17 上游邊界條件麻善大橋流量歷線 ... 69 圖 3.18 上游邊界條件麻善大橋入砂量歷線 ... 69 圖 3.19 下游邊界條件國聖大橋水位歷線 ... 70 圖 4.1 DHL 彎道定床模擬結果流速分佈圖 ... 85 圖 4.2 DLFM 彎道定床模擬結果流速分佈圖 ... 85 圖 4.3 DHL 彎道動床模擬結果流速分佈圖 ... 86 圖 4.4 DLFM 彎道動床模擬結果流速分佈圖 ... 86 圖 4.5 DHL 彎道動床模擬結果底床變動(△Z/h0)地形圖 ... 87 圖 4.6 DHL 彎道動床模擬結果底床變動地形圖(張，2005) ... 87 圖 4.7 DLFM 彎道動床模擬結果底床變動(△Z/h0)地形圖 ... 88 圖 4.8 DLFM 彎道動床實測底床變動(△Z/h0)地形圖 ... 88 圖 4.9 DHL 底床變動(△Z/h0)縱向剖面圖 ... 89 圖 4.10 DLFM 底床變動(△Z/h0)縱向剖面圖 ... 89 圖 4.11 Run 1 模擬結果格網延展趨勢圖 ... 90 圖 4.12 Run 1 案例格網延展結果圖 ... 90 圖 4.13 Run 1 模擬結果流速分佈圖 ... 90 圖 4.14 Run 2 模擬結果格網延展趨勢圖 ... 91 圖 4.15 Run 2 案例格網延展結果圖 ... 91 圖 4.16 Run 2 模擬結果流速分佈圖 ... 91 圖 4.17 Run 1 岸壁變形比較圖上：Rahman (1997)實驗結果 ... 92 中：Nagata (2000) 模擬結果下：CCHE2D 模擬結果 ... 92 圖 4.18 Run 2 岸壁變形比較圖上：Rahman (1997)實驗結果 ... 93 中：Nagata (2000) 模擬結果下：CCHE2D 模擬結果 ... 93 圖 4.19 CASE III 定床水理模擬結果各斷面水面線圖... 94 圖 4.20 CASE I 模擬結果平面流場分佈圖 ... 95 圖 4.21 CASE II 模擬結果平面流場分佈圖 ... 96 圖 4.22 CASE III 模擬結果平面流場分佈圖 ... 97 圖 4.23 CASE I 各斷面不同入砂條件動床模擬結果底床地形圖 ... 98 圖 4.24 CASE II 各斷面不同入砂條件動床模擬結果底床地形圖 ... 99 圖 4.25 CASE III 各斷面不同入砂條件動床模擬結果底床地形圖 ... 100 圖 4.26 CASE III-1 各斷面不同時段動床模擬結果底床地形圖 ... 101 圖 4.27 不同流量上游無來砂各斷面動床模擬結果底床地形圖 ... 102 圖 4.28 CCHE2D 模擬 95 年案例終了模擬結果地形圖 ... 103 圖 4.29 CCHE2D 模擬 95 年案例終了模擬地形變化圖 ... 103 圖 4.30 曾文溪下游 95 年案例深槽最低點與實測資料比較圖 ... 104 圖 4.31 曾文溪下游 95 年案例結果與斷面實測資料代表斷面比較圖 ... 105 圖 4.32 曾文溪下游 95 年案例結果與大斷面實測資料比較圖(1/4) ... 106 圖 4.33 曾文溪下游 95 年案例結果與大斷面實測資料比較圖(2/4) ... 107

(11)

X 圖 4.34 曾文溪下游 95 年案例結果與大斷面實測資料比較圖(3/4) ... 108 圖 4.35 曾文溪下游 95 年案例結果與大斷面實測資料比較圖(4/4) ... 109 圖 4.36 CCHE2D 模擬 98 年案例終了模擬結果地形圖 ... 110 圖 4.37 CCHE2D 模擬 98 年案例終了模擬地形變化圖 ... 110 圖 4.38 CCHE2D 模擬 98 年案例洪峰時段流速分佈圖 ... 111 圖 4.39 曾文溪下游 98 年案例結果與大斷面實測資料代表斷面比較圖 ... 112 圖 4.40 曾文溪下游 98 年案例結果與大斷面實測資料比較圖(1/4) ... 113 圖 4.41 曾文溪下游 98 年案例結果與大斷面實測資料比較圖(2/4) ... 114 圖 4.42 曾文溪下游 98 年案例結果與大斷面實測資料比較圖(3/4) ... 115 圖 4.43 曾文溪下游 98 年案例結果與大斷面實測資料比較圖(4/4) ... 116 圖 4.44 曾文溪下游 98 年案例模擬深槽最低點與實測資料比較圖 ... 117

(12)

III

符號說明

b：渠寬 coh C ：土壤黏性係數；：中值粒徑 F：福祿數 fs：安全係數； Cor f ：科氏力參數；：初始水位 H：岸壁高度； L：泥砂非平衡調適長度； ：彎道段長度 R：水力半徑：渠道中心線平均曲率半徑：雷諾數 * t k q _{：單位寬度河床質載輸砂量；} ：初始渠道坡度：渠道坡度 sk  ：泥砂粒徑 k 之沈降速度； ：馮卡門係數； b



：底床剪應力； τc：岸壁臨界剪應力；

(13)

1

第一章

緒論

1.1 前言 台灣河川豐、枯水期流量差異甚大，枯水期時河川裸露之高灘地，常有民眾植栽，於颱洪時期，洪水歷程急促，致河道侵蝕嚴重，河道深槽變遷，垂直向侵蝕導致河槽刷深，側向岸壁沖蝕改變河道邊界，危及河防建造物之安全及兩岸土地之安定。河道之變遷常造成可用土地之流失，並可能間接污染水質。因此，不論就河川治理規劃方面，如水道治理計畫線之劃設；或河川管理方面，如高灘地管理，均需對河川變遷之可能型態、形成機制等進一步研究。 1.2 研究動機與目的 坡陡流急、富含大量砂石為台灣河川之特色，且蜿蜒曲折，造成河川之流場更為複雜。由於彎道所造成之離心力與底床剪力無法平衡，而在徑向(radial direction)形成二次流(secondary flow)。於凹岸處，接近水面之二次流沖向凹岸，受阻後轉而向下沖擊底床，而後沿底床向凸岸流動。因此，凹岸處之河床受到淘刷而崩塌退縮；而這些被淘刷之河床質則被二次流帶向凸岸，導致凸岸之淤積成長。隨著凹凸岸沖淤之變化，除了垂向侵蝕導致河床嚴重刷深外，亦有側向侵蝕造成岸壁崩坍。河岸退縮使河川蜿蜒幅度可能逐漸加大並有往下游擺盪之趨勢，岸壁沖刷為探討河道變遷的重要

(14)

2 因素之一。蜿蜒型的沖積河道，大多由較深的主槽以及兩側的洪水平原構成。洪水平原多有植被覆蓋，以及各種土地利用、經濟活動等人類社會行為。為保障洪水期間，洪水平原上居民的生命財產安全，在洪水平原建造堤防以限制洪水寬度，因而形成複式斷面河槽。然而台灣本島地形狹長，河川均發源於高山，多源短坡陡，於枯水期時複式斷面的流量較小，水流多集中於主深槽(main channel)，此時河川所呈現出的流況幾乎皆為簡單之蜿蜒單槽緩變速流；然而每當夏秋之際颱風豪雨來臨，河川流量隨降雨暴增，當洪水漫淹至洪水平原時，此時河川之流況則轉變為蜿蜒複式斷面 (meandering compound channel)之水流。複式斷面之水理原乃屬三維之複雜流場，當洪水漫過高灘地時，此種流況之水流結構與單槽渠流相當不同，不僅主深槽與洪水平原之流速及水深明顯不同，且深槽和洪水平原之間有複雜的三維性混合過程。一般河川之水深大多小於水面寬度，實際應用上仍是以二維模式演算橫斷面的水理及沖淤變化為較經濟常用的方式。以往蜿蜒複式河槽的研究主要針對固定的洪水平原探討水流結構，有關蜿蜒複式河槽主深槽變遷及河床沖淤方面需要更詳細的分析。因此本研究將探討水位較低、流量僅在深槽的單槽蜿蜒渠流，及洪水漫淹至洪水平原時之蜿蜒複式河槽溢岸流分別進行動床模擬。

(15)

3 蜿蜒河川可視為由數個不同曲率半徑及彎曲角度之彎道(bend)夾接只有數倍河寬長之近似直線段所構成。本研究先收集 CCHE2D 模式模擬彎道及蜿蜒河道相關之模擬案例，並選擇實驗室之 140°及 180°單一彎道試驗進行動床模擬，以了解 CCHE2D 模式在彎道二次流機制的功能及特色。而在單槽蜿蜒渠流方面，除了使用原有的二次流機制及動床輸砂模組，本研究亦使用能考慮岸壁因水流沖刷造成岸壁基腳破壞、河岸退縮的岸壁沖刷模組(bank erosion)，選擇正弦連續彎道岸壁破壞過程之實驗室案例，模擬連續彎道渠流岸壁後退河道側向變形的過程。在蜿蜒複式河槽中流量漸漸增加，水流漫淹至洪水平原時即成為溢岸流(overbank flow)。圖 1.1 為 Wormleaton et al. (2004)實驗三種類型蜿蜒複式河槽示意圖，圖中箭頭指示的粗線為決定洪水平原型態的河岸邊界，細線為蜿蜒之主深槽。此三種類型可由粗線條的河岸邊界區分，(A)為邊界順直型，(B)為邊界蜿蜒型，(C)為邊界平行於主深槽的狹窄型；本研究選擇常見之邊界順直、洪水平原寬廣的類型(A)為研究對象。因 CCHE2D 模式為水深平均之二維模式，並未考慮垂直方向的變化。在蜿蜒複式河槽溢岸流的研究中，先分析定床水理模擬的結果，了解二維模式模擬的複式河槽的水理特性，後針對三組不同流量及不同上游來砂條件進行動床模擬，探討不同流況下底床的變動情形，並應用於曾文溪下游之複式斷面河道，進行現場案例的模擬。

(16)

4 1.3 文獻回顧 1.3.1 彎道動床數模 在沖積河川中彎曲河段之床形變化，是一項困難且極為重要之課題， Engelund (1974)為探討渦流理論對彎道底床變動情形之影響，發展了一適用於穩態流場之彎道動床模式，並於模式中加入擬似二次流概念之徑向流速演算；Kikkawa et al. (1976)發展一適用於穩態流場之二維彎道動床模式，展示了以非耦合理論(分開求解水理及底床高程)模擬完全發展流況下之彎道底床變動情形；Struiksma et al. (1985)發展一適用穩態流之二維數值模式，並於底床載源運移演算時，加入曲率半徑、縱向及徑向之流速和底床坡降的影響； Shimizu and Itakura (1989)所發展適用於穩態流之二維數值模式，在沉滓運移演算中，此模式並無考量縱向底床坡降和懸浮載源之影響；Yeh and Kennedy (1993)發展一個適用於穩態流之二維數值模式，該模式乃著重於應用 MOM (moment of momentum)方程式，將曲率半徑、橫向底床坡降及縱向和橫向的流速等影響因子，導入沉滓運移之演算中。連(1999)則是考慮水理二次流效應對彎道流場之影響，並探討延散剪應力在彎道段之分布情形；Kassem and Chaudhry (2002)所發展一適用於非穩態流之二維數值模式，將曲率半徑、縱方向及橫方向之流速和底床坡降引入沉滓運移演算中，將有效應力與水理方程式結合；許(2002)發展一個適用於非穩態流之二維數值模式，於水理演算中導入曲率半徑之效應，於沉滓演算中；Duc et al. (2004)

(17)

5 發展一適用於非穩態流之二維數值模式，並將二次流傳輸效應反應於水深平均之ε-k 模式中的擴散係數，在底床載源的演算，包含曲率半徑、縱方向及橫方向之流速和底床坡降之影響。隨著電子計算機的功能快速成長，許多學者開始發展三維模式。Shimizu et al. (1990)發表考慮底床載及懸浮載的三維模式，模擬蜿蜒河道底床演變情形。Gessler et al. (1999)將 CH3D-SED 動床模式應用於密西西比河蜿蜒區段。Wu et al. (2000)使用 FAST3D 之非平衡輸砂公式模擬 180°彎道渠槽試驗。Ruther and Olsen (2005)將側向底床坡度修正項加入 SSIIM 模式，模擬 90°彎道的底床變化及流速分佈。Zeng et al. (2008)發表非靜力三維模式模擬非均值粒徑之實驗室彎道及直線道。然而，三維模式能直接提供水深方向的資訊但數值演算方法較為複雜，對於長時間的底床演變模擬，亦須耗費相當大量的計算時間，演算時所需要的邊界條件亦不易取得，因此目前尚較不具實用性。. 1.3.2 岸壁沖刷理論 關於岸壁之沖刷機制，歷來研究亦將非凝聚性沈滓與凝聚性沈滓分開考量。對於非凝聚性沈滓而言，一般有兩種主要方法用以計算沖刷岸壁所致淤積物之分布情形。Pizzuto (1990)、Li and Wang (1993, 1994)假設當岸壁坡度大於其安息角時，則破壞面上之土體，將下滑至岸壁基腳處堆積。Wiele

(18)

6

(1992)、Kovacs and Parker (1994)則是以泥砂連續方程式，將岸壁沖刷所致之堆積以側向沈滓通量處理。其中後者較符合沈滓顆粒間作用之機制；而前者則較符合邊坡穩定之理論(Wiele, 1992)。

對於凝聚性沈滓岸壁沖刷機制而言，Simon et al. (1991)、Darby and Thorne (1996)假設由岸壁所沖刷下來之物質將迅速堆積至岸壁基腳，並均勻地分布於近岸壁區域。而 Mosselman (1992)則假設這些沖刷物質均勻分布於近岸之計算格點。Osman (1985)、Alonso and Combs (1986)並未提及這些堆積物確切分布之方法，僅提及這些物質於很靠近岸壁基腳處堆積。Borah and Bordoloi (1989)則是根據局部之泥砂輸送能力，使用一線性分布假設。

近年關於岸壁沖刷機制之研究有 Sun et al. (2001a, b)將岸壁沖刷機制納入蜿蜒河道變遷之考量，Wu et al. (2004)引用 Osman and Thorne (1988)平面崩坍機制於其一維不平衡輸砂模式中模擬網狀河系，Jang and Shimizu (2005) 應用移動邊界調適坐標系統模擬寬淺渠道，Lauer and Parker (2005)以洪水平原至深水河槽沈滓之淨轉移量估算岸壁變化情形。

數模研發與應用方面，Nagata et al. (2000)假設岸壁遵守一簡易質量守恆之破壞行為，將該機制建置於一水平二維非平衡輸砂模式，模擬連續彎道之岸壁沖刷過程，及直線渠道底床交替砂洲(alternative bar)現象。Darby et al. (2002)以凝聚性沈滓之岸壁沖刷機制為研究對象，模擬崩坍後土堆坍落於

(19)

7

岸壁基腳(bank toe)情形，以兩組實驗案例及一組現場案例進行模式驗證。 Jang and Shimizu (2005)模擬直線渠道之交替砂洲過程，並比較不同岸壁土壤強度之模擬成果差異。Wang et al. (2008)模擬黃河下游縱向及側向河道變遷過程，同時比較岸壁可沖蝕係數、岸壁土壤臨界剪應力參數之敏感度，得知增加土壤臨界剪應力將明顯減緩岸壁沖刷情形。Jia et al. (2010)提出河岸坍塌的三維水沙數值模型，模擬實驗室河岸侵蝕的過程及長江中游河岸石首河彎雙層沉積結構的變遷趨勢。姜(2011)發展並應用河岸退縮預測數值計算模式，模式中採用未飽和地下水流理論，利用數值計算獲致河岸土體內部孔隙水壓隨時間變化之分布，藉以建構河岸邊坡塊體破壞機制，但未考慮河川底床沖淤過程及河寬變化對河岸退縮的影響。 1.3.3 複式河槽 蜿蜒複式斷面水流研究大多於試驗室水槽進行，至今除英國 SERC (Science and Engineering Research Council)之水槽模型為中尺度外，其他模型皆為小尺度，SERC 試驗分為三階段：A 階段試驗著重於直線及斜向 (Skewed)定床複式斷面；B 階段試驗為平面蜿蜒的定床形態；C 階段則為直線和平面蜿蜒的動床試驗。

Smith (1978)發現於高流量時較直線複式渠流有較高之水位，據此推論蜿蜒複式渠道產生額外阻力。Mckeogh & Kiely (1989)研究直線及蜿蜒複式

(20)

8

斷面水流特性，發現蜿蜒複式斷面水流紊流強度(turbulence intensity)較直線複式斷面水流大；最大紊流強度發生於連接槽(crossing reach)起訖點斷面之主深槽內側；蜿蜒複式斷面之洪水平原可傳輸較多的水流。

Lorena (1992)比較蜿蜒渠道及直線複式渠道溢岸流(overbank flow)之橫向流速和縱向流速之比值，認為渠道之幾何型態會影響二次流的強度。 Toebes & Sooky (1967)、Shiono & Muto (1993)、Willetts & Hardwick (1993) 等發現蜿蜒複式斷面渠流於溢岸和未溢岸(inbank)時，彎道頂(apex)點處之二次流發生機制不同，前者由剪力造成，後者則由離心力所導致而成。

Sellin et al. (1993)發現水流從洪水平原流入主深槽時，其水流將ㄧ分為二，如圖 1-1 所示；其中溢漫於主深槽上方之水流稱為溢漫流(overtopping flow)，捲入主深槽下方之水流稱為捲入流(plunging flow)，且從主深槽處亦有強大脫離流(vigorous expulsion flow)脫離至洪水平原上。Ervine et al. (1993) 認為影響蜿蜒複式斷面渠流能量損失因素有：(1)底床阻力；(2)底床形狀； (3)河道平面形狀；(4)螺旋流(水流之螺旋運動)(spiral motion flow)；(5)橫向剪應力；(6)蜿蜒度；(7)洪水平原植披等。

Willetts & Hardwick (1993)發現水流阻力，隨主深槽蜿蜒度增加而增大；主深槽與洪水平原交界面處之相互作用所產生阻力，與主深槽之幾何條件有關。Shiono et al. (1994)研究蜿蜒洪水平原邊壁型態對於水流結構的影響，

(21)

9

指出二次流發展於溢岸時，其旋轉方向與蜿蜒單槽之方向相同，由以上結果可知洪水平原邊壁之相位，對二次流之水流結構有所影響。Shiono & Muto (1998)等量測蜿蜒複式渠流流速，發現在主深槽滿岸水深下方之主流，沿著主深槽而流動；在主深槽滿岸水深上方之主流，則平行於洪水平原邊壁之方向而流動。並由可視化之觀測，發現主深槽於低流量時，表面水流之路徑隨著主深槽而流動；而主深槽於高流量時，其表面水流之流動路徑則平行於洪水平原之邊壁。Wormleaton et al. (2004)探討定岸洪水平原不同流量下水流結構及主深槽的底床變化情形，實驗使用的三種類型蜿蜒複式河槽如圖 1.1 所示。發現在溢岸流的情況下決定流速分佈及底床型態的兩個主要因素為(1)主深槽與洪水平原在直線段區因水流方向不同產生的橫向剪應力； (2)深槽蜿蜒產生的離心力。此兩個力為相反方向，會互相抵消。而圖 1.1 中(C)種類的蜿蜒複式河槽因主深槽與洪水平原水流方向相同，直線段區的橫向剪應力減弱，深槽蜿蜒的離心力的作用相對增加。李(2006)與姜(2008)分別施行光滑及粗糙洪水平原之水工模型試驗，並藉此探討蜿蜒複式渠流之主流流速剖面、主流等流速線及二次流之變化；並藉由可視化之觀測以瞭解表面水流之變化。賀(2009)研究粗糙洪水平原蜿蜒複式渠流流速分佈並與光滑洪水平原條件之試驗結果比較，其結果表示出洪水平原糙率提升僅影響低相對水深(洪水平原水深/主深槽水深)洪水平原流速分佈。對於主深槽以及高相對水深洪水平原之流速分佈，洪水平原

(22)

10 糙率提升不具影響力。近年來也有許多學者利用實驗成果發展數值模擬。李(2007)以水工模型試驗及二維水深平均模式分析深槽蜿蜒之複式斷面直渠的水理特性，由量測結果之分析可看出橫向底床水平的蜿蜒深槽，最大流速係在彎道處接近凸岸，在直線段則約在中間。Zhao (2007)將非靜力三維模式應用於實驗室蜿蜒複式河槽，模擬結果之流場圖與實驗結果相符，並展示二次流及底部剪應力的分佈。Ghani (2010)使用三維模式 FLUENT 模擬蜿蜒複式河槽，比較標準 k-ε、RNG k-ε和 RNG k-ω紊流模式，模擬結果顯示標準 k-ε模式能模擬完整溢岸流的流場結構。

(23)

11 圖 1.1 三種類型蜿蜒複式河槽示意圖(Wormleaton, 2004)

(A)

(B)

(C)

LF Floodplain Edge RF Floodplain Edge LF Floodplain Edge RF Floodplain Edge LF Floodplain Edge RF Floodplain Edge

Downstream

(24)

12

(25)

13

第二章

CCHE2D 模式理論及數值方法

CCHE2D 模式是一套水平二維明渠流動床模式，由美國國家計算水科學及工程中心(NCCHE)所開發，適用於變量流(unsteady)與紊流(turbulent flow)流況，此模式亦具有處理均勻和非平衡輸砂條件之能力。此外，模式包含三種紊流封閉演算法(turbulent closure schemes)用以解決不同流況和不同精度，可以處理河道之乾溼點問題，並可供使用者在非流動邊界設定全滑移、無滑動和指數型滑移之條件。以下摘錄由 NCCHE 所發表之 Wu (2001)、 Jia (2001a)、Jia (2008)等相關報告，說明模式水理、輸砂模式理論、岸壁沖刷理論、模式功能限制及代表性檢定驗證案例。 CCHE2D 模式水理模組採用交錯格網求解連續方程式之水位，特別之速度校正法來求解系統方程式；變量流模擬採用隱式時間推進法求解，乾濕點問題之處理採用移動邊界法；紊流模式部分，除了兩種零方程渦黏滯模式，還可考慮 k-ε 紊流模式；懸浮載傳輸為求解對流傳輸方程式，推移載傳輸採用非平衡輸砂方程式，且考慮底床坡度與彎道二次流效應。 2.1 水理模式 許多明渠流都為淺水波問題，垂直方向之運動通常可忽略，對於明渠水力問題而言，水深積分之二維控制方程式通常都允許合理及有效地被採用，在卡式座標下水深積分之動量方程式可表示如下：

(26)

14 v f h y h x h h x g y u v x u u t u Cor bx xy xx                                 1 (2-1) u f h y h x h h y g y v v x v u t v Cor by yy yx _ _                               1 (2-2) 式中，u 與 v 為 x 與 y 方向水深積分後之速度分量，t 為時間，g 為重力加速度，為水面高程，為水之密度，h 為局部水深，fCor為科氏力參數，xx、 xy  _、_yx_、_yy_{為水深積分之雷諾剪應力，}_bx_、by為底床剪應力，由於此版本模式不考慮風力，因此忽略水表面之剪應力，自由水面之高程藉由求解水深積分之連續方程式，假設為水面高程， 為底床高程，連續方程式可表示如下： 0             y vh x uh t t   (2-3) 由於底床改變量遠小於水理之變化，此方程式廣泛被應用在計算二維模式之自由表面高程，當底床高程受到沖淤影響急速改變時， /t_項不為零且應考慮之。式(2-1)、(2-2)中之雷諾剪應力取決於 Bousinesq 之假設，與水深積分流場之主要應力變化渦黏滯係數有關： ij i j j i t j i ij k x u x u v u u      3 2 ) (            (2-4)

(27)

15 渦黏滯係數為流體之函數，且與流體性質有關，在 CCHE2D 模式中採用兩種方法計算渦黏滯係數，第一種方式為使用水深積分拋物線渦黏滯係數公式： h u C A v_t  _xy _s _* (2-5) 式中u*為剪力速度，為馮卡門係數(0.41)，h 為水深，Axy為渦黏滯係數之校正係數，預設值為 1 且其值介於 1-10 間。除此之外，可採用深度積分之混合長度模式計算渦黏滯係數： 2 2 2 2 2 2 2 _                                       z U x v y u y v x u l v_t (2-6) 式中， h d h dz h z z h l 1  1  1 (1 )  0.267 0           



 (2-7) 沿著垂直方向之水深積分速度梯度可解釋為底床表面之紊流效應，(2-6) 式中決定之渦黏滯係數在均勻流況下可忽略為零，在缺乏其他項之情況下，渦黏滯係數可視為與均勻流相同，假設流速沿著水深方向為對數分佈，其垂直梯度可表示為： z u z U  *    (2-8) 將式(2-8)積分，可得水深垂直梯度之積分，可表為：

(28)

16    h u C dz z h u dz z U h z U m zo * * 1 1 _ _     



(2-9) 式中 U 為總速度，Cm為係數，Cm不直接計算，此係數可以某種方式給定，使得式(2-8)能夠取代式(2-12)中之所有水平方向速度梯度，此係數給定為 2.34375。除了上述兩種零方程渦黏滯模式，還可考慮採用 k-ε 紊流模式： kV k t k t _P _P y k y x k x y k v x k u t k                      _     ] [ ] [ (2-10) V t t P k c P k c y y x x y v x u t _ _               _ _ _                   2 2 1 ] [ ] [ _(2-11) ] ) ( ) ( 2 ) ( 2 [ 2 2 2 x v y u y v x u v x u u u P _t j i j i _                  (2-12) 2 4 3 , h U C P h U C PkV k V   _  _ _ _(2-13)    c c c C c C v u c U f f k f 3/4 2 2 2 6 . 3 , 1 , ) (      (2-14) 式中c_=0.09、k=1.0、=1.3、c1=1.44、c2=1.92，渦黏滯係數可由下式計算：   2 k c v_t  (2-15)

(29)

17 kV P 、PV增加於二維模式中，用來模擬由於底床摩擦所造成之紊流能量與消散。 2.2 輸砂模式 對輸砂模式來講，通常可分離為推移載與懸浮載，或依泥砂來源角度，可區分為河床質載與沖洗載。因此，CCHE2D 模式使用兩種概念計算總載之輸砂行為，其一為分離推移載與懸浮載，其二合併為河床質載(一般地形變遷動床模式模擬時忽略沖洗載對床形影響)，前者之控制方程式如下所示。懸浮載之對流傳輸方程式：



 





k k



sk d syk k s sxk k s k k k C C D y C h y D x C h x y VhC x UhC t hC                                                    ) ( (2-16) 式中Ck為某一粒徑 k 懸浮載之水深平均濃度；Ck為某一粒徑 k 平衡條件下懸浮載之水深平均濃度，s為泥砂之擴散係數；d為懸浮載之延散係數

(non-equilibrium adaptation coefficient)；Dsxk、Dsyk為水深方向非均勻流速與

泥砂濃度所造成之消散通量 (the dispersion terms due to the effect of

(30)

18 非平衡推移載輸砂方程式可表示如下：  





_ _ 0 1 ) ( _ _ _         k b bk bk by bk bx bk _q _q L y q x q t c    _(2-17) 式中為推移載區域之厚度；cbk為推移載區域之平均濃度；bx、 by 為推移載運動之餘弦函數，通常假設沿著底床剪應力方向，但當考慮彎道二次流效應與陡坡之重力影響時需做調整；qbk為某一粒徑 k 推移載之實際輸砂率；qbk為粒徑 k 推移載輸砂能力；L 為泥砂非平衡調適長度。 底床改變可由下式計算：   sk k k  bk b k k b m q q L C C t z p      _          1 1  _* (2-18) 式中pm 為河床質孔隙率；



zb t



k為某一粒徑 k 之底床改變率。 CCHE2D 模式中有四種輸砂公式可計算Ck與qbk，此四種輸砂公式之

挑選經由數十個實驗與現場資料檢定驗證而得(Wu and Wang, 2003)，可考慮非均勻輸砂之各種效應，此四種輸砂公式又可寫成以下一般化形式：  k  pbkCk C ; qb_k  pbkqbk (2-19) 式中pbk為混合層之床質粒徑分佈，  bk q 與C_k分別為某一粒徑 k 之潛勢推移載、懸浮載輸砂量。

(31)

19 床質粒徑分佈通常沿著垂直方向變化，因此不沖刷層以上之河床可再切割為數層，最上層定義為混合層，此層與移動之泥砂直接參與交換，在混合層中之粒徑分佈可以下式計算：                        t z t p t z t p _m _b bk k b bk m   ) * ( (2-20) 式中，m為混合層厚度；zb /t為底床總變形率，由  



_



 



N k b k b t z t z 1 累加計算，N 為代表粒徑總數；當_m/tz_b /t 0，式中 * bk p 可由混合層之床質分佈pbk替代之，反之當m/tzb/t0， * bk p 為混合層下層粒徑 k 所佔 之百分比。

此外，CCHE2D 模式可模擬凝聚性沈滓(Wu and Wang, 2004b)，可考慮植生影響之水流與輸砂作用(Wu et al., 2005)，考慮彎道之二次流效應輸砂 (Wu and Wang, 2004a)，並且能模擬岸壁基腳沖刷與破壞過程。

2.2.1 二次流機制 在自然的河道中經常出現彎道段及蜿蜒的型態，當水流流經不同曲率的彎道時，水流受到離心力的影響產生了二次流。二次流在水深上半部指向外岸，在水深下半部指向內岸，造成泥砂由外岸向內岸的橫向運移。底床改變後水流狀態隨之改變，兩者之間互相影響，二次流主導著彎曲段河道的底床變化，因此在模擬彎曲河道時不能忽略二次流造成的底床變化。

(32)

20 但水深平均模式並沒有考慮垂直方向的變化特性，因此二維模式大都使用經驗式來推估二次流效應，希望能合理地展現泥砂在彎道河段的運移。而這個部分最重要的參數即是主要水流方向與底床附近剪應力的角度。目前的 CCHE2D 模式中使用 Engelund (1974)提出的公式推估： (2-21) 式中的 r 為縱向格網的曲率半徑。根據 Engelund (1974)的研究成果，此推估公式的誤差為 3%。圖 2.1 展示了泥砂在底床運移的情形。泥砂在底床橫向坡度β的底床因重力作用往下游移動，而二次流則是趨向與底床橫向坡降不同方向。當兩者互相抵消達到平衡之後泥砂隨著主流移動，決定縱向的運移路徑。 2.2.2 輸砂公式 目前版本之 CCHE2D 模式中，共有四種輸砂公式可供使用者依不同需求來選取使用，分別為 SEDTRA module、Wu et al’s formula、Modified Ackers and White’s formula、Modified Engelund and Hansen’s formula，以下對各公式加以介紹。

1、SEDTRA module (Garbrecht, Kuhnle and Alonson, 1995)

本模組可依不同底床粒徑，來選取合適之輸砂公式求取河床質載輸砂量，當床質粒徑介於 0.01mm~0.25mm，使用 Laursen’s (1958) formula；介於

(33)

21

0.25mm~2.0mm，使用 Yang’s (1973) formula；介於 2.0mm~50.0mm，使用 Meyer-Peter and Mueller’s (1948) formula。

河川底床質載輸砂量可由  



k k k t p C C _* 式來計算，其中C_*t為河床質載輸砂能力(ppmw)，C_*k為某 k 粒徑之輸砂能力(ppmw)，

p

_k為某一粒徑 k 佔總 粒徑的百分比。

同時為了模擬非均勻輸砂中底床床質之隱顯效應(hiding and exposure)，

應用(2-20)式計算出某 k 粒徑之臨界水流強度，其中d_k為粒徑 k 之平均值，

d

_m 為底床質之平均粒徑，而 x 為一經驗參數。 x m k k ek d d d d         (2-22) 本輸砂模組藉由 Brownlie 之實驗數據(1981)來進行檢定驗證，共有 5,263 組實驗室資料及 1,764 組現地觀測資料，其中現地資料皆是在相同或極相似狀況下進行測量。

2、Wu, Wang and Jia’s (2000) Formula

此輸砂公式應用隱顯參數(hiding and exposure factor)來建立非均勻輸砂之啟動門檻、推移載及懸浮載，將推移載及懸浮載加總即為所需之河床質載輸砂量，其輸砂公式表示如下：

(34)

22





2.2 3/ 2 , 3 0.0053 1 / 1 bk b ck s bk k q _n n p gd     _ _  _ _    _ _         (2-23)





1.74 3 0.0000262 1 / 1 sk ck k s bk k q _U p gd      _ _       _ _       (2-24)





0.03 _k _s ck  dk     (2-25) m ek k hk p p  _{ } _   (2-26) 1 N j hk bj j k j d p p d d   



(2-27) 1 N k ek bj j k j d p p d d   



(2-28) 其中

q

_bk及

q

_sk分別為某粒徑 k 之推移載與懸浮載單位寬度輸砂量(m2 /s)， bk

p

為某粒徑 k 所佔總粒徑之百分比， n為糙度係數曼寧 n 值，_n,_為相對糙 度曼寧 n 值，



_b為底床剪應力，



_ck為某粒徑 k 啟動之臨界剪應力，k為某 粒徑 k 的隱顯參數，



_k為某 k 粒徑之沉降速度，m 為實驗參數(模式設定為 0.6)，本輸砂公式之粒徑適用範圍介於 0.01mm~250mm。 上述公式在均勻輸砂之檢定驗證上，是從 Brownlie (1981)之實驗數據挑出 1,859 組河床質載資料來進行比較；對於推移載輸砂量，又再從 1,859 組

(35)

23

資料中挑出 1,345 組推移載資料來進行比較調整。

在非均勻輸砂之檢定驗證上，挑選自 Toffaleti (1968)之實驗數據，分別是來自三組不同實驗團隊：Nomicos、Einstein-Chien、Vanoni-Brooks 於 Rio Grande River、Middle Loup River、Niobrara River 及 Mississippi River 的實驗結果。

3、Modified Ackers & White’s formula (Proffitt & Sutherland,1983)

本輸砂公式應用 Bagnold 之水流功率概念及相似性理論推導而成，但由於原始公式會超量預測輸砂量，HR Wallingford (1990)將原始公式加以修正，使粒徑適用範圍更為廣泛。









1 * , 1/ 2 32 log 10 / / 1 n n gr k k k s k U V F h d gd             _ _   (2-29) k  為校正係數(correction factor)，表示如下：其中dk為某粒徑 k 之平均值，du為相對粒徑。 , , 1 m gr k gr k F G C A    _  _   (2-30)





1.30, 0.53log / 1.0, 0.40, k dk du    _    / 3.7 0.075 / 3.7 / 0.075 k u k u k u d d d d d d    

(36)

24 * , / n k gr k bk k s C h U G p d  V    _ _   (2-31) 其中Ck即為總輸砂濃度(ppmw)。

本公式經與 Proffitt and Sutherland 之實驗數據比較，其輸砂量跟實驗數據相似度較高，但跟現地案例結果比較，其精準度較差。

4、Modified Engelund & Hansen’s 1967 Formula (with Wu et al’s correction factor, 2000) 本公式應用 Bagnold 之水流功率概念及相似性理論推導而成，並導入 Wu et al., 2000)輸砂公式之校正係數來調整隱顯作用在非均勻輸砂中之影響，公式如下：





5/ 2 ' * 0.1 k k k f    (2-32) ' 2 2gSd f U  ；_*_k ₀/_



 _s/



d_k_； 0.45 ek k hk p p  _{ } _   (2-33)





3 * / / 1 k qt k pbk s gdk   _    _ (2-34) 其中*k為無因次底床剪應力，k為隱顯效果校正係數，k為無因次輸砂量，qt k* 為單位寬度河床質載輸砂量(m2_/s)。

(37)

25 2.2.3 輸砂參數特性說明 以下針對 CCHE2D 模式中之相關參數進行其物理特性說明，瞭解各參數之物理意義做為模式參數檢定驗證的基礎。 1. 調適長度 L (adaptation length) 調適長度為河床上平衡至非平衡輸砂過程之特徵長度，L 通常與案例之推移載輸砂過程相關，由推移載過程至形成床形(bed forms)之距離通常即為調適長度大小。若河床多由砂漣(sand ripples)所構成，L 可為泥砂平均跳躍長度或砂漣長度；若河床多由砂丘(sand dunes)所構成，L 可為 5~10 倍之水深；若河床多由交替砂洲(alternate bars)所構成，L 可為交替砂洲長度或 6.3 倍之河道寬度；若研究河段為蜿蜒變遷型河道，L 可為蜿蜒波長或 12 倍河道寬度；若調適長度小於格網大小，則可能產生極大計算誤差。在天然河川之工程應用上，平均河寬為調適長度檢定之建議初始值，使用者可慢慢增加調適長度至模擬結果達到最佳。 2. 調適係數α (adaptation coefficient) 在水深平均模式中，靠近底床之泥砂濃度 Cb 與水深平均濃度 C 存在某種簡易比例關係(Cb = α*C)，式中α即為調適係數。而α與輸砂過程中許多因子相關，包含沈降速度、紊流作用、床形、斷面形狀等。當調適係數較小時，代表懸浮載與底床之作用關係較小，在一維模式中，α=1 時，代

(38)

26 表嚴重沖刷；α=0.5 時，代表介於沖刷與淤積之間；α=0.25 時，代表嚴重淤積。在二維模式中，上述理論仍可適用，但需注意α仍為檢定參數，且與案例輸砂條件相關，其為不大於 1 之正值，模式檢定時，建議可從α=0.1 開始進行測試，若推移載與懸浮載同時計算時，使用者應先檢定調適長度 L，再進行調適係數α之檢定。 3. 施密特數 (Schmidt number) 施密特數為無因次係數，為動黏滯係數與擴散係數之比值，用來描述同時有動量擴散與質量擴散之流體，可校正懸浮載泥砂之擴散程度。根據一般實驗數據，模式中施密特數預設值為 0.5，對於複雜流況或高含砂水流 問題，施密特數建議範圍約在 0.5~0.01 間。 2.3 岸壁沖刷理論 一般來說岸壁沖刷模擬包含岸壁面沖刷、岸壁基腳沖刷 (bank toe erosion)、岸壁材料破壞等過程，不同破壞過程皆與以下幾種參數相關，如泥砂粒徑組成、岸壁材料凝聚性、岸壁坡度、泥砂同質性、土壤濕度、地下水水位及岸壁高度等，且岸壁沖刷速率又與河道水流之剪應力、水深、曲率等因子有關。雖然數值模式為目前進行岸壁沖刷研究較有力之工具，但岸壁沖刷預測仍受到以下之限制：(1)難以完整收集分析所需之基本資料， (2)並非所有岸壁沖刷過程皆能正確瞭解，(3)水流模擬之正確性受到水理模

(39)

27 式方法論與計算機計算容量之限制。 CCHE2D 模式具有計算岸壁沖刷過程之功能，模式可計算由水流剪應力與岸壁塌陷所造成之岸壁土壤破壞機制，採用格網延展(mesh stretching) 技巧處理岸壁沖刷後之格網變形過程，模式計算流程如圖 2.2 所示。 2.3.1 岸壁破壞模組

引用 Osman and Thorne (1988)之理論，假設岸壁破壞是由近岸區域之力平衡所控制，因其簡易之假設，模式只需少數參數即可應用於實際案例。由於 Osman and Thorne (1988)提出之岸壁破壞機制僅適用於一維渠道模式，使用水平二維模式計算岸壁沖刷須計算靠近岸邊之剪應力，同時也須考慮計算格網點之變形，因二維模式計算側向剪力分布時須合理計算格網點分布，才能適當反映流場與地形變化。平面破壞機制之土堤示意圖如圖 2.3 所示，方程式如下： d r s F F f  _(2-35)      sin tan tan 2 sin 2 ' 2 2            W H y H F d t d (2-36)     

 2 tan tan cos tan

sin ) ( 2 2 '2              H y C H y H F_r d coh d ₍_2-37₎ 式中，fs為安全係數；Fd與 Fr為岸壁土堤崩坍下滑力與抵抗摩擦力；Wt 崩

(40)

28 坍土體重量； 崩坍平面傾角； 岸壁土體比重；H 岸壁高度；yd 張力裂縫深度； '  H 未發生沖刷之岸壁高度； 岸壁傾斜角度；Ccoh 土壤黏性係數；





摩擦角。當安全係數之值小於 1 時，表示土堤抵抗摩擦力小於崩坍下滑力，土堤即沿著傾角之破壞平面發生崩坍。崩坍平面傾角可由下式計算：





           _    _ _  tan 1 tan 2 1 1 2 K H H (2-38) 式中，K 為張力裂縫指數，為張力裂縫與岸壁高度之比值。 一般岸壁崩坍後產生之土質初期會堆積至岸壁基腳位置，若水流強度夠大，則土質將隨著水流作用而帶往下游，對大型河川來說，崩坍土質對水流流況影響不大，但對小型河川則可能會造成水流流況之擾動。在現階段 CCHE2D 模式中，崩坍土質假設為推移載之額外增量，此近岸區域之額外輸砂增量會使推移載濃度增大，由於岸壁崩坍發生之時間間距相較一般輸砂模組大許多，此崩坍之額外增量可視為均勻輸入。

在 Osman and Thorne (1988)之模式中，初始岸壁坡度在發生崩坍後，則產生新的岸壁坡度，並維持此岸壁坡度之穩定條件。若研究之岸壁為已發生崩坍且維持穩定一段時間，因此可合理假設其岸壁坡度為已知，且僅需計算其平行退縮過程，在此假設前提下，岸壁側向退縮距離 CD 為：

(41)

29  tan H H CD   (2-39) 而新舊岸壁高度H與 ' H 之比值可由下式計算： ] 4 ) ( [ 2 1 2 1 2 1 2         H H (2-40)    

1 cos sin cos2 tan (2-41)

H C K coh     2 2(1 ) (2-42) 利用上述理論計算岸壁崩坍時，應進行相關模式驗證，以確認模式適用性。表 2.1 表 2.2 為一般現場案例岸壁相關參數之參考範圍。 2.3.2 岸壁面沖刷模組 本研究採用水平二維模式計算河道流場，因此岸壁面之沖刷處理必須與 Osman and Thorne (1988)之一維概念有所不同，CCHE2D 模式由水理模式計算靠近岸壁之剪應力，其剪應力與局部區域流場中之流速、水深、曼寧糙度值、水面坡降等相關，而上述變數沿著河道側向之分布會因位置而有所不同，因此水平二維模式之岸壁沖刷計算上，與沿著河道斷面剪應力維持定值之一維模式有所區隔。由於岸壁面上之糙度會影響水流，靠近岸壁之流速會呈現小於河道中心之流速。

(42)

30 與臨界剪應力差異之比值，如下式： c c dB dW      (2-43) 式中，dW 為岸壁側向沖刷率，dB 為初始岸壁沖刷率，可由下式計算： c e dB223104_c 0.13 (2-44) 其中，τ 為斷面之剪應力，τc為岸壁臨界剪應力，可由岸壁本身材質套用經驗公式與查表計算，如圖 2.4 所示，其中 SAR 為土樣鈉吸附比、CONC 為孔隙流體鹽類濃度、∆ε0為電界質擴散。岸壁臨界剪應力範圍約介於 0.2~100 dyne/cm2，。 Partheniades (1965)使用下式計算黏性材料之沖刷率： a c k(₀  )    (2-45) 式中， 為沖刷率，k 為可沖刷係數，τ0為斷面之剪應力，a 為常數。Hanson and Simon (2001)採用約 200 組之現場資料分析出 k 的函數形式如下： ] / [ 10 * 2 . 0 6 0.5 m3 NS k  _c (2-46) 從上式中發現，黏性材料的沖刷率與臨界剪應力之減少成某種比例關係。 2.3.3 格網延展模組 CCHE2D 模式假設岸壁線為沿著外圍格網之邊界線，當岸壁沖刷發生

(43)

31 時，岸壁可沿著格網進行延展。一般有兩種不同概念計算岸壁沖刷後之底床幾何改變情形，概述如下： 1、變形格網當河道岸壁沿線受到岸壁沖刷而移動，水平計算格網點會變形以擴張原河道範圍。此概念對於水流、輸砂及岸壁沖刷計算上較有效率，但須花費較多工作量於格網變形進行隨時間變化之處理。 2、固定格網水平計算格網點在模擬過程中固定不變，計算區域設定足夠範圍使高灘地區域模擬時為乾區，岸壁沖刷發生時將沖刷岸壁區域與高灘地，使原乾區成為河道中之濕區。此概念不需要重新分配格網，但需要較大計算區域以調適計算格網之空間，此將會增加計算機使用量。

CCHE2D 模式採用第一種概念，圖 2.5、圖 2.6 為 Nagata et al. (2000) 實驗渠道格網延展前後之變形過程，其格網變形並非在一次計算疊代下就能完成，而是分成數次微小岸壁沖刷疊代來進行延展。當格網延展時，計算範圍即會改變，包含所有距離、方向與格網位置等，因此計算範圍內之所有變數皆須重新進行記算。使用者僅需定義研究範圍中左右邊界位置，並給定相關岸壁參數，在計算過程中模式可根據水流與輸砂條件自動校正邊界之格網，以計算岸壁沖刷過程。

(44)

32 舉例來說，岸壁沖刷發生之時距相較一般水流與輸砂過程大，即水流與輸砂演算數個步階後，才考慮岸壁沖刷情形，若該時距內岸壁沖刷所造成之格網延展非常顯著，在重新計算水流條件前，必須適當內插前一時刻與延展後之模擬結果，因此水流、輸砂與岸壁沖刷間之時距可設定不同。 2.3.4 岸壁沖刷模組功能限制 由於 CCHE2D 模式為水平二維模式，因此模式不適用於非水深平均二維流場之問題，如海洋流場、水庫流場等；同時對於突變區域(sharp changing zone)之模擬精度較差，如有極大落差之大峽谷等；岸壁沖刷部分，不適用懸臂土體岸壁崩塌問題(見圖 2.7)。 2.4 模式功能與限制 CCHE2D 模式可模擬穩態及非穩態明渠流、亞臨界、超臨界、混合流，大尺度之天然河川與小尺度之實驗渠道都經過檢定驗證程序，可應用於實驗室或天然河川之定量流、變量流、紊流流場等，亦包含複雜天然環境之自由液面流、輸砂、水質、污染傳輸等機制。表 2.3 表 2.4 將 CCHE2D 模式功能與限制進行彙整。 2.5 代表性檢定驗證案例 本節彙整目前國內外與彎道、蜿蜒河道相關之模擬案例，以便了解

(45)

33

CCHE2D 模式在實驗室條件及現場狀況下的表現，彙整於。以下整理自 Wu

and Wang (2004a)、Jia and Wang (2001)、「美國國家計算水科學及工程中心

河道變遷模式之引進及應用研究」(2009)。 1、180 °U 型彎道水理模擬此案例使用 De Vriend (1979)之實驗佈置模擬 180 °U 型彎道，渠道寬 1.7m，彎道中心線曲率半徑=4.25m，正弦彎道水理模擬參數表(θ0=30°)列於表 2.5。圖 2.8 為 U 型彎道水理模擬格網圖，邊界部分略做加密，為 23× 49 之不均勻格網。圖 2.9 為 180 °U 型彎道水理模擬與實測流速比較圖，此流速為水深平均流速，圖中可以看出 CCHE2D 模式可以把此彎道水深平均流速的特徵完整的展現出來。圖 2.10 為彎道水理模擬與實測水位比較圖，在彎道段，外岸水位較內岸來的高，此圖可以看出模擬與實測資料趨勢相同，雖然存在些許誤差，推測誤差原因可能為水深平均模擬沒辦法表現二次流的現象。 2、正弦彎道流場模擬(θ0=30°) 圖 2.11 為正弦彎道水理模擬格網圖(θ0=30°)，模擬之主流與側向方向格網為結構性規則格網，格網點數為(17x42)，實驗渠槽之流速測量依循渠道之正弦函數形狀，正弦彎道水理模擬參數表(θ0=30°)列於表 2.7。圖 2.12 為模擬與實測流速比較圖，模擬結果靠近凸岸之流速有些許誤差，但其他

(46)

34 位置幾乎吻合，當水流靠近凸岸時有加速作用，此加速導致水流改變方向。 3、正弦彎道流場模擬(θ0=110°) 在河道蜿蜒發展的過程中，彎道曲率會隨這蜿蜒度(河流長度與河谷長度之比值)增加而增加，此模擬為θ0=110°的正弦彎道水理模擬，上游為定量流下游為固定水深。圖 2.13 為正弦彎道水理模擬格網圖(θ0=110°)，使用 17×106 之均勻網格，水理模擬參數表列於表 2.8。圖 2.14 為模擬與實測流速比較圖(θ0=110°)，模擬結果大致與實測資料符合。唯在接近彎道的內岸及遠離彎道之後的外岸處模擬結果略低於實測值，在彎道段較準確，在直線過度段略有不同，可能因為蜿蜒度較大所造成。 4、180 °彎道變量流動床模擬

此模擬係針對 Yen and Lee (1995)之彎道動床模型試驗。模型渠道為一 180°急彎，渠寬 1m，曲率半徑為 4m，初始底床坡度為 0.002。為變量流流況，上游入流歷線為三角形歷線，基流為 0.02m3 /s，基流水深為 0.054m，於 Run4 模擬案例中，洪峰流量為 0.053m3 /s，延時為 240 分鐘，洪峰發生時間為總延時之前 1/3 處，泥砂條件為非均勻砂，中值粒徑為 1.0mm，標準偏差為 2.5，計算格網點為(91x31)，時間間距為 1 分鐘，水理模擬參數表列於表 2.8。圖 2.15 為模擬與實測彎道部分底床等高線比較圖，圖 2.16 為四個斷面之底床變化剖面圖，計算之床形變化與實測值趨勢相同。

(47)

35 5、蜿蜒河道底床地形變化模擬圖 2.17 為 CCHE2D 模擬蜿蜒渠道底床變化過程圖，為θ0=110°之正弦連續彎道，初始地形為平坦底床。本案例使用 CCHE2D 模式中的二次流推估公式(Englund, 1974))進行動床模擬。模擬前期時，渠道凸岸淤積形成沙洲，凹岸沖刷形成深潭；模擬後期時，淤積及沖刷的變化增強，且下游部分發展較上游劇烈。另一特色是，沙洲和深潭的位置並不對稱於彎道頂點，而是稍微往下游平移。隨著砂洲的淤積漸漸增加後，出現部分流速為零的乾點，流路逐漸窄縮。此種河道變化的行為為合理的過程，可在實驗室及自然河川中觀察到。 6、德國內卡河南德水庫模擬南德水庫位於德國斯圖加特(Stuttgart)之內卡河上，為一個河槽型水庫，全長約 12km，由上游比迪堪(Bietigheim)至下游南德，在河心距 129.0km ~ 132.3km 處經由兩段 180°之急彎，且於主槽兩岸有大片之高灘地。自 1950 年起，由於河川上游建置一連串 13 個壩，使得水庫上游懸浮載入砂量減少。南德流量站之流量變化約介於 14.1m3 /s~1,650m3/s 間，平均流量約為 88.5m3/s，Enz 河匯流進入內卡河，平均流量約為 20m3/s。圖 2.18 為此區段不規則床形之變量流模擬，洪峰流量約為 1,644m3 /s，圖中藍色區塊代表岸壁高程高於水位之區域，即所謂之乾點區，從河心距 130km~131.2km 與

(48)

36 130.55km~130.8km 可看出主槽流場與高灘地區域交界之二次環流效應；此外，主槽中靠近內側岸壁之流速大於外側岸壁流速。 7、濁水溪名竹大橋至中沙大橋岸壁沖刷案例模擬範圍為彰雲大橋至中沙大橋。以 87 年 DEM 資料為初始地形，模擬民國 87 至 96 年串接之洪水事件，因此挑選流量大於 4,000cms 者進行模擬。由模擬結果可知，河道左岸之岸壁沖刷情形與實測相比有良好模擬結果，沖刷位置點可準確預測模擬，河道側向變遷沖刷之距離可達數百公尺以上，最大側向沖刷距離約有 800 m，亦與實測相符。圖 2.19 為岸壁沖刷實際發生位置配合斷面圖之比較，可見斷面 70 位置之模擬與實測良好，無論是實測 DEM 地形或大斷面之比較，都能準確模擬。

(49)

37

表 2.1 Ackers & White’s formula 公式係數(1973, 1990)

1973 1990 60 gr d  n=0.0 A=0.17 m=1.50 C=0.025 n=0.0 A=0.17 m=1.78 C=0.025 1dgr 60 1.00 0.56log _gr n  d 0.5 0.23 _gr 0.14 A d   1 9.66 _gr 1.34 m d  





2

logC 3.53 2.86 log d_gr logd_gr

1.00 0.56log _gr n  d 0.5 0.23 _gr 0.14 A d   1 6.83 _gr 1.67 m d  





2

logC 3.462.79 logd_gr 0.98 logd_gr

表 2.2 岸壁相關參數建議範圍(Simon et al., 2003)

參數名稱範圍平均值

γ (kN/m3₎ _20.2~23.0 _21.34

C_coh(kPa) 0~78.9 16.0

(50)

38 表 2.3 CCHE2D 模式功能彙整表 2.4 CCHE2D 模式功能與限制表 2.5 CCHE2D 模式代表性檢定驗證案例彙整表 CCHE2D模式功能與特點彙整 1. 考慮植生影響 2. 考慮河川彎道影響與岸壁沖刷 3. 乾濕點處理 4. 水質(磷、氮、泥砂吸附/ 5. 去吸附) 模擬 6. 污染傳輸模擬 7. 圖形化使用者介面 8. 穩態、非穩態之自由液面流於 9. 複雜天然環境模擬 10. 隱式時間推進有效元素法

11. (Efficient element method)

12. 亞臨界、混合流、超臨界流模擬 13. 非平衡、非均勻輸砂模擬 14. 考慮凝聚性沈滓 15. 考慮水工結構物：壩、堰、閘門 CCHE2D模式之限制 1. 不適用非水深平均二維流場之問題

2. 突變區域(sharp changing zone)之模擬準確度較差

3. 不適用懸臂土體岸壁基腳掏刷與崩塌問題 CCHE2D模式代表性檢定驗證案例 1. 180 °U型彎道水理模擬 2. 正弦彎道流場模擬(θ₀=30°) 3. 正弦彎道流場模擬(θ₀=110°) 4. 180 °彎道變量流動床模擬 5. 蜿蜒河道底床地形變化模擬 6. 德國內卡河南德水庫模擬 7. 濁水溪名竹大橋至中沙大橋岸壁沖刷案例

(51)

39 表 2.6 180 °U 型彎道水理模擬參數表表 2.7 正弦彎道水理模擬參數表(θ0=30°) 表 2.8 正弦彎道水理模擬參數表(θ0=110°) 表 2.9 180 ° U 型彎道動床模擬模擬參數表 Q (m3/s) B (m) hm (m) um (m/s) Re* Fr B/hm 0.18 1.7 0.1953 0.542 1513 0.392 5.1 Q (l/s) B (m) hm(m) Sb um (m/s) Re* Fr B/hm 2.1 0.4 3.2 1/1000 6.4 5250 0.086 12.5 Q (l/s) B (m) hm(m) Sb um (m/s) Re* Fr B/hm 2.1 0.4 3 1/1120 16.7 5000 0.095 13.3 B (m) Rc (m) Sb h0(m) D50 (mm) mesh △t (sec) 1 4 0.002 0.054 1 91x31 60

(52)

40 圖 2.1 泥砂在底床的運移(Engelund, 1974) 圖 2.2 CCHE2D 模式岸壁沖刷計算流程計算二維水理流場計算懸浮載與推移載輸砂計算岸壁沖刷與岸壁移動計算格網點延展與分布重新計算新格網下之初始與邊界條件計算初始與邊界條件計算底床變化

(53)

41

圖 2.3 岸壁土堤破壞示意圖(Osman and Thorne, 1988)

圖 2.4 岸壁臨界剪應力

d W

(54)

42

圖 2.5 CCHE2D 模式未延展前之格網

(55)

43

圖 2.7 懸臂土體破壞形式(模式之限制)

(56)

44

圖 2.9 180 °U 型彎道水理模擬與實測流速比較圖(De Vriend, 1979)

(57)

45

圖 2.11 正弦彎道水理模擬格網圖(θ0=30°)

圖 2.12 正弦彎道模擬與實測流速比較圖(θ0=30°)

(58)

46

圖 2.14 正弦彎道模擬與實測流速比較圖(θ0=110°)

(59)

47

圖 2.16 180°彎道變量流動床模擬底床變化斷面剖面圖 (Yen and Lee, 1995))

(60)

48 圖 2.17 CCHE2D 模擬蜿蜒渠道底床變化過程圖圖 2.18 內卡河之變量流模擬 3.3 2.9 2.5 2.2 1.8 1.4 1.1 0.7 0.4 0.0 130.0 km 131.2 km Velocity [ m/s ] 131.2 km 130.0 km 131.2 km Island 130.9 km 129.5 km

(61)

49 圖 2.19 CCHE2D 模擬彰雲大橋下游岸壁沖刷比較結果 A B C D CS70 CS78

(62)

50

第三章

模擬案例

前一章節介紹了 CCHE2D 模式理論基礎及模式限制，並彙整與彎道、蜿蜒河道相關的代表性檢定驗證案例。二次流是彎道中決定泥砂運移方向及底床型態的重要因素，但以往的案例較多是探討模式在單槽彎道的流速及水位，對河道底床變化敘述較少。故本研究在模擬蜿蜒複式河槽溢岸流之前，先模擬兩個單一彎道動床案例，了解 CCHE2D 模式的二次流機制在不同曲率彎道底床變動情形。而河道側向變形方面，「美國國家計算水科學及工程中心河道變遷模式之引進及應用研究」(2009)已有應用於濁水溪之現地案例，故本研究選擇條件相對清楚且有實測資料比對的實驗室案例為模擬對象。以下說明模式模擬各項案例的條件及參數，包含實驗室設計案例及現場案例模擬。實驗室案例分為三部分，140°及 180°單一彎道動床、正弦連續彎道及蜿蜒複式河槽溢岸流動床，現場案例則應用於曾文溪下游麻善大橋至國聖大橋河段。 3.1 實驗室設計案例 3.1.1 單一彎道動床 140°彎道案例採用 Struiksma (1983)於荷蘭臺夫特科技大學水利實驗室 (Delft Hydraulic Laboratory)所作之室內水槽試驗組別 T2，簡稱 DHL 實驗。