河道三維高含砂水流沉滓運移模式發展與應用

國

立

交

通

大

學

土木工程學系

博士論文

河道三維高含砂水流沉滓運移模式發展與應用

Development and Application of a Three-Dimensional Model for

Hyper-Concentrated Flow and Sediment Transport in Alluvial Channels

研 究 生 ： 鍾浩榮

指導教授 ： 楊錦釧博士

謝德勇博士

河道三維高含砂水流沉滓運移模式發展與應用

Development and Application of a Three-Dimensional Model for

Hyper-Concentrated Flow and Sediment Transport in Alluvial Channels

研 究 生 ：鍾浩榮

Student： Hau-Rong Chung

指導教授 ：楊錦釧

Advisor： Jinn-Chuang Yang

：

謝德勇

Te-Yung Hsieh

國 立 交 通 大 學

土木工程學系

博 士 論 文

A Dissertation

Submitted to Department of Civil Engineering College of Engineering

National Chiao Tung University in partial Fulfillment of the Requirements

for the Degree of Doctor of Philosophy

in

Civil Engineering July 2012

Hsinchu, Taiwan, Republic of China

河道三維高含砂水流沉滓運移模式發展與應用

河道三維高含砂水流沉滓運移模式發展與應用

河道三維高含砂水流沉滓運移模式發展與應用

河道三維高含砂水流沉滓運移模式發展與應用

研究生：鍾浩榮 指導教授： 楊錦釧 謝德勇 國立交通大學土木工程學系

摘要

摘要

摘要

摘要

為能探討高含砂水流運移現象並提供實際規劃應用參考，本研究採水平正交 曲線、垂向 sigma 座標系統之靜水壓淺水波方程式，發展一三維高含砂水流沉滓 運移數值模式。控制方程式以水平垂直分離概念，分為水深平均二維方程式與流 速差異量方程式，搭配三維連續方程式求解三維流場。沉滓運移控制方程式分為 三維質量傳輸方程式、作用層連續方程式與底床連續方程式。高含砂效應考量二 次式流變關係反應非牛頓流體特性，以狀態函數反應含砂濃度對密度之影響，並 採用考慮高含砂效應之懸浮載與底床載經驗式，反應高含砂水流沉滓運移與濃度 變化。 高含砂水流及沉滓運移模式應用上需率定之相關參數甚多，因此首先進行模 式參數敏感度分析，探討影響水理、濃度與底床沖淤模擬結果重要參數之權重。 另研選數組高含砂水流實驗案例進行模式測試，探討水理及沉滓運移特性。探討 重點包含：流變關係之阻力使潰壩湧波移動一段距離後發生停止運移之現象；彎 道水位超高變化受流變關係之影響；流變關係與密度對底床載運移量之影響；並 釐清落淤速度與剪力速度之比值，可判斷兩種紊流擴散係數分布於模擬三維懸浮 載濃度分布之適用性。最後以實際應用角度，建立高含砂水流效應下，定床水深、 流速與底床剪力相對於清水流之增量及以動床觀點探討底床載對底床沖淤之影 響，以簡易關係式之型式供工程規劃設計應用之參考。 關鍵 關鍵 關鍵 關鍵詞詞：詞詞：：： 高含砂水流、高含砂水流高含砂水流高含砂水流、、三維模式、三維模式三維模式三維模式、、、、流變關係流變關係、流變關係流變關係、、沉滓運移、沉滓運移沉滓運移沉滓運移、、懸浮載、、懸浮載懸浮載、懸浮載、、、底床載底床載底床載底床載

Development and Application of a Three-Dimensional Model for

Hyper-Concentrated Flow and Sediment Transport in Alluvial Channels

Student：Hau-Rong Chung Advisor： Jinn-Chuang Yang

Te-Yung Hsieh Department of Civil Engineering National Chiao Tung University

Abstract

To investigate the transport behavior of hyper-concentrated flow and as a result to provide as a reference for engineering planning of practical cases, a hydrostatic three-dimensional model for hyper-concentrated flow and sediment transport in alluvial channels was developed in this study. By following the vertical and horizontal splitting concept (VHS), the shallow water flow governing equations were split into two parts including the depth-averaged two-dimensional equations and velocity defect equation in vertical direction. The former one was transformed into orthogonal curvilinear coordinate system; the latter one was derived as the form of sigma coordinate. Incorporated with continuity equation, the three-dimensional velocity field can therefore be solved. Sediment transport governing equations include three-dimensional mass transport equation, active-layer continuity equation, and bed-layer continuity equation. The effects of hyper-concentrated flow were treated as follows: a quadratic rheological relation was used to reflect the characteristics of non-Newtonian fluid; a state function was used to reflect the influence of concentration to density; the empirical suspended- and bed-load formulae with hyper-concentrated flow effect were used for the sediment transport computation.

Sensitivity analysis was performed first to identify the weighting of parameters to be calibrated in the model. The influence extent induced by the parameters on water

flow, suspended sediment concentration, and bed evolution, thereafter were examined and justified. To further investigate the characteristics of hyper-concentrated flow and sediment transport, several sets of experimental cases collected from the literatures were simulated. The case of dam-break wave propagation of non-Newtonian fluid was simulated to demonstrate the limit of traveling distance of hyper-concentrated flow. A channel bend flow experiment with 40% volume concentration of sediment was studied to investigate the effect of hyper-concentrated flow on super-elevation of water surface. A criterion of ratio of falling velocity and shear velocity was numerically examined and justified based on experimental data as a baseline for choosing the proper distribution type of turbulent diffusivity for simulation of suspended-load movement. At last, simple and concise regression relations for the increments of water depth, velocity, bed shear stress, and bed change caused by the hyper-concentrated flow were established for application of engineering planning and design.

Key words：：：： hyper-concentrated flow, three-dimensional model, rheological relation,

目錄

目錄

目錄

目錄

摘要 ... I Abstract ... II 目錄 ... IV 表目錄 ... VII 圖目錄 ... VIII 符號表 ... XI 第一章 緒論 ... 1 1.1 研究動機 ... 1 1.2 文獻回顧 ... 1 1.2.1 理論機制部分 ... 2 1.2.2 數值模式部分 ... 4 1.3 研究方法 ... 6 第二章 理論基礎 ... 11 2.1 水理部分 ... 11 2.1.1 三維方程式 ... 11 2.1.2 二維水深平均方程式 ... 13 2.1.3 流速差異量方程式 ... 15 2.1.4 高含砂水流與紊流輔助關係式 ... 17 2.1.5 邊界條件 ... 20 2.2 沉滓運移部分 ... 21 2.2.1 基本控制方程式 ... 21 2.2.2 高含砂水流與紊流輔助關係式 ... 22 2.2.3 邊界條件 ... 26 2.3 假設條件 ... 27 第三章 數值架構 ... 33

3.1 水理部分 ... 33 3.1.1 求解架構 ... 33 3.1.2 數值差分 ... 37 3.2 沉滓運移部分 ... 40 3.2.1 求解架構 ... 40 3.2.2 數值差分 ... 42 3.3 模式演算流程 ... 44 第四章 模式參數敏感度分析 ... 48 4.1 參數彙整 ... 48 4.2 案例設計 ... 49 4.3 水理部分 ... 50 4.4 濃度部分 ... 51 4.5 底床沖淤部分 ... 52 第五章 模式功能測試與驗證 ... 59 5.1 水理部分 ... 59 5.1.1 湧波傳遞案例 ... 59 5.1.2 彎道案例 ... 61 5.2 沉滓運移部分 ... 63 5.2.1 底床載案例 ... 63 5.2.2 懸浮載案例 ... 65 5.2.3.1 紊流黏滯係數檢定與驗證 ... 66 5.2.3.2 紊流擴散係數檢定與驗證 ... 67 5.2.3 彎道流場對懸浮載濃度分布之影響 ... 71 第六章 高含砂水流效應影響程度分析與應用 ... 96 6.1 前言 ... 96 6.2 定床情況高含砂水流之影響程度分析 ... 97 6.2.1 無因次參數影響程度分析 ... 97

6.2.2 重要參數之關聯性分析與應用 ... 100 6.3 底床載效應對底床沖淤影響程度分析 ... 104 6.3.1 無因次參數影響程度分析 ... 104 6.3.2 重要參數之關聯性分析與應用 ... 106 第七章 結論與建議 ... 123 7.1 結論 ... 123 7.2 建議 ... 125 參考文獻 ... 127

表

表

表

表目錄

目錄

目錄

目錄

