以數值模擬探討高含砂效應對渠流流況之影響

國立交通大學

土木工程學系

碩士論文

以數值模擬探討高含砂效應對渠流流況之影響

Numerical Investigation of Hyperconcentrated Flow

Behavior in Open Channel

研 究 生：羅冠顯

指導教授：楊錦釧 博士

謝德勇 博士

以數值模擬探討高含砂效應對渠流流況之影響

Numerical Investigation of Hyperconcentrated Flow

Behavior in Open Channel

研 究 生：羅冠顯 Student: Guan-Sian Luo

指導教授：楊錦釧 Advisor: Jinn-Chuang Yang

謝德勇

Te-Yung Hsieh

國 立 交 通 大 學

土 木 工 程 研 究 所

碩 士 論 文

A Thesis Submitted to Civil Engineering

College of Engineering

Nation Chiao Tung University

in Partial Fulfillment of the Requirements

for the Degree of Master

in

Civil Engineering

July 2008

誌 謝

承蒙恩師 楊教授錦釧與謝博士德勇於論文研究期間之悉心指 導，使得本論文得以順利完成，在此謹致由衷的敬意與感謝。在論文 審定期間，感謝口試委員詹教授錢登與陳教授樹群之細心匡正與寶貴 建議，使本論文更臻完善；此外，亦感謝系上老師於課業上的啟發， 讓學生在兩年交大生活獲益匪淺。 感謝就學期間研究室祥禎學長、夢祺學長、胤隆學長、昇學學長、 世偉學長、浩榮學長、建華學長、弘恩學長、宥達學長、欣瑜學姊、 歆婷學姊於課業與生活上的提攜與照顧，亦感謝同學俊哲、鏡如、仙 蕓、誠達、仁凱、思廷、佑民的相互砥礪與扶持，及學弟妹全謚、振 家、俊宏、歆淳、冠曄、宇翔、柏傑平時的幫忙協助。 同時，感謝摯友仁傑、杰昇、上源、國晟、峻榮、昆廷、家賓、 鳳玉、展毅、作欣、承龍、晉星、濡岳、其均、國棟、兆成於此段求 學時間傾聽於壓力上之排解與精神上之鼓勵，在我苦悶時陪我ㄧ起歡 度快樂時光，這堅定的友情是一輩子都不會變的。 最後僅以本論文獻給我親愛的父母、妹妹芊琇、弟弟冠麟，感謝 你們在精神上的支持與鼓勵，而女友奐文在背後默默的支持更是我前 進的動力，沒有你們的體諒、包容，相信這兩年的生活將是很不一樣 的光景，謹以學位論文與你們分享我成長的喜悅，感恩之心無以言喻。

以數值模擬探討高含砂效應對渠流流況之影響

學生：羅冠顯 指導教授：楊錦釧 謝德勇 國立交通大學土木工程學系

中文摘要

楊與謝(2003)發展一水深平均二維水理模式，鍾(2007)進一步考 量非牛頓流體之流變關係，將模式擴充為高含砂水流模式。本研究以 上述兩模式為基礎，藉此探討高含砂水流對明渠水流水深、流速與剪 力之效應，並分析其影響程度，其結果可提供與模式使用者作為選用 合適模式之參考。

關鍵字：水深平均二維水理模式、非牛頓流體、流變關係、高含砂水

流

Numerical Investigation of Hyperconcentrated Flow

Behavior in Open Channel

Student：Guan-Sian Luo Advisors：Jinn-Chuang Yang Te-Yung Hsieh Department of Civil Engineering

National Chiao-Tung University

ABSTRACT

A numerical experiment is carried out to study the hyperconcentrated flow behavior in open channel. Two kinds of 2D depth-averaged model are considered in this study of which the hyperconcentrated flow model (Chung 2007) includes the Newtonian and non-Newtonian fluid rheological relations in bottom shear stress, and the clear flow model (Yang and Hsieh 2003) only includes the effect of Newtonian fluid in bottom shear stress.

The maximum relative discrepancy of the depth, velocity, and shear stress obtained from the comparison of these two models, are use as the criterion to judge the hyperconcentrated effect on open channel flow. The results can provide a guideline for model users to determine the proper approach to simulate the hyperconcentrated flow problem by either using the clear flow model or the hyperconcentrated flow model.

Key words: depth-averaged 2-D flow model

、non

-Newtonian fluid

、

rheological relation

、

hyperconcentrated flow

目錄

誌謝... I 中文摘要...II 英文摘要... III 目錄... IV 表目錄...VII 圖目錄...VIII 符號表...XII 第一章 緒論... 1 1.1 研究動機與方向 ... 1 1.2 文獻回顧 ... 2 1.3 研究目的與方法 ... 2 1.4 章節介紹 ... 3 第二章 理論基礎... 6 2.1 控制方程式 ... 6 2.2 輔助關係式 ... 7 2.3 邊界條件 ...11 第三章 數值架構... 15 3.1 雙階分割操作趨近法 ... 15

3.2 數值差分式 ... 17 第四章 高含砂水流模式功能測試 ... 20 4.1 穩態流案例測試 ... 20 4.1.1 亞臨界流案例測試 ... 20 4.1.2 超臨界流案例測試 ... 22 4.2 湧波傳遞案例測試 ... 23 第五章 高含砂效應影響程度之模擬分析 ... 31 5.1 因次分析 ... 31 5.2 模擬案例設定 ... 33 5.3 參數重要性分析 ... 34

5.4 高含砂效應對_MaxH*_、MaxU*_及Maxτ*_{影響分析... 35}

5.4.1 福祿數效應影響分析 ... 35 5.4.2 直線道案例影響分析 ... 37 5.4.3 彎道案例影響分析 ... 38 5.5 高含砂效應重要性分析 ... 39 5.5.1 直線道案例 ... 39 5.5.1.1 _MaxH*_{影響分析... 39} 5.5.1.2 _MaxU*_{影響分析... 40} 5.5.1.3 _Maxτ*_{影響分析... 41} 5.5.2 彎道案例 ... 43

5.5.2.1 _MaxH*_{影響分析... 43} 5.5.2.2 _MaxU*_{影響分析... 44} 5.5.2.3 _Maxτ*_{影響分析... 45} 第六章 結論與建議... 89 6.1 結論 ... 89 6.2 建議... 91 參考文獻... 92 附錄 FLO-2D 數值模式介紹... 95

表目錄

表1.1 FLO-2D 與 RESED2D 模式比較表... 5 表2.1 開放邊界處理原則... 13 表4.1 層流阻力係數參考表... 25 表4.2 亞臨界流測試案例邊界條件 ... 25 表4.3 超臨界流測試案例邊界條件 ... 25 表5.1 輸入變數統計特性表... 47 表5.2 測試案例一覽表... 48 表5.3 流變參數一覽表... 49 表5.4 彎道_MaxH*_不同 b θ 之E_rms與ρ_c表... 50 表5.5 彎道_MaxU*_不同 b θ 之E_rms與ρ_c表... 50 表5.6 彎道_Maxτ*_不同 b θ 之E_rms與ρ_c表... 50

圖目錄

圖2.1 彎道二次流示意圖... 14 圖3.1 控制體積法示意圖(a)實際區域；(b)計算區域... 19 圖4.1 亞臨界流平均水深與平均流速展示圖... 26 圖4.2 亞臨界流底床剪應力項與體積濃度關係圖 ... 27 圖4.3 超臨界流平均水深與平均流速展示圖... 28 圖4.4 超臨界流底床剪應力項與體積濃度關係圖 ... 28 圖4.5 高含砂水流測試案例邊界條件 ... 29 圖4.6 高含砂水流模式測試比較-體積濃度 25% ... 29 圖4.7 高含砂水流模式測試比較-體積濃度 30% ... 30 圖4.8 高含砂水流模式測試比較-體積濃度 40% ... 30 圖5.1 設計案例渠道示意圖... 51 圖5.2 直線道_MaxH*_{與無因次參數關係圖... 52} 圖5.3 直線道_MaxU*_{與無因次參數關係圖... 53} 圖5.4 直線道_Maxτ*_{與無因次參數關係圖... 54} 圖5.5 彎道_MaxH*_{與無因次參數關係圖... 56} 圖5.6 彎道_MaxU*_{與無因次參數關係圖... 58} 圖5.7 彎道_Maxτ*_{與無因次參數關係圖... 60}

