顯式有限解析法模式於動床主支流之研究

(1)

顯式有限解析法模式於動床主支流之研究

學生：郭峰志指導教授：葉克家國立交通大學土木工程學系

摘要

一般情況下，河川水系常有主支流交匯之情形。因此對於交匯後河川流量、水深、底床之變化，應予以探討，以供河川整治及工程規劃之需。本研究係延續Hsu and Yeh’s (1996)一維有限解析法動床模式 (Explicit Finite Analytic method, EFA)，發展擴充至具有模擬主支流渠道之功能。在交匯區內部邊界條件採用水位相等之概念來處理；在上、下游邊界處則根據流況給予流量或水位歷線。

由於動床主支流實驗案例甚為缺乏，故本模式採用Suryanarayana (1969)之均勻與非均勻沉滓動床實驗資料，假設另一條相同渠道作為支流匯入，進行含主支流渠道之沖刷與淤積模擬分析，以進行模式之測試。在天然河川水系方面，收集大漢溪及三峽河之斷面資料，以大漢溪為主流、三峽河為支流進行颱洪事件模擬，並與新海橋水位站之實測水位值比較，以檢驗本模式於天然河川之適用性。經由數模比較，顯示成果尚屬合理。

(2)

Study On Mobile-Bed Confluence Flows Using Explicit Finite Analytic Model

Student:：Feng-Chih Kuo Advisor：Keh-Chia Yeh Institute of Civil Engineering

National Chiao Tung University

Abstract

Under normal conditions, river systems often include junction and division flows. In order to rebuild and map out engineering plans for the river systems, one would have to probe into the research of the discharges, the depths and the changes of river channels after confluence. This study extends Hsu & Yeh’s (1996) 1-D finite analytic mobile-bed model (EFA), and expands it to incorporate the function of simulating the confluence flows. Equal stage at the region of confluence is adopted into the internal boundary condition; with the discharge or stage hydrograph specified at the upstream or downstream boundary according to the given flow characteristics.

Due to the lack of experimental cases for mobile-bed confluence

flows, Suryanarayana’s (1969) experimental data of uniform and

non-uniform bed material are adopted. This model simulates the

scour and deposition of the tributary, under the assumption that the

simulated channel and the tributary are the same. In regard to the

natural river system, the section data of Tahan Creek and Sanxia

River are collected in which they are considered as the mainstream

and the tributary, respectively. To test the suitability of this model in

the natural river system, the simulation of the typhoon flood incident

is carried out and the results are compared with the measuring value

of the Hsinhai Bridge water level station. By way of this model

simulates, the simulation results are satisfactory.

(3)

誌謝

本文承蒙吾師葉克家教授於研究期間不厭其煩地細心指導與諄諄教誨，得以順利完成，在此致上最誠摯的謝意。亦感謝楊錦釧教授、

張良正教授、蔡武廷教授與許勝田教授於課業上之教導與啟發，讓學生獲益匪淺。誠摯感謝口試委員國立成功大學蔡長泰教授、國立台灣大學水工試驗所賴進松博士之細心匡正與建議，使得本論文得以更趨完整。

感謝璁哥、添哥、昇佑、哲宏、昇學、仕哲、士偉、小宇學長以及曉萍、秀容學姊在研究與生活上之指導與照顧；感謝學長姊勇哥、

東哥、阿伯、D 哥、Dr.張、大槍、老王、俊傑、绣真、秋蓉及好友璨仲、俊毅、于軒、益家、弘恩的砥礪與扶持；感謝學弟妹力瑋、仲達、宣汝、欣瑜、浩榮、偉國、雅婷、宥達的幫助與關心。因為有了大家的參與，使我的研究生涯多采多姿，永生難忘。

最後，謹將此拙作獻給我親愛的父母、兩位妹妹及陪伴我的怡婷，感謝你們不斷地支持與鼓勵，使我求學的路上能無後顧之憂，因為有了你們，我才能成功跨出這一步，謝謝你們。

(4)

目錄

中文摘要... I 英文摘要...II 致謝...III 目錄... IV 表目錄... VII 圖目錄...VIII

第一章緒論...1

1.1 前言...1

1.2 文獻回顧...3

1.2.1 定床交匯渠流...3

1.2.2 動床及動床交匯渠流...6

1.3 論文架構...11

第二章理論基礎 ...13

2.1 模式概述...13

2.2 水流控制方程式...14

2.3 輸砂控制方程式...15

2.4 無黏性沉滓輸砂輔助關係式 ...17

2.4.1 河床載通量(bed load flux, Qbk)...17

(5)

2.4.2 作用層厚度 (

E

_m) ...18

2.4.3 懸浮載源 (Sk )...19

2.4.4 作用層源 (Sak) ...20

2.5 黏性沉滓輸砂輔助關係式 ...20

2.6 顯式有限解析法之簡介 ...22

第三章數值方法 ...24

3.1 求解水理方程式之數值方法 ...24

3.2 求解輸砂方程式之數值方法 ...25

3.3 輸砂數值解析法(Newton-Raphson法)...28

3.4 水理邊界條件(外部邊界條件)...29

3.5 交匯區邊界條件(內部邊界條件)...31

3.6 穩定性分析...34

3.7 交錯格網...34

3.8 模式演算之流程...35

第四章模擬成果與分析 ...38

4.1 主支流實驗渠道之模擬 ...38

4.1.1 模式演算參數設定...39

4.1.2 率定案例主支流模擬...42

4.1.3 應用案例主支流模擬...47

(6)

4.1.4 小結...50

4.2 主支流河川案例模擬...50

4.2.1 流域概述...50

4.2.2 模擬範圍概述...51

4.2.3 模式演算所需資料...52

4.2.4 模擬結果與分析...54

4.2.5 小結...55

第五章結論與建議 ...58

5.1 結論...59

5.2 建議...60

參考文獻 ...62

(7)

表目錄

表4-1 各模擬案例之上游入流量、下游水位條件與支流之入流量..70 表4-2 大漢溪斷面樁號 ...71 表4-3 大漢溪檢討河段各斷面河床質粒徑分組重量累計百分比…..72 表4-4 大漢溪模擬河段於艾莉颱洪之上下游邊界條件(1/2) ...73 表4-4 大漢溪模擬河段於艾莉颱洪之上下游邊界條件(2/2) ...74

(8)

圖目錄

圖3-1 懸浮沈滓之移流特性軌跡 ...75

圖3-2 一般合流交匯區隔點示意圖 ...32