表 1.1 高含砂水流之流變關係彙整表 ... 9 表 1.2 高含砂水流流變關係實驗數據蒐集整理列表 ... 10 表 4.1 流變關係參數係數統計特性列表 ... 53 表 4.2 參數敏感度分析案例設計範圍表 ... 53 表 5.1 定床湧波扇形擴展案例停止運移時前端位置比較表 ... 76 表 5.2 彎道定床高含砂水流模擬參數表 ... 76 表 5.3 彎道定床高含砂水流模擬結果比較表 ... 76 表 5.4 高含砂水流底床載案例實驗參數表 ... 77 表 5.5 高含砂水流懸浮載案例實驗參數表 ... 78 表 5.6 懸浮載案例之流變關係參數列表 ... 78 表 5.7 紊流擴散係數分布差異分析之參數範圍表 ... 78 表 5.8 寬深比>5 動床實驗案例之ω_f /U_*參數表 ... 79 表 5.9 彎道設計案例條件列表 ... 79 表 6.1 定床影響程度分析參數統計特性 ... 109 表 6.2 定床影響程度分析設計案例表(1/2) ... 109 表 6.3 定床影響程度分析設計案例表(2/2) ... 110 表 6.4 定床影響程度分析應用案例參數表 ... 110 表 6.5 底床載影響程度分析沉滓與底床糙度統計特性 ... 110 表 6.6 底床載影響程度分析設計案例參數範圍列表 ... 111 表 6.7 無因次參數對 之影響程度迴歸係數表 ... 112 表 6.8 福祿數與流體比重對 之影響程度迴歸係數表 ... 112 * / b Max Z∆ H * / b Max Z∆ H

圖

圖

圖

圖目錄

目錄

目錄

目錄

圖 2.1 正交曲線座標轉換示意圖 ... 29 圖 2.2 σ 座標轉換示意圖 ... 29 圖 2.3 水深方向流速剖面示意圖 ... 30 圖 2.4 動床模式之懸浮載、作用層變數配置示意圖 ... 31 圖 2.5 落淤速度修正係數與顆粒雷諾數關係圖 ... 32 圖 3.1 水平格網控制體積法示意圖 ... 45 圖 3.2 交錯格網示意圖 ... 46 圖 3.3 三維流場與濃度離散之控制體積法示意圖(計算域) ... 46 圖 3.4 數值模式演算流程圖 ... 47 圖 4.1 模式參數對水深平均速度之敏感度係數柱狀圖 ... 54 圖 4.2 不同案例之無因次縱向流速剖面圖與基值案例 ... 54 圖 4.3 模式參數對水面速度之敏感度係數柱狀圖 ... 55 圖 4.4 模式參數對近底床速度之敏感度係數柱狀圖 ... 55 圖 4.5 模式參數對水深敏感度係數柱狀圖 ... 56 圖 4.6 模式參數對深度平均濃度之敏感度係數柱狀圖 ... 56 圖 4.7 不同案例之無因次濃度剖面圖與基值案例 ... 57 圖 4.8 模式參數對最大濃度垂直分布差異之敏感度係數柱狀圖 ... 57 圖 4.9 模式參數對底床沖淤量敏感度係數柱狀圖 ... 58 圖 5.1 直線道湧波傳遞之潰壩實驗案例初始水深與模擬結果圖 ... 80 圖 5.2 扇形擴展之潰壩案例水深模擬結果立體圖 ... 81 圖 5.3 扇形擴展之潰壩案例渠道中心線水深沿程變化圖 ... 82 圖 5.4 彎道定床高含砂水流案例水槽幾何條件與格網配置圖 ... 82 圖 5.5 彎道定床以側向水深檢定流變關係參數比較圖 ... 83 圖 5.6 曲率半徑 0.4 m 彎道定床高含砂水流水位超高比較圖 ... 83 圖 5.7 曲率半徑 0.6 m 彎道定床高含砂水流水位超高比較圖 ... 84

圖 5.8 曲率半徑 1.0 m 彎道定床高含砂水流水位超高比較圖 ... 84 圖 5.9 底床載案例之底床載通量比較圖 ... 85 圖 5.10 底床載案例流變關係效應對底床剪力之影響比較圖 ... 85 圖 5.11 底床載案例流變關係影響底床載通量比較圖 ... 86 圖 5.12 不同分布之紊流黏滯係數沿水深分布比較圖 ... 86 圖 5.13 垂直紊流黏滯係數對流速剖面影響比較圖 ... 86 圖 5.14 懸浮載案例(水利署 2011)流速剖面模擬結果(寬深比>5) ... 87 圖 5.15 懸浮載案例(水利署 2011)流速剖面模擬結果(寬深比<5) ... 87 圖 5.16 懸浮載案例(水利署 2008)流速剖面模擬結果 ... 88 圖 5.17 不同型態之紊流擴散係數沿水深分布比較圖 ... 88 圖 5.18 不同紊流擴散係數分布型態懸浮載濃度剖面之影響 ... 89 圖 5.19 不同紊流擴散係數分布型態對濃度影響之敏感度分析柱狀圖 ... 89 圖 5.20 懸浮載案例之無因次紊流擴散係數沿無因次水深分布圖 ... 90 圖 5.21 Coleman(1970)之無因次紊流擴散係數沿無因次水深分布圖 ... 90 圖 5.22 應用w_f /U_*判斷不同紊流擴散係數分布之濃度剖面模擬結果... 91 圖 5.23 拋物線分布紊流擴散係數對彎道濃度之影響 ... 92 圖 5.24 拋物線-常數分布紊流擴散係數對彎道濃度之影響 ... 92 圖 5.25 彎道入口處不同紊流擴散係數分布對彎道濃度影響比較 ... 93 圖 5.26 彎道中點處不同紊流擴散係數分布對彎道濃度影響比較 ... 93 圖 5.27 彎道出口下游不同紊流擴散係數分布對彎道濃度影響比較 ... 94 圖 5.28 彎道設計案例 BF1 底床沖淤等值圖 ... 94 圖 5.29 彎道設計案例 BF2 底床沖淤等值圖 ... 95 圖 5.30 彎道設計案例 BF3 底床沖淤等值圖 ... 95 圖 6.1 水深差異量與不同無因次參數關係圖 ... 113 圖 6.2 流速差異量與不同無因次參數關係圖 ... 114 圖 6.3 剪力差異量與不同無因次參數關係圖 ... 115 圖 6.4 T_m與ρm*對 相對差異分析圖 ... 116 * Maxτ

圖 6.5 高含砂水流增量與無因次賓漢剪應力 關係圖 ... 116 圖 6.6 高含砂水流剪力增量與 關係圖 ... 117 圖 6.7 高含砂水流定床案例應用水深與流速增量結果比較圖 ... 117 圖 6.8 以 計算高含砂水流剪力增量結果比較圖(案例 Af2、Af3 與 Af4) ... 118 圖 6.9 以 計算高含砂水流剪力增量結果 (案例 Af5) ... 118 圖 6.10 定床案例下高含砂水流效應影響程度之水深比較 ... 119 圖 6.11 定床案例下高含砂水流效應影響程度之剪力差異比較圖 ... 120 圖 6.12 坡度對 影響程度比較圖 ... 120 圖 6.13 與主控因子迴歸關係圖 ... 121 圖 6.14 不同案例之主控因子對 影響程度分類圖 ... 121 圖 6.15 底床載影響程度之驗證(Rickenmann, 1991 案例 12) ... 122 圖 6.16 底床載影響程度之驗證(水利署, 2008; 水利署, 2011) ... 122 m T * m ρ * m ρ m T * / b Max Z∆ H * b Max Z∆ * b Max Z∆

符號表

符號表

符號表

符號表

B=渠道寬度； C=懸浮沉滓濃度； C_{=懸浮沉滓水深平均濃度；} , a e C =近底床沉滓平衡濃度； a C =近底床沉滓濃度； c c =Chezy 糙度係數； 2 / f c C =g c = 摩擦係數； 50 d = 懸浮沉滓中值粒徑； D=水深； b D =河床沉滓代表粒徑； m D =編號第 m 組沉滓之粒徑； s D = 懸浮載向下通量； ms D =懸浮沉滓代表粒徑； m D_∗ =無因次顆粒粒徑； p E =作用層厚度； rms E =均方根誤差； s E = 懸浮載向上通量； F_r = 福祿數； f_ξρɶ、 f_ξρɶ _{= 主流與側方向受密度差異量影響之移流作用力；} g = 重力加速度； H =正常水深； H_m=高含砂效應下，最大差異發生處之水深； H_w =清水流下，最大差異發生處之水深； 1 h 、h₂ =ξ、η方向轉換係數； l K = 多項式流變關係之水深平均層流阻力係數； k_s = 粗糙高度；

κ =卡門係數(Von Kármán’s 係數)； M = 沉滓代表粒徑總數目； m =第 m 組沉滓編號 n=曼寧糙度係數； Pɺ= 動量方程式之壓力項； p=孔隙率； Q = 總流量； q = 單寬流量； m bi q =i 方向某一粒徑 m 之河床載通量； R=彎道任一點之曲率半徑； L R =彎道左岸之曲率半徑； h R =水力半徑； c r =彎道中心線曲率半徑； 0 S =渠道坡降； m s S =某一粒徑 m 之懸浮載源； g s =砂比重； E S =作用層源； w S =水面坡降； *m S =某一粒徑 m 水流挾砂能力； k T =輸送參數； 11 T 、T₁₂、T₂₂ =有效剪應力項； t=時間； u =主流方向深度平均速度； u=主流方向速度； U= 斷面平均速度； U_∗ =剪力速度； U_m =高含砂效應下，最大差異發生處之水深平均速度； U_w =清水流下，最大差異發生處之水深平均速度；