圖5.8 福祿數效應對直線道_MaxH*_{影響分析圖... 61} 圖5.9 福祿數效應對直線道_MaxU*_{影響分析圖... 61} 圖5.10 福祿數效應對直線道_Maxτ*_{影響分析圖... 62} 圖5.11 福祿數效應對彎道_MaxH*_{影響分析圖... 64} 圖5.12 福祿數效應對彎道_MaxU*_{影響分析圖... 66} 圖5.13 福祿數效應對彎道_Maxτ*_{影響分析圖... 68} 圖5.14 直線道 2 8 B B w K U U H μ τ ρ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ 對 * MaxH 、_MaxU*_與Maxτ*_{關係圖... 70}

圖5.15 不同θ_b之 2 8 B B w K U U H μ τ ρ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ 與 * MaxH 關係圖(1)... 72 圖5.16 不同θ_b之 2 8 B B w K U U H μ τ ρ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ 與 * MaxU 關係圖(1)... 74 圖5.17 不同θ_b之 2 8 B B w K U U H μ τ + ρ ⎝ ⎠ * ⎛ ⎞ ⎜ ⎟ 與Maxτ 關係圖(1)... 76 圖5.18 不同θ_b之 2 8 B B w K U U H μ τ ρ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ 與 * MaxH 關係圖(2)... 77 圖5.19 不同θ_b之 2 8 B B w K U U H μ τ ρ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ 與 * MaxU 關係圖(2)... 77 圖5.20 不同θ_b之 2 8 B B w K U U H μ τ + ρ ⎝ ⎠ * ⎛ ⎞ ⎜ ⎟ 與Maxτ 關係圖(2)... 78 圖5.21 直線道高含砂效應對_MaxH*_{影響分類圖... 78} 圖5.22 直線道 * 8 B B K U 2 w MaxH U H μ τ ⎛ ⎞ −_⎜ + _⎟ ⎝ ⎠ ρ 關係曲線... 79 圖5.23 直線道高含砂效應對_MaxU*_{影響分類圖... 79}

圖5.24 直線道 * 8 B B K U 2 w MaxU U H μ τ ⎛ ⎞ −_⎜ + _⎟ ⎝ ⎠ ρ 關係曲線... 80 圖5.25 直線道 * m ρ 與_Maxτ*_{關係圖... 80} 圖5.26 直線道 2 8 B B w K U U H μ τ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ ρ 與 * m ρ 對_Maxτ*_{相對誤差分析圖... 81} 圖5.27 直線道 2 8 B B w K U U H μ τ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ ρ 與 * m ρ 對_Maxτ*_{影響分類圖... 81} 圖5.28 直線道 * 8 B B K U 2 w Max U H μ τ −⎛_⎜τ + ⎞_⎟ ρ ⎝ ⎠ 關係曲線... 82 圖5.29 直線道 * m Maxτ −ρ *關係曲線... 82 圖5.30 彎道 2 8 B B w K U U H μ τ ρ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ 與 * MaxH 關係圖... 83 圖5.31 彎道高含砂效應對_MaxH*_{影響分類圖... 83} 圖5.32 彎道 * 8 B B K U 2 w MaxH U H μ τ ⎛ ⎞ −_⎜ + _⎟ ⎝ ⎠ ρ 關係曲線... 84 圖5.33 彎道 2 8 B B w K U U H μ τ ρ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ 與 * MaxU 關係圖... 84 圖5.34 彎道高含砂效應對_MaxU*_{影響分類圖... 85} 圖5.35 彎道 * 8 B B K U 2 w MaxU U H μ τ ⎛ ⎞ −_⎜ + _⎟ ⎝ ⎠ ρ 關係曲線... 85 圖5.36 彎道 2 8 B B w K U U H μ τ + ρ ⎝ ⎠ * ⎛ ⎞ ⎜ ⎟ 與Maxτ 關係圖... 86 圖5.37 彎道 * m ρ 與_Maxτ*_{關係圖... 86}

圖5.39 彎道 * 8 B B K U 2 w Max U H μ τ −⎛_⎜τ + ⎞_⎟ ρ ⎝ ⎠ 關係曲線... 87 圖5.40 彎道 * m Maxτ −ρ *關係曲線... 88

符號表

B = 渠道寬度； v C = 體積濃度； Cr = Courant number； f C = 摩擦因子； c = Chezy 係數； d = 水深； rms E = 均方根差； Fr = 福祿數； s G = 土砂比重； g = 重力加速度； H _{= 平均水深；} 1 h 、h₂ = ξ、η 方向轉換係數； K = 層流阻力係數； k = von Karman’s 係數； L = 渠道長度； * MaxH = 水深之最大相對差異； * MaxU _{= 流速之最大相對差異；} * Maxτ _{= 剪應力之最大相對差異；} total N = 比較值的數量； n = 曼寧糙度係數； q_ξ = ξ方向單寬流量； q_η = η方向單寬流量；

R = 水力半徑； c r = 渠道中心線曲率半徑； SI = 二次流強度因子； 0 S = 渠道坡度； 11 T 、T₁₂、T₂₂ = 有效剪應力項； t = 時間； U = ξ 方向平均速度； u = ξ 方向速度； * u = 剪力速度； V = η 方向平均速度； v = η 方向速度； X = 表平均數； i x = 某比較點回歸線值； Y = 表平均數； i y = 相對應比較點回歸線值； b z = 底床高程； 0 b z = 原始底床高程； s z = 水面高程； H B = 寬深比； c ρ _{= 相關係數；} B μ = 賓漢黏滯係數； B τ = 賓漢降伏應力； b θ = 彎道長度因子； * m ρ = 無因次含砂水流密度； m ρ = 含砂水流密度；

s ρ = 泥砂密度； 1 b τ 、τb₂ = 底床剪應力在 ξ 與 η 方向之分量； model τ = 模式模擬所得之底床剪應力； regression τ = 迴歸函數計算所得之底床剪應力； 2 8 B B K U U H μ τ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ ρw = 無因次非牛頓流變剪應力項； 1 α = 賓漢黏滯係數係數項； 1 β = 賓漢黏滯係數指數項； 2 α = 賓漢降伏應力係數項； 2 β = 賓漢降伏應力指數項； ξ、η = 平面上兩正交曲線座標方向； t Δ = 時間間距。 上標 n = n tΔ _{時刻之已知變數；} 1 n+ = (n+ Δt1) _{時刻之未知變數；} 1 2 n+ = (n+ Δ1) t與n tΔ _{間之未知變數；} ( ) = 時間平均； ( ) = 水深平均； ( ) = 時間平均瞬時擾動量。 ' 下標 s = 變數在水面的值； b = 變數在底床的值； m = 高含砂水流變數； w = 清水流變數；

第一章 緒論

1.1 研究動機與方向 台灣本島由於地形上較為特殊，因此河流特徵多為河身短、坡度 陡峻、水流湍急，也因此使河川蘊藏較高的沖刷潛勢。再加上台灣地 區地質脆弱、地震頻繁且雨量豐沛，這些因素組合起來，使得台灣地 區由於自然因素所造成之土砂災害很多。且台灣地區三分之二為山 區，加上人口密度不斷增高，因此只好轉向山區發展，由於山區過度 開發的結果，使得山坡地極為不穩定。先天自然因素不佳加上人為因 素破壞，使台灣山坡地災害頻傳，甚至引發高含砂水流災害造成人民 生命財產之威脅。 因此，學者們便常藉由數學模式來模擬分析高含砂水體的流動行 為及流動特性，進而評估可能造成之災害損失。但高含砂水流模式相 較於清水流模式除所需模式計算時間較長外，數值穩定性亦較差，因 此若能釐清高含砂效應對於兩模式的影響程度，提出兩模式選用之標 準，則可提供與模式使用者作為選用合適模式之參考依據。 綜觀上述，制定本研究之研究方向為釐清高含砂效應對渠流流況 之影響程度，提出在模擬分析高含砂水流問題時，何時需考量高含砂 效應的影響，何時可忽略高含砂效應的影響，希冀能提供與模式使用 者作為選用合適模式之參考依據。

1.2 文獻回顧 高含砂水流為水與泥砂顆粒之混合體，且流變特性隨著泥砂含量 與泥砂粒徑的改變而有很大的變化；如Xu(2004)與 Zhao(1996)表示， 當含砂水流懸浮載濃度超過某一極限值，此時含砂水流之物理特性與 力學性質便可能產生變化，使含砂水流脫離牛頓流體範疇進而轉變為 賓漢流體；張德茹等(2000)在研究大陸洛惠渠之高含砂水流現象時指 出，當含砂水流之含砂濃度達到一定程度時，具有非牛頓流體之特性。 隨著電腦科技蓬勃發展，專家學者常藉由數值模式來模擬分析高