圖3-2 一般分流交匯區隔點示意圖 ...32

圖3-4 一維特性曲線交錯格點示意圖 ...76

圖3-5 程式流程圖 ...77

圖4-1 模型佈置圖(摘錄自Suryanarayana 1969) ...78

圖4-2 非均質粒徑分布曲線圖(摘錄自Suryanarayana 1969) ...79

圖4-3 Run21 淤積案例 (單一渠道)...80

圖4-4 Run21 懸浮質濃度沿程變化 (單一渠道)...81

圖4-5 Run21 懸浮載源沿程變化 (單一渠道)...82

圖4-6 Run21-A淤積案例 (主支流渠道) ...83

圖4-11 Run21-A 懸浮質濃度沿程變化 (主支流渠道)...88

圖4-12 Run21-A 懸浮載源沿程變化 (主支流渠道)...89 圖4-13 Run21-B淤積案例 (主支流渠道，支流流量＝0.1 主流流量)90

(9)

圖4-14 Run21-B淤積案例 (主支流渠道，支流流量＝0.1 主流流量)91 圖4-15 Run21-B淤積案例 (主支流渠道，支流流量＝0.1 主流流量)92 圖4-16 Run21-B淤積案例 (主支流渠道，支流流量＝0.1 主流流量)93 圖4-17 Run21-B淤積案例 (主支流渠道，支流流量＝0.1 主流流量)94 圖4-18 Run21-B 懸浮質濃度沿程變化 (主支流渠道，支流流量＝0.1

主流流量) ...95

圖4-19 Run21-B懸浮載源沿程變化(主支流渠道，支流流量＝0.1 主流流量) ...96

圖4-20 Run22 沖刷案例 (單一渠道)...97

圖4-23 Run22-A沖刷案例 (主支流渠道) ...100

圖4-26 Run22-A 懸浮質濃度沿程變化 (主支流渠道)...103

圖4-27 Run22-A 懸浮載源沿程變化 (主支流渠道)...104

圖4-28 Run22-B沖刷案例 (主支流渠道，匯流點延後)...105

(10)

圖4-31 Run22-B 懸浮質濃度沿程變化 (主支流渠道)...108

圖4-32 Run22-B 懸浮載源沿程變化 (主支流渠道)...109

圖4-33 Run22-C沖刷案例 (主支流渠道，支流流量＝0.2 主流流量) ...110

圖4-36 Run22-C 懸浮質濃度沿程變化 (主支流渠道，支流流量＝0.2 主流流量) ...113

圖4-37 Run22-C 懸浮載源沿程變化 (主支流渠道，支流流量＝0.2 主流流量) ...114

圖4-38 Run25 淤積案例 (單一渠道)...115

圖4-41 Run25 淤積案例 (主支流渠道)...118

圖4-46 Run25 懸浮質濃度沿程變化 (主支流渠道)...123

(11)

圖4-47 Run25 懸浮載源沿程變化(主支流渠道)...124

圖4-48 Run14 沖刷案例 (單一渠道)...125

圖4-50 Run14 懸浮載源沿程變化(單一渠道)...127

圖4-51 Run14 粒徑組成百分比與中值粒徑剖面圖(單一渠道) ...128

圖4-52 Run14 沖刷案例 (主支流渠道)...129

圖4-55 Run14 懸浮質濃度沿程變化 (主支流渠道)...132

圖4-56 Run14 懸浮載源沿程變化(主支流渠道)...133

圖4-57 Run14 粒徑組成百分比與中值粒徑剖面圖(主支流渠道) ....134

圖4-57(續) Run14 粒徑組成百分比與中值粒徑剖面圖(主支流渠道) ...135

圖4-58 模擬範圍之流域示意圖 ...136

圖4-59 艾莉颱洪上游流量歷線 ...137

圖4-60 艾莉颱洪下游水位歷線 ...137

圖4-61 艾莉颱洪新海橋水位歷線比較 ...138

圖4-62 艾莉颱洪斷面淡 52、斷面淡 51 及斷面三 01 之流量歷線.138 圖4-63 艾莉颱洪大漢溪水位及底床高程沿程之變化...139

圖4-64 艾莉颱洪大漢溪懸浮質濃度沿程變化...140

(12)

圖4-65 艾莉颱洪大漢溪懸浮載源沿程變化...140

圖4-66 艾莉颱洪大漢溪之粒徑組成百分比與中值粒徑剖面圖...141

圖4-66(續) 艾莉颱洪大漢溪之粒徑組成百分比與中值粒徑剖面圖 ...142

圖4-67 艾莉颱洪三峽溪水位及底床高程沿程之變化...143

圖4-68 艾莉颱洪三峽溪懸浮質濃度沿程變化...144

圖4-69 艾莉颱洪三峽溪懸浮載源沿程變化...144

圖4-70 艾莉颱洪三峽溪之粒徑組成百分比與中值粒徑剖面圖...145

圖4-70(續) 艾莉颱洪三峽溪之粒徑組成百分比與中值粒徑剖面圖 ...146

(13)

第一章緒論

1.1 前言

天然河川水系多屬網路型河道，係由許多大小支流匯入到主流，

或由主流分離出支流組合而成，而都市排水系統、下水道系統及農田水利灌溉系統也與天然河川水系類似，都有渠道交匯或分流之情形。

河川或渠道交匯處由於流量增加，動量交換與能量型態轉換，導致流況變為複雜，往往造成水位劇烈變化，且河道輸砂能力也受到水流變化影響，導致河床過度的沖刷或淤積，近而危及附近橋墩及水工結構物之安全。因此，河川之整治、排水系統之規劃設計等，均應考慮到匯流問題。溪流系統動床水理演算為流域經營及防洪工程規劃設計之重要依據，若演算結果不佳，可能導致工程費偏差或洪災損失過大，

甚或危及居民性命之後果。

河川分流現象常使用於人工渠道和地下水排水系統，分流區之流況因流量改變與能量、動量之傳遞轉移為一複雜的水理現象，當水流流經分離區時可能發生亞臨界流、超臨界流以及水躍現象。而在輸砂現象中，河川中流動之泥沙，依其移動型態可大致區分為推移載及懸浮載兩大類，前者導因於作用河床面泥沙顆粒在水流方向之阻抗力；

而後者則導因於水流亂流之擴散現象，兩者之物理或力學機制不一，

(14)