U_∗c =臨界剪力速度； v =側方向深度平均速度； v=側方向速度； w=z 方向速度 f w =沉滓沉降速度； fh w =濃度影響下沉滓沉降速度； z = 卡式座標系統之垂直方向座標 b z = 底床高程； s z = 水面高程； l z = 最接近底床之垂直格網厚度； Z_u = 懸浮沉滓傳輸係數； 1 α _、α2= 流變關係參數常數項係數； m β _{=第 m 組代表粒徑在作用層中所佔百分比；} 1 β _、β2 = 流變關係參數指數項係數； sH ε 、ε_sV= 懸浮載連續方程式中水平、垂直方向之紊流擴散係數； l υ = 層流黏滯係數； t υ = 紊流黏滯係數； H υ = 水平方向渦流黏滯係數； V υ = 垂直方向渦流黏滯係數； t µ = 動力黏滯係數 B µ =賓漢黏滯係數； m θ = 第 m 組代表粒徑承受之無因次剪力； cm θ = 第 m 組代表粒徑之無因次臨界剪力； ρ = _{含砂水流流體密度；} s ρ = 乾砂密度； B τ =賓漢降伏應力； m τ =高含砂效應下，最大差異發生處之底床剪力； w τ =清水流下，最大差異發生處之底床剪力；

1 b τ 、 2 b τ =底床剪應力在ξ、η方向之分量； 1 s τ 、τ_s2=風剪應力在ξ、η方向之分量； c τ = 泥砂顆粒啟動之臨界剪力； γ = 含砂水流比重量； ξ、η =平面上兩正交座標方向； σ = 水深分層之垂直方向座標，無因次水深； ω= σ _{座標方向之垂向速度；} t ∆ =時間間距； η ∆ 、∆ =ξ ξ、η方向兩相鄰格網點之格網間距。

上標

上標

上標

上標

n= ∆n t時刻之已知變數； 1 ( 1) n+ = + ∆n t時刻之未知變數； 1 2 ( 1) n+ = + ∆n t與n t∆ 間之未之變數；

( )

時間平均；

( )

水深平均；

( )

ɶ 差異量（流速為與深度平均流速之差，密度為與清水密度差）；

( )

′ 紊流擾動量。

下

下

下

下標

標

標

標

i = ξ方向格點編號； j = η方向格點編號； k = σ 方向格點編號； 1、2= ξ、η方向變數； p、e、w、n、s =通量所在格點位置； a = 近底床高度 5%水深處變數； b = 近底床變數。

第一章

第一章

第一章

第一章 緒論

緒論

緒論

緒論

1.1 研究動機

研究動機

研究動機

研究動機

台灣河流特徵多陡坡且流急，使水流蘊藏較高的沖刷潛勢，部分山區地質鬆 軟且地震頻繁，大量崩塌覆蓋各集水區，產生大量砂源。而近年來氣候變遷現象 明顯，極端降雨事件頻傳，時有山洪挾帶土砂而下，對河道屢生衝擊。各地不時 傳出土砂災害，且越趨嚴重，如民國 85 年賀伯颱風在台東太麻里溪沿岸造成大 量土石崩落；民國 90 年桃芝颱風在南投東埔蚋溪造成堤防毀損、民房沖毀；民 國 98 年莫拉克颱風在中南部、東部造成了大量土砂災難。 廣義而言，由於砂土在水流中的含量過高，導致密度、黏滯係數等物理特性 異於清水流情況，使水體的流動、沉滓運移特性受到改變，皆可稱為高含砂水流。 Takahashi et al. (2006) 提到，水體含砂量達體積濃度 2%以上時，流體特性開始受 到影響，因此本研究將主要針對體積濃度 2%以上之含砂水流進行探討。面對高 含砂水流可能逐漸形成河川輸砂的常態現象，本研究建構適用之高含砂水流動床 數值模式，RESED3D，配合實驗數據探討高含砂水流與沉滓運移之現象。

1.2 文獻回顧

文獻回顧

文獻回顧

文獻回顧

一般而言，高含砂水流之研究範圍多為河川中上游區域，含砂水流有較高比 例之粗顆粒泥砂，土石碰撞、摩擦效應明顯，且流動一段距離後粗顆粒泥砂會有 明顯落淤現象，Wan and Wang (1994)定義此區域為土石流(debris flow)運動區。在 土石流運動區下游，河道坡度較緩，泥砂顆粒較細，也就是顆粒間碰撞效應可以 忽略，較鬆散、細小之泥砂容易被水流帶往下游，Wan and Wang (1994)定義此為 高含砂水流(hyper-concentrated flow)運動區，Takahashi et al. (2006)稱為非成熟土 石流(immature debris flow)運動區。此區域既可能有凝聚性沉滓產生非牛頓流體 (non-Newtonian fluid)特性，也有非凝聚性沉滓以移流擴散效應傳輸，此區域為本

研究主要探討重點。就高含砂水流運動區而言，相關文獻可分為理論機制與數值 模式發展兩個部分，以下分別說明。 1.2.1 理論機制部分理論機制部分理論機制部分 理論機制部分 描述流體應力與應變關係的函數稱為流變關係(rheological relation)，高含砂 水流在水體含砂量高到某種程度時，水流流動特性可能由牛頓流體(Newtonian fluid)的流變關係轉變為非牛頓流體之流變關係。文獻中已有學者探討適用於高含 砂水流的非牛頓流體流變關係，Wan and Wang(1994)整理後作一比較，其成果如 表 1.1 所示。表中 88%的高含砂水流研究係以賓漢流(Bingham fluid)作為描述非 牛頓流體之流變關係。

賓漢流體之賓漢降伏應力(Bingham yield stress)、賓漢黏滯係數(Bingham viscosity) 為 高 含 砂 水 流 中 描 述 流 變 關 係 之 常 見 參 數 ， 稱 為 流 變 關 係 參 數 (rheological relation parameters )，通常與體積濃度有某種函數關係，如下列式(1.1) 與式(1.2)的指數函數。表 1.2 為根據文獻所蒐集、彙整之現場土樣流變關係參數， 其中α1、β1、α2、β2等係數可對應式(1.1)與式(1.2)。 1 1 C B e β α

τ

= (1.1) 2 2 C B e β µ =α (1.2) 其中 τ_B、µ_B分別為賓漢降伏應力與賓漢黏滯係數，C 為含砂體積濃度。 除流體特性外，沉滓運移機制也是高含砂水流的研究課題，可分為懸浮載 (suspended load)與底床載(bed load)兩部分做探討。高含砂水流之泥砂運移機制， 普遍以挾砂能力(sediment carrying capacity) (曹如軒，1987；張紅武與張清，1992； 劉峰與李義天，1997；費祥俊與舒安平，1998；Shu and Fei, 2008)作為主控水體 中含砂量的參數。挾砂能力係在對應的水流條件下，維持沉滓懸浮所耗費之水流

紊流能量，可視為平衡狀態(equilibrium condition)的水深平均濃度。前述文獻中用 以迴歸、建構挾砂能力經驗式之數據均係以單一代表粒徑(如佔粒徑百分比 50% 之沉滓粒徑)為主，又均採用現場資料，黏性沉滓含量對整體沉滓運移機制之影響 無法忽略，其造成的絮團(flocculation)結構則增加經驗公式的不確定性。水利署 (2008)、水利署(2011)與莊巧巧等(2011)採用不同粒徑、摒除黏性沉滓之沉滓，利 用循環水槽研析適用於台灣河川粒徑特性的挾砂能力經驗式。另一方面，也可藉 由解析解之濃度垂直剖面直接積分(Rouse, 1937; Zhang et al., 2001)，量化懸浮載 運移量，並且以經驗公式推估近底床濃度，作為濃度剖面在近底床處之邊界條 件。直接積分方法較能考慮濃度在深度方向不均勻分布之特性，但是僅考慮懸浮 參數(suspension parameter, Rouse, 1937)之解析解，難以反應濃度剖面受三維複雜 流場之影響。張紅武等(1994)推導濃度剖面解析解並探討高含砂水流懸浮載沿水 深之分布，但仍是在解析解流速剖面適用、恆定流(steady flow)與平衡狀態之假設 下，簡化三維濃度質量傳輸方程式而得到。 由前述懸浮載之經驗公式的發展回顧，可知相關研究頗豐。然而至目前為 止 ， 探 討 高 含 砂 水 流 密 度 、 流 體 特 性 影 響 的 底 床 載 運 移 機 制 不 多 。 Rickenmann(1991)以水槽實驗解析高含砂水流之底床載運移機制，考慮高含砂水 流效應之密度、底床剪力對底床載運移之影響，指出水體密度為影響底床載運移 之主控參數。然而，真正考量底床載效應並分析高含砂水流動床沖淤的研究在文 獻上少見，以黃河流域為例，由於坡度緩且粒徑較細，多為沖洗載，相對較粗的 沉滓也均屬於懸浮載的範圍。Yang et al.(1996)蒐集黃河流域現場沉滓數據並分析 適用於高含砂水流情況下的總載(total load)經驗公式，其數據中已經將沉滓區分 為水體的懸浮沉滓與在河床表面之沉滓，水利署(2011)將兩個區間的沉滓分別計 算懸浮參數(Rouse number)後，發現無論是懸浮沉滓或是底床沉滓之懸浮參數均 小於 2.0。根據錢寧與萬兆惠(1983)，在懸浮參數小於 2.0 可歸類為懸浮載。由此