含砂水流流動行為，如Wei (1990)、Zhang et al.(2001)、Ni et al.(2004)、

Cao (2006)針對大陸黃河所發展的高含砂水流模式；黃與萬(2001)針

對中國大陸北方之多砂河川，假設高含砂水流為賓漢流體(Bingham

fluid)所研發之高含砂水流模式；Liu and Huang (2006)參考Julien and

Lan (1991)所發展的流變模式將高含砂水流視為賓漢流體進行模擬。 1.3 研究目的與方法 本研究旨在釐清高含砂效應對渠流流況之影響，藉由比較高含砂 水流與清水流模式之模擬差異性，探討高含砂效應對渠流流況之影響 程度，提供與模式使用者作為研選合適模式之參考。表1.1 為 FLO-2D 模式與RESED2D 模式的比較表，分別比較兩模式後可發現，FLO-2D

含砂水流流經彎道時，水深與流速所產生的變化情形，造成模擬分析 上之誤差。因此，為有效模擬分析高含砂水流之流動行為，本研究在 此便藉由鍾(2007)延續楊與謝(2003)發展之水深平均二維模式，結合 高含砂水流流變關係，所發展的高含砂水流定床模式RESED2D 來進 行模擬分析。 為探討高含砂效應對渠流流況之影響程度，本研究首先藉由因次 分析決定相關無因次參數，再經由參數重要性分析決定影響渠流流況 之重要參數，最後並建立高含砂效應與重要參數之關係圖，藉此釐清 高含砂效應對渠流流況之影響程度，並可提供與模式使用者作為研選 合適模式之參考。模式驗證方面，清水流模式部份謝(2003)已分別針 對迴水演算、彎道、環流流場及潰壩等案例進行模式功能測試，且均 得到良好的結果，故於本研究中便不多作贅述；高含砂水流模式部 份，則採用流變關係與 RESED2D 相同之商用 FLO-2D 模式來驗證模 式的合理性。 1.4 章節介紹 前三節已闡述本研究之研究動機和方向、文獻回顧、研究目的與 方法，以下將扼要說明本論文各章節之內容。 第一章緒論，針對本研究之緣起和方向作說明，回顧相關文獻 後，再提出本研究之目的，並於章末做論文架構說明。

第二章理論基礎，分別闡述本研究之二維正交曲線座標系統水理 控制方程式、輔助方程式及相關之邊界條件。 第三章數值架構，說明水理控制方程式所採用數值方法和差分型 式。 第四章模式功能測試，藉由與FLO-2D 商業模式之比較，來測試 模式模擬高含砂水流之功能。 第五章高含砂效應影響程度之模擬分析，藉由分析高含砂水流模 式與清水流模式水深、流速與剪力模擬差異性，來探討高含砂效應對 渠流流況之影響程度，提供與模式使用者作為選用合適模式之參考依 據。 第六章結論與建議，除對本研究成果作綜合性的論述外，並針對 內容不盡完善與日後可改善之處提出建議。

表 1.1 FLO-2D 與 RESED2D 模式比較表 FLO-2D模式 RESED2D模式 空間維度 二維 ● ● 亞臨界流 ● ● 超臨界流 ● ● 定量流 ● ● 變量流 ● ● 回流流場 ● 支流入流 ● 水理條件 彎道(二次流效應) ● 河床質載 ● ● 懸浮載與河床載分離 ● 河床質篩分甲護 ● ● 凝聚性沉滓 ● 支流入流 ● ● 彎道 ● 輸砂公式選擇 ● 輸砂條件 高含砂水流動床模式耦合計算 ● 流變關係 由輸砂模式計算體積濃度 ● 數值方法 有限差分 ● ● 其他 障礙物或水工建造物 ● ●

第二章 理論基礎

透過座標系統轉換將控制方程式轉換為正交曲線座標系統方程 式，再將此控制方程式作時間平均及水深平均後，即可推得水深平均 二維正交曲線座標模式所需之控制方程式。茲將模式採用的理論基礎 敘述如下： 2.1 控制方程式 為適度簡化複雜的控制方程式，需對數學模式作若干假設，分別 為(1)不可壓縮牛頓流體(incompressible Newtonian fluid)；(2)靜水壓分 布；(3)忽略風剪力；(4)忽略科氏力。則水深平均二維正交曲線座標 水理控制方程式可表示為 (1)水流連續方程式 0 1 2 2 1 d h h + (h ud)+ (h vd) t ξ η ∂ ∂ ∂ ₌ ∂ ∂ ∂ _(2.1) (2)水流動量方程式

ξ

方向： 2 1 2 1 2 1 2 1 2 1 h 1 h u u u v u uv v t h h h h h h ∂ ∂ ∂ ∂ ∂ ∂ + ∂ξ + ∂η + ∂η − ∂ξ 2 11 1 12 1 1 2 1 2 1 1 ( _b ) ( ) ( m m g z d h T h T h h h d h h d ) ∂ ∂ ∂ ∂ξ ρ ∂ξ ρ ∂η = − + + + 1 1 2 12 22 1 2 1 2 1 1 b m m h h T T h h d h h d _md τ ∂ ∂ ρ ∂η ρ ∂ξ ρ + − − 2 11 2 11 1 12 1 12 1 2 1 ( ) s ( ) b ( ) s ( ) s b s b m z z z h h h h h h d τ τ τ τ b z ρ ξ ξ η ⎡ ∂ ∂ ∂ ⎤ + _⎢− + − + _⎥ ∂ ∂ ∂ ⎣ η⎦ ∂ ∂ (2.2)

η

方向：

2 2 1 1 2 1 2 1 2 1 h 1 h v u v v v uv u t h h h h h h ∂ ∂ ∂ ∂ ∂ ∂ + ∂ξ + ∂η + ∂ξ − ∂η 2 12 1 22 2 1 2 1 2 1 1 ( _b ) ( ) ( m m g z d h T h T h h h d h h d ) ∂ ∂ ∂ ∂η ρ ∂ξ ρ ∂η = − + + + 2 1 2 11 12 1 2 1 2 1 1 b m m h h T T h h d h h d _md τ ∂ ∂ ρ ∂η ρ ∂ξ ρ − + − 2 12 2 12 1 22 1 22 1 2 1 ( ) s ( ) b ( ) s ( ) s b s b m z z z h h h h h h d τ τ τ τ b z ρ ξ ξ η ⎡ ∂ ∂ ∂ ⎤ + _⎢− + − + _⎥ ∂ ∂ ∂ ⎣ η⎦ ∂ ∂ (2.3) 式中 2 2 11 [ 11 ' ( ) s b z m m z T =

∫

τ −ρ u −ρ u−u ]dz (2.4) 2 22 [ 22 ' ( ) s b z m m z T =

_∫

τ −ρ v −ρ v−_v 2_{]dz (2.5)} 12 21 [ 12 ' ' ( )( )] s b z m m z T =T =

∫

τ −ρ u v −ρ u u v v d− − z _(2.6) 以上諸式中， ξ、_η = 平面上兩正交曲線座標方向，其中ξ為縱 方向，_η為側方向；h₁ = ξ 方向轉換係數； = h₂ η 方向轉換係數； = u

ξ

方向速度；v = η 方向速度；ρ_m = 含砂水流密度；d = 水 深；g = 重力加速度； = 時間； = 底床高程；t z_b z = 水面高程；_s i b

τ

=底床剪應力在

_ξ

與

_η

方向之分量; ( )= 時間平均; ( ) =水深平 均；( ) = 時間平均瞬時擾動量；下標 、b 分別代表變數在水面與

底床的值； , _{,T = 有效剪應力項(effective stress term)，包含層} 流剪應力、亂流剪應力與延散剪應力(dispersion stresses)。 ' s 11 T T_{12 22} 2.2 輔助關係式 (1)含砂水流密度

本模式採用狀態函數(State Function)來反映體積濃度 對含砂水 流密度 v C m ρ 之影響，關係式如下所示： ρm =ρsCv+ρw(1−Cv) (2.7) 式中，_ρs=乾砂密度；_ρw=清水密度；C_v=沉滓體積濃度。 (2)底床剪應力 在模擬高含砂水流時，所需考慮之應力應包含凝聚降伏應力 c