其區分亦無明顯界限；推移載受到床面亂流的作用，部分較細顆粒被帶走而形成懸浮質，懸浮質又因重力作用而變成推移質。而同一種粒徑顆粒由於局部的流況不同，可同時成為推移質或懸浮質。故如何釐定兩者之分界，成為近年來學者專家們所關心之課題。

雖然有許多學者針對輸砂模式特定的物理現象做過深入研究，如泥沙起動條件、跳躍特性、沉降速度等，配合水工模型試驗，相關的成果也相當豐碩。利用水工模型試驗雖可以直接觀測流場之變化，但由於模型尺度難以拿捏，經費龐大且無再利用價值，因此較適於局部現象之探討。反觀現今之計算機速度日趨進步，數值模式則具有較省時省力，且適用於大區域流場及長期現象模擬及再利用性等優點。然而，迄今仍無法充分瞭解水流與泥沙間之力學機制，尤其對非均勻沉滓更是如此，因此在模式的建立過程中常需藉助經驗公式來計算相關物理量，而這些經驗公式本身的適用性便影響模式的可靠度，使得數值模式之適用性與實用範圍常受到限制，這也是數值模式歷經二、三十年的研究發展至今仍不間斷的原因。

本研究即利用一維 EFA 動床模式，考慮河川合分流的現象與河川中泥沙河床載與懸浮載間之交換，並以實驗案例及實際案例來驗證模式之正確性。

(15)

1.2 文獻回顧

1.2.1 定床交匯渠流

早期學者對於交匯渠流之探討，常假設匯流處無能量損失，如 Preissmann type scheme (Cunge,1980) 或匯流處水位相同，如 Abbott-Ionescu type scheme (Cunge,1980)。Preissmann type scheme 利用交匯區三格點能量相等方法，作為內部邊界條件；Abbott-Ionescu type scheme 則利用交錯格點的配置方法，於主支流交匯點計算水位。

上述兩種方法，皆利用雙掃法求解矩陣且皆屬於隱式解法。吳(1995) 模擬主支渠交匯流況時，以單側法模擬主渠流之側向對流效應，故以消散界面有限差分運算子將時間項類比為前項差分法，但於前一時段之值則依變數之方向性以本身和鄰近各點之權重表示擴散至相鄰各點之程度，並視水深為無方向性之擴散，而速度則具方向性之擴散。

Taylor (1944)以動量方程式推導主渠上、下游參數之匯流公式，

其預測值在45∘匯流角時結果尚佳，但於 135∘時差異較大。同時他亦指出此項偏差是由於大的交匯角度於渠流匯合後下游速度分佈發生扭曲，與他所作的渠流與渠壁平行的假設不符。Webber (1965)則假設非旋性與無滯性流體，導出一上下游水深，主支渠交角與下游福祿數之關係式，並以實驗修正其關係式。Abbott (1979)利用消散介面概

(16)

念（dissipative interface）於二維明渠直交匯流之模擬，但是其模式僅可模擬直交匯流的明渠流況，對於非直交匯流的流況則無法模擬。

Atkin (1980)等人則認為高雷諾數二維分離流之數值預測，使用單側差分運算子之模擬結果較使用中央差分法準確。Hager (1984)以能量率和動量法則解析匯流流況；他採用Gardel (1957)之資料認為在匯流角為90∘時，側流流向角為匯流角之 8/9 倍。Pavlovic et al. (1984)、

清水與板倉(1986)曾量測匯流段之流場，其主要目的均在驗證所建立之二維數值模式。在Pavlovic et al.(1984)之研究中，匯流渠道均為一寬度漸變之渠道，下游渠寬為上游之兩倍；在清水與板倉(1986)之實驗中，支流寬度僅為主流寬度之1/16，除有分流區外，均無明顯之水面洩降或水躍現象發生。

Ramamurthy et al. (1988)量測支流邊牆壓力差，於流量比為 0.4~0.77 時，顯示支流傳遞至主流向之動量與左、右岸壓力差之比值為流量比之正切函數，且下游能量修正係數和流量比成線性關係，但其結果僅適用於主、支渠流量比 0.5 左右。Ramamurthy and Satish (1988)對於 90∘主、支流交匯所產生的分流現象提出理論性的模式解，此模式利用質量守恆、能量守恆及動量守恆的觀念而得，並把主支流渠道的寬度比率納入考量、且提出主支流流量比跟主流下游的福祿數有一定的相關性。Hager (1989)考慮壓力修正項而求得上游水深

(17)

和臨界水深比值之函數，並推估因迴流區造成之寬度收縮係數；經由其實驗值之比較，發現於 90∘時之水深比預測值有高估現象。其表示式於側流對主流動量傳遞純由匯流角度控制，而在 90∘時無側流之動量傳遞，此和 Ramamurthy et al. 之結果有明顯之差異。 Ramamurthy et al. (1990)則提出主、支流流量比跟主流上、下游福祿數的相關性。Shettar and Murthy (1996)利用數值方法求解二維水深平均方程式的方法，採取控制體積的方式，求解直交主、支流交匯問題，

對於紊流應力影響採用深度平均k-ε 方式設計，模式則利用實驗資料與之比對，在水面高程和流速剖面的結果比對上有不錯的結果。Hsu et al. (1998)應用質量守恆、動量守恆與能量守恆方程式，綜合成一無因次多項式，透過實驗資料所得之主支流各項矯正係數，透過迴歸分析的方法，進行直角匯流實驗無因次水深比的預測。

國內學者王與何(1984)應用模型實驗將交匯渠流損失水頭分為摩擦損失及混流損失，並將損失表示為各相關水理參數和幾何參數之迴歸關係式，其結果用於實際河段之模擬且十分吻合。伍(1985)、徐(1987) 等均利用動量守恆原理並配合實驗資料分析，建立匯流段上、下游水深間的關係式。

國外McGuirk and Rodi (1978)及許與巫(1990)利用深度平均二維模式計算直交側入射流和主流交匯所造成之迴流長度；其結果和主流

(18)

與側流平均流速之動量比呈對數上升，然其主流與側流之寬度差異相當大，故並不適用於一般明渠匯流之模擬。吳與許(1992)應用質量、

能量和動量守恆理論推導匯流處附近亞臨界流流場之一維穩態解析解，並進行可變交匯角匯流實驗加以分析。王等(1997)應用第二種多方式特徵法為核心數值解法，以矩陣解法或試誤法解決所有的外界與內部的邊界條件問題，將外部、內部邊界點與渠道節點分開處理。前者以矩陣疊代法求解，後者以顯式法求解。