可知，黃河流域之沉滓特性應可僅以懸浮載機制探討底床沖淤。但是，在坡度相 對較陡、空間與時間降雨均相當集中、流況變化劇烈的河川時，例如台灣部分河 段，尤其是在河川中上游，底床載運移量則可能貢獻較高的底床沖淤量(郭朝雄 等，1984-1988；水規所，2007)。 1.2.2 數值模式部分數值模式部分數值模式部分 數值模式部分 為探討河道維持底床穩定的條件，泥砂輸送(sediment transport)的計算一直是 數值模式探討之對象。由於沉滓運移之量化係依靠對應之水理條件、沉滓特性， 流場計算成為數值模式發展的首要重點。就現有電子計算機的發展而言，水深平 均二維水理模式在現場應用有快速、簡便的優點，且諸多應用在近年來已經獲得 許多研究之驗證，例如近幾年之 Molls and Chaudhry (1995)，Ye and McCorquodale (1997)，Lien et al.(1999)與 Hsieh and Yang(2003)等，均是應用在複雜幾何之渠道、 模擬諸多複雜流場，並得到良好的結果。為能夠得到更精確的流場分布，應用三 維模式於現場之數值模擬也逐漸發展，例如 Leschziner and Rodi (1979)，Sinha et al. (1998)，Wu et al. (2000)，Meselhe and Sotiropoulos (2000)與 Zeng et al. (2008)等， 均是考慮完整的三維流體動力(fully three-dimensional hydrodynamic)方程式，可以 解析局部流場之現象。但是若以時間成本考量，且所關心的現象著重於大範圍流 場、垂向流速遠小於水平方向流速之應用時，則水深方向靜水壓(hydrostatic)假設 之數值模式亦可提供合理之三維流速資訊。例如，在海洋與湖泊的研究當中，由 於溫度與鹽分之影響，導致水體密度沿深度產生分層(stratified)之現象，在研究範 圍廣情況下，又需要節省計算時間，漸多研究提出水平垂直分離演算的概念發展 三維數值模式，如 Lardner and Cekirge (1988)，Jin and Kranenburg (1993)，Wang (1994)與 Lin and Huang (2008)。水平垂直分離演算的概念主要是先以深度平均模 式求解水深與深度平均流速，再以三維子模式求解垂直方向分布的流速剖面。洪 聖翔(2011)發展一水平垂直分離演算概念之數值模式，並應用於彎道明渠流解析

二次流(secondary current)現象。

水理條件得到較精確流場資訊後，結合底床載、懸浮載或總載經驗公式作為 輔助關係式，可進一步發展動床數值模式。近年來由於數值計算能力大幅進步， 漸多研究直接以三維數值模式探討河道動床機制，如 Gessler et al.(1999)、Fang and Wang (2000)、Wu et al.(2000)與 Zeng et al.(2010)等。主要係因天然河道幾何形狀 多彎曲、漸縮漸擴等不規則側壁邊界，僅採用平面二維數值模式難免會由於水深 積分，而無法考慮垂向流場資訊，或是構造物附近的局部流場，難以合理計算二 次流效應或是垂向流速資訊，如此則無法合理的計算相關效應對動床沖淤之影 響。前述三維動床模式研究均著重在底床載(bed load)佔較高比例的沖刷案例，含 砂濃度較低。其中 Zeng et al.(2010)係以懸浮載為主控之案例(比例約佔總載之 50% 至 60%)進行三維水理與動床數值模擬，其所模擬的案例中，懸浮載比例較高， 認為複雜流場顯著影響懸浮載源通量沿水深之變化，相對於僅考慮近底床之懸浮 載源通量，較能合理反應懸浮載對底床沖淤之貢獻。在含砂水流較高之河川，計 算較完整、合理的懸浮載三維濃度資訊，除了有助於解析底床沖淤外，引水(intake) 防砂的水資源問題(Ruether et al. 2005)也可藉由三維濃度的分布資訊，做為引水設 施置放高程與平面位置之參考。 如前所述，水體含砂量增加後，水體的流動特性受到改變，O'Brien et al. (1993)、 Naef et al. (2006)、Liu and Huang (2006)與 Cetina et al. (2006)發展數值模 式，視高含砂水流為連續流體，假設運移過程為定床，以流變關係解析非牛頓流 體運動特性。描述不同流體的流變關係中，如牛頓流體的曼寧(Manning’s)公式(紊 流條件下)、土體內摩擦角流變關係(Voellmy,1955)與多項式(quadratic)流變關係 (O'Brien et al., 1993)，均是反應於淺水波(shallow water)方程式之底床剪力。Naef et al. (2006)以實驗數據比較不同流變關係之間差異，在考量曼寧公式(紊流效應) 時，如多項式流變關係，數值模擬結果較能符合高含砂水流流動情況。目前流變

關係理論發展漸豐，仍鮮少有研究加入動床機制共同探討，主要係因為在應用面 上，如 Canuti et al.(2002)、Liu and Huang (2006)與 Cetina et al. (2006)、Sosio et al. (2007)與 Armento et al.(2008)等，均主要假設沉滓與流體為完整混合且濃度與密度 均固定的均勻連續流體，高含砂水流之影響範圍為其研究重點。假設濃度不隨時 間空間改變是較為簡化的方法，越向下游則由於濃度不考慮改變機制，同時流體 流變關係也僅為濃度的函數時，無法藉由濃度改變反應合理的流場機制。在所關 心的課題若為影響範圍時，會得到保守之運移距離，若是考慮剪力對底床、結構 物之影響，則可能高估剪力。而改變濃度的合理方法，就是結合高含砂水流之沉 滓運移機制，考慮沉滓沿程交換，計算濃度改變量。

高含砂水流之動床模式如黃遠東等(2003)、Ni et al.(2004)與 Cao et at. (2006) 等，考量修正的沉滓落淤速度，以挾砂能力概念解析高含砂水流懸浮載通量，探 討黃河流域的高含砂水流洪水對底床沖淤之影響。但是 Ni et al.(2004)與 Cao et at. (2006)等係一維模式，而黃遠東等(2003)係水平二維模式，考慮高含砂水流之懸浮 載常主控整體沉滓運移機制，僅以一維或二維模式進行分析難免忽略流場較複雜 處之濃度資訊。van Maren et al.(2009)採用 Delft3D 靜水壓三維模式模擬黃河的高 含砂水流洪水事件，但是僅採用低濃度下適用之經驗式與黏性沉滓的沖淤經驗式 進行模擬，並未採用適用之懸浮載經驗式，同時現場沉滓特性十分複雜，因此難 以良好模擬現場底床沖淤趨勢。目前高含砂水流動床數值模式均尚未考慮流變關 係造成的底床剪力，及底床剪力對沖淤機制之影響。同時，底床載經驗式除了理 論發展外，鮮少整合至高含砂水流數值模式，針對底床載機制進行分析與應用。

1.3 研究方法

研究方法

研究方法

研究方法

本研究擬發展一數值模式，基於謝德勇(2002)之水深平均水理與沉滓運移模 式(RESED2D)及洪聖翔(2011)之靜水壓三維水理模式，建構一考慮高含砂水流效 應之三維動床模式(RESED3D)。針對高含砂水流非牛頓流體特性，採用多項式流

變關係，計算水深平均二維模式之底床剪力。密度採用狀態函數反應含砂濃度對 密度之影響，且考慮濃度在三維空間之不均勻分布特性，計算密度差造成之壓力 梯度(baroclinic pressure gradient)。模式以沉滓運移連續方程式計算高含砂水流的 濃度變化與沉滓運移量，並將底床載與懸浮載分開計算。選用文獻中考量高含砂 水流效應而發展之經驗公式，計算底床載通量與懸浮載源，以解析動床機制。 本研究所探討之高含砂水流，係針對體積濃度 2%以上之含砂水流，主要是 因為 Takahashi et al. (2006) 認為此濃度含砂水流已經對水流特性產生影響。體積 濃度更高的情況，費祥俊與舒安平(2004)提出一水砂重量的判斷式：

(

)