τ (Cohesive yield stress)、莫爾庫倫剪應力τ_mc(Mohr-Coulomb shear stress)、粘滯剪應力τ_v(Viscous shear stress)、紊流剪應力τ_t(Turbulent shear stress)和離散剪應力τ_d(Dispersive shear stress)等五項應力，關係 式可表示如下： τ τ τ= +_{c mc}+ + +τ τ τ_{v t d} (2.8) O´Brien 和 Julien(1985)將式(2.8)改寫為 _{τ ( ) ( )}2 B B du du C dy dy τ = +μ + (2.9) 其中 _{τ τ} B c mc τ = + (2.10) 1 3 2 * ₁ 2 m i m v C C l a d C ρ ρ ⎡⎢⎛ ⎞ ⎤⎥ = + _{⎜ ⎟} − ⎢_{⎝ ⎠} ⎥ ⎣ ⎦ s (2.11) 式中，_τ =賓漢降伏應力； B μB=賓漢粘滯係數；l=Prandtl混合長度； =

經 驗 係 數(=0.01) ； = 最 大 靜 體 積 濃 度(Maximum static volume concentration)； i a * C s d =沉滓粒徑。

4 3 2 2 2 τ γ 8γ B B y v f td md md d K u n u S =S +S + S = + μ + (2.12) τ 2 2 (2.13) Cv B =α eβ 1 1 (2.14) Cv B e β μ =α 其中，_Sy=降伏坡降；_Sv=黏滯坡降；_Std=紊流坡降；_γm=水砂混合比 重； =層流阻力係數；_{K n}=曼寧糙度係數；_αi、_βi=經驗參數，由流變 試驗獲得。 應用(2.12)式修正(2.2)及(2.3)兩式中之底床剪力項，可得 1 2 2 2 2 τ 8 b B B Cf u u v u v u K u d τ = μ + ρ_m + + + (2.15) 2 2 2 2 2 τ 8 b B B Cf v u v u v v K v d τ = μ + ρ_m + + + (2.16) 式中，_Cf =摩擦係數。 (3)層流與亂流剪應力 採用Boussinesq之渦流黏性理論，層流與亂流剪應力可合併表示 為 2 11 1 1 1 2 1 ' 2 m h u v u h h h τ _υ ρ ξ ⎡ ⎤ η ∂ ∂ − = _{⎢ ∂} + ∂ ⎣ ⎦⎥ (2.17) 2 22 2 2 1 2 1 ' 2 m h v u v h h h τ _υ ρ η ⎡ ∂ ∂ ⎤ − = _⎢ + ∂ ∂ ⎣ ξ ⎦⎥ (2.18) 12 2 1 1 2 2 1 ' ' 2 m h v h u v h h h h τ _υ ρ ξ ⎡ ∂ ⎛ ⎞ ∂ ⎛ ⎞⎤ − = _{⎢ ⎜ ⎟}+ _{⎜ ⎟} ∂ _{⎝ ⎠} ∂ _{⎝ ⎠} ⎣ ⎦ u η ⎥ (2.19)

式中，_{υ υ υ= + ；l t υl} = 層流黏滯係數；_υt = 亂流黏滯係數 = (Falcon 1979)； * / 6 ku d 1/ 2 * ( / )b u = τ ρ = 剪力速度； = von Karman’s 係數(約等於0.4)。 k (4)延散剪應力 為積分水深平均所產生之延散剪應力項，須對流速剖面作一適當 假設，本模式在延散剪應力的處理方面，則僅考量二次流的影響。 由於水流在進入彎道後，流場隨水流而彎曲，致使流線因彎曲而 產生徑向慣性力，水面因而形成超高以產生徑向靜水壓差，得以與徑 向慣性力取得平衡。在這兩種力之作用下，水流除了以縱向方向流動 外，在徑向尚產生兩層水流，上層水流之外岸慣性力大於靜水壓差， 下層水流則反之，因此造成上層水流流動方向朝向外岸，下層水流則 為朝內岸流動，稱之二次流，如圖2.1所示。

本模式中二次流速度剖面參照Hsieh and Yang(2003)所採用 de Vriend (1977)之假設： 1 g g ln _m u u uf kc kc ζ ( )ζ ⎡ ⎤ = ⎢ + + ⎥= ⎢ ⎥ ⎣ ⎦ (2.20) 1 2 2 ( ) 2 ( ) ( ) 2(1 ) ( ) m g g ud v vf F F f k r kc kc m ζ ⎡ ζ ζ ζ ⎤ = + ⎢ + − − ⎥ ⎢⎣ ⎥⎦ (2.21) 式中， 1 ₀1 ln ( ) 1 F ζ ζ dζ ζ = −

∫

； 2 1 2 ₀ ln ( ) 1 F ζ ζdζ ζ = −

∫