1.2.2 動床及動床交匯渠流

目前常見之動床數值模式多採用分離演算法，缺少描述懸浮載、

河床載兩者間之運移及交互作用。Simons et al. (1979) 發展之模式中並沒有考慮懸浮載之運移，且在單一粒徑限制下，無法描述混合砂礫層之變化行為。Pavlovic et al. (1985) 發展模式中，懸浮載只利用經驗公式求解而忽略移流及擴散特性，且河床載與懸浮載間無交換作用，

亦無法描述非均勻床質之變化。Usseglio-Polatera & Cunge (1985)發展之數值模式，只注重於懸浮載之運移，對於河床載及河床之發展則缺乏描述，而且僅限於均勻質條件。TABS2 模式(Thomas & McAnally, 1985)係利用有限元素法求解污染傳輸及河床變動問題，但該模式只能處理單一粒徑均質輸砂問題，無法考慮非均勻沉滓之差異沉淤特性，且模式採用完全分離(uncoupled)演算法，對於河床與水流之相互

(19)

作用將無法正確模擬。HEC-6 (Thomas,1977)、CHARIMA (Yang, Holly and Cunge,1986)也為分離演算法模式之屬。而運用結合演算之模式有 Holly 及 Rahuel(1990) 之一維有限差分法結合演算模式 (SEDI-COUP)，其考量懸移載、河床載之運移及兩者間之作用，並考慮期間之交換機制(mechanism)、空間延遲效應(spatial-delay effects) 之不平衡推移載及不均勻河床質之篩分及護甲；Correia et al. (1992) 發表之完全變量流結合模式；劉氏(1992)發展之非均勻質結合演算模式，其後吳氏(1993)之研究則加入護甲效應；陳氏(2002)則採用一維顯式有限解析法(Explicit Finite Analytic method, EFA, 葉等，1996、

1997)，考慮懸浮載與河床載不同之運移機制，將兩者予以分開計算。

本模式延續陳氏(2002)採用一維顯式有限解析法，將懸浮載與河床載予以分開計算。當水體流動對底床所產生之作用力大於河床質之起動條件時，便會形成河床載傳輸，此種傳輸型態受到重力之限制僅侷限於底床附近所謂的侷限層範圍內，此範圍之定義可參考 Bagnold (1973)及 Einstein (1950)。

關於河床載的計算方式，由於侷限層的高度較一般水深小很多，

加上受重力的侷限，可以直接以傳輸通量來表示，無須考慮其在垂直方向上的變化。文獻中計算河床載傳輸通量的經驗公式頗多，常見者如Meyer-Peter & Muller (1948)、Einstein (1950)及 van Rijn (1984a)等，

(20)

這些輸砂公式通常以某些代表粒徑(最常見者為中值粒徑)及粒徑分佈標準偏差來表示河床質粒徑特性，並依此計算整體河床載傳輸通量，

再依照河床質的質量守恆定律求解底床高程的改變量。

懸浮載傳輸的範圍即在參考高度(reference level)上方至水面的空間，可藉由求解移流-擴散方程式來加以計算。為完整的呈現出懸浮質濃度分佈，引用描述濃度剖面型態之經驗式，並根據局部流況與沉滓條件求解近似的「平衡濃度剖面」與「水體承載之沉滓濃度剖面」

(flow-laden sediment concentration profiles)。

關於描述平衡濃度剖面型態之經驗式，從早期適用於細砂的 Rouse (1937)方程式，進而針對粗砂修正以及考慮懸浮沉滓對 Von Karman 係數的影響等(Einstein & Chien, 1955)，有關平衡濃度剖面的研究已漸趨成熟且適用範圍更廣泛。Chiu et al. (2000)利用二相流 (two-phase flow)方程式推導定量均勻流情況下沉滓濃度剖面，其與擴散-沉降方程式所表示之平衡濃度剖面的差異，在於此法考慮了沉滓顆粒間相互作用與沉滓的慣性作用，為較一般化的表示式，而在忽略以上兩種作用的情況下，則退化成如同Rouse 方程式的形式。

而在具備有模擬支流功能之動床模式方面，前人也有相關之研究，就前人發展之動床主支流模式簡述如下：

(21)

1. CHARIMA 模式

CHARIMA模式為一維變量水流及沉滓演算模式，可用於辮狀沖積河槽水面及河床演變之長期模擬。該模式係由愛荷華大學之愛荷華水理研究所所發展出來，係採用Presissmann's有限差分近似法及 Newton-Raphson疊代法來模擬一維河床質載之輸砂現象。其模擬之河川渠道可為規則與不規則；另外，此模式能夠處理由單一渠道至網路渠道之河床演變，但本模式只適用於亞臨界流況之水理現象，對於超臨界流況發生時，尚未有處理能力。

2. HEC-2SR 模式

HEC-2SR（HEC 2 with sediment routing）模式為李日明於 1980 年所發展完成，其目的在模擬集水區沉滓之產量及河川淤積與沖刷現象。採用分離演算法就水流及輸砂部份各別演算；亦即當進行每一時段水流部份演算時，假設河床之變動微小可忽略不計算，係借自美國工兵團所發展之HEC2 模式，該模式為目前被廣泛應用之迴水演算模式，惟不考慮動床及輸砂現象，模式中另包括一輸砂演算模式，銜接 HEC2 模式演算之結果，而建構一完整之沖積河川數值模擬模式。

3. KUWASER 模式

KUWASER（known discharge , uncoupled , water and sediment routing）模式為 Colorado State University 之 Simons et al.於 1979 年所

(22)

發展完成，其應用僅限於亞臨界流，模式無法推測河槽粗粒化或 2D flow 之效應，但可反映主槽及支流之流況變化，亦可反映具有沙洲、

沙灘、島嶼之合分流河槽，適用於長期河川沖淤分析之需。

4. NCCHE 模式

NCCHE 一系列模式係由王書益教授所帶領美國國家計算水科學及工程中心(National Center for Computational Hydroscience and Engineering， The University of Mississippi )所發展之模式，此模式具有計算水理與輸砂功能，可模擬水利結構物如丁壩、圍堰、溢流堰、

河道縮收與突闊等結構物，亦可考慮到河道植生與主支流之狀況，更成功地應用於海岸建築物周邊地形變化的驗證，波浪所導致的河口地形變化的預測，以及港口防砂堤的方案設計。

5. NETSTARS 模式

