國立交通大學
土木工程學系
碩士論文
二維彎道動床模式之發展研究
Development of A 2-D Depth-Averaged Mobile-Bed Model
for Channel Bends
研 究 生：張益家
指導教授：楊錦釧 教授
謝德勇 博士
二維彎道動床模式之發展研究
Development of A 2-D Depth-Averaged Mobile-Bed Model
for Channel Bends
研 究 生：張益家 Student: Yi-Jia Chang
指導教授：楊錦釧 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 2005
誌 謝
承蒙吾師楊教授錦釧及謝博士德勇於論文研究期間之悉心指 導、諄諄教誨，得以使此研究順利完成。在學期間，得吾師、葉教授 克家與許博士勝田於課業上之傳授和為人處事之啟發，讓學生在兩年 交大生活獲益匪淺。 於論文審定期間，感謝口試委員台大生工所許教授銘熙、美國休 士頓大學土木與環境工程系王教授克漢及台大水工所賴博士進松之 細心匡正、會賜卓見，使本論文更臻完善。 此外，感謝實驗室東霖學長、夢祺學長、恩添學長、祥禎學長、 胤隆學長、世偉學長及曉萍學姊於課業和生活上之指導，亦感謝同學 弘恩、俊毅、峰志、于軒及璨仲於此兩年裡陪伴我一同成長努力，及 學弟妹欣瑜、浩榮、仲達、宣汝、力瑋等之陪伴。 另外感謝好友哲一及嘉俊不斷予我鼓勵和經驗分享；而交大登山 社所有伙伴於這兩年裡與和我一同享受徜徉在山林的懷抱中，體驗最 真實的生命，於此特別提出致意。 最後感謝我生活在彰化田中的爸媽、大哥、大嫂、姊姊、姪女汶 琳及姪子廷瑋，你們是我努力成長的最大支柱，由衷感謝你們長久以 來的包容和關心，讓我在求學的路上無後顧之憂，謹以學位論文與你 們分享我成長的喜悅，感恩之心無以言喻。二維彎道動床模式之發展研究
學生：張益家 指導教授：楊錦釧 謝德勇 國立交通大學土木工程研究所摘 要
本研究旨在發展一水深平均水平二維之彎道動床模式。模式在 發展過程中，通盤考量二次流之效應，水理演算不忽略延散剪應力 項；懸浮載之運算包含環流傳輸項；底床載的演算，則將其運移之修 正角度納入考量，包含流速、底床坡降及曲率半徑之影響。在數值處 理上，水理採用雙階分割操作趨近法，將水理控制方程分割成延散步 驟和傳播步驟求解，數值差分則是採用隱式法，使模式可採用較大之 演算間距；沉滓運移則是採用耦合演算法(同時求解所有沉滓運移控 制 方 程 式 ) ， 利 用 Newton-Raphson 疊 代 聯 立 求 解 。 藉 由 模 擬 Struiksma(1983) DHL-T2 之彎道實驗，發現在固定參數下，模式並無 法適切模擬出彎道入口段過度沖淤(overshoot)之現象。經由模式彎道 二次流效應之敏感度分析，發現底床載二次流效應對底床沖淤之影響 最大。因此，透過調整底床載二次流效應之程度，可檢定彎道模擬時 之合適底床載參數。最後，將此參數檢定架構應用於 DHL-T1 及 DLFM 實驗之模擬，驗證此架構之合宜性。關鍵字：水深平均、動床、二次流、過度沖淤、彎道
Development of A 2-D Depth-Averaged Mobile-Bed Model
for Channel Bends
Student：Yi-Jia Chang Advisors：Jinn-Chuang Yang Te-Yung Hsieh
Department of Civil Engineering National Chiao-Tung University
ABSTRACT
The purpose of this study is to develop a 2-D depth-averaged mobile-bed model for channel bends. To the model, the secondary current effect is thoroughly considered for the water and sediment movements, which include the dispersion stress term in the flow momentum equation, the circulatory transport term in the suspension transport equation, and the angle of sediment transport in the bed load transport equation. The angle of sediment transport includes the effects of the flow velocity, the bed slope, and the radius of channel bend. As for the numerical solution procedure, the split-operator approach is used to solve the flow equations. Based on the decomposition of the momentum equations, the split-operator procedure involves two steps in which one is the dispersion step, and the other one is propagation step includes the pressure and bed friction terms. Implicit discretized method is adopted to relax the time step restriction to allow large time steps. The direct coupling, i.e. simultaneous solution of the sediment equations and Newton-Raphson methods, is used to solve the sediment governing equations. By using the
uniform parameters in the DHL-T2 case of Struiksma’s(1983) experiment, the simulated results can’t demonstrate the overshoot phenomenon near the entrance of channel bend. The sensitive analysis indicates that the secondary-current effect in the bed-load transport equation has significant effect on the bed evolution. Hence, the nonuniform parameters of the secondary-current effect in the bed load transport equation are the key for the model calibration. Finally, the calibration process is applied to DHL-T1and DLFM cases to assess its rationality.
Key words: depth averaged, mobile bed, secondary current, overshoot, channel bend
目錄
誌謝... I 中文摘要... II 英文摘要...III 目錄...V 表目錄... VII 圖目錄...VIII 符號表...X 第一章 緒論...1 1.1 研究動機與方向 ...1 1.2 文獻回顧 ...1 1.3 研究目的與方法 ...4 1.4 章節介紹 ...5 第二章 理論基礎...6 2.1 水理部分 ...6 2.1.1 控制方程式 ...6 2.1.2 輔助關係式 ...8 2.1.3 邊界條件 ...9 2.2 沉滓部分 ...10 2.2.1 控制方程式 ...10 2.2.2 輔助關係式 ...11 2.2.3 邊界條件 ...15 第三章 數值架構...18 3.1 水理部分 ...183.1.1 隱式雙階分割操作趨近法 ...18 3.1.2 數值差分式 ...19 3.2 沉滓運移部分 ...21 3.2.1 耦合演算法 ...21 3.2.2 數值差分式 ...22 3.2.3 數值解析法 ...25 第四章 模式適用性及二次流效應之敏感度測試分析 ...29 4.1 模式適用性測試分析及問題解析...29 4.2 二次流效應之敏感度測試分析...31 4.2.1 水理二次流效應分析 ...32 4.2.2 滑移邊界效應分析 ...34 4.2.3 懸浮載二次流效應分析 ...35 4.2.4 底床載彎道效應分析 ...37 第五章 底床載源參數檢定架構探討及應用分析 ...49 5.1 參數檢定架構探討 ...49 5.1 參數檢定應用分析 ...52 第六章 結論與建議...64 6.1 結論 ...64 6.2 建議 ...65 參考文獻...66
表目錄
表 4.1 不同θb之 * MaxU 與SI之回歸函數一覽表(謝，2002) ...39 表 4.2 二次流效應對彎道動床之影響分析一覽表 ...39 表 4.3 二次流彎道效應對底床變動影響分析一覽表 ...40 表 5.1 角度修正因子之參數變動對彎道動床影響一覽表_1( fs =2)....54 表 5.2 角度修正因子之參數變動對彎道動床影響一覽表_2(A* =9.1) .54圖目錄
圖 2.1 彎道二次流示意圖...17 圖 2.2 作用層示意圖...17 圖 3.1 控制體積法示意圖 (a)實際區域；(b)計算區域 ...27 圖 3.2 質點運移軌跡示意圖...28 圖 3.3 質點運移軌跡在計算平面之示意圖 ...28 圖 4.1 Nico Struiksma (1983) DHL-T2 實驗水槽幾何形狀圖...41 圖 4.2 模擬 DHL-T2 實驗之徑向斷面非均勻格網示意圖 ...41 圖 4.3 模擬 DHL-T2 實驗彎道段之水面高程示意圖 ...42 圖 4.4 模擬 DHL-T2 實驗之彎道縱向底床變化圖 ...42 圖 4.5 模擬 DHL-T2 實驗之彎道縱向水深示意圖 ...43 圖 4.6 不同θb之 * MaxU 與SI關係圖(謝，2002) ...44 圖 4.7 水理二次流效應對縱向水深影響示意圖_1(DHL-T2 實驗)...44 圖 4.8 水理二次流效應對縱向水深影響示意圖_2(DHL-T2 實驗)...45 圖 4.9 考慮滑移邊界與否之縱向水深比較圖(DHL-T2 實驗)...45 圖 4.10 懸浮載二次流效應對縱向水深影響示意圖_1(DHL-T2 實驗).46 圖 4.11 懸浮載二次流效應對縱向水深影響示意圖_2(DHL-T2 實驗).46 圖 4.12 懸浮載二次流效應對縱向水深影響示意圖_3(DHL-T2 實驗).47 圖 4.13 模擬不同粒徑對沖淤影響示意圖(DHL-T2 實驗)...47 圖 4.14 考慮底床載彎道效應與否之縱向水深比較圖(DHL-T2 實驗).48 圖 5.1 參數修正對縱向水深影響示意圖(DHL-T2 實驗)...55 圖 5.2 參數調整對縱向水深變動量之影響示意圖(DHL-T2 實驗)...55 圖 5.3 參數調整對彎道縱向平均水深差異示意圖(DHL-T2 實驗)...56 圖 5.4 參數修正對縱向水深影響示意圖_ for A(DHL-T2 實驗) ...56圖 5.5 參數修正後演算結果之縱向底床變化比較圖(DHL-T2 實驗)...57 圖 5.6 參數修正後演算結果之縱向水深比較圖(DHL-T2 實驗)...57 圖 5.7 參數修正後底床變動(ΔZb/ h0)地形圖(DHL-T2 實驗) ...58 圖 5.8 參數修正後彎道水深示意圖(DHL-T2 實驗)...58 圖 5.9 參數修正後彎道流速分佈圖(DHL-T2 實驗)...59 圖 5.10 參數修正後演算結果之縱向水深比較圖(DHL-T1 實驗)...59 圖 5.11 參數修正後縱向底床變化之時變圖(DHL-T2 實驗)...60 圖 5.12 參數修正後之縱向水深比較圖(DHL-T1 實驗)...60 圖 5.13 參數修正後底床變動(ΔZb/ h0)地形圖(DHL-T1 實驗) ...61
圖 5.14 Sutmuller and Glerum (1980)實驗水槽幾何形狀 ...61
圖 5.15 參數未修正前演算結果之縱向水深示意圖(DLFM 實驗)...62
圖 5.16 參數修正後演算結果之縱向水深示意圖(DLFM 實驗)...62
符號表
A = 渦流係數； a = 沙丘高度之一半； B = 渠道寬度； C = 濃度； Cr = Courant number； f C = 摩擦係數； 0 C = 常數； c = Chezy 係數； 1 c = 顆粒蔡司係數； k D = 顆粒 k 之粒徑； m D = 不產生移動的最小顆粒粒徑； *k D = 無因次顆粒粒徑； d = 水深； 0 d = 平均水深； E = 糙度因子； m E = 作用層厚度； mzp E = 底床變動指標； rms E = 均方根差； 2 e = 環流傳輸係數； d F = 可沖刷的深度； s f = 沉滓顆粒之形狀因子； g = 重力加速度； H = 平均水深；1 h 、h2 = ξ、η 方向轉換係數； k = von Karman’s 係數； L = 渠道長度； r L = 環流流場的環流長度； ' M = 沖刷係數； * MaxU = ξ方向速度之最大相對差異； total N = 模擬案例總計算格點數目； n = 曼寧糙度係數； P = 沉降機率； p = 孔隙率； bi q = i 方向某一粒徑之底床載通量； R = 水力半徑； c R = 渠道中心線平均曲率半徑； r = 曲率半徑； c r = 渠道中心線曲率半徑； S = 懸浮載源； SI = 二次流強度因子； s = 砂比重； f S = 作用層源； 0 S = 渠道坡度； k T = 輸送參數； 11 T 、T12、T22 = 有效剪應力項； t = 時間； U = ξ 方向平均速度； u = ξ 方向速度；
w u = 近固體邊界的水深平均速度； * u = 剪力速度； *c u = 臨界剪力速度； V = η 方向平均速度； v = η 方向速度； lk w = 顆粒 k 之躍起速度； fk w = 顆粒 k 之沉降速度； , i j X = 某計算格點於考慮指定模擬條件時之模擬結果； , i j Y = Xi j, 相對應計算格點於未考慮指定模擬條件時之模擬結果； w y = 固體邊界與鄰近固體邊界格點的距離； b z = 底床高程； 0 b z = 原始底床高程； bip z = 某計算格點於考慮指定模擬條件時之底床高程； bnp z = 某計算格點於未考慮指定模擬條件時之底床高程； s z = 水面高程； sm z = 平均水面高程； α = 底床載源運移角度； β = 粒徑百分比； δ = 底床剪應力之方向； s ε = 沉滓運移模擬之收斂容許誤差值； 1 ε 、ε2 = ξ、η方向之亂流傳輸係數； μ = 流體動力黏滯係數； l υ = 層流黏滯係數； t υ = 亂流黏滯係數；
b θ = 彎道長度因子； θ = Shields paramerter； ρ = 流體密度； L ρ = 濕密度； s ρ = 泥砂密度； 1 b τ 、 2 b τ = 底床剪應力在 ξ 與 η 方向之分量； cd τ = 臨界沉降剪應力； cem τ = 臨界塊狀剝蝕剪應力； ces τ = 臨界表層沖刷剪應力； ξ、η = 平面上兩正交曲線座標方向； ς = 距離底床之高度與水深之比值； h ζ = 隱藏因子； t Δ = 時間間距； x Δ 、Δy = ξ、η方向之格網間距。 上標 b = BF 模式之模擬結果； m = 疊代次數； n = n tΔ 時刻之已知變數； 1 n+ = (n+ Δ1) t時刻之未知變數； 1 2 n+ = (n+ Δ1) t與n tΔ 間之未知變數； ( )= 時間平均; ( ) =水深平均； (') = 時間平均瞬時擾動量。
下標 A、D =沉滓到達和離去點； c = 模擬值； M = A、D之中點值； m = 實驗值； s = 變數在水面的值； b = 變數在底床的值。
第一章 緒論
1.1 研究動機與方向 在近十幾年來，地狹人稠的台灣飽受風災、水災、旱災等天然災 害的影響，輕則公共設施毀損，交通受阻，重則人員傷亡、重大公共 建設頹傾、民生停滯、傷及國本；台灣在 1999 年九二一集集大地震 之後，土壤鬆動情形加劇，致使每當遇豪雨或颱風時，往往造成人民 之身家安全遭受水流、土石等侵襲傷害。 坡陡流急、富含大量砂石為一向是為台灣河川之特色，並且特別 蜿蜒曲折，造成河川之流場更為複雜，因此河川流況及沉滓運移乃為 台灣水利規劃中相當重要之課題。然而當水流進入彎曲渠道時，其複 雜流場更易造成彎道外岸土石崩落，甚或掏刷堤防基礎，因此彎道水 理與輸砂之研究更顯其重要性。 本研究旨在發展一彎道動床模式，希冀能提供日後水利防災、水 利工程設計所用。 1.2 文獻回顧 對水利工程師而言，在沖積河川中彎曲河段之床形變化，是一項 困難且極為重要之課題，包含其水流特性、輸砂情形及岸壁穩定性， 均需加以探討其物理現象；而彎道二次流效應乃為影響整體彎道底床 沖淤之重點，其中包含水理控制方程之延散剪應力項、懸浮載控制方 程之環流傳輸項及底床載通量之彎道曲率半徑效應。 Rozovskii(1961)和 Yen(1965)依據理論和實驗來研究平床之彎道 特性；Yen(1967，1970)利用一蜿蜒實驗渠槽，探討床形之平衡狀況 和其對固定渠道邊界之水流發展的影響；de Vriend(1977, 1978)依據實驗探討非平床之定床彎道水理狀況，並針對二次流的現象稍作討論； Struiksma(1983)利用一彎道動床實驗，探討整體彎道沖淤情形；Ikeda and Nishimura (1986)的研究成果顯示，在砂粒-粉土為主控河床質載 之河川，懸浮載源影響了最大沖刷深度，於其所模擬之案例可影響最 大沖刷深度達 8%之程度，並認為二次流效應具有影響彎道床形相位 延遲之作用。在電子計算機發達的時代，在最近三十年來，已有眾多 學者陸續發展彎道動床數值模式，將於下一一簡述之。 早期模式多為定床模式，僅探討固定床形對水理部份之發展影 響，如 Huang 等(1967)、de Vriend (1976)、Smith and McLean (1984) 和 Ali (1985)等人研究發展；而 Leschziner and Rodi (1979)和 Tamai and
Ikeya (1985)則發展了k−ε模式。 然而，亦有許多學者從事探討水理和底床變動之交互作用的研 究，Engelund (1974)為探討渦流理論對彎道底床變動情形之影響，發 展了一適用於穩態流場之彎道動床模式，並於模式中加入擬似二次流 概念之徑向流速演算，並求算徑向底床載源之通量，期望能反應彎道 底床沖淤之現象，但在底床載演算過程中並無包含縱向流速和縱向底 床坡降之影響，且在沉滓運移之演算過程中並無考量懸浮載源之影 響；Kikkawa 等 (1976)發展一適用於穩態流場之二維彎道動床模式， 展示了以非耦合理論(分開求解水理及底床高程)模擬完全發展流況 下之彎道底床變動情形，於沉滓運移的演算過程中無考慮縱方向及橫 方向底床坡降的影響，並無考量懸浮載源之效應；Struiksma 等(1985) 發展一適用穩態流之二維數值模式，並於底床載源運移演算時，加入 曲率半徑、縱向及徑向之流速和底床坡降的影響，以數組實驗資料加 以驗證該模式演算彎道動床之能力，但此模式在水理部份，並無加入
發展一個適用穩態流之三維數值模式，認為彎道二次流和床形變化皆 導因於橫向之流體移動，並利用質量-流量平衡(mass-flux balance)之 概念連結動量方程式與彎道床形平衡方程式，彎道徑向底床坡降之演 算則僅導入曲率半徑和縱方向流速之效應，並且無考慮懸浮載源之效 應，此模式亦僅適用於同一渠寬且彎道曲率半徑遠大於渠寬之河道； Shimizu and Itakura(1989)所發展適用於穩態流之二維數值模式，在沉 滓運移演算中，此模式並無考量縱向底床坡降和懸浮載源之影響；Yen and Ho (1990)發展一個適用於穩態流之二維數值模式，該模式在沉滓 運移演算中，並不包含縱向及橫向之流速影響，且彎道處亦無納入曲 率半徑之影響，整體沉滓演算亦無包含懸浮載源之效應；Yeh and Kennedy(1993)發展一個適用於穩態流之二維數值模式，該模式乃著 重於應用 MOM(moment of momentum)方程式，將曲率半徑、橫向底 床坡降及縱向和橫向的流速等影響因子，導入沉滓運移之演算中，但 並無演算懸浮載源之效應；連氏(1999)則是考慮水理二次流效應對彎 道流場之影響，並探討延散剪應力在彎道段之分布情形，但於底床載 源之演算中則未考慮曲率半徑、底床坡降和流速之影響，於懸浮載源 之演算，則無考慮環流傳輸項之影響；Kassem and Chaudhry(2002)所 發展一適用於非穩態流之二維數值模式，將曲率半徑、縱方向及橫方 向之流速和底床坡降引入沉滓運移演算中，將有效應力與水理方程式 結合，但模式中並未包含水理二次流和懸浮載源之演算；許氏(2002) 發展一個適用於非穩態流之二維數值模式，於水理演算中導入曲率半 徑之效應，於沉滓演算中，忽略懸浮載之擴散及延散通量，於底床載 部份則未考慮曲率半徑、底床坡降和流速之影響；Duc 等(2004)發展 一適用於非穩態流之二維數值模式，並將二次流傳輸效應反應於水深 平均之k−ε 模式中的擴散係數，在底床載源的演算裡，則是包含曲率
半徑、縱方向及橫方向之流速和底床坡降之影響，並演算懸浮載源效 應，但無包含環流傳輸項之演算，且未處理其水理控制方程式之延散 剪應力項。 1.3 研究目的與方法 完整之彎道動床模式應包含：(1)水理二次流效應；(2)懸浮載二次 流效應；(3) 徑向與縱向之流速和底床坡降，及彎道中曲率半徑對底 床載源通量之影響；(4)適用於穩態流和非穩態流流場之模擬。 由 1.2 節相關文獻之分析探討中，可以發現迄今尚沒有一完整考 慮上列各項因素之彎道動床數值模式，因此本研究將發展一完整考量 前述彎道各項影響因子之動床模式。 本研究之研究重點旨在延續國立交通大學楊錦釧教授和謝德勇 博士(2002)所發展之 RESED-2D 模式，利用此模式在水理及懸浮載二 次流之完備性(包含水理部份延散剪應力項和懸浮載部份環流傳輸 項)，於沉滓運移演算的過程中植入各項彎道動床之影響因子，使此 模式在沉滓運移部份之理論基礎趨於完備，進而探討蜿蜒河道之底床 沖淤情形。此模式為一水深平均二維非穩態之水理、沉滓運移數值模 式，座標系統採用正交曲線座標，水理控制方程採用連續及動量方程 式求解，為使水理及輸砂之二次流理論更為完善，模式在發展過程中 不忽略水理延散剪應力及懸浮載源之環流傳輸項，及底床載源考慮流 速、底床坡降和彎道曲率半徑之影響；沉滓傳輸型態則是採用懸浮載 及底床載分開演算之方式。在數值方法之選用上，水理部份採用隱式 雙階分割操作趨近法，將水理控制方程分割成延散步驟和傳播步驟； 沉滓運移部份則採用耦合演算法(同時求解所有沉滓運移控制方程 式)，利用 Newton-Raphson 疊代求解。而為驗證模式之適用性及合理
性，以採用實際實驗案例之資料，予以模擬並檢定相關參數和結果可 靠度。 1.4 章節介紹 前三節已闡述本研究之研究動機和方向、文獻回顧、研究目的與 方法，以下將簡要敘述本論文各章節之內容。 第一章為緒論，針對本研究之緣起和方向作說明後，回顧相關之 文研究文獻後，再提出本研究之目的與研究方法，並於章末作本論文 架構說明。 第二章為理論基礎，分別闡述本研究之二維正交座標曲線系統所 需水理部份及輸砂部份的控制方程式、輔助方程式及相關之邊界條 件。 第三章為數值架構，說明水理部份之隱式雙階分割操作法和數值 差分式，並於此章節闡述輸砂部份之結果演算法、數值差分式及數值 解析法。 第四章為模式測試與應用，針對本研究之重點「彎道動床」部份， 進行模式在沉滓運移部份之適用性分析，並針對初步之演算結果進行 問題解析，以進行模式各項二次流效應之敏感度分析 第五章為檢定底床載參數在彎道動床模式之適用性，並擬定出合 適之底床載參數修正方式，最後再作驗證其正確性。 第六章為結論與建議，針對本研究成果作綜合性之歸納說明，並 對不盡完備或日後能改進之處提出建議。
第二章 理論基礎
透過座標系統轉換將控制方程式轉換為正交曲線座標系統方程 式，再將此控制方程式作時間平均及水深平均後，即可推得水深平 均二維正交曲線座標模式所需之控制方程式。其中水理控制方程包 含水理連續及動量方程式；沉滓運移部份則將輸砂通量分離為非均 勻質之懸浮載與底床載，所以其控制方程包含某一粒徑之懸浮載質 量守恆方程式、某一粒徑於作用層(active layer)之質量守恆方程式及 整體河床沉滓運移之質量守恆方程式。茲將水理及沉滓運移部分之 理論基礎敘述如下： 2.1 水理部分 2.1.1 控制方程式 為適度簡化複雜的控制方程式，需對數學模式作若干假設，分 別為(1)不可壓縮牛頓流體(incompressible Newtonian fluid)；(2)靜水 壓分布；(3)忽略風剪力；(4)忽略科氏力。則水深平均二維正交曲線 座標水理控制方程式可表示為 (1)水流連續方程式 1 2 ( 2 ) ( 1 ) 0 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 ) ( ) ( ) g z d h T h T h h h d h h d ∂ ∂ ∂ ∂ξ ρ ∂ξ ρ ∂η = − + + + 1 1 2 12 22 1 2 1 2 1 h 1 h b T T h h d h h d d τ ∂ ∂ ρ ∂η ρ ∂ξ ρ + − − 2 11 2 11 1 12 1 12 1 2 1 ( ) s ( ) b ( ) s ( ) b s b s b z z z z h h h h h h d τ τ τ τ ρ ξ ξ η η ⎡ ∂ ∂ ∂ ∂ ⎤ + ⎢− + − + ⎥ ∂ ∂ ∂ ∂ ⎣ ⎦ (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 ) ( ) ( ) g z d h T h T h h h d h h d ∂ ∂ ∂ ∂η ρ ∂ξ ρ ∂η = − + + + 2 1 2 11 12 1 2 1 2 1 h 1 h b T T h h d h h d d τ ∂ ∂ ρ ∂η ρ ∂ξ ρ − + − 2 12 2 12 1 22 1 22 1 2 1 ( ) s ( ) b ( ) s ( ) b s b s b z z z z h h h h h h d τ τ τ τ ρ ξ ξ η η ⎡ ∂ ∂ ∂ ∂ ⎤ + ⎢− + − + ⎥ ∂ ∂ ∂ ∂ ⎣ ⎦ (2.3) 式中， 2 2 11 [ 11 ' ( ) ] s b z z T =
∫
τ −ρu −ρ u−u dz (2.4) 2 2 22 [ 22 ' ( ) ] s b z z T =∫
τ −ρv −ρ v−v dz (2.5) 12 21 [ 12 ' ' ( )( )] s b z z T =T =∫
τ −ρu v −ρ u−u v−v dz (2.6) 以上諸式中， ξ、η = 平面上兩正交曲線座標方向；h1 = ξ 方向轉換 係數；h2 = η 方向轉換係數；u = ξ 方向速度；v = η 方向速度；ρ = 流體密度；d = 水深；g = 重力加速度； t = 時間；zb = 底床高 程；zs = 水面高程； i b τ =底床剪應力在 ξ 與 η 方向之分量; ( )= 時間 平均; ( ) =水深平均；(') = 時間平均瞬時擾動量；下標s、b 分別代 表變數在水面與底床的值； 11 T ,T12,T22 = 有效剪應力項(effective stressterm) ， 包 含 層 流 剪 應 力 、 亂 流 剪 應 力 與 延 散 剪 應 力 (dispersion stresses)。
2.1.2 輔助關係式
(1)底床剪應力
底床剪應力採用Rastogi and Rodi (1978)之經驗式 1 2 2 1/ 2 ( ) b Cf u u v τ = ρ + (2.7) 2 2 2 1/ 2 ( ) b Cf v u v τ = ρ + (2.8) 式中， 2 c g Cf = = 摩擦係數；c = Chezy 係數。 (2)層流與亂流剪應力 採用Boussinesq之渦流黏性理論，層流與亂流剪應力可合併表 示為 2 11 1 1 1 2 1 ' 2 u v h u h h h τ υ ρ ξ η ⎡ ∂ ∂ ⎤ − = ⎢ + ⎥ ∂ ∂ ⎣ ⎦ (2.9) 2 22 2 2 1 2 1 ' 2 v u h v h h h τ υ ρ η ξ ⎡ ∂ ∂ ⎤ − = ⎢ + ⎥ ∂ ∂ ⎣ ⎦ (2.10) 12 2 1 1 2 2 1 ' ' 2 h v h u u v h h h h τ υ ρ ξ η ⎡ ∂ ⎛ ⎞ ∂ ⎛ ⎞⎤ − = ⎢ ⎜ ⎟+ ⎜ ⎟⎥ ∂ ⎝ ⎠ ∂ ⎝ ⎠ ⎣ ⎦ (2.11) 式中，υ υ υ= +l t；υl = 層流黏滯係數；υt = 亂流黏滯係數 = ku d* / 6 (Falcon 1979)； 1/ 2 * ( b/ ) u = τ ρ = 剪力速度；k = von Karman’s 係數(約 等於0.4)。 (3)延散剪應力
由於水流在進入彎道後，流場隨之而彎曲，致使流線因彎曲而 產生徑向慣性力，水面因而形成超高以產生徑向靜水壓差，得以與 徑向慣性力取得平衡。在這兩種力之作用下，水流除了以縱向方向 流動外，在徑向尚產生兩層水流，上層水流之外岸慣性力大於靜水 壓差，下層水流則反之，因此造成上層水流流動方向為朝外岸，下 層水流則為朝內岸流動，稱之二次流(secondary current)，如圖2.1所 示。 為積分水深平均所產生之延散剪應力項，須對流速剖面作一適 當假設。本模式在延散剪應力的處理方面，則是僅考量二次流的影 響，並採用 de Vriend (1977)對二次流速度剖面之假設： 1 ln m( ) g g u u uf kc kc ζ ζ ⎡ ⎤ = ⎢ + + ⎥= ⎢ ⎥ ⎣ ⎦ (2.12) 1 2 2 ( ) 2 ( ) ( ) 2(1 ) ( ) m m g g ud v vf F F f k r kc kc ζ ⎡ ζ ζ ζ ⎤ = + ⎢ + − − ⎥ ⎢ ⎥ ⎣ ⎦ (2.13) 式中， 1 1 0 ln ( ) 1 F ζ ζ dζ ζ = −
∫
； 2 1 2 0 ln ( ) 1 F ζ ζdζ ζ = −∫
；ζ =(z−zb)/d = 距離底 床之高度與水深之比值；r = 曲率半徑；。 de Vriend (1977)二次流速度剖面之適用範圍為(1)水深遠小於渠 道寬度；(2)渠道寬度遠小於渠道之曲率半徑；(3)單一二次流(single secondary eddy only)；(4)完全發展流況(developed flow)。2.1.3 邊界條件
本模式目前考量三種邊界條件設定，分別為渠道入流、渠道出 流與固體邊界。一般而言，渠道入流邊界條件設定為單位寬度入流
量，渠道出流邊界條件則採用水位高程設定。在固體邊界處，應用 側壁理論(Law of the wall)，設定靠近固體邊界的邊界條件為：
* 1 ln( ) w u Ey U k + = (2.14) 式中，uw=近固體邊界的水深平均速度；E=糙度因子=9.0 (Lien 等 1999b)；y y Uw */υ + = ， w y =固體邊界與鄰近固體邊界格點的距離。 2.2 沉滓運移部分 2.2.1 控制方程式 輸砂通量可分離為非均勻質之懸浮載與底床載，因此，沉滓運 移部分控制方程式包含某一粒徑之懸浮載質量守恆方程式、某一粒 徑於作用層(active layer)之質量守恆方程式及整體河床沉滓運移之質 量守恆方程式。 (1)某一粒徑懸浮載之質量守恆方程式 1 2 C u C v C t h ξ h η ∂ ∂ ∂ + + ∂ ∂ ∂ 1 1 1 2 2 1 2 2 1 2 2 1 1 ( h C) ( e )h C S h h d ξ ε h ξ h h d η ε h η ρd ⎡ ⎤ ∂ ∂ ∂ ∂ = + ⎢ + ⎥+ ∂ ∂ ∂ ⎢⎣ ∂ ⎥⎦ (2.15) (2)某一粒徑於作用層之質量守恆方程式 1 2 1 2 2 1 ( ) (1 ) m ( ) ( ) 0 s b b f E p h h h q h q S S t ∂ β ∂ ∂ ρ ∂ ∂ξ ∂η − + + + − = (2.16) (3)整體河床沉滓運移之質量守恆方程式 0 ) ( ) ( ) 1 ( 2 1 1 2 2 1 ⎥= ⎦ ⎤ ⎢ ⎣ ⎡ + ∂ ∂ + ∂ ∂ + ∂ ∂ −
∑
h q hq S t z h h p b b b s ξ η ρ (2.17) 以上諸式中，C = 濃度；ε 、ε = ξ、η方向之亂流傳輸係數；e =環流傳輸係數；ρs = 泥砂密度；β = 粒徑百分比； p = 孔隙率； m
E = 作用層厚度；qbi = i 方向某一粒徑之底床載通量；S = 懸浮載
源；Sf = 作用層源(active-layer floor source)。
2.2.2 輔助關係式 (1)亂流傳輸 亂流傳輸係數採用Elder (1979)之經驗公式，可表示為 1 5.93u d* ε = , ε2 =0.23u d* (2.18) (2)環流傳輸 η方向之環流傳輸係數採用Fischer等(1979)之經驗公式，可 表示為 2 2 * 25 ud e u r ⎛ ⎞ = ⎜ ⎟ ⎝ ⎠ (2.19) Fischer等(1979)經驗式之適用範圍為水深遠小於渠道寬度之 流況。 (3)底床載通量 qb
非凝聚性沉滓之底床載通量則採用Van Rijn (1984a)之經驗式
(以ξ方向為例)： 1 1 2.1 0.3 * ( ) 0.053 ( 1) k k b b k s k k T q q D s gD D D ρ = = − (2.20) 上式中， 1 3 * 2 ( 1) k k s g D D ν − ⎡ ⎤ = ⎢ ⎥ ⎣ ⎦ = 無因次顆粒粒徑；Dk = 顆粒 k 之粒徑；
2 2 * * 2 * ( ) ( ) c k k c k u u T u − = = 輸 送 參 數 ； u*c = 1 u g c = 臨 界 剪 力 速 度 ； 1 90 12 18 log( ) 3 d c D = = 顆粒蔡司係數；s ρs ρ = =砂比重。 假設底床載運移僅發生在作用層內，並考慮較小粒徑在水體中 會形成懸浮載及較細顆粒可能被隱藏在較粗顆粒之間較不易被水流 帶動之機制，則某一粒徑之底床載通量可進一步修正為： ( ) i k i b h k b k q =ζ β q D (2.21)
式中，ζh = 隱藏因子(hiding factor) (Karim等 1987)。
由於水流在渠道中運行時，底床載源之變化將受到(1)縱向及橫 向之底床坡降S0ξ、S0η；(2)縱向及橫向之流速u、v的影響；若為彎 道處，將再受(3)彎道曲率半徑r之影響；本研究採用Struiksma (1985) 所引用之公式，因此某一粒徑之底床載通量可再進一步修正為： ξ方向： ( ) cos k b h k b k q q D ξ =ζ β α (2.22) η方向： ( ) sin k b h k b k q q D η =ζ β α (2.23) 上兩式中， 1 1 sin tan 1 cos b s b s z f z f δ θ η α δ θ ξ − ∂ ⎛ − ⎞ ⎜ ∂ ⎟ ⎜ ⎟ = ∂ ⎜ − ⎟ ⎜ ∂ ⎟ ⎝ ⎠ (2.24) 1 1 tan v tan Ad u r δ = − ⎛ ⎞− − ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ (2.25) 2 2 1 g A k kc ⎛ ⎞ = ⎜⎜ − ⎟⎟ ⎝ ⎠ (2.26)
其中，α = 沉滓運移角度；δ = 底床剪應力之方向； fs = 沉滓之形 狀因子，1 2 s f ≤ ≤ ； θ = 50 2 2 2 ) 1 ( ) ( D s c v u − + = Shields paramerter； A = 渦流係 數(Jansen 1979)。 (4)懸浮載源 S 本研究著重於探討非凝聚性沉滓對彎道動床之影響，其中懸浮 載源係屬於河床質載之一部份，是由懸浮質向下之通量與底床亂流 剪力作用產生河床質向上之通量交互作用之結果。懸浮質下移到河 床表面，主要是受到重力的影響。對某一粒徑 k 之懸浮質而言，其 向下之通量可表為： k k d fk d S = −ρw C (2.27) 式中， * [3.25 0.55ln( )] k fk d k w C C u κ = + (Lin 1984)；Ck：顆粒 k 之深度平均 濃度； fk w = 顆粒 k 之沉降速度。 其中wfk則在不同之粒徑下採用不同之經驗式(Van Rijn,1984b)， 如下三式所示： v gD s w k fk ) 1 ( 18 1 − = ；{Dk <100μm} (2.28) ⎪⎭ ⎪ ⎬ ⎫ ⎪⎩ ⎪ ⎨ ⎧ − ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + − =10 1 0.01( 1) 1 5 . 0 2 3 v gD s D v w k k fk ；{100μm≤Dk ≤1000μm} (2.29)
(
)
[
]
0.5 1 1 . 1 k fk s gD w = − ；{Dk ≥1000μm} (2.30) 另一方面，底床載成為懸浮質，主要受到底床之亂流作用所造 成。對某一粒徑 k 而言，底床載向上之通量可表為：k k e lk k e S =ρ βw C (2.31) 式中， 3 . 0 * 5 . 1 015 . 0 k k D T a D C k k e = (Van Rijn 1984b)；a：沙丘高度之一半； lk w = 顆粒 k 之躍起速度。 顆粒躍起速度定義為河床質發生跳躍(saltation)時，離開底床之 瞬間垂直速度。本研究採用Hu and Hui (1996)提出之經驗公式
* 3.2 4.5 log 12 3.1 1.2 lk w u − Θ Θ < ⎧ = ⎨ Θ > ⎩ (2.32) 式中， ( ) b s gDk τ ρ ρ Θ = − 。 故由(2.27)及(2.31)式知，對某一粒徑 k 之懸浮載源可表為： ( ) k k k lk k e fk d S =ρ w β C −w C (2.33) (5)作用層源 Sf 作用層示意圖如圖2.2所示。作用層源係由於母層(active stratum) 頂面之升降而產生，當其下降時， (1 ) [( ) ( )] f s s k b m S p z E t ∂ ρ β ∂ = − − − (2.34) 式中，(βs k) = 母層中某一顆粒 k 之粒徑百分比。 如母層之厚度增加，及其頂面上升時，(2.34)式中之(βs k) 則改為 k β 。 (6)作用層厚度 Em
以下式表示： 1 ( n n) m em b b E = −C z + −z (2.35) 式中，Cem為數值參數(本模式暫取為 20)。 當河床表面接近護甲條件時，作用層厚度接近零，在這種情況 下，可用 Borah 等 (1982) 所提出護甲層之厚度，予以修正： 1 1 ( ) 1 n n m m b b K k m k D E C z z p β + = = − − + −
∑
(2.36) 式中，Dm：不產生移動的最小顆粒粒徑。 另外，作用層在淤積期間可定義為： 1 1 ( ) n n n n m m b b E + =E + z + −z (2.37) 2.2.3 邊界條件 沉滓運移求解的未知變數為深度平均懸浮載濃度 C、粒徑百分 比 β及河床高程 zb。 C 邊界條件的處理，將渠道入流邊界條件設定為入流濃度分 佈，渠道出流處邊界條件為∂ ∂ =C ξ 0，在固體邊界處，邊界條件設 定為∂ ∂ =C η 0。β與zb部分，入流邊界條件利用已知值設定，下游 邊界處之底床坡度通常較為平緩，其邊界條件可設定為∂ ∂ =β ξ 0與 0 b z ξ ∂ ∂ = ， 在 固 體 邊 界 處 ， 其 邊 界 條 件 則 設 定 為 ∂ ∂ =β η 0 與 0 b z η ∂ ∂ = 。圖 2.1 彎道二次流示意圖
圖 2.2 作用層示意圖 (摘自 Spasojevic 1988)
第三章 數值架構
3.1 水理部份 3.1.1 隱式雙階分割操作趨近法 本研究基於分割操作之觀念，將動量方程式分割成二個步驟(延散步 驟及傳播步驟)，並利用隱式數值方法求解。延散步驟求解移流項和擴散 項，傳播步驟求解壓力項、底床剪應力項和連續方程式。據此，水理控制 方程可改寫成： 延散步驟 1 1 1 2 2 1 2 ( ) n n n n V V V T t ρ + + + ∂ ⎛ ⎞ = − ⋅∇ + ∇ ⋅ ⎜ ∂ ⎟ ⎝ ⎠ (3.1) 傳播步驟 1 1 2 1 ( ) n n n b b V V g z d t t d τ ρ + + + ∂ ∂ ⎛ ⎞ −⎛ ⎞ = − ∇ + − ⎜ ∂ ⎟ ⎜ ∂ ⎟ ⎝ ⎠ ⎝ ⎠ (3.2) 1 0 n V + ∇ ⋅ = (3.3) 式中，V 表示速度向量；T表示擴散及延散項；n+1表示(n+ Δ1) t時刻之未 知變數； n 1 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 h T 1 h T 1 h 1 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 ( ) b s b s b z z z z h h h h h h d τ τ τ τ ρ ξ ξ η η ⎡ ∂ ∂ ∂ ∂ ⎤ + ⎢− + − + ⎥ ∂ ∂ ∂ ∂ ⎣ ⎦ (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 h T 1 h T 1 h 1 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 ( ) b s b s b z z z z h h h h h h d τ τ τ τ ρ ξ ξ η η ⎡ ∂ ∂ ∂ ∂ ⎤ + ⎢− + − + ⎥ ∂ ∂ ∂ ∂ ⎣ ⎦ (3.5) 傳播步驟 2 2 1 (zb d) C u uf v u g t h ξ d + ⎛∂ + ⎞ ∂ = − − ⎜ ⎟ ∂ ⎝ ∂ ⎠ (3.6) 2 2 2 (zb d) C v uf v v g t h η d + ⎛∂ + ⎞ ∂ = − − ⎜ ⎟ ∂ ⎝ ∂ ⎠ (3.7) 和 2 1 1 2 ( ) ( ) 0 h ud h vd d h h t ξ η ∂ ∂ ∂ + + = ∂ ∂ ∂ (3.8) 針對n+1時刻的水深值( n 1 d + )做線性化處理，且僅保留一階項，(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τ α = − Δ ； 1 2 1 2 2 1 1 [ ] n n n zb h h g t d u Cτ C hτ ∂ ∂ β ∂ξ ∂ξ + + Δ = − + ; 1 1 n d γ =β ； 1 2 2 n h g t d C hτ α = − Δ ; 12 1 1 1 2 2 [ ] n n n zb h h g t d v Cτ C hτ ∂ ∂ β ∂η ∂η + + Δ = − + ; 2 2 n d γ =β ； 1 1 2 2 2 2 ( ) ( ) 1 n n f n C u v C t d τ + + + = + Δ ； n 1 n d d + d Δ = − 。 3.1.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 w p ξ ξ + + + ⎛∂Ψ⎞ Ψ − Ψ = ⎜∂ ⎟ Δ ⎝ ⎠ (3.10) 1 1 1 n n n n s p η η + + + ⎛∂Ψ⎞ =Ψ − Ψ ⎜∂ ⎟ Δ ⎝ ⎠ (3.11) 式中， 1 1 1 1 1 1, , 0.5 ( ) 0.5 ( ) n n n n n e E P i j i j + + + + + + Ψ = ⋅ Ψ + Ψ = ⋅ Ψ + Ψ ； 1 1 1 1 1 , 1, 0.5 ( ) 0.5 ( ) n n n n n w P W i j i j + + + + + − Ψ = ⋅ Ψ + Ψ = ⋅ Ψ + Ψ ； 1 1 1 1 1 , 1 , 0.5 ( ) 0.5 ( ) n n n n n n N P i j i j + + + + + + Ψ = ⋅ Ψ + Ψ = ⋅ Ψ + Ψ ； 1 1 1 1 1 , , 1 0.5 ( ) 0.5 ( ) n n n n n s P S i j i j + + + + + − Ψ = ⋅ Ψ + Ψ = ⋅ Ψ + Ψ ； Ψ可表為 u, v, h1, h2, d, zs 和 zb。 混合型上風法為上風法(upwind scheme)與中央差分法組合而成，當移 流效應重要時，採用上風法；移流效應不重要時，則採用中央差分法。至 於 移 流 效 應 重 要 性 的 判 斷 ， 則 採 用 格 網 雷 諾 數 (mesh Reynolds number)Rx、Ry作為判斷的因子，當 Rx 或 Ry 大於2時，代表移流效應重要，
差分方法採用能反映方向性的上風法； x R 或 Ry 小於等於2時，移流效應 可視為不重要，差分方法採用中央差分法。 混合型上風法應用於本研究移流項的處理可表成 , 1 1 1 1 1 , 1, , , 1, 1 1 0.5 (1 ) (1 ) i j n n n n n n n i j i j i j i j i j x x u u h ξ h α ξ α ξ + + + + + + − ⎡ ⎛Φ − Φ ⎞ ⎛Φ − Φ ⎞⎤ ⎛∂Φ ⎞= − + + ⎢ ⎜ ⎟ ⎜ ⎟⎥ ⎜ ∂ ⎟ ⎜ Δ ⎟ ⎜ Δ ⎟ ⎢ ⎥ ⎝ ⎠ ⎣ ⎝ ⎠ ⎝ ⎠⎦ (3.12) , 1 1 1 1 1 , , 1 , , , 1 2 2 0.5 (1 ) (1 ) i j n n n n n n n i j i j i j i j i j y y v v h η h α η α η + + + + + + − ⎡ ⎛Φ − Φ ⎞ ⎛Φ − Φ ⎞⎤ ⎛∂Φ ⎞ = ⎢ − ⎜ ⎟+ + ⎜ ⎟⎥ ⎜ ∂ ⎟ ⎜ Δ ⎟ ⎜ Δ ⎟ ⎢ ⎥ ⎝ ⎠ ⎣ ⎝ ⎠ ⎝ ⎠⎦ (3.13) 其中 0 2 1 2 1 2 x x x x R R R α ⎧ ≤ ⎪ ⎪ =⎨ > ⎪ ⎪− < − ⎩ ; 0 2 1 2 1 2 y y y y R R R α ⎧ ≤ ⎪ ⎪ =⎨ > ⎪ ⎪− < − ⎩ (3.14) 上列諸式中， , 1 , n i j i j x u h R ξ μ ρ Δ = ； , 2 , n i j i j y v h R η μ ρ Δ = ；μ = 流體動力黏滯係數 (dynamic viscosity)；Φ 可表成 u 或 v 。 3.2 沉滓運移部分 3.2.1 耦合演算法 假定底床是由 k 種粒徑所組成，則與沉滓運移有關之控制方程式包 含 k 個深度平均懸浮載之質量守恆方程式((2.15)式)、k 個作用層質量之 守恆方程式((2.16)式)及一個整體河床沉滓運移之質量守恆方程式((2.17) 式)。求解的未知變數為 k 個深度平均懸浮載濃度 C、k 個粒徑百分比 β 及河床高程 zb。
由於各物理量具有高度之相關性，如河床質與懸浮質藉著懸浮載源 S 連結，河床高程 zb 及作用層內之粒徑百分組成 β 又因懸浮載源 S 而變動等，故有賴耦合演算方式以求解此三條方程式。(2.16)及(2.17)式為 雙曲線偏微分方程式，但(2.15)式為雙曲線-橢圓形混合之偏微分方程式， 因此聯立解此三方程式時，須先對(2.15)式作一適當處理。 (2.15)式為延散方程式，解此方程如同解水流運動方程式，可分成下 列二個步驟： (1)移流及反應(advection-reaction)步驟 1 2 C u C v C S t h h d ∂ ∂ ∂ ∂ + ∂ξ + ∂η = ρ (3.15) 式中 C t ∂ ∂ 為移流及懸浮載源 S 所引起的濃度變化，可以( ) a C t ∂ ∂ 表示 之。 (2)擴散(diffusion)步驟 1 1 1 2 2 1 2 2 1 1 ( C)d ( C)a ( h C) ( e )h C t t h h d ξ h h h d η η h ∂ ∂ ε ε ∂ ∂ η η η η ⎡ ⎤ ∂ ∂ ∂ ∂ − = + ⎢ + ⎥ ∂ ∂ ∂ ⎢⎣ ∂ ⎥⎦ (3.16) 經上述方法分解後之懸浮載方程式，在移流及反應步驟為一雙曲線偏 微分方程，而於擴散步驟則為一橢圓型偏微分方程。因此，先聯立求解 (2.16)、(2.17)及(3.15)式；然後，所求得各變數之中間值再與(3.16)式反覆 疊代至收斂為止。 3.2.2 數值差分式
採用與水理部分相同的差分概念，以圖 3.1 中之主格點 P 為例，(2.16) 及(2.17)式之離散化方程式可分別表為： 1 1 1 1 1 2 2 1 2 (1 ) 0.5 [( ) ( ) ] [ ( ) ( ) ] e p n n n n s m p m p b e w b w p p E E h q h q t h h ρ − β + − β + + − + Δ 1 1 2 2 1 1 2 2 1 1 1 2 1 2 0.5 0.5 [ ( ) ( ) ] [ ( ) ( ) ] e w n s p p p p n n n n b e b w b n b s h q h q h q h q h h h h + + + − + − 2 2 1 1 2 1 2 0.5 [ ( ) ( ) ] 0.5 0.5 ( ) 0 n w p p n n n n b n b s p p f p h q h q S S S h h + + − + + − = (3.17) 1 1 1 1 2 2 1 1 2 (1 ) 1 [( ) ( ) ] {0.5[ ( ) ( ) ] e p n n n n s b p b p b e w b w p p z z h q h q t h h ρ − + − + + − + Δ
∑
1 1 2 2 1 1 2 2 1 1 0.5[ ( ) ( ) ] 0.5[ ( ) ( ) ] e w n s n n n n b e b w b n b s h q h q h q + h q + + − + − 2 2 1 1 1 0.5[ ( ) ( ) ]} [0.5 0.5 ] 0 n s n n n n b n b s p p h q h q S + S + − +∑
+ = (3.18) 懸浮載質量守恆方程式之移流及反應步驟可採用特性法求解，(3.15) 式之等式左邊各項，可以全微分表示，即 1 2 C u C v C DC t h h Dt ∂ ∂ ∂ ∂ + ∂ξ + ∂η = (3.19) 則(3.15)式變為： DC S Dt = ρd (3.20) 在正交曲線座標上，其方向軌跡為： 1 h d u dt ξ = 和 h d2 v dt η = (3.21) 圖 3.2 為質點運動軌跡之示意圖，A 點為 1 ( 1) n t + = n+ Δt時之到達點，D 點為 n t 時之離開點。要解(3.20)式之全微分項，必須先解(3.21)式，即求出其積分路徑。由於模式採用固定格點，故離開點不能保證剛好落在格點 上，而必須以鄰近格點採內差法求之。 由於底床載方程式與懸浮載方程式在時間尺度上相差很大，即底床載 移動量遠小於懸浮質移動量。在必須聯立求解前提下，移流及反應步驟之 懸浮載方程式((3.19)式)必須使用較大之可蘭數(Courant number)，但大可 蘭數將導至軌跡線必須穿過若干個格網方能落在 n t 平面上，為求出質點 之移動軌跡，路徑求解採分段方式處理，如圖 3.3 所示。 (1)運移軌跡之計算 對某一粒徑，如以從 Dl 積至 Al 之軌跡為例，則： ( ) ( 1) ( 1) ( ) 1 ( ) ( ) 2 m m A D m D A m A D M u u t t h ξ + =ξ − + − + (3.22) ( ) ( 1) ( 1) ( ) 2 ( ) ( ) 2 M m m A D m D A m A D v v t t h η + =η − + − + (3.23) ( ) ( ) ( ) ( ) 1 2 ( 1) ( ) ( ) 1 1 2 2 ( ) ( ) min[ ; ] ( ) ( ) M M m m m m D A D A m D A m m A D A D h h t t u u v v ξ ξ η η + = − − − + + (3.24) 式中，下標A、D表示在到達和離去點的軌跡，M 為A、D之中點值；上 標m表疊代次數。此疊代過程係為求得離去點Dl之座標值。 (2)移流及反應步驟之離散方程式 當運動軌跡已知後，可積分(3.20)式，以圖 3.3 之軌跡為例，其離散 化之方程式為： 1 1 1 1 1 1 2 [( ) ] ( ) 2 2 l l l l L l l D D D D D A A D A D l L D D D A S S t t S S t t C C d d d d ρ ρ −− − ρ ρ = − − − =
∑
+ + + (3.25)另外，懸浮載質量守恆之擴散項之數值方法，同水理部分之延散步 驟，時間項採用前項差分方法，擴散項採用中央差分法。 3.2.3 數值解析法 在主格點 P(如圖 3.1 所示)上之離散化方程式共有 k 個(3.17)式，1 個 (3.18)式及 k 個(3.25)式。而在 P 點上之未知量，可以如下向量表式之： 1 1 1 1 ( , , ,..., , ,..., , ) n n b k k K K sG + = z + C β C β C β (3.26) 或更簡潔的記作 1 1 2 2 1 ( , , ) 1, n k k s + = s s s + k = K G (3.27) 式中 k 為粒徑 k 之代號。則(3.17)、(3.18)及(3.25)式可分別寫成： 1 1( ) 0 n F sG + = (3.28) 1 2 ( ) 0 1, n k F sG + = k= K (3.29) 1 2 1( ) 0 1, n k F + sG + = k= K (3.30) (3.28) 至 (3.30) 式 為 一 非 線 性 代 數 式 ， 可 加 以 線 性 化 後 ， 利 用 Newton-Raphson 法疊代求解： 1 1 1 [ F ] s F(msn ) s ∂ ∂ + Δ = −G G G (3.31) 1 2 2 [ k] (m n ) 1, k F s F s k K s ∂ ∂ + Δ = −G G = G (3.32) 1 2 1 2 1 [ k ] (m n ) 1, k F s F s k K s ∂ ∂ + + + Δ = −G G = G (3.33)
式中，[ F] s ∂ ∂G 為 Jacobian 係數矩陣中之一列向量； 1 m n sG + 為前一次疊代未知 向量；ΔsG為疊代修正向量，可表為Δ = Δ ΔsG ( s1, s2k,Δs2k+1) 。解得修正向量ΔsG 後，可得新的m 1 n 1 s + G + 向量： 1 1 1 m n m n s s s + G + = G + + ΔG (3.34) 當達到Δ ≤sG εs(εs為收斂容許誤差值)之收斂條件時，疊代得以結束。
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)計算區域
圖 3.2 質點運移軌跡示意圖 (摘自 Spasojevic 1988)
圖 3.3 質點運移軌跡在計算平面之示意圖 (摘自 Spasojevic 1988)
第四章 模式適用性及二次流效應之敏感度測試分析
經由第二章與第三章在模式理論基礎及數值方法之說明介紹 後，將於本章節進行模式在沉滓運移部份之理論基礎的適用性分析， 並進行模式各項二次流效應之敏感度分析，以期檢定各項影響因子 在彎道動床模式之影響程度及提供模式參數修正之方向。 4.1 模式適用性測試分析及問題解析根據 Yeh and Kennedy(1993)所描述之一般在具有動床特性之彎 曲渠道會有四項特徵：(1)有徑向沉滓運移之情形；(2)傾斜之底床形 狀；(3)偏斜之主流位置，一般位於彎道外岸；(4)具有二次流之流場。 因此模式演算以能達到理論應展現之彎道動床特性為目的。
測試案例係採用 Struiksma(1983)於荷蘭臺夫特科技大學水利實 驗室(Delft Hydraulic Laboratory) 所作之室內水槽試驗組別 T2，此後 簡稱 DHL-T2 實驗。 DHL-T2 實驗之水槽佈置如圖 4.1 所示，從斷面 AD至斷面 BD及 斷面 CD至斷面 DD為直線道，長度均為 15m，斷面 BD為彎道起點， 從此處沿著渠道至斷面 CD處為一蜿蜒 140°之固定曲率半徑之動床矩 形渠道，渠道寬度為 1.5m，渠道中心線之曲率半徑為 12m，彎道全 長 29.32m，全渠道之縱向底床坡降為 0.00203，Chezy 係數為 28.8 m1/2/s，底床作用層採用粒徑為 0.45mm 之均勻砂。上游邊界條件為固 定入流量 0.062 m3 /s，下游邊界條件則為實測水位 0.1m。採用 121×37 之非均勻計算格網，其中側方向格網分佈情形如圖 4.2 所示，此案例 之總模擬時間為 400 min，程式演算
Δ
t
=1.0 sec。 由於以往之學者在數值模式中，針對彎道底床載二次流所應用之參數均為固定值，並不因不同斷面或不同渠道位置而有所差異，其中
包括沉滓形狀因子 fs與渦流係數A，故本研究於初步模擬時，仍是應
用固定值之參數，採用 Kassem and Chaudhry(2002)之設定值，其中 fs
為 2，A則為 9.10。 圖 4.3 為模擬 DHL 實驗 T2 之彎道段各橫斷面無因次水面高程示 意圖，分別取0°(彎道入口)、23.3°、46.6°、70°、93.3°、116.6°、140°(彎 道出口)、彎道出口 5m 和彎道出口 10m 等九個斷面作說明，其中縱 軸為該斷面各點實際水位(Zs)與平均水位(Zsm)之比較值，由圖可明確 觀察到在彎道段，外岸有明顯水面超高之情形；水流流出彎道後，水 面超高情形隨之減小，且徑向慣性力和徑向靜水壓之差異逐漸消失， 而二次流效應亦隨之衰退。 圖 4.4 和圖 4.5 分別為模擬 DHL-T2 實驗之縱向底床變化圖及縱 向水深示意圖，此二圖中實線段為實驗值，方塊點狀線為本研究所模 擬之結果。另外，圖 4.5 虛線段為 Struiksma(1985)所模擬之結果。由 圖 4.4 中可清楚看到無論是實驗結果或是本研究模擬結果，皆有外岸 沖刷、內岸淤積之定性上的共通點；由圖 4.5 中可明顯看出彎道內岸 水深較淺，彎道外岸水深值則明顯高於內岸水深值。 惟在定量上觀察圖 4.4 和圖 4.5，無論本模式之模擬結果或 Struiksma(1985)模擬結果皆與實驗結果有一定之差距，尤其是以彎道 入口段之水深或底床高程特別明顯。由此二圖中實驗段之量測值來觀 察，在彎道入口段的外岸區域有沖刷較深和內岸區域有淤積較高之情 形，稱之為「過度沖淤」(overshoot)現象，數值模擬結果與實驗結果 明顯有所不符。Struiksma(1985)評斷是為影響彎道動床模式之因素相 當複雜，可能會受到模式之參數或相關演算條件之設定而致使無法反
應實際現象之情形，並針對此結論予以分析探討，提出三項可能影響 結果之分析：(1)上游給定之邊界條件，可能影響沖淤深度；(2)求解 水理和輸砂之過程中，可能對於其運動形式之定義有所不足，於此提 高其渦度源項之梯度，亦即增加整體水流徑向流速之坡降，彎道外岸 流速加大，內岸流速減低；(3)模式中之各項參數可能針對不同之案 例會有不同之值，然而這大多為經驗公式，含有相當高之不確定性， 尤其是在沉滓運移角度修正因子中之 fs和A值，須加以探討其值。 4.2 二次流效應之敏感度測試分析 影響彎道動床之因素可概括為水理二次流、懸浮載二次流及底床 載二次流效應，另外 Lien 等(1999)發現在彎道處之水理滑移邊界效 應，影響了水理二次流對彎道流速之分佈情形，而渠道之沖淤現象往 往與流場息息相關，故本研究在探討水理二次流時，並同時探討滑移 邊界效應對彎道底床沖淤之影響。由 4.1 節可得知，在已引用之理論 下，數值模式並無法直接模擬出符合實際彎道之床形，因此需透過測 試上述因素對底床變動之敏感度，予以探討並修正模式演算之各項設 定值，以期模擬出符合實驗結果之彎道床形。以下章節將針對各項彎 道二次流效應進行分析，逐一探討其對彎道底床沖淤之敏感度。
為展示各影響因子之敏感度，選用均方根差(root mean square
error)Erms來作為評斷該因子之敏感程度，此指數敘述如下：
(
)
2 1/2 , , ⎟⎟ ⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎜ ⎝ ⎛ − =∑
total N j i j i rms N Y X E total (4.1) 其中，Xi,j表某計算格點於考慮指定模擬條件時之模擬結果；Yi,j表Xi j, 相對應計算格點於未考慮指定模擬條件時之模擬結果；N 表模擬案例總計算格點數目。 另外，為適切量化各影響因子對彎道底床變動程度，乃定義一指 標Emzp來代表，該式如下所示： 0 b bip bnp bip mzp z z z z E − − = (4.2) 式中，zbip = 某計算格點於考慮指定模擬條件時之底床高程；zbnp = 某計算格點於未考慮指定模擬條件時之底床高程；zb0 = 原始底床高 程。當Emzp之值愈趨近於 0 時，表示該模擬條件對底床變動之影響甚 小，亦即反應其沖淤現象較微弱，當Emzp之值愈大，表示此模擬條件 對底床變動之影響較大，亦即反應其沖淤現象較為強烈。 本節將以Erms及Emzp等兩指標來評斷各項彎道動床影響因子之敏 感程度，並作為模式參數修正之依據。 4.2.1 水理二次流效應分析 根據謝(2002)對彎道水理二次流所作的分析中，水深平均二維模 式可分為兩種型態，一種為傳統模式(Conventional Model；CN 模式)， 一種為彎道模式(Bend-Flow Model；BF 模式)，這兩種模式最主要之 差異在於動力方程式中對延散剪應力的處理。CN 模式忽略彎道二次 流之效應，即假設垂直方向之流速為均勻分佈，因此忽略延散剪應力 項；BF 模式則為引進深度方向二次流流速剖面來處理延散剪應力 項，以適切反映彎道二次流效應之影響。透過謝(2002)之分析，得知 當二次流強度因子(SI)小於 0.02 時，二次流效應對流場之影響可忽略 不計，因此本研究採用SI值大於 0.02 之案例作為探討對象，以完整 了解水理二次流效應對沉滓運移之影響程度，並對其敏感度加以分析
探討。二次流強度因子計算式如下所示： f c C r H SI = (4.3) 此部份分析採 DHL-T2 實驗作為模擬對象，該實驗之SI 值為 0.0766，此SI值比對於圖 4.6 中屬較小之值，其中圖 4.6(謝，2002)為 不同θb之 * MaxU 與SI關係圖，θb為彎道長度因子=L/(2πrC)，MaxU*為ξ 方向速度之最大相對差異= ( / ) b u u Max Δ ，此案例之SI值大於忽略水理 二次流效應之SI值 0.02，故模擬此案例理應考慮水理二次流效應；表 4.3(謝，2002)為不同θb之 * MaxU 與SI之回歸函數一覽表，以彎道中心斷 面70°為例，其θb等於 0.1944，利用表 4.1 中之回歸函數可得其MaxU*等 於 0.531，代表水理二次流影響ξ方向之流速達 50%以上，因此必須考 量水理二次流效應。 圖 4.7 為模擬 DHL-T2 實驗時，水理二次流效應對縱向水深影響 示意圖，並參照表 4.2 可查得模式在計算水理二次流與否之水深Erms值 並不高，在距離彎道內岸 0.375m 之縱斷面為 0.001674，距離彎道外 岸 0.375m 之 縱 斷 面 則 為 0.001361 ， 若 以 全 彎 道 網 格 來 看 則 為 0.001484；若以Emzp之值來觀察，由表 4.3 可得知模擬此案例時，模式 針 對 計 算 水 理 二 次 流 與 否 影 響 整 體 彎 道 之 平 均 底 床 變 動 量 約 10.2%，在距離彎道內岸 0.375m 之縱斷面為 12.95%，距離彎道外岸 0.375m 之縱斷面則為 7.35%，因此可發現水理二次流在此案例對底床 沖淤仍有所影響。 但由於本案例之SI值實屬較小之值，因此本研究為探討當模擬案 例在較大SI值的時候，水理二次流效應對底床沉滓運移之影響程度， 以 DHL 實驗同一渠道形狀之案例，採用同一粒徑砂( = mm)，
設計較大SI之流況，而設計案例與原 DHL -T2 實驗差異之相關資訊 為：(1)設計流量=0.1163 cms；(2)平均水深=0.15 m；(3)平均流速=0.51 m/s；(4)二次流強度因子SI=0.1149。 圖 4.8 為當SI值為 0.1149 時，水理二次流效應對縱向水深影響示 意圖，由圖 4.7 與圖 4.8 可以明顯看出模式在開啟和未開啟水理二次 流的條件下，在較大SI值的時候(約為 DHL-T2 實驗之 1.5 倍)，其水 深之差異明顯大於小SI值時之差異，若以水深之Erms值來觀察，距離 彎道內岸 0.375m 之縱斷面為 0.007097，距離彎道外岸 0.375m 之縱斷 面則為 0.004544，若以全彎道網格來看則為 0.005914；若以Emzp之值 來觀察，模式在此案針對計算水理二次流與否影響整體彎道之平均底 床變動量約 12.26%，在距離彎道內岸 0.375m 之縱斷面為 9.0%，距 離彎道外岸 0.375m 之縱斷面則為 15.0%，整體高於模擬 DHL-T2 實 驗之水深Emzp值；由此可得知SI值愈大，水理二次流效應對底床載沖 淤之影響則較為明顯。 4.2.2 滑移邊界效應分析 在模擬天然渠道或實驗渠道時，理論上邊界之流速和流量應設定 為 0，亦稱之為不滑移邊界(non-slip boundary)，但程式在處理邊界處 常會因為鄰近格網之流速不連續而造成邊界數值演算之不穩定性，因 此偶爾會以滑移邊界(slip boundary)來模擬案例。 圖 4.9 為模擬 DHL-T2 實驗考慮滑移邊界與否之縱向水深比較 圖，縱軸為模擬水深，橫軸為彎道內之縱向距離，虛線為設定不滑移 邊界之水深，實線為設定滑移邊界之水深，由圖中可觀察得不滑移邊 界和滑移邊界所模擬出之水深差異極小，由表 4.2 可見得其水深Erms值 在內外岸 0.375m 處分別為 0.000337 和 0.000220，若考慮全彎道網格