；ζ=(z-z )/db = 距離底 床之高度與水深之比值；r = 曲率半徑。 de Vriend (1977)二次流速度剖面之適用範圍為(1)水深遠小於渠 道寬度；(2)渠道寬度遠小於渠道之曲率半徑；(3)單一二次流(single

2.3 邊界條件 邊界條件為數值模式中相當重要的一環，所設定的個數必須符合 物理意義。邊界條件一般可分為開放式邊界條件與固體邊界條件，茲 分述如下： 開放式邊界條件主要設定在渠道上游與下游邊界處，依上下游不 同的流況可分為 一、 上游 (1) 超臨界流：此時上游需給定三個邊界條件，即

ξ

、

η

方向 單寬流量q_ξ、q_η與水深d 。 (2) 亞臨界流：此時上游僅需給定兩個邊界條件，一般給定

ξ

、

η

方向的單寬流量q_ξ、q_η。 二、 下游 (1) 超臨界流：下游不需任何邊界條件。 (2) 亞臨界流：需設定一個下游邊界條件，通常為下游水深 d 。 詳細的開放邊界條件處理原則如表2.1 所示。 固體邊界在沿法線方向(

η

方向)為非貫穿條件，所以該處的流速 為零(vwall = )；沿固體邊界切線方向(0

ξ

方向)而言，可分為滑移條件 與非滑移條件，處理規則如下： (1) 固體邊界為非滑移條件時，則固體邊界處

ξ

方向的流速為 零(uwall = ) 0

(2) 固體邊界為滑移條件時，則固體邊界處

ξ

方向的流速等於 相鄰格網點的流速(uwall =uwall−1)。 其中vwall為固體邊界處

η

方向水深平均流速，uwall為固體邊界處ξ 方向水深平均流速，uwall−1為相鄰固體邊界格點ξ方向水深平均流速。 此外，在模式體積濃度設定方面，本研究目前所模擬的案例，並不考 慮濃度延層變化之機制，因此在模擬過程中，均將濃度視為一恆定值 來進行模擬。

表2.1 開放邊界處理原則

位置 上游 下游

流況 超臨界流 亞臨界流 超臨界流 亞臨界流

邊界條件 q_ξ、q_η、d q_ξ、q_{η 無 d}

內岸 外岸

第三章 數值架構

3.1 雙階分割操作趨近法 本研究基於分割操作之觀念，將動量方程式分割成二個步驟(延 散步驟及傳播步驟)，並利用隱式數值方法求解。延散步驟求解移流 項和擴散項，傳播步驟求解壓力項、底床剪應力項和連續方程式。據 此，水理控制方程可改寫成： 延散步驟 1 1 1 2 2 1 2 ( ) n n n n m V V V T t ρ + + + ∂ ⎛ ⎞ _{= − ⋅∇ + ∇ ⋅} ⎜ _∂ ⎟ ⎝ ⎠ _(3.1) 傳播步驟 1 1 2 ₁ ( ) n n n b b m V V g z d t t d τ ρ + + + ∂ ∂ ⎛ ⎞ ₋⎛ ⎞ _{= − ∇ + −} ⎜ _∂ ⎟ ⎜ _∂ ⎟ ⎝ ⎠ ⎝ ⎠ _(3.2) 1 ₀ n V + ∇ ⋅ = _(3.3) 式中，V 表示速度向量；T 表示擴散及延散項；n+1表示 時 刻之未知變數； ； 表示 (n+ Δ1) t 1 n n t t + t Δ = − n n tΔ 時刻之已知變數； 2 1 + n 表示在 與 間之未知變數。 (n+ Δ1) t n tΔ (3.1) ~ (3.3)的一般式可表示成： 延散步驟 2 1 2 1 2 1 2 1 h h u u u v u uv v t h ξ h η h h η ξ ⎡∂ ∂ ⎤ ∂ _{= −} ∂ ₋ ∂ _{− −} ⎢ ⎥ ∂ ∂ ∂ _⎣∂ ∂ _⎦ 2 11 1 12 1 2 12 22 1 2 1 2 1 2 1 2 ( ) ( ) 1 1 1 1 m m m m h T h T h h T T h h d h h d h h d h h d ρ ξ ρ η ρ η ρ ∂ ∂ ∂ + + + − ∂ ∂ ∂ ξ ∂ ∂

2 11 2 11 1 12 1 12 1 2 1 ( ) s ( ) b ( ) s ( ) s b s b m z z z h h h h h h d τ τ τ τ b z ρ ξ ξ η ⎡ ∂ ∂ ∂ ⎤ + _⎢− + − + _⎥ ∂ ∂ ∂ ⎣ η⎦ ∂ ∂ (3.4) 2 2 1 1 2 1 2 1 h h v u v v v uv u t h ξ h η h h ξ η ⎡∂ ∂ ⎤ ∂ _{= −} ∂ ₋ ∂ _{− +} ⎢ ⎥ ∂ ∂ ∂ _⎣∂ ∂ _⎦ 2 12 1 22 1 2 11 12 1 2 1 2 1 2 1 2 ( ) ( ) 1 1 1 1 m m m m h T h T h h T T h h d h h d h h d h h d ρ ξ ρ η ρ η ρ ∂ ∂ ∂ + + − + ∂ ∂ ∂ ξ ∂ ∂ 2 12 2 12 1 22 1 22 1 2 1 ( ) s ( ) b ( ) s ( ) s b s b m z z z h h h h h h d τ τ τ τ b z ρ ξ ξ η ⎡ ∂ ∂ ∂ ⎤ + _⎢− + − + _⎥ ∂ ∂ ∂ ⎣ η⎦ ∂ ∂ (3.5) 傳播步驟 2 2 2 2 2 1 ( ) τ 8 f b m m B B C u u v z d u g t h _{d u v} d u K u d ξ _ρ ρ μ ⎛ ₊ ⎞ ⎛∂ + ⎞ ∂ _{= − −⎜ + ⎟} ⎜ _{⎟ ⎜ ⎟} ∂ _⎝ ∂ _{⎠ ⎝ +} + _⎠ (3.6) 2 2 2 2 2 2 ( ) τ 8 f b m m B B C v u v z d v g t h _{d u v} d v K v d η _ρ ρ μ ⎛ ₊ ⎞ ⎛∂ + ⎞ ∂ _{= − −⎜ + ⎟} ⎜ _{⎟ ⎜ ⎟} ∂ _⎝ ∂ _{⎠ ⎝ +} + _⎠ (3.7) 和 2 1 1 2 ( ) ( ) 0 h ud h vd d h h t ξ η ∂ ∂ ∂ _{+ +} ∂ ∂ ∂ = (3.8) 針對n+1時刻的水深值(_dn+1_{)做線性化處理，且僅保留一階項，(3.8)} 式可改寫成 1 2 1 1 1 2 2 2 ( ) ( ) 0 d d d h h d d t ξ α ξ β γ η α η β γ ⎛ ⎞ ⎛ ∂ ∂ ∂ Δ ∂ ∂ Δ + _⎜ + Δ + _⎟+ _⎜ + Δ + ∂ ∂ _⎝ ∂ _⎠ ∂ _⎝ ∂ ⎞ = ⎟ ⎠ (3.9) 式中 2 1 1 n h g t d C hτ α = − Δ ； 12 1 2 2 1 1 [ ] n n n z_b h h g t d u C_τ C h_τ ∂ ∂ β ∂ξ ∂ξ + + Δ = − + ；_{γ1 =β1d}n_；

1 2 2 n h g t d C hτ α = − Δ ； 1 2 1 1 1 2 2 [ ] n n n zb h h g t d v C_τ C h_τ ∂ ∂ β ∂η ∂η + + Δ = − + ；_{γ2 =β2d}n； 1 1 2 2 2 2 ( ) ( ) 1 n n f n C u v C t d τ + ₊ + = + Δ ；_{Δ =d d}n+1_−dn_。 3.2 數值差分式 在數值差分方法選用考量上，利用顯示數值方法求解時，演算時 間間隔受到很大的限制，在模擬天然明渠水流問題時將耗費冗長的演 算時間與龐大的電腦計算量，在應用上有其困難存在，因此，本研究 採用隱式數值方法求解。 本模式採用控制體積(control volume)法的觀念來離散控制方程 式，控制體積法的基本概念如圖3.1所示，其中(a)圖為實際區域，(b) 圖為計算區域，E、W、N、S 表相鄰格點，e、w、n、s 表控制面。 在水理控制方程式中，除了移流項採用一階精度混合型上風法(hybrid scheme)(Spalding 1972)差分外，所有空間差分均採用二階精度的中央 差分法。另外，時間項則採用簡單的前向差分方法。 中央差分法可表示成 1 _{1 1} n _{n n} e p ξ + _{+ +} ⎛∂Ψ⎞ Ψ − Ψ = ⎜ _∂ ⎟ _Δ ⎝ ⎠ w ξ (3.10) 1 _{1 1} n _{n n} n s p η η + _{+ +} ⎛∂Ψ⎞ Ψ − Ψ = ⎜_∂ ⎟ _Δ ⎝ ⎠ (3.11) 式中 1 1 1 1 1, , 0.5 ( ) 0.5 ( ) n n n n e E P i j + + + + + Ψ = ⋅ Ψ + Ψ = ⋅ Ψ + Ψn1 i j + 1 , n i j + ； 1 1 1 1 , 1 0.5 ( ) 0.5 ( ) n n n n w P W i j + + + + − Ψ = ⋅ Ψ + Ψ = ⋅ Ψ + Ψ ；

1 1 1 1 , 1 , 0.5 ( ) 0.5 ( ) n n n n n N P i j + + + + + Ψ = ⋅ Ψ + Ψ = ⋅ Ψ + Ψn1 i j + 1 1 n i j + ； 1 1 1 1 , , 0.5 ( ) 0.5 ( ) n n n n s P S i j + + + + − Ψ = ⋅ Ψ + Ψ = ⋅ Ψ + Ψ ； Ψ可表為 u , v , h₁, h₂, , d z_s 和 z_b。 混合型上風法為上風法(upwind scheme)與中央差分法組合而 成，當移流效應重要時，採用上風法；移流效應不重要時，則採用中 央差分法。至於移流效應重要性的判斷，則採用格網雷諾數(mesh Reynolds number)R_x、R_y作為判斷的因子，當 Rx 或 Ry 大於2時，代表 移流效應重要，差分方法採用能反映方向性的上風法；_Rx 或 _Ry 小於 等於2時，移流效應可視為不重要，差分方法採用中央差分法。 混合型上風法應用於本研究移流項的處理可表示成 , 1 , 1, , , 1 1 0.5 (1 ) (1 ) i j n n n n n n i j i j i j i j i j x x u u h ξ h α ξ α + + − ⎡ ⎛Φ − Φ ⎞ ⎛Φ − Φ ⎤ ⎛∂Φ ⎞ = ⎢ − ⎜ ⎟+ + ⎜ ⎥ ⎜ _∂ ⎟ _{⎜ Δ ⎟ ⎜} ⎢ ⎥ ⎝ ⎠ _{⎣ ⎝ ⎠ ⎝ ⎦} 1, n ξ ⎞ ⎟⎟ Δ _⎠ (3.12) , 1 , , 1 , , 2 2 0.5 (1 ) (1 ) i j n n n n n n i j i j i j i j i j y y v v h η h α η α + + − ⎡ ⎛Φ − Φ ⎞ ⎛Φ − Φ ⎤ ⎛∂Φ ⎞_{= − + +} ⎢ ⎜ ⎟ ⎜ ⎥ ⎜ _∂ ⎟ _{⎜ Δ ⎟ ⎜} ⎢ ⎥ ⎝ ⎠ ⎣ ⎝ ⎠ ⎝ ⎦ , 1 n η ⎞ ⎟⎟ Δ _⎠ (3.13) 其中 0 2 1 2 1 2 x x x x R R R α ⎧ _≤ ⎪ ⎪ =_⎨ > ⎪ ⎪− < ⎩ − ; 0 2 1 1 2 y y y R R R α ⎧ 2 y ≤ ⎪ ⎪ = ⎨ ⎪ ⎪ > − < − ⎩ (3.14) 上列諸式中， i jn, 1 ,i j x m u h R ξ μ ρ Δ = ； i j,n 2 ,i j y m v h R η μ ρ Δ = ；μ = 流體動力黏滯 係數(dynamic viscosity)； 可表成 Φ u 或 v 。

P i,j e E i+1,j w n s W i-1,j i,j+1 N S i,j-1 P i,j e _E i+1,j s n W i-1,j w S i,j-1 N i,j+1 ξ η X Y (a) (b) 圖 3.1 控制體積法示意圖 (a)實際區域；(b)計算區域

第四章 高含砂水流模式功能測試

本章節將針對高含砂水流模式與清水流模式進行模式功能測 試。清水流模式部份謝(2003)已分別針對不同流場特性之案例，如迴 水演算、彎道、環流流場及潰壩等案例進行模式功能測試，均得到良 好的結果。因此，針對清水流模式驗證部分，本研究在此便不多加著 墨。高含砂水流模式部份，目前由於尚無合適的實驗數據可供模式驗 證，因此本研究在此採用流變關係與 RESED2D 相同之商用 FLO-2D 模式(相關理論置於附錄)來進行高含砂水流行為模擬分析比較，測試 模式模擬高含砂水流之功能。 4.1 穩態流案例測試 在此分別比較分析兩模式於不同含砂濃度中，水體之平均水深 與平均流速模擬結果，來測試模式於超、亞臨界流況中模擬高含砂水 流之功能。 4.1.1 亞臨界流案例測試 案例設計條件為矩形渠道長1000 m，渠寬 100 m，河床坡度設為 0.01，曼寧係數設定為 0.03，流變關係設定方面，分別設定賓漢黏滯 係數 1 與賓漢降伏應力 之係數項與指數項為α 1 Cv B e β μ =α 2 2 τ Cv y =α eβ 1= 0.06287、β1= 18.712、α2= 0.70421、β2= 16.121，土壤比重設定為 2.58，

示)，因此本研究在此初步將層流阻力係數K假設為 100 來進行模擬。 藉由比較分析兩模式於不同含砂濃度中，水體之平均水深與平均流速 模擬結果，來測試模式於亞臨界流況中模擬高含砂水流之功能。邊界 條件設定方面，由於FLO-2D與RESED2D兩模式在邊界條件設定上之 不同(FLO-2D下游邊界為均勻流設定，僅需設定上游邊界便可進行模 擬)，為考慮兩模式邊界條件設定之一致性，RESED2D下游邊界設定 部份便參考FLO-2D下游模擬結果來設定，如表 4.2 所示。 圖4.1(a)、(b)分別為兩模式收斂後，含砂水體之平均水深與平均 流速展示圖，可發現兩模式所模擬之結果展現高度相似性；並由圖 4.1(a)可發現兩模式於體積濃度 30% ~ 50%間水位皆有急遽增加之趨 勢；且圖4.1(b)中兩模式皆反應出因水位急遽增加而使流速大幅降低 之結果，進而展現模擬結果之合理性。將上述案例模擬所得到之流 速、水深引進(2.15)式底床剪力中，分別比較三項之大小，分析結果 如圖4.2 所示。由圖中可發現，體積濃度 20%以下之區域，賓漢流變 特性並不明顯，底床剪應力主要受到紊流項影響；體積濃度大於30% 以上的區域，賓漢流變特性明顯，底床剪應力的增加主要由賓漢降伏 應力主導。分別比較圖4.1 與圖 4.2 可發現，兩模式於體積濃度 30% ~ 50%間，水位與流速急遽變化之情形，乃是由於賓漢流變特性的影 響，致使底剪應力急遽增加，所產生之結果。