NETSTARS （ network of stream tube model for alluvial river simulation）模式為李鴻源教授等人於 1993 年所發展完成，整合 BRALLUVIAL、GSTARS 與 CHARIMA 三個模式之模擬功能後，陸續發展修正而成。為一個擬似二維沖淤模式，可模擬主支流、複雜河系、定量流及變量流之水理及相對應底床沖淤特性，同時藉由流管之重新分配進行輸砂演算，可以反應河床橫斷面之不規則變化。適用於一般網路型河川及水工模型試驗沖淤問題的模擬。

(23)

6. UUWSR 模式

UUWSR（uncoupled , unsteady water and sediment routing）模式為Colorado State University 之 Tucci et al. 於 1979 年所發展完成，具有模擬島嶼、分歧、合流等中度複雜之河流系統，亦可模擬堤防、碼頭、船塢及水壩等水利構造物之效應。

1.3 論文架構

近年來主支流之處理利用水深相等及交匯區無摩擦損失(Cunge, 1980)之觀念，有了一維模式之處理方法，本模式乃採顯式法處理斷面流量與面積。陳氏(2003)已利用此法應用於自然河川水系之研究，

但僅止於定床數值模式。

本模式延續陳氏(2002)，採用一維顯式有限解析法(explicit finite analytic method, EFA, 葉等，1996、1997)之動床數值模式，並結合陳氏(2003)一維定床模式應用於主支流交匯之研究，擴充至含主支流之動床數值模式。以顯式有限解析法求解水理連續方程式及動量方程式，輸砂數值則採用半隱式法差分聯立求解懸移載質量守恆方程式、

作用層質量守恆方程式及整體河床輸砂質量守恆方程式。本研究於處理主支流交匯之方法應用水位相等法，此法除陳氏(2003)應用於自然河川水系研究上之外，於李氏(2002)之定床合分流研究，應用於淡江

(24)

大學水工試驗所之實驗資料，由模式驗證之結果與實測值相當接近。

在模式應用方面，有鑑於動床主支流實驗案例之缺乏，故延續陳氏 (2002)採用 Suryanarayana (1969)之水槽動床實驗案例之單一渠道模擬，假設相同條件渠道作為支流匯入之模擬；而於天然河川之應用，

收集大漢溪及三峽河流域之斷面資料，以大漢溪為主流，三峽河為支流，進行颱洪事件模擬，評估模式於天然河川之適用性，並期望能廣泛應用於洪水預報、河川治理規劃等方面。

(25)

第二章理論基礎

2.1 模式概述

本動床數值模式可分成水理及輸砂兩大部分。在水理計算方面，

在滿足de Saint Venant 之基本假設下，採用顯式有限解析法（EFA）

模式直接求解水流之連續方程式及動量方程式。EFA 法有推導容易及精度良好之優點，且此法在計算流力與水力計算領域之應用已證明成果相當不錯。

在沉滓運移計算方面，因天然河床係由多種不同粒徑之沉滓所組成，且由於水庫迴水區甚長，能流淤在壩址前之沉滓一般均為細顆粒者，故經由水庫防淤操作排放至下游河道之沉滓多屬細小之黏土或粉土，其特性將與其下游河床上之底質有所差異，因此所採用之模式須具有模擬非均勻沉滓之特點，以反映水庫排淤部份以及河道原有沉滓部份之不同特性。此外，為考慮懸浮載與河床載不同之運移機制，故將兩者予以分開計算，並考慮沉滓在渠道底床附近發生沉淤與再懸浮之情形，因此模式引入懸浮載與河床載間之交換機制，藉以推估水體中各懸浮沉滓之濃度變化，以及河床上床質粒徑之組成。本研究採用適於雙曲線型方程式之特性法求解懸浮載質量守衡方程式，並與河床載質量守衡方程式及整體河床輸砂之質量守衡方程式進行結合演

(26)

算，利用 Newton-Raphson 疊代聯立求解。茲分述水流及沉滓運動之控制方程式與模式所採用之數值方法如下:

2.2 水流控制方程式

水理演算係根據de Saint Venant 所推導之一維緩變非穩流控制方程式計算之，而de Saint Venant 之基本假設如下：

1.流速均勻分佈:

流速均勻分佈在通水面積上，即每一個通水斷面積僅存在一個流速，此即一維水流。

2.靜水壓分佈:

假設渠道中水流之垂向流線曲率很小而且忽略其垂直加速度，因此水深方向速度梯度為零，可忽略垂向加速度，則假設成立。

3.渠道定量流摩擦損失估計:

渠底摩擦與紊流效應對水流所造成的損失，可以定量流摩擦律估算。

4.底床坡度甚小：

當假設成立時，重力沿渠道所造成的分力將會很小，甚至可忽略不計，亦即水深可以垂向水面與渠底高程差表示。

(27)

5.忽略柯氏力及風力的影響:

以重力、靜水壓、摩擦力為主。

對於不可壓縮水流之控制方程式，包括水流連續方程式與水流動量方程式，為如下形式。

水流連續方程式 − =0

∂ +∂

∂

ql

x Q t

A (2-1)

水流運動方程式

⁰

2

=

−

∂ + + ∂

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛

∂ + ∂

∂

l l f

w gAS qu

x gA Z A Q x t

Q β (2-2)

式中，A =通水斷面積；Q =流量；t =時間；x =沿渠道中心線之距離；

g =重力加速度；Z =水位；

β

w_{=動量校正係數；}_q

l =單位渠長之支流側流量，q _l為正屬合流之處理，q _l為負屬分流之處理；u_l =支流在主流方向的速度分量；

3 4 2

2

R A

n Q

S_f =Q =摩擦坡降，其中 R =水力半徑，n = 曼寧值。

2.3 輸砂控制方程式

輸砂控制方程式將河道輸砂通量分離為非均勻之懸浮載與河床載兩部份，同時求解某一粒徑之懸浮載、河床載質量守衡，及整體河床質載之質量守衡等控制方程式，分別表示如下：

(28)

k k k

k S

x AK C x x

Q C t

A

C +

∂

= ∂

∂ +∂

∂

∂( ) ( ) ( ) k=1,2,...,n (2-3)

( ) 0 )

1