0 1 sC C ρ =ρ − (1.3) 其中ρs與ρ0分別為乾砂密度與清水流密度，在單位體積的高含砂水流中，若乾砂 比重為 2650 kg/m3_{，水流所挾帶的最大砂量則為 27.4%。這也是後續所模擬高含} 砂水流案例的體積濃度上限，超過此範圍則應考慮土石流之運移機制作探討。 綜合前述水理部分的流變關係、密度效應與高含砂水流沉滓運移機制，模式 發展考量了眾多參數之影響，為釐清模式參數之重要性，以參數敏感度分析的方 法，權衡數值模式中各參數對流場、濃度、底床沖淤之影響程度。 模式功能測試與驗證方面，首先解析高含砂水流流變關係對流場之影響，以 潰壩案例解析流變關係造成的阻力與流體自重平衡後，產生的停止流動現象，另 外也以彎道水位超高案例，解析超臨界流下，高含砂水流水位在平面空間不均勻 分布的特性。考量高含砂水流情況下適用之底床載經驗式，以實驗案例作測試， 驗證高含砂水流流變關係與密度影響下，模式模擬底床沖淤之能力。過去高含砂 水流動床模式研究多以現場案例測試數值模式之正確性，如 Ni et al.(2004)與黃遠 東等(2003)，但由於現場沉滓特性複雜，難以排除沖洗載效應並單獨探討懸浮載 受高含砂水流效應之影響。本研究以水利署(2008)、水利署(2011)與莊巧巧(2011)

之實驗案例，檢定影響參數並模擬三維濃度分布，探討影響三維濃度分布之沉滓 運移機制與紊流效應。最後為能展現三維模式模擬複雜流場之模擬能力，並解析 複雜流場對高含砂水流之影響，以假設之彎道案例解析二次流效應對濃度分布、 底床沖淤之影響。 本研究數值模式可同時提供二維與三維之流場、濃度資訊，考量實務應用上 仍以二維資訊較為普遍，例如規劃設計所需資訊或是沉滓運移經驗式之流速資 訊，均以水深平均或斷面平均速度為主要參數，又由於高含砂水流水理、動床之 分析工具仍屬少見，因此在研究的最後，以二維流場、沖淤資訊為基礎，提出一 套簡單的圖表關係，讓高含砂水流效應可由清水流角度或是在進行模式分析前， 得知高含砂水流效應之影響程度。定床部分，本研究為釐清高含砂水流效應對明 渠流之影響程度，以清水流模式為對照，分析高含砂水流與清水流差異的主控參 數，以期進一步瞭解高含砂水流效應對明渠流水深、流速與剪力之影響。動床部 分，考量底床載對河床穩定的重要性，解析底床載效應考慮與否對高含砂水流情 況下底床沖淤影響程度之主控參數。

1.

表 1.1 高含砂水流之流變關係彙整表

Material 流變關係 作者 年份

Loess deposite/water Bingham Y.Sha 1947 Deposits in reservoir at

estuary, in river/water

Bingham N.Qian & H.Ma 1958 Deposite in river/water Pseudo-plastic Y.Zhou 1963 Loess deposits/water Bingham Z.Wan et al. 1979 Loess deposits/water Bingham H.Zhang et al. 1980

Clay/water Bingham Pseudo-plastic

J.Dai et al. 1980 Loess deposits/water Bingham M.Wang et al. 1983

Bentonite/water, kaoline/water

Bingham Z.Wan 1982

Loess deposits/water Bingham X.Fei 1982, 1983 Slurry of debris flow Bingham Y.Wang 1982 Slurry of debris flow Bingham S.Sheng & S.Xie 1983

Deposits in bay, reservoir, kaoline,

limestone

Bingham C.Migniot 1968

Clay/water Bingham H.E.Babbit & D.H.Caldwell 1940 Kaoline/water Bingham D.G.Thomas 1961, 1963

Clay/water Bingham G.W.Govier & M.D.Winning 1948 Clay/water Bingham N.I.Heywood & J.F.Richardson 1978 Clay/water Bingham,

Crowley-Kitzes

D.G.Thomas 1962 Clay/water Bingham B.A.Firth & R.J.Hunter 1976 Clay/water Bingham K.Emeya 1970 Clay/water Bingham H.D.Weymann &

M.C.Chuang,R.A.Ross

1973 Clay/water Bingham H.D.Weymann 1965 Clay/water Bingham G.F.Brooks & R.L.Whitmore 1968 Clay/water Bingham V.M.Dobrychenko et al. 1975 Clay/water Bingham E.C.Bingham &T.C.Durham 1911 Clay/water Bingham H.Pazwash & J.M.Robertson 1975

表 1.2 高含砂水流流變關係實驗數據蒐集整理列表

數據出處 α1 β1 α2 β2

Dai et al. (1980) 0.26 0.1748 0.00075 0.1439 O'Brien and Julien(1988):1 0.0348 0.225 0.00373 0.22 O'Brien and Julien(1988) :2 0.0228 0.176 0.0002 0.268 O'Brien and Julien(1988) :3 0.0052 0.19 0.000505 0.195 O'Brien and Julien(1988) :4 0.0044 0.21 0.000502 0.231 O'Brien and Julien(1988) :5 0.0043 0.2 0.000441 0.218 O'Brien and Julien(1988) :6 0.0149 0.16 0.01017 0.128 O'Brien and Julien(1988) :7 0.001 0.228 0.003619 0.15 O'Brien and Julien(1988) :8 0.0009 0.232 0.000831 0.25 Coussot and Piau(1994) :1 0.031 0.244 0.00103 0.17 Coussot and Piau(1994) :2 0.043 0.22 0.00535 0.17 Coussot(1995) 0.0658 0.215 0.002856 0.175 Chen(1998): 1 0.0924 0.236 0.00052 0.301 Chen(1998): 2 0.0057 0.234 0.000315 0.155 Malet et al.(2005) :1 0.0066 0.21 0.020937 0.169 Malet et al.(2005) :2 0.0136 0.204 0.090642 0.161 Malet et al.(2005):3 0.0421 0.175 0.150151 0.153 Malet et al.(2005):4 0.0519 0.19 0.101 0.113 Sosio et al.(2007):1 0.0098 0.194 0.000143 0.21 Sosio et al.(2007):2 0.0015 0.207 0.000157 0.204 Martino(2003) 0.0059 0.17 0.005424 0.12 Malet et al.(2003):1 0.0138 0.207 0.025797 0.18 Malet et al.(2003):2 0.0083 0.212 0.030805 0.166 Malet et al.(2003):3 0.003 0.221 0.004702 0.171 Malet et al.(2003):4 0.0029 0.211 0.008558 0.157 Malet et al.(2003):5 0.0516 0.204 0.056075 0.141 Malet et al.(2003):6 0.0065 0.216 0.001119 0.198 Malet et al.(2003):7 0.0113 0.221 0.014546 0.174 Schatzmann et al.(2003 ) 0.1486 0.223 0.008385 0.185 神木村(徐昌益 1996) 0.0811 0.1372 0.000462 0.1124 東埔蚋溪(水利署 2008) 0.0157 0.17 0.000047 0.221 陳有蘭溪(水利署 2011) 0.00254 0.198 0.000824 0.176 濁水溪(水利署 2011) 0.0041 0.177 0.0032 0.263 太麻里溪(水利署 2011) 0.00004 0.225 0.247 0.095 1 1 C B e β α

τ

= 、 2 2 C B e β µ =α 修改自王志賢(2007)

第二章

第二章

第二章

第二章 理論基礎

理論基礎

理論基礎

理論基礎

2.1 水理

水理

水理部分

水理

部分

部分

部分

本研究以靜水壓為主要假設之三維動量方程式為基礎，採水平垂直分離演算 概念，分為水平二維方程式與流速差異量方程式(洪聖翔，2011)，並考量高含砂 水流造成的密度差與流變關係，建構高含砂水流水理模式。為能適當表達天然河 道不規則的幾何形狀，模式在水平方向採用正交曲線座標(orthogonal curvilinear coordinate)系統，如圖 2.1 所示，在垂直方向則採用 σ 座標系統(Blumberg and Mellor, 1983)，如圖 2.2 所示，如此能將不規則的計算區域轉換至矩形計算區域 求解，且能便利地處理渠道的側壁、自由液面及底床邊界。 2.1.1 三維三維三維方程式三維方程式方程式方程式 基於不可壓縮流之假設下，對那威爾-史托克司(Navier-Stokes)方程式取時間 平均後，控制方程式水平方向以正交曲線座標系統、垂直方向為卡式座標表示如 下： 連續方程式 連續方程式 連續方程式 連續方程式 2 1 1 2 (h u) (h v) (h h w) 0 z ξ η ∂ ₊ ∂ ₊ ∂ ₌ ∂ ∂ ∂ (2.1) 動量方程式 動量方程式 動量方程式 動量方程式 ξ 方向：