4.1.2 超臨界流案例測試 案例設計條件除將河床坡度設定為 0.3 外，其餘參數設定皆與 4.1.1 節相同。藉由比較分析兩模式於不同含砂濃度中，水體之平均 水深與平均流速模擬結果，來測試模式於超臨界流況中模擬高含砂水 流之功能。邊界條件設定方面，由於 FLO-2D 與 RESED2D 兩模式在 邊界條件設定上之不同，為考慮兩模式邊界條件設定之一致性， RESED2D 上游邊界設定部份參考 FLO-2D 上游模擬結果來設定，如 表4.3 所示。 圖4.3(a)、(b)分別為兩模式收斂後，含砂水體之平均水深與平均 流速展示圖，由圖中可發現兩模式模擬結果展現高度相似性。且由圖 4.3(a)可發現兩模式於體積濃度 30% ~ 50%間水位皆有急遽增加之趨 勢；圖4.3(b)中兩模式皆反應出因水位急遽增加而使流速大幅降低之 結果，展現模擬結果之合理性。將上述案例模擬所得到之流速、水深 引進(2.15)式底床剪力中，分別比較三項之大小(如圖 4.4 所示)。由圖 中亦可發現，體積濃度20%以下之區域，底床剪應力主要由紊流項主 導；體積濃度大於30%以上的區域，賓漢流變特性明顯，底床剪應力 的增加主要由賓漢黏滯項主導。分別比較圖 4.3 與圖 4.4 可發現，可 發現兩模式於體積濃度 30% ~ 50%間，水位與流速急遽變化之情形，

果。 4.2 湧波傳遞案例測試 為比較分析兩模式於不同體積濃度中，湧波傳遞之模擬結果，本 研究在此設計一長2000 m，寬 100 m矩形渠道，底床曼寧係數假設為 0.03，流變關係設定方面，分別設定賓漢黏滯係數μ_B 賓漢 降伏應力τy 係數項與指數項為α 1 1 Cv eβ α = 與 Cv eβ α = 2 之 2 1= 0.00047、β1= 22.1、α2= 0.157、β2= 17.1(東埔蚋溪流變關係，2007)，土壤比重設定為 2.58， 層流阻力係數K設定方面，參考表 4.1 將其假設為 100 來進行模擬， 初始乾床，邊界條件為假設之變量流，如圖4.5 所示。底床坡度設定 方面，因濃度較低的案例其坡度效應會大於底床剪應力，為比較水體 在靜止時之深度分佈，因此將河床坡度設定為平床；在濃度較高的案 例中，將河床坡度設定為 0.003，以藉由坡度效應來分析坡度驅動對 底床剪應力之影響。 圖 4.6 至圖 4.7 分別為 FLO-2D 與 RESED2D 兩模式於體積濃度 25%與 30%之模擬結果，圖中分別顯示兩模式在第 0.5 及 3 小時的水 位與距離分佈。比較兩模式在第 3 小時水位穩定後之結果可發現，體 積濃度 25%之模擬結果較體積濃度 30%更為相近；體積濃度 30%時 兩模式前端行進距離約相差 7.7%；在上游水位部份兩模式水位壅生 高度亦非常相似。此外，比較兩模式在第0.5 小時水深分佈，可發現

不論在體積濃度 25%或是體積濃度 30%時，水深分佈趨勢皆有所差 異。然而，就水體分佈型態與湧昇波前形狀來看，高含砂水流在前端 的運移行為有堆高的現象實屬正常，顯示 RESED2D 所模擬之結果較 具合理性。 圖4.8 為加入底床坡降驅動後所模擬之結果，由圖中可發現兩模 式在收斂過程中水深之分佈亦會有所差異，就水體分佈型態與湧昇波 前形狀來看，顯示RESED2D 所模擬之結果較合理。由圖中可發現兩 模式在第3 小時之水深分佈與行進距離都呈現非常高之相似性，進而 顯示RESED2D 模式與 FLO-2D 模式之相似性。 藉由與FLO-2D 模式進行模擬比較分析後，發現兩模式雖在收斂 過程水深分佈有所差異，但在水位趨於穩定後兩模式模擬結果則非常 相似。

表4.1 層流阻力係數參考表（FLO-2D 使用手冊，2006） 層流阻力係數K 地面條件 K 值範圍 混凝土/瀝青 24-108 裸露的砂土 30-120 表層級配 90-400 被侵蝕的裸露黏土-被侵蝕的壤土 100-500 稀疏植被 1000-4000 矮草原地 3000-10000 早熟禾屬植物草地 7000-50000 表 4.2 亞臨界流測試案例邊界條件 案例 體積濃度 Cv 上游邊界₂ ( ) q m s 下游邊界 ( ) d m 福祿數 Fr 1 0.0 0.8 0.43 0.914 2 0.1 0.8 0.43 0.914 3 0.2 0.8 0.46 0.816 4 0.3 0.8 0.60 0.552 5 0.4 0.8 1.03 0.242 6 0.5 0.8 2.24 0.077 表 4.3 超臨界流測試案例邊界條件 案例 體積濃度 Cv 上游邊界₂ ( ) q m s 上游邊界 ( ) d m 福祿數 Fr 1 0.0 2.0 0.295 4.608 2 0.1 2.0 0.295 4.608 3 0.2 2.0 0.275 4.438 4 0.3 2.0 0.300 3.796 5 0.4 2.0 0.430 2.285 6 0.5 2.0 0.730 1.017

Subcritical flow Concentration (%) 0 10 20 30 40 50 60 A v e rag e flo w d epth (m) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 RESED2D FLO-2D (a) Subcritical flow Concentration (%) 0 10 20 30 40 50 60 Aver