( + − =

∂ +∂

∂

− ∂ ^k ^m ^bk S_k S_ak x

Q t

p β BE

k=1,2,...,n (2-4) ⁽¹ ⁾ ⁽ ⁾ ⁽ ⁾ ⁰

1

=

∂ + + ∂

∂

− ∂

∑

ⁿ₌

k

k b

b S

x Q t

p BZ (2-5)

上三式中， =某一代表粒徑 k 之懸浮載濃度；u =水體流速；K=懸浮載之擴散通量(flux)；

Ck

ρs=沉滓密度； =孔隙率；p βk=作用層內粒徑 k之百分組成； =作用層厚度； =粒徑 k 之河床載通量； =粒徑 k之懸浮載資源項； =粒徑 k 於作用層底部資源項； =底床高程；

B=渠道寬；n =非均勻沉滓之代表粒徑數。

Em Q_bk S_k

Sak Z_b

在本模式主支流輸砂控制方程式處理上，由於程式將河道輸砂通量分離為懸浮載與河床載兩部分，故支流之側向輸砂量也分兩部分給定之。在懸浮載方面，將支流匯流點處之懸浮濃度給定主流，再將主流匯流點處之懸浮濃度代回支流匯流點處，反覆運算之；而河床載處理上，將支流匯流點處之河床載通量給定主流，此時主流之匯流點處之河床載通量即為主流河床載通量加支流河床載通量，如此反覆疊代運算而成。

為求解上述三式，尚需一些補助關係式，用以決定、、及

。當水庫排放此種具凝聚性或黏性之沉滓時，由於濃度極高，其運移行為不同於無黏性沉滓者。為求模式之一般性，除考慮無黏性沉滓

Qb E_m S

Sa

(29)

之運移機制外，亦將黏性沉滓運移之機制加入模式中，茲將其相關理論分述如下。

2.4 無黏性沉滓輸砂輔助關係式

2.4.1 河床載通量(bed load flux,

Qbk

)

Van Rijn (1984a)係以中值粒徑為代表粒徑來計算非均勻河床載之體積通量。其後，Spasojevic (1990)針對每一粒徑(D_k)，將 Van Rijn 公式稍加修正。假設河床載運移僅發生在作用層內，其中某粒徑之百分組成表示為βk，且考慮較小粒徑在水體中會形成懸浮載，可引入懸浮載與河床質載之比值γ (Van Rijn, 1984b)，對河床載通量作修正。此外，在一般非均勻之河床質中，較細顆粒可能被隱藏在較粗顆粒之間，而不易被水流帶動，故Karim et al. (1987)提出一經驗因子，稱之為隱藏因子(hiding factor，^ζ )，對河床載通量予以修正。綜合上述之影響因子可得粒徑為D_k之河床載通量如下：

^B

D D T gD s B

D q Q

k k k k j

k k

t b k k

bk 0.3

* 1 . 2

) 1 ( ) 053 . 0 ( ) 1 ( ) ( )

1

( − = − −

= γ ζ β γ ζ β (2-6)

上式中，

− =

= ₂ ¹³

* ( 1) ]

[ ν g D s

D_k _k 無因次顆粒粒徑； − =

= ₂

* 2

*

) (

c k c k

k u

u

T u 輸送參數；

c g

u_* = u 有效河床剪力速度; = ⋅ )= 3 log(12

18 D

c d

90

顆粒蔡司係數；

(30)

=

= ρ ρ_s

s 砂比重； ν =運動滯度；u_*_c = 臨界剪應力； ⁰^.⁸⁵

50

) (D

D_k

k = ζ

⎪⎪

⎩

⎪⎪⎪

⎪

⎨

⎧

≤

<

+

≥

=

4 . 0 0

10 4

. 0 )

ln(

325 . 0 25 . 0

10 1

*

k k k

k

w u

w u w

u

w u

γ

上式中， wk=粒徑 k 之沉降速度。

2.4.2 作用層厚度 ( ^E

m

)

沖刷現象發生時，根據 Bennet 與 Nordin (1977)之研究，可以下式表示：

Em

⁽ ¹ bⁿ⁾ (2-7)

n b

m C z z

E =− ⁺ −

式中，C 為數值參數。當河床表面接近護甲條件時 (armored condition)，作用層厚度接近零，在這種情況下，可用 Borah et al. (1982) 所提出護甲層之厚度(armored-layer thickness)，予以修正：

^E ^C ^Z ^Z ^K ^D^m_p

m k

k n

b n b

m =− − + −

∑

= +

1 ) 1

( ¹

β (2-8)

式中，為不動之最小顆粒粒徑。另外，作用層在淤積期間可定義為：

D

_m

¹ ⁽ ¹ bⁿ⁾ (2-9)

n b n m n

m E Z Z

E ⁺ = + ⁺ −

上式中之上標代表n 及(n+1)時刻。

(31)

2.4.3 懸浮載源 (

Sk

)

懸浮載源係由懸浮質向下之通量與底床亂流剪力作用產生河床質向上之通量交互作用之結果。使懸浮質下移到河床表面，主要是受到重力之影響。懸浮質之沈降速度決定向下沈淤之通量，所以對某一粒徑k 之懸浮質而言，其向下之通量可表為：

Qd_k =BwkCd_k (2-10) 式中，

_d ^k C_k u C w

k [3.25 0.55ln( )]

κ *

+

= (Lin,1984)；

其中，C_k=顆粒 k 之平均濃度；

κ

=von Karman 係數。

另一方面，床面沈滓成為懸浮質，主要受到底床之亂流作用所造成。Bennett and Nordin (1977)認為對某一粒徑 k 而言，河床質向上之通量可表為：

q_e_k =

ρ

_sw_lkC_a_k (2-11) 式中，

Θ

−

=3.2 4.5log u*

w_lk

2 .

< 1

Θ

( Hu & Hui, 1996 ) ； =3.1

Θ > 1 . 2

s gD

b

) (ρ ρ

τ

= − Θ

(32)

3 . 0

* 5 . 1

015 . 0

Dk

T a C_ak = D^k ^k

( Van Rijn, 1984b ) ；

wlk =顆粒向上躍起之速度；

τ

b=底床剪應力； =砂丘高度之一半。

由(2-10)及(2-11)式知，對某一粒徑 k 之懸浮載源可表為：

a

S_k = B(w_lkC_a_k −w_kC_d_k ) (2-12)

而針對沈滓交換速率可能會高估的問題，陳氏(2002)已引用許氏(2002) 之積分法修正之，在此不加詳述。

2.4.4 作用層源 (

Sak

)

作用層源係由於母層(active stratum)頂面之升降而產生，當其下降時，

ak ⁽¹ ⁾ ^[ sk⁽BZb BEm^)]

p t

S −

∂

− ∂

−