2 2 2 1 2 2 1 2 1 2 1 2 1 2 2 2 2 1 2 2 1 2 1 2 1 2 1 2 2 11 1 1 2 1 1 2( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 1 2( ) ( ) ( ) ( ) ( ) ( ) 1 1 ( ) h h h uv u v u u uv uw t h h z h h h h h h h h h u v u v u u v u w h h z h h h h h h P h h h h ρ ρ ρ ρ ρ ρ ρ ξ η η ξ ξ ρ ρ ρ ρ ρ ρ ξ η η ξ ξ τ ξ ξ ∂ ∂ ∂ ∂ ₊ ∂ ₊ ∂ ₊ ∂ _{+ + −} ∂ ∂ ∂ ∂ ∂ ∂ ∂ ′ ′ ∂ ′ ∂ ′ ∂ ∂ _′ ∂ _{′ ′} ∂ _{′ ′} + + + + + − ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = − + ∂ ∂ ɺ 12 1 22 2 1 12 1 2 13 1 2 1 2 (h ) (h h ) h h z h h h h τ τ τ τ η η ξ  ₊ ∂ ₊ ∂ ₊ ∂ ₋ ∂  _{∂ ∂}  _{∂ ∂}   (2.2) η 方向： 2 2 2 2 1 1 2 1 1 2 1 2 1 2 2 2 2 2 1 1 2 1 1 2 1 2 1 2 1 22 2 1 2 1 1 2( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 1 2( ) ( ) ( ) ( ) ( ) ( ) 1 1 ( ) h h h uv v u v v uv vw t h h z h h h h h h h h h u v v u v u v v w h h z h h h h h h P h h h h ρ ρ ρ ρ ρ ρ ρ η ξ ξ η η ρ ρ ρ ρ ρ ρ η ξ ξ η η τ η η ∂ ∂ ∂ ∂ ₊ ∂ ₊ ∂ ₊ ∂ _{+ + −} ∂ ∂ ∂ ∂ ∂ ∂ ∂ ′ ′ ∂ ′ ∂ ′ ∂ ∂ _′ ∂ _{′ ′} ∂ _{′ ′} + + + + + − ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = − + ∂ ∂ ɺ 12 2 11 1 2 12 1 2 23 1 2 1 2 (h ) (h h ) h h z h h h h τ τ τ τ ξ ξ η  ₊ ∂ ₊ ∂ ₊ ∂ ₋ ∂  _{∂ ∂}  _{∂ ∂}   (2.3) z 方向： 在一般淺水的天然河道，垂向流速遠小於水平方向流速時，垂直方向之動量 方程式可用靜水壓分布來簡化， 0 P g z ρ ∂ + = ∂ ɺ (2.4) ，採用連續方程式(2.1)可求解垂直方向流速w。 以上諸式中，ξ、η、z為三維正交曲線座標方向，其中ξ、η為水平方向，z 為垂直方向；下標 1、2、3 分別代表物理量在ξ、η、z方向代號；h₁、h₂分別為 ξ、η方向之轉換係數；u、v、w分別為ξ、η、z方向流速；g為重力加速度； t為時間；τ 為層流剪應力；( )表時間平均；( )′ 表時間平均瞬時擾動量；ρ 為密 度；Pɺ為水體內任一點壓力。 為能考慮高含砂水流與清水流之密度差異ρɶ，空間中各點需考量垂直方向的 密度差異量積分所造成的壓力增加量，將式(2.4)之Pɺ靜水壓與壓力增量合併為：

0 ( ) ( , , , ) s z s z Pɺ=ρ g z − +z g

∫

ρ ξ ηɶ z t dz′ _(2.5) 2.1.2 二維二維二維水深平均方程式二維水深平均方程式水深平均方程式水深平均方程式 將 2.1.1 節之三維水理控制方程式(2.1)、(2.2)、(2.3)轉換為σ座標系統，該系 統之轉換方式為： b s b z z z z σ = − − (2.6) 再利用萊布尼茲法則(Leibuniz rule)對深度方向積分，加上水面與底床運動邊 界條件(kinematic boundary condition)，並取深度平均，可得水平二維水理控制方 程式： 連續方程式 連續方程式 連續方程式 連續方程式 1 2 ( 2 ) ( 1 ) 0 D h h h uD h vD t ξ η ∂ ₊ ∂ ₊ ∂ ₌ ∂ ∂ ∂ (2.7) 動量方程式 動量方程式 動量方程式 動量方程式 ξ 方向： 2 1 2 0 1 2 1 2 1 2 1 1 0 0 1 1 1 1 2 1 2 11 12 11 12 22 1 2 1 2 ( ) 1 1 1 2 s s b h h u u u v u uv v f t h h h h h h g z g D D d d h h h h h T T T T T Dh h Dh Dh D ρ ξ σ ρ ξ η η ξ ρ ρ _{ρ σ σ σ} ξ ξ σ ξ τ τ ξ η ξ ξ η _{∂ ∂ ∂ ∂ ∂}  + + + − + =   ∂ ∂ ∂ ∂ ∂   ∂  _{∂ ∂ ∂}  − −  +  ∂  ∂ ∂ ∂   _{∂ ∂ ∂}  _{∂ ∂ −} +  _∂ + _∂ − _∂ + _∂ + _∂ +  

∫ ∫

ɶ ɶ ɶ (2.8) η 方向：

2 2 1 0 2 1 1 2 1 2 1 1 0 0 2 1 2 2 1 2 1 22 12 22 12 11 1 2 2 1 ( ) 1 1 1 2 s s b h h v v v u v uv u f t h h h h h h g z g D D d d h h h h h T T T T T Dh h Dh Dh D ρ η σ ρ η ξ ξ η ρ ρ _{ρ σ σ σ} η η σ η τ τ η ξ η η ξ _{∂ ∂ ∂ ∂ ∂}  + + + − + =   ∂ ∂ ∂ ∂ ∂   ∂  _{∂ ∂ ∂}  − −  +  ∂  ∂ ∂ ∂   _{∂ ∂ ∂}  _{∂ ∂ −} +  + − + + + ∂ ∂ ∂ ∂ ∂  

∫ ∫

ɶ ɶ ɶ (2.9) 式中， 1 1 1 2 2 0 0 0 1 2 1 1 1 2 2 2 2 1 2 2 0 0 0 1 2 1 1 ( ) ( ) 2 1 ( ) ( ) ( ) f u d u u d uv uv d t h h h h h uv uv d u u d v v d h h ρ ξ ρ σ _ξ ρ σ _η ρ σ ρ σ ρ σ ρ σ η ξ ξ ∂ ∂ ∂ = + − + − ∂ ∂ ∂  ∂ ∂ ∂  + _ − + − − − _ ∂ ∂ ∂  

∫

∫

∫

∫

∫

∫

ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ (2.10) 1 1 1 2 2 0 0 0 2 1 1 1 1 2 2 2 2 2 1 1 0 0 0 1 2 1 1 ( ) ( ) 2 1 ( ) ( ) ( ) f v d v v d uv uv d t h h h h h uv uv d v v d u u d h h ρ η ρ σ _η ρ σ _ξ ρ σ ρ σ ρ σ ρ σ ξ η η ∂ ∂ ∂ = + − + − ∂ ∂ ∂  ∂ ∂ ∂  + _ − + − − − _ ∂ ∂ ∂  

∫

∫

∫

∫

∫

∫

ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ (2.11) 1 2 2 11 ₀ ( 11 ) T =

∫

Dτ −ρuɶ −ρu′ dσ _(2.12) 1 2 2 22 22 0 ( ) T =

∫

Dτ −ρvɶ −ρv′ dσ (2.13) 1 12 ₀ ( 12 ) T =

∫

Dτ −ρuvɶ ɶ−ρu v d′ ′ σ _(2.14) 以上諸式中，D 為水深；z_s為水面高程；z_b為底床高程；τ_si為 i 方向水面剪 應力；τbi為 i 方向底床剪應力；( )表水深平均；( )ɶ 表物理量與平均值之差異量 (流速為三維流速與水深平均流速之差值:uɶ= −u u、 _vɶ= −_{v v}，密度為含砂水流與 清水密度之差值:ρ ρ ρɶ = − 0)；下標 s、b分別為水面及底床代號；τb1、τb2為主流 與側方向之底床剪力； f_ξρɶ 與 f_ηρɶ為主流與側方向受密度差異量影響之移流作用 力；T_ij為有效剪應力項(effective shear stress term)，包含層流剪應力(laminar shear stress)、紊流擴散(turbulent diffusion)與延散效應(dispersion)。

2.1.3 流速差異量方程式流速差異量方程式流速差異量方程式 流速差異量方程式 將 2.1.1 節 之 三 維 水 理 控 制 方 程 式 (2.2) 、 (2.3) 代 入 靜 水 壓 假 設 ， 並 令 u = +uɶ u、_v = +_vɶ _{v (如圖 2.3)，使三維流速成為流速差異量(uɶ}、_{vɶ)與水深平均流} 速(u 、v )之函數，對所得到之三維控制方程式減去與水平二維水理控制方程式 (2.8)、(2.9)重複之項，即可得到流速差異量方程式如下： ξ 方向方向方向方向：：：： 1 0 1 1 2 2 1 2 1 2 1 2 2 2 2 2 1 2 1 2 1 2 1 2 1 2 1 2 2 1 2 h u u u u u v u v u u uv uv uv t h h h h h h h h h h h vv h v u u u v u uv h v h h h h h t h h h h h h u uv u t h ρ ω ξ ξ η η σ η ρ ξ ξ ξ η η ξ ρ ρ ξ   ∂ ∂ ∂ ∂ ∂ ∂ ∂  _{+ + + + + + + +}    _{∂ ∂ ∂ ∂ ∂ ∂ ∂}        ∂ ∂ ∂ ∂ ∂ ∂ ∂ − _∂  _∂ + _+  _∂ + _∂ + _∂ + _∂ − _∂       ∂ ∂ + + + ∂ ∂ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ 2 2 1 2 1 1 1 0 1 1 13 ( 1 1) 1 ( n ) s b u uv uw f u h t h h D g g D D d D D d d h h Horizontal Diffusion i D D ρ ξ σ σ ρ ρ σ σ σ η σ ξ η ρ ρ σ _σ ρ ρ σ _{σ σ} ξ σ ξ ξ σ ξ τ τ τ _ξ σ  ∂ _{− +}∂ ∂ ₊ ∂ ₊ ∂ ₊    ∂ ∂ _ ∂ ∂ ∂   _{∂ ∂ ∂}   _{∂ ∂ ∂}  = − _ + _ + _ + _ ∂ ∂ ∂ ∂ ∂ ∂     ∂ − + − + ∂

∫

∫ ∫

ɶ ɶ ɶ ɶ ɶ ɶ ɶ (2.15) η 方向方向方向：方向：： ： 2 0 2 2 1 1 1 2 1 2 1 2 2 2 1 2 1 1 2 1 2 2 1 1 2 1 2 2 2 1 2 h v v v v v u v u v v uv uv uv t h h h h h h h h h h h uu u v v v u v uv h u h h h h h t h h h h h h v uv v t h h ρ ω η η ξ ξ σ ξ ρ η η ξ ξ η ρ ρ ρ η   ∂ ∂ ∂ ∂ ∂ ∂ ∂  + + + + + +  + +  _{∂ ∂ ∂ ∂ ∂ ∂ ∂}  _{ }      ∂ ∂ ∂ ∂ ∂ ∂ −  + _+  + + + −  ∂   ∂ ∂ ∂ ∂ ∂  ∂ ∂ ∂ + + + ∂ ∂ ∂ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ 2 2 1 1 1 1 0 2 2 23 ( 2 2) 1 ( n ) s b v uv vw f v t h h D g g D D d D D d d h h Horizontal Diffusion i D D ρ η σ σ ρ σ σ σ ξ σ η ξ ρ ρ σ _σ ρ ρ σ _{σ σ} η σ η η σ η τ τ τ _η σ ∂  ∂ ∂ ∂  − +_∂  _∂ + _∂ + _∂ +     _{∂ ∂ ∂}   _{∂ ∂ ∂}  = −  _∂ + _{∂ ∂}  +  _∂ + _{∂ ∂}      ∂ − + − + ∂

∫

∫ ∫

ɶ ɶ ɶ ɶ ɶ ɶ (2.16) 式中，

1 1 2 2 u u v v w t h h h h z σ σ σ σ σ σ ω ξ ξ η η ∂ ∂ ∂ ∂ ∂ ∂ = + + + + + ∂ ∂ ∂ ∂ ∂ ∂ ɶ ɶ (2.17) 1 D D t t σ −_σ ∂ _{= −} ∂ ∂ ∂ (2.18) 1 zb 1 D D D σ _σ ξ − ξ − ξ ∂ ∂ _{= − −} ∂ ∂ ∂ ∂ (2.19) 1 zb 1 D D D σ _σ η − η − η ∂ ∂ _{= − −} ∂ ∂ ∂ ∂ (2.20) 1 D z σ − ∂ = ∂ (2.21) 2 1 2 11 12 22 1 2 2 12 11 2 2 2 2 2 1 1 1 1 2 1 1 2 1 1 2 2 1 2 1 ( ) 2 2 2 2 1 1 1 2 H H H H H h h h Horizontal Diffusion in T T T Dh h T T u h u Dh Dh h h h h h h h v v h h h h h h ξ ξ η ξ υ υ ρυ ρ η ξ ξ ξ ξ ξ ρυ ρυ ξ η η ξ  _{∂ ∂ ∂}  = − _ + − _ ∂ ∂ ∂       ∂ ∂ ∂ ∂ ∂ ∂ − − +   + + ∂ ∂ ∂  ∂   ∂  ∂  _∂   ∂ ∂ ∂ + _∂  _∂ + _∂  _∂     1 2 2 1 2 3 1 1 2 2 2 2 2 2 2 1 1 2 2 1 1 2 2 1 1 2 H H H H u h u h h h h u v u h v h R h h h h h h h h ρυ η η ρυ _ρυ ρυ η ξ η η η ξ   _∂   _∂   +   _∂  _∂                 _{∂ ∂} ∂ ∂ ∂ ∂ +  +   +   − + ∂    ∂   ∂  ∂ ∂ ∂ (2.22) 其中 3 u R 為 sigma 轉換之殘餘項，以式(2.22)等號右邊第 8 項為例，做座標轉換可 型如： 2 3 2 ... ... ... u u u u R u u u u σ ξ ξ ξ ξ σ ξ σ σ σ ξ ξ σ ξ σ ξ σ ξ ξ   ∂ ∂ ∂ ∂ ∂ ∂ = _{∂ ∂} + =_∂ _∂ +_{∂ ∂} +       ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = _ + _+ _ + _ + ∂ ∂ ∂ ∂  ∂ ∂ ∂ ∂ ∂ 下式也以相同方式表示為 3 v R 。

1 2 1 22 12 11 1 2 2 22 12 1 2 2 2 2 1 2 2 1 2 2 2 2 2 1 2 1 1 2 1 ( n ) 2 2 2 2 1 1 1 2 H H H H H h h h Horizontal Diffusion i T T T Dh h T T v h v Dh Dh h h h h h h h u v h h h h h h η η ξ η υ υ ρυ ρ η ξ η η η η ρυ ρυ η ξ ξ ξ  _{∂ ∂ ∂}  = − _ + − _ ∂ ∂ ∂       ∂ ∂ ∂ ∂ ∂ ∂ − _∂ − _∂ + _∂  _∂  + _∂ + _∂      _∂   ∂ ∂ ∂ +   +   ∂  ∂ ∂  ∂ 1 1 2 1 2 3 2 2 1 2 2 2 2 2 1 2 1 2 2 1 1 2 2 1 1 2 H H H H v h u h h h h v v u h u h R h h h h h h h h ρυ ξ η ρυ _ρυ ρυ ξ ξ η ξ ξ η   _∂   _∂   +       ∂ ∂               _{∂ ∂} ∂ ∂ ∂ ∂ + _∂  +  _∂  + _∂   _∂ − _{∂ ∂} +        (2.23) 其中υ_H為水平方向渦流黏滯係數= +υ υ_{l t}；υ_l為層流黏滯係數，約為 1.12 6 10− × ；υ_t 為紊流黏滯係數。式(2.22)與式(2.23)中，本研究以 3 u R 與R3v表示為水平紊流剪力 中經 sigma 轉換後的高次項，假設其數值計算結果可忽略，在模式計算時並未考 慮。 2.1.4 高含砂高含砂高含砂水流高含砂水流水流水流與紊流與紊流輔助與紊流與紊流輔助輔助輔助關係關係關係關係式式式式 密度變化 密度變化 密度變化 密度變化 高含砂水流造成在空間與時間上濃度的不均勻分布，而在含砂濃度主控的情 況下，可由狀態函數表示密度變化： 0 (1 C) C _s ρ= − ρ + ρ _(2.24) 藉由上式，可反應濃度對流體密度ρ 之影響，其中 C 為含砂水流濃度，ρ₀為 清水流密度，ρs為乾砂密度。本研究在水流方程式當中考慮密度變化，以及其對 壓力梯度在空間三個維度方向之影響，以反應密度效應對高含砂水流流動特性之 影響。 卡門係數 卡門係數 卡門係數 卡門係數 高含砂水流情況下之卡門係數並非常數，引用張紅武等(1994)：

(

)

0.4 1 4.2 C 0.365 C κ = _ − − _   (2.25) ，其中C為水深平均濃度。 底床剪應力 底床剪應力 底床剪應力 底床剪應力 高含砂水流情況下，高比例的含砂量造成水流流動特性改變，剪力與應力應 變關係可能由牛頓流體轉變為非牛頓流體特性。本研究以含砂濃度反應牛頓與非 牛頓兩種流體特性之主控參數，採用 O'Brien et al.(1993)提出之多項式流變關係計 算高含砂水流情況下之流變關係參數。流變關係主要作用於於二維水深平均動量 方程式之底床剪力，並且亦在式(2.15)與式(2.16)等號右邊第 4 項當中： 2 2 1 _{2 2} τ 8 l b f B B _{C u u v} u v K u u D τ = µ + ρ + + + (2.26) 2 2 2 _{2 2} τ 8 l b f B B _{C v u v} u v K v v D τ = µ + ρ + + + (2.27) τ_B=賓漢降伏應力；µ_B=賓漢黏滯係數；C_f = 2 / _c g c (c_c= Chezy 糙度係數)；K_l =多項式流變關係之層流阻力係數，為底床糙度函數(O'Brien et al.,1993)。 紊流黏滯係數 紊流黏滯係數 紊流黏滯係數 紊流黏滯係數 本研究採用零方程(zero-equation)紊流模式求解紊流黏滯係數，並選用文獻中 較常見之形式作為後續研究探討，如下分為兩種： Elder(1959) (type-I) ： * / 6 t U D υ κ= (2.28)

Jobson and Sayre(1970) (type-II) ：

* (1 )

t U D

上列υt即為υH = +υ υl t與υV = +υ υl t中之υt。兩式所依據的假設，是水深遠小於渠 寬，流速剖面依循對數分布。 層流與紊流剪應力 層流與紊流剪應力 層流與紊流剪應力 層流與紊流剪應力 式(2.2)與式(2.3)之層流剪應力與紊流擴散效應採用 Boussinesq 之渦流黏性理 論簡化為： 2 11 1 1 1 2 1 2 _H u v h u h h h τ _υ ρ ξ η  _{∂ ∂}  ′ − =  +  ∂ ∂   (2.30) 2 22 2 2 1 2 1 2 _H v u h v h h h τ _υ ρ η ξ  _{∂ ∂}  ′ − =  +  ∂ ∂   (2.31) 12 2 1 1 2 2 1 2 _H h v h u u v h h h h τ _υ ρ ξ η  _∂   _∂   ′ ′ − =   +   ∂   ∂     (2.32) 垂向層流剪應力以流速梯度表示時，若考慮正交曲線座標以及σ 座標時，可 表示如下： 13 1 1 1 V u w D h τ _υ ρ σ ξ  _{∂ ∂}  =  +  ∂ ∂   (2.33) 23 2 1 1 V v w D h τ _υ ρ σ η  _{∂ ∂}  =  +  ∂ ∂   (2.34) 式中，υ_V = +υ υ_{l t}為垂直方向渦流黏滯係數，在本研究中由於應用案例之格網 尺度差異均不大，因此假設υV 與υH可採用相同之數值，υt如式(2.29)、式(2.29) 所示。

2.1.5 邊界條件邊界條件邊界條件 邊界條件 水平 水平 水平 水平二維部分二維部分二維部分二維部分 水平二維方程式考量三種邊界條件設定，分別為渠道入流、渠道出流與固體 邊界。一般而言，渠道入流邊界條件設定為單位寬度入流量，渠道出流邊界條件 則採用水位高程設定。在固體邊界處，沿固體邊界法線方向採不透水邊界條件， 而沿固體邊界切線方向可分為滑移與非滑移邊界條件。另外，上游為亞臨界流時 需給予流量，若為超臨界流則另需提供水位資訊；下游若為亞臨界流需給予水 位，若為超臨界流則無須設定。 流速差異量部分 流速差異量部分 流速差異量部分 流速差異量部分 流速差異量方程式考量渠道入流、渠道出流、自由液面及底床邊界條件。在 渠道入流條件為流速，出流處採用主流方向零通量設定。自由液面採風剪力邊界 條件：

(

)

s1 V u u Dτ σ υ ρ ∂ + = ∂ ɶ (2.35)

(

)

s2 V v v Dτ σ υ ρ ∂ + = ∂ ɶ (2.36)

而底床則採用 French (1986)之壁函數(wall function)經驗式：

(

)

2 30 2.5 ln 2.72 l b b V s u u D z u u k σ υ −   ∂ +   =    ∂    ɶ (2.37)

(

)

2 30 2.5 ln 2.72 l b b V s v v _{D z} v v k σ υ −   ∂ +   =    ∂ _{  } ɶ (2.38) 式中，u_b、v_b分別為ξ、η方向之近底床流速；z_l為邊界層厚度，模式中採用近底

床格網之厚度；ks為粗糙高度。

2.2 沉滓運移

沉滓運移

沉滓運移部分

沉滓運移

部分

部分

部分

沉滓運移部分在懸浮載、底床載分開求解概念下，求解懸浮載之質量傳輸 (mass transport)方程式、作用層連續方程式(active-layer continuity equation)、底床 連續方程式(bed-layer continuity equation)(謝德勇，2002)，其中懸浮載之質量傳輸 為三維控制方程式。由於本研究模式各項變數之間的空間關係較為複雜，包含懸 浮載(suspended load)濃度、懸浮載底床濃度、懸浮載向上與向下通量、底床載(bed load)、作用層(active layer)等，將各變數於深度方向之關係表示於圖 2.4。本研究 將沉滓運移計算分為懸浮載與底床載兩部分，以底床高程z_b為分界。懸浮載運移 在zb之上，底床載運移在zb之下，且底床載運移發生於作用層內。作用層主要是 一理論上的假設物理量，模式假設底床載運移與粒徑變化均僅發生於作用層當 中，頂面位置為底床zb所在高程。懸浮載之底床濃度則是位於zb，交換通量也發 生於z_b。以下逐一介紹圖 2.4 中的各變數與計算方式。 2.2.1 基本控制方程式基本控制方程式基本控制方程式 基本控制方程式 沉滓運移控制方程式將沉滓粒徑分為數種代表粒徑，分開計算懸浮載與底床 載，計算各代表粒徑在懸浮載中的濃度以及在作用層當中之比例，並計算整體沉 滓在交換後所造成之底床沖淤量，將方程式說明如下 : 質量傳輸 質量傳輸 質量傳輸 質量傳輸方程式方程式方程式方程式 2 2 1 2 2 2 1 2 2 2 2 2 2 1 2 1 1 1 2 1 1 1 2 2 2 2 1 1 1 ( ) ( ) fh fh sV sV sH sH sH sH sH sH w w C C C C C C C C uh vh t h h D D D D C h h C h C h h h h h C h h C h C h h h ε ε _ω ξ η σ σ σ σ σ σ ε ε ε ξ ξ ξ ξ ξ ε ε ε η η η η η ∂   _∂ ∂ ₊ ∂ ₊ ∂ ₋ ∂ ₋ ∂ ₊ ∂ ₋ ∂ ₋   ∂  ∂ ∂  ∂ ∂ ∂ ∂ ∂ ∂  _{∂ ∂ ∂ ∂ ∂} =  + + ∂ ∂ ∂ ∂ ∂   ∂ ∂ ∂ ∂ ∂ + + +  ∂ ∂ ∂ ∂ ∂  (2.39)

上式移流項部分藉由 sigma 座標轉換後，產生ω速度項，係以式(2.17)求解。 水平擴散項由於 sigma 座標轉換後產生多項二次微分之乘積，本研究中假設其值 可予以忽略。 作用層連續方程式 作用層連續方程式 作用層連續方程式 作用層連續方程式 2 1 1 2 1 2 1 2 ( ) (1 ) ( ) ( ) 0 m m m m p s b b s E E p h q h q S S t h h h h ∂ β ∂ ∂ ρ ∂ ∂ξ ∂η − + + + − = (2.40) 底床連續方程式 底床連續方程式 底床連續方程式 底床連續方程式 2 1 1 2 1 1 2 1 2 (1 ) ( ) ( ) 0 m m m M b s b b s m z p h q h q S t h h h h ∂ ∂ ∂ ρ ∂ = ∂ξ ∂η   − +  + + =  

∑

(2.41) 上列諸式中，C = 某一代表粒徑的懸浮載體積濃度；εsV、εsH=垂直、水平方 向紊流擴散係數;w_fh= 泥砂落淤速度；ρ = 含砂水流密度；ρ_s = 乾砂密度；β_m= 第 m 組代表粒徑在作用層中所佔百分比；p = 孔隙率；E_p= 作用層厚度，模式中 假設為 0.2 - 0.5 (m)； SE = 作用層源(source of active layer)；zb = 底床高程；qb1m、

2m b q =水流方向、水流側方向方向第 m 組代表粒徑之底床載通量； m s S 為第 m 組代 表粒徑懸浮載源(source of suspended load)。

2.2.2 高含砂水流高含砂水流高含砂水流與紊流高含砂水流與紊流與紊流輔助與紊流輔助輔助輔助關係式關係式關係式關係式 落淤速度 落淤速度 落淤速度 落淤速度 落淤速度(falling velocity)是沉滓在底床交換、流體中垂直方向移流之重要物 理機制，在分析高含砂水流懸浮載運移機制上，是本研究不可或缺之參數，常以 經驗公式推估之。一般而言，清水流中的沉滓落淤速度主要變數為粒徑 Dm，van Rijn(1984)蒐集整理不同粒徑區間(Dm)之落淤速度(wf)計算式如下：