age flow velocity

( m /s) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 RESED2D FLO-2D (b)

Subcritical flow Concentration (%) 0 10 20 30 40 50 60 Bed shear st ress t e rm 0 20 40 60 80 100

Yield stress term (RESED2D) Bingham viscosity term (RESED2D) Turbulence term (RESED2D) Yield stress term (FLO-2D) Bingham viscosity term (FLO-2D) Turbulence term (FLO-2D)

(

)

2 Nm 圖 4.2 亞臨界流底床剪應力項與體積濃度關係圖 Supercritical flow Concentration (%) 0 10 20 30 40 50 60 A v e rag e flo w d epth (m) 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 RESED2D FLO-2D (a)

Supercritical flow

Concentration (%)

0 10 20 30 40 50 60

Aver

age flow velocity

( m /s) 0 5 10 15 20 RESED2D FLO-2D (b) 圖 4.3 超臨界流平均水深與平均流速展示圖 Supercritical flow Concentration (%) 0 10 20 30 40 50 60 B e d sh e a r st re ss t e rm 0 200 400 600 800 1000

1200 Yield stress term (RESED2D) Bingham viscosity term (RESED2D) Turbulence term (RESED2D) Yield stress term (FLO-2D) Bingham viscosity term (FLO-2D) Turbulence term (FLO-2D)

(

)

2

Flow Rate Time(hr) 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Q(CM S ) 0 10 20 30 40 50 Inflow hydrograph 圖 4.5 高含砂水流測試案例邊界條件 Cv=25% distance(m) 0 500 1000 1500 2000 D epth(m) 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 RESED2D 0.5Hrs RESED2D 3.0Hrs FLO-2D 0.5hrs FLO-2D 3.0hrs RESED2D 0.5hr RESED2D 3hr FLO-2D 0.5hr FLO-2D 3hr Distance (m) 圖 4.6 高含砂水流模式測試比較-體積濃度 25%

Cv=30% distance(m) 0 500 1000 1500 2000 De pt h(m ) 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 _{RESED2D 0.5hrs} RESED2D 3.0Hrs FLO-2D 0.5hrs FLO-2D 3hrs RESED2D 0.5hr RESED2D 3hr FLO-2D 0.5hr FLO-2D 3hr Distance (m) 圖 4.7 高含砂水流模式測試比較-體積濃度 30% Cv=40% Distance (m) 0 500 1000 1500 2000 De pth (m) 0 1 2 3 4 5 6 7 RESED2D 0.5hr RESED2D 3hr FLO-2D 0.5hr FLO-2D 3hr Channel bed

第五章 高含砂效應影響程度之模擬分析

為探討高含砂效應對渠流流況之影響程度，本研究首先藉由因次 分析決定相關無因次參數，再經由參數重要性分析決定影響渠流流況 之重要參數，最後並建立高含砂效應與重要參數之關係圖，藉此釐清 高含砂效應對渠流流況之影響程度，並可提供與模式使用者作為研選 合適模式之參考。 5.1 因次分析 藉由比較高含砂水流模式與清水流模式間水深、流速與剪力差異 性，可探討高含砂效應對渠流流況之影響程度，並可提供與模式使用 者作為選用合適模式之參考依據。為探討高含砂效應對渠流流況之影 響程度，本研究在此定義三個無因次參數進行比較分析，分別為水深

增量_MaxH*_{、流速增量MaxU}*_{與剪力增量Maxτ}*_{，並將上述無因參數定}

義如下： * ( m w) m H H MaxH H − = (5.1) * ( m w) m U U MaxU U − = (5.2) * ( m w) m Maxτ τ τ τ − = (5.3) 式中，H_m、U_m與τ_m分別為為高含砂水流模式模擬結果之水深、流速 與剪力，H_w、U_w與τ_w則為清水流模式模擬結果之水深、流速與剪力。 為決定影響高含砂效應的相關物理參數，本研究將高含砂效應影

響因子分為四類，分別為流體性質、土壤特性、水力特性與幾何型態。 流體性質包括流體密度ρ、黏滯係數μ、體積濃度Cv、賓漢黏滯係數 係數項α₁、賓漢黏滯係數指數項β₁、賓漢降伏應力係數項α₂、賓漢降 伏應力指數項β₂等；土壤特性包括土砂比重G_s、層流阻力係數K等； 水力特性包括平均速度U及V 、平均水深H、渠道寬度B、渠道長度 L、渠道中心線曲率半徑 、渠道坡度 、重力加速度 等；幾何形 態包含底床型態之糙度係數 。彙整上述影響因子後可將與 c r S₀ g c MaxH*、 * MaxU 及_Maxτ*_{相關的因子整理為} * * * 1 1 1 2 2 0 , , ( , , , , , , , _s, , , , , , , , ,_c , )

MaxH MaxU Maxτ = f ρ μ Cvα β α β G K U V H B L r S g c (5.4)

式中，c=Chezy 糙度係數=R1 6 n ；n=曼寧糙度係數；R=水力半徑。 在本研究採用正交曲線座標的情況下，V 之效應可忽略不計。因 此，(5.4)式在剩下 17 個獨立變數的情況下，利用柏金漢(Vaschy - Buckingham)理論可得 14 個無因次參數，可表示成 1 1 2 2 Re=ρUH μ,Cv, , ,α β α β, ,Gs K, ,H _B,θ_b =L (2πr_c) 2 * 0 ( _c) (_{c f} ), , , _f SI UH u r= =H r C S Fr U= gH C =g c (5.5) 經由柏金漢理論雖可獲得(5.5)式之 14 個無因次參數，但尚應考 量各參數無因次參數本身所代表之物理意義，因此本研究在此將(5.5) 式改寫為

Re , 2, * , , (2 ) 8 m B B w m b w H K U UH U L r H B ρ μ c ρ μ τ ρ ρ θ π ρ ⎛ ⎞ = _⎜ + _⎟ = = ⎝ ⎠ 2 * 0 ( _c) (_{c f} ), , , _f SI UH u r= =H r C S Fr U= gH C =g c (5.6) 式中，Re=雷諾數(Reynolds number)； 2 8 B B K U U H μ τ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ ρw =無因次非牛 頓流變剪應力項； * m ρ =無因次含砂水流密度；H B =寬深比(depth-width ratio)；θ_b=彎道長度因子(relative length of channel bend)； =二次流

強度因子(relative strength of secondary current)； =福祿數(Froude

number)； SI Fr f C =摩擦因子(friction factor)；u_*= C_fU=剪力速度。 因此，可將無因次分析之結果表示為 * * * 2 * 2 0 , , (Re, , , , , , , , ) 8 B B w m b H K U f

MaxH MaxU Max f U SI S Fr C

H B μ τ = ⎛_⎜τ + ⎞_⎟ ρ ρ θ ⎝ ⎠ (5.7) 其中，(5.7)式若用於分析直線道案例時則不將 與SI θ_b兩無因次參數納 入影響分析，如下所示 * * * 2 * 3 0 , , (Re, , , , , , ) 8 B B w m H K U f

MaxH MaxU Max f U S Fr C

H B

μ

τ = ⎛_⎜τ + ⎞_⎟ ρ ρ

⎝ ⎠ (5.8)

_MaxH*_、MaxU*_及Maxτ*_{可用來作為判斷高含砂效應的影響程度，}

並可作為高含砂水流與清水流模式間模式選用之判斷依據。 5.2 模擬案例設定 為探討上述(5.7)、(5.8)兩式中，各無因次參數對_MaxH*_、MaxU*_及 * Maxτ 影響程度，本研究根據模擬參數之統計特性(如表 5.1)，經由蒙 地卡羅法衍生 928 組測試案例進行模擬分析工作(如表 5.2 所示)，其 中所衍生出資料型態假設為常態分佈(normal distribution)。表 5.1 中

1 α 、β₁、α₂、β₂與K之設定為參考表 5.3 與表 4.1 來設定參數模擬範 圍。在每個案例中，渠道幾何均為矩形斷面之單一彎道案例，且在此 單一彎道前後各連接一長10 m 之直線渠道(如圖 5.1 所示)，以避免受 下游邊界之影響；邊界條件設定方法如表2.1 所示。 本研究測試案例所涵蓋的範圍：Re範圍為3.05×103~ 2.65×106、 / H B範圍為2.46×10-3 ~ 1.567、θ_b範圍為0 ~ 0.25、 範圍為 0.002 ~ 0.645、Fr範圍為0.093 ~ 5.260、 SI f C 範圍為2.11×10-3 ~ 9.69×10-2、Cv範 圍為 1% ~ 50%、α₁範圍為 4×10-5 ~ 1.254×10-2、β₁範圍為 8.35016 ~35.99127、α₂範圍為0.02568 ~ 2.58988、β₂範圍為4.10035 ~ 25.871、 K範圍為31.963 ~ 476.056、G_s範圍為2.6 ~ 2.9。由設定範圍可知，本 研究所設計之案例範圍相當廣，可概括性描述各種水理發生情形，以 貼近真實河川。 5.3 參數重要性分析