= β (2-13) 其中，βsk為母層內某一粒徑之百分組成。如母層之厚度增加，即其頂面上升時，(2-10)式中之βsk則改為βk。

2.5 黏性沉滓輸砂輔助關係式

1. 沉降速度

沉降速度( )影響其淤積速度以及懸浮載在垂直方向之濃度剖面。細小非黏性圓球體在靜止水體中之沉降遵循Stokes law：

ws

(33)

^ρ _ρ ^ρ _ν

18 gD2

w_s ^s −

= (2-14) 對黏性沉滓而言，由於顆粒表面分子間之吸力，將凝聚成團，加速沉降之速度，根據Migiot (1989)之研究，可表為：

ws_,agg =F ×ws_,part (2-15) 式中，，＝凝聚體與單其一顆粒之沉降速度；F＝凝聚因子，

與粒徑之關係為，其中D的單位為微米。

agg

ws_, w_s_,_part

8 .

250 ⁻1

= D

F

2. 沉淤速率

根據Krone (1962)之研究指出，單位面積及時間之沉淤量可表為：

Qd =wsPdCk (2-16) 上式中，為粒徑 k 之濃度；為黏性沉滓黏附於床面之機率，可表為

Ck Pd

] 1

[ _*d

Pd = −τ τ ，當τ 大於τ_*_d時P_d為零；τ 為底床剪應力；τ_*_d為發生沉淤之臨界剪應力，Shrestha 及 Orlob (1996)根據試驗資料率定τ_*_d為 0.06 Nm^-2。

1. 沖刷速率

根據Cormault (1971)之實驗資料，單位面積及時間之沖刷量為

_⎥

⎦

⎢ ⎤

⎣

⎡ −

= 1

*e

e M

Q τ

τ (2-17)

上式中，τ_*_e為發生沖刷之臨界剪應力，且當τ小於τ_*_e時，沖刷量為零；M 為沖刷係數，採 0.0015。根據 Shrestha 及 Orlob (1996)文獻，

Qe

(34)

對於剛落淤於底床上之砂層τ_*_e採 0.06 Nm^-2，而對於較底部之砂層因有壓密因素，其沖刷臨界剪應力變大，故採0.082 Nm^-2。

2.6 顯式有限解析法之簡介

本研究採用之有限解析法為美國佛羅里達大學陳景仁(Chen C. J.) 教授所創，可分為隱式及顯式有限解析法兩種。隱式有限解析法已廣泛運用於計算流體力學而得到極佳的成果(Chen & Chen，1984，Chen

& Choi，1990)，此方法的特色為:

(1) 對於不規則的幾何邊界，仍可在卡式座標使用結構性(structured) 格網求解。

(2) 由於在個別元素體的數值離散化(discretization)係採用局部解析 (local analytic solution)來近似，因此可將數值演算的捨入誤差降到最低。

(3) 數值解之穩定性極佳，屬於無條件穩定之數值模式。

由於隱式有限解析法之解法係針對橢圓或拋物線型(elliptic or parabolic)之偏微分方程式求解，如用於雙曲線型(hyperbolic)偏微分方程(明渠水流即為一例)並不太適合，因此乃有顯式有限解析法的發展及提出(Dai，1994)，但僅針對無自由表面之 Navier-Stokes 方程式。

在求解對流傳輸方程式中，對流項以特性法觀念解得式中變量之局部

(35)

解析解，且使依時變量再透過適當給定之初始條件而求得，此即為顯式法之求解精神。顯式有限解析法之穩定性，和一般常見的顯式數值方法一樣，受可蘭(Courant)數須小或等於 1 的限制，雖不如隱式法般屬於無條件穩定，但解法較簡單，故在應用上仍有其吸引力。

(36)

第三章數值方法

河道水理計算之控制方程式為非線性聯立方程組，本數值模式沿用葉等(1996、1997)成果進行河道水理演算。水理模式採用顯式有限解析法(EFA)進行水理控制方程式之離散化，此數值方法主要係用以求解雙曲線型偏微分方程式，符合移流項之數學形式。輸砂模式方面採用與水理分離演算(uncoupled)的計算方式，即在每一計算時段內先求解水理條件，再以此推估輸砂量與底床沖淤量等，反之水理條件受到輸砂行為的影響則在累進時間的過程中反應。

3.1 求解水理方程式之數值方法

連續方程式保存保守型方程式之特性，並以控制體積的觀念來差分之，用以求得水位變化量。動量方程式則因其具有雙曲線型方程式之特性，故針對移流項之部份採用顯式有限解析法予以處理。經離散後之連續控制方程式如下：

(

1

)

¹ ¹ 0

1 1 1 1 1

1 ⎥⎤=

⎢⎡

⎟⎟⎞

⎜⎜⎛

Δ

± Ψ −

−

⎟⎟+

⎜⎜ ⎞

⎛

Δ

± Ψ −

Δ +

− ₊⁺ ₋⁺ ⁺ ₊ ₋

+ A Q Q q Q Q q

A _iⁿ _iⁿ _liⁿ

c n

li n i n i c n i n i

2

2 ⎥⎦

⎢⎣ ⎝ x ⎠ ⎝ x ⎠

t (3-1)

式中，為未知數，上標為(n+1)者，係先給定 n 時刻之量測值，經反覆疊代後，再將(n+1)時刻所計算之值帶入；

+1 n

Ai

ψ

c 與以下的

ψ

m_分別

為連續方程式與動量方程式之時間加權因子(time weighting

(37)

factor) ，其範圍在［0,1］之間。為支流流量，合流時為負，分流時為正。

qli q_li

qli

經離散後之動量控制方程式如下：

( ) ( )

0

1

1 1 1

1 1

1 1 1 1

=

− +

⎥⎦

⎢ ⎤

⎣

⎡

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛

Δ Ψ −

−

⎟⎟+

⎠

⎜⎜ ⎞

⎝

⎛

Δ Ψ −

+

⎥⎦

⎢ ⎤

⎣

⎡ ⎟⎟

⎠

⎜⎜ ⎞

⎝

⎛ Δ Ψ −

−

⎟⎟+

⎠

⎜⎜ ⎞

⎝

⎛ Δ Ψ −

Δ +

−

+ + +

+ +

+ + +

+

n li n li n f n i

d n l n n r m d

n l n n r i m

d n l n n r m d

n l n n r i m n n i

u q S gA

x n

Z gA Z

x n

Z gA Z

x n

v Q v

x n

v Q v

t Q Q

ξ

ξ β ξβ

i

(3-2)

式中，為未知數，而結合特性線與有限解析法之觀念，依水流方向及流況採用不同的差分方式。當流況為亞臨界流(-1<Fr<1)時，則

=i+1， =i-1， =2，代表中央差分；當流況為正（向下游）之超臨界流(Fr＞1)，則 =i，l =i-1， =1，代表後項差分；當流況為負（往上游）之超臨界流時(Fr＜-1)， =i+1，l=i， =1，代表前項差分，i 為計算點位置。上標者為時間點，下標者為空間位置，為時間間距，為二斷面之間距，下標符號

+1 n

Qi

r

l nd

r

n_d

r

n_d

Δ t

Δ x

^ξ 者為 n 時段上之特性線位置，

該特性線係由n+1 時段上之計算點向後（backward）投射，此為顯式有限解析法求解特色之一。

3.2 求解輸砂方程式之數值方法

由於各物理量具有高度之相關性，如河床質與懸浮質藉著懸浮載源Sk連結，河床高程及作用層內之粒徑百分組成又因懸浮載源而變動

(38)

等，故有賴結合演算方式求解此三條方程式，所採用之方法為半隱式法差分聯立求解。因懸浮載之質量守恆方程式依物理特性可分割為移流及反應項(advection-reaction term)與擴散項(diffusion term)二部份，

故首先將懸浮載質量守恆之移流及反應項與作用層質量守恆方程式、整體河床質載守恆方程式，利用 Newton-Raphson 法疊代聯立求解；然後，根據所獲得之變數值，再與懸浮載質量守恆方程式之擴散項反覆疊代至收斂為止。各控制方程式離散後之形式如下：

(i)作用層質量守恆差分式：

] ) ( ) 2 [(

] ) (

) )[(

1

( ₁

1 1

1

1 +

− +

+

+ −

Δ + ⋅ Δ −

− _n

i bk n

i m k n

i m

k Q Q

BE x t BE

p β β θ

0 ) )(

1 ( ) (

] ) ( ) 2 [(

) 1

( ₁

1

1− + + + − + =

Δ

⋅

+ − Q_bk _iⁿ₊ Q_bk ⁿ_i₋ S_k S_a ⁿ_i⁺ S_k S_a ⁿ_i

x

θ θ

θ

(3-3)

(ii)整體河床質載守恆差分式：

[ ]

∑ {

=

+

− +

+

+ −

Δ + ⋅ Δ −

− ⁿ

k

n i bk n

i bk n

i b n

i

b Q Q

BZ x t BZ

p

1

1 1 1

1

1 ( ) ( )

2 ] 1 ) ( ) )[(

1

( θ

[

⁽¹ ⁾

]

⁰

]}

) ( ) )[(

1 (

1

1 1

1 − + + − =

−

+

∑

=

+

− +

n

k

n k n

k n

i bk n

i

bk Q S S

Q θ θ

θ (3-4)

圖 3.1 以一維空間為例示意沈滓之移流軌跡。就計算時刻之計算點 A 而言，移流軌跡存在於與時刻之間，定義時刻之端點 A 為到達點 (arrival point)，時刻之端點 D 為

+1

tn

+1

tn tⁿ

+1

tn tⁿ

(39)

離開點(departure point)。離開點 D 之懸浮沈滓濃度可由初始條件求得，但由於模式採用固定格點，故離開點 D 並不能保證剛好落在格點上，因此該點之濃度必須藉由鄰近格點濃度以內插的方式求得。

一般情況下，河床質移動量遠小於懸浮質移動量，因而兩者在時間尺度上相差甚大，在必須聯立求解的前提下，懸浮載方程式必須使用較大之可蘭數(Courant number)，此將導致移流軌跡穿越若干個計算格點空間，因而必須採分段處理以求得較正確之移流軌跡。假設移流軌跡從離開點 D 至到達點 A 共跨越 LNS 個計算格點空間，將該軌跡進入及離開各計算 格點空間之座標依序編號為(LNS+1)個節點，各節點間的相對 位置可以表示為：

) 2 (

) (

1 1

1 k k

k k

k

k u u t t

x

x = ₊ − ⁺ + ₊ − l=1,2…LNS (3-5)

其中，下標為節點編號， =LNS+1 代表到達點 A， =1 代表 離開點 D。利用式 (3-5)推求各節點位置，必須要先知道各節點上的移流速度，但移流速度又與節點位置有關，可利用疊代收斂的方式來推求一正確的移流軌跡。

l l l

(iii) 懸浮載質量守恆特性方程式：

(40)

當運動軌跡已知後，可積分(2-3)式，其離散化之方程式為：

2 ] ) (

[ )

( )

( ₁ ¹

1

k k k k LNS

k dep k arr k

t S t

S A

C A

C −

+

=

− ₊ ⁺

∑

= ⁺

t C C

x C t AK C

C x C

AK n

i s k i k n

i s n

i s k i k n

i

s − + Δ

Δ + − Δ +

Δ − ⁺ ⁻

+

− +

+

+ [( ) 2( ) ( ) ]

) (

) 1 ] ( ) ( ) ( 2 ) ) [(

( ² ¹ ¹

1 1 1

1 2 1

θ

θ -

x u Au C_i ⁱ ⁱ

Δ

− ₋

+

2

1

1 (3-6)

由上式可清楚瞭解到，在水深平均模式中沈滓交換速率所造成之影響會直接反應在水深平均濃度的改變上，應被視為一源項，其正確性將左右模式的表現，為影響沈滓濃度分佈的重要物理量。

S

3.3 輸砂數值解析法(Newton-Raphson 法)

經由以上數值離散處理後，在非均勻沈滓共區分為TK個粒徑區間的情況下(TK ≥2)，任一計算點共可得(2TK +1)條代數關係式，包括TK

條懸浮載質量守恆離散式、TK條作用層質量守恆離散式及 1 條整體河床輸砂質量守恆離散式。但在考慮均勻沈滓的情況下，則僅存懸浮載與整體河床輸砂質量守恆離散式各1 條，至於作用層質量守恆離散式則退化成β =1的恆等式，符合單一粒徑時之情況。

為方便說明起見，在計算點上之未知量可以向量形式表式如下：

⁺¹ =( b, ₁, ₁,..., k, k,..., TK, _TK)ⁿ⁺¹

n z c c c

s^v β β β (3-7)

或更簡潔地寫成：

顯 式 有 限 解 析 法 模 式 於 動 床 主 支 流 之 研 究