圖 5.2 至圖 5.4 分別為直線道_MaxH*_、MaxU*_及Maxτ*_{對(5.8)式中}

各無因次參數之關係圖，由圖中分別可發現 2 8 B B w K U U H μ τ ρ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ 對 *

MaxH 、_MaxU*_與Maxτ*_{有較明顯的關係趨勢。此外，由圖5.4 中可發}

現除 2 8 B B K U U H μ τ ρ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ w 外， * m ρ 對_Maxτ*_{亦有明顯關係趨勢存在，針對} * m ρ 與 2 8 B B w K U U H μ τ ρ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ 間對 * Maxτ 之影響分析將於後續章節詳加探

圖 5.5 至圖 5.7 分別為彎道_MaxH*_、MaxU*_及Maxτ*_{對(5.7)式中各} 無因次參數之關係圖，由圖亦可發現 2 8 B B w K U U H μ τ ρ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ 對 * MaxH 、 * MaxU 與_Maxτ*_{有較明顯之關係趨勢。此外，由圖 5.7 亦可發現除} 2 8 B B w K U U H μ τ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ ρ 外 * m ρ 對_Maxτ*_{亦有明顯關係趨勢存在，關於兩無因} 次參數對_Maxτ*_{之影響分析同樣將於後續章節詳加探討。}

5.4 高含砂效應對MaxH*、_MaxU*_及Maxτ*_影響分析

5.4.1 福祿數效應影響分析

為探討福祿數效應對_MaxH*_、MaxU*_與Maxτ*_{之影響，本研究在此}

根據表5.2 所設計之 928 組測試案例，來模擬分析福祿數效應對水深

增量_MaxH*_{、流速增量MaxU}*_{及剪力增量Maxτ}*_之影響。

圖5.8 為福祿數效應對直線道_MaxH*_{影響分析圖，由圖中可發現}

Fr＜1 與 Fr＞1 兩關係曲線隨著高含砂效應變小而略顯差異，Fr＜1

之曲線皆高於Fr＞1 的曲線，但整體仍呈現相同之趨勢。為判斷上述

Fr＜1 與 Fr＞1 兩曲線差異性，本研究在此藉由均方根差(root mean

square error)E_rms與相關係數(correlation coefficient)ρ_c來作為判斷兩曲

線相似性的指標。 1/ 2 2 ( ) total i i N rms total x y E N ⎛ − ⎞ ⎜ = ⎜ ⎜ ⎟ ⎝ ⎠

∑

_⎟ ⎟ (5.9) 2 ( )( ) ( ) ( i i c i i x X y Y 2 ) x X y Y ρ = − − − −

∑

∑

∑

(5.10)

其中，x_i表某比較點回歸曲線值，y_i表相對應比較點回歸曲線值，X、 表平均數， 表比較值的數量。 可視為兩曲線的平均誤差， 其值越小表示兩曲線性差異越低； Y N_total E_rms c ρ 越趨近於1 則表示兩曲線相關性 越高。 經由上述定義發現，圖 5.8 中兩曲線E_rms及ρ_c值約為 0.045 與 0.994，顯示兩曲線相似性甚高，亦表示兩曲線間差異性極為微小。 圖5.9 為福祿數效應對直線道_MaxU*_{影響分析圖，由圖中可發現Fr＜} 1 與 Fr＞1 兩關係曲線呈現些微之差異，兩關係曲線E_rms及ρ_c值約為 2.597 與 0.996。圖 5.10 為福祿數效應對直線道_Maxτ*_{影響分析圖，由} 圖中可發現Fr＜1 與 Fr＞1 兩關係曲線亦隨高含砂變效應變小而有所 差異，且Fr＞1 之曲線皆高於 Fr＜1 的曲線，與_MaxH*_{呈相反之趨勢，} 兩關係曲線E_rms及ρ_c值約為 0.063 與 0.998，顯示兩關係曲線差異性甚 低。

圖5.11 至圖 5.13 分別為福祿數效應對彎道_MaxH*_、MaxU*_與Maxτ*

影響分析圖，由圖中分別可發現Fr＜1 與 Fr＞1 兩曲線趨勢皆有明顯

的差異性，不同關係曲線之E_rms及ρ_c值如表 5.4 至表 5.6 所示。探究

造成此差異性之緣由，乃是因為模式在延散剪應力項中，流速剖面之

面適用性恐會有所問題。因為超臨界流流速較快，易受邊牆幾何因素 的作用力影響，而產生交波(cross wave)之現象。這樣的現象恐會使流 況無法維持完全發展流狀態，造成de Vriend 所假設之流速剖面產生 不適用之情況，因而使模擬上產生誤差，致使兩關係曲線呈現較大之 差異。 綜觀上述，Fr＜1 與 Fr＞1 兩關係曲線於直線道時雖有些微之差 異，但整體趨勢仍非常相似，顯示超、亞臨界流況於直線道之差異性 甚低。反觀彎道時，由於超臨界流在流經彎道所產生交波現象之影 響，造成兩關係曲線有較大的差異性，為避免此現象所產生之模擬誤 差，因此在接續章節之分析中，便不將超臨界流納入彎道水深增量 *

MaxH 、流速增量_MaxU*_{及剪力增量Maxτ}*_{比較分析中。}

5.4.2 直線道案例影響分析 圖 5.14(a)～(c)分別為直線道 2 8 B B w K U U H μ τ ρ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ 對 * MaxH 、_MaxU*

及_Maxτ*_{關係圖。由圖中分別可發現MaxH}*_、MaxU*_及Maxτ*_均會隨 2 8 B B K U U H μ τ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ ρw 增大而變大，各參數與 2 8 B B K U U H μ τ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ ρw 間為正比 關係。上述關係圖亦可用為判斷在不同流變參數下，高含砂水流相較 於清水所產生的水深、流速及剪力增量值，以提供與模式使用者在分 析直線道高含砂水流問題時作為參考。此外，在考慮關係曲線本身所 代表之物理意義下，本關係曲線適用範圍為_{Ln MaxH(} *_{)及Ln Max( τ}*_)小

於0 之部份，大於 0 之範圍則不適用。 5.4.3 彎道案例影響分析 圖5.15 至圖 5.17 分別為彎道不同θ_b之 2 8 B B w K U U H μ τ ρ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ 對 * MaxH 、 *

MaxU 與_Maxτ*_{關係圖，其中MaxH}*_、MaxU*_及Maxτ*_{之值係取自與長}

度因子θ_b相對應之斷面處。由圖中可看出各θ_b之_MaxH*_、MaxU*_及 * Maxτ 均會隨 2 8 B B K U U H μ τ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ ρw 增大而變大， *

MaxH 、_MaxU*_及Maxτ*_各

參數與 2 8 B B K U U H μ τ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ ρw 間為正比關係；圖 5.15 至圖 5.17 中之關係 曲線圖可作為判斷在彎道中不同θ_b處，於不同流變參數下，高含砂水 流相較於清水所產生的水深、流速及剪力增量值，以提供與模式使用 作為參考。然而，在考慮圖中各關係曲線本身所代表之物理意義下， 關係曲線圖適用範圍為_{Ln MaxH(} *_{)及Ln Max( τ}*_{)小於 0 之部份，大於 0} 之範圍則不適用。 圖5.18 為將圖 5.15 中不同θ_b處之 2 8 B B w K U U H μ τ ρ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ 與 * MaxH 關係 曲線繪於同一張圖上，在固定 2 8 B B K U U H μ τ ⎛ ₊ ⎞ ⎜ ⎟ ⎝ ⎠ ρw 的情況下，由圖中可 發現_MaxH*_會隨著 b θ 增大而變小。藉由E_rms及ρ_c判斷不同θ_b關係曲線 後發現，E_rms與ρ_c約為0.172 及 0.987，各關係曲線間雖有些許差異， 但整體仍呈現相同之趨勢。依據上述方法同樣可將圖5.16 與圖 5.17 繪製如圖5.19 與圖 5.20 之關係曲線圖，由圖 5.19 中亦可發現在固定 Kμ U ⎛ ⎞