輸砂公式

本研究有關輸砂量之推估，使用以下之輸砂公式：

(1)Engelund-Hansen公式

^s

(

_g

) (

_s

)

^fd ^s

(3) 楊氏公式(Yang, 1973)，此一公式適用於砂質河床

上列各式中之各變數定義：

u ：平均流速；

γ ：水之單位重；

γs ：砂之單位重；

S ：底床坡降；

Sg ：沈滓之比重；

d50 ：底床質之中值粒徑；

BTi ：各粒徑所佔輸砂量之比例；

dmi ：底床質載之平均粒徑；

R ：水力半徑；

ucj ：粒徑之臨界起動流速；

ucmin ：最小粒徑dsj之臨界起動流速；

uct ：平均臨界起動速度；

2 1,α

α ：率定係數；

ω

：沈降速度；

U* ：剪力速度；

v ：運動黏滯係數。

2.5 不平衡輸砂量

由於上述輸砂公式所算得者為平衡輸砂量，但清水沖刷、減載及 超載等為不平衡之輸砂行為，故應考慮不平衡輸砂時空間稽延效應

（spatial-delay effects）之影響。由Bell and Sutherland (1983)之負載率

(loading law)，積分可得不平衡輸砂量公式如下：

∗

−

−

∗ ⎥

⎦

⎢ ⎤

⎣

⎡ + −

= ^{k x x} _s

s s

s e q

q

q q ^{l( 1)}

0

0 1)

(

1 ( 2.39 )

式中，

qs ：單位寬度之不平衡輸砂量；

qs

＊ ：單位寬度之平衡輸砂量；

kl ：負載率係數(loading law coefficient)，本文參考

張氏(1991)之經驗，取模擬河段長之半之倒數；

x1 ：分析河段之上游端位置；

qso ：x1位置時之不平衡輸砂量；

qso

＊ ：x1位置時之平衡輸砂量。

取上游端為原點，即 x1 = 0 ，若為清水沖刷時， qso= 0， 則(2.39)式成為：

q_s = −1 e^{−k xl( )} q_s∗

( 2.40 )

2.6 非均勻沉滓輸砂量之計算

為計算非均勻沉滓各粒徑之輸砂量，需考慮在不同沉滓級配下，

粗細顆粒相互作用之效應，如：較細顆粒可能被隱藏較粗顆粒之間，

不易被水流帶動，故須考慮細顆粒所受之遮蔽效應，以及粗顆粒較容 易顯露於河床表面，故須考慮粗顆粒之暴露效應對輸砂量之影響；然 而此種物理機制至今仍難以有效掌握，一般以遮蔽因子(hiding factor) 或暴露因子(exposure factor)作為修正依據；本研究使用二種修正方法 供選擇，第一種方法乃根據 IALLUVIAL 模式中使用之加權方法，第

先由輸砂公式計算該斷面輸砂量q_s，再考慮各粒徑所佔輸砂量之 比例

( )

BTi ，而各粒徑在傳輸質中所佔之比例

( )

BTi 是以起動機率與原交 換層內之比例來計算。由（2.32）式求得各粒徑之臨界起動流速，且

Gessler (1967)認為粒徑之起動機率，為高斯常態機率分佈，由下述積

分式可求得各粒徑之起動機率：

上式中，

Hd ：砂丘高度；

h ：水深

θ ：無因次河床剪應力 =

( )

50 0

1 gd S_g − ρ

τ ；

τ0 ：底床剪應力；

b0 ：0.079865； b1 ：2.23897； b2 ：-18.1264； b3 ：70.9001； b4 ：-88.3293。

則交換層厚度 Tm 假設如下：

d Tm

m C H

T = ( 2.46 )

在 IALLUVIAL模式中，C_Tm為一比例常數，其值為0.5。

2. Yalin (1964)之砂丘高度法

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛ −

=

0

6 1 τ τ_cr

m

T h ( 2.47 )

上式中，h ：水深 ； τ_cr：底床臨界剪應力 ； τ₀：底床剪應力。

2.8 水力篩選與河床分層

天然河川中河床質多由不同粒徑之沈滓組成，又因水流對不同粒 徑之沈滓之輸送能力不同，故河床質組成會改變，除了在水流方向產

生變化外，在垂直方向亦會發生變化，此種現象即水力篩選作用 (hydraulic sorting)。

為模擬此種現象，假設存在一交換層厚度(如圖 2.1所示)，在演 算時距(Δt)內，底床沖淤之變化量皆小於此厚度，沈滓之交換皆在此 層內進行，根據質量守恆之觀念，可推導出粒徑連續條件式[(2.9)]； 在淤積過程中，部份交換層內之河床質將離開交換層而進入不動層；

在沖刷過程中，部份不動層內之河床質將進入交換層。

為反應底床在沖淤交替時，交換層與不動層間之粒徑組成變化關 係，必須對河床粒徑予以分層記錄。因為在淤積之過程中，部份交換 層內之河床質將進入不動層，如此將改變交換層下之粒徑組成；而在 沖刷之過程中，部份不動層內之河床質將進入交換層，故除了定出交 換層中粒徑分佈外，尚須定出交換層下若干深度之不動層粒徑組成。

藉著河床分層之詳細記錄，可正確的定出各層沈滓隨時間的變化。

本模式在處理底床分層時，自河床表面往下分層，由上而下分別 為，護甲層、第零層、第一層、第二層、...、第M層，故共有M+2 層(如圖2.2，其中 DZ 為單位時間 Δt 內之刷深量)。

2.9 河床面之護甲作用

底床經水流沖刷時，由於水力篩選的作用，細顆粒沈滓較粗顆粒 沈滓容易被水流所帶動，導致粗顆粒停留在河床表面上，經若干時間

後，形成一粗化床層，稱之為護甲層(armor coat)；若河床完全護甲化，

底床將不再沖刷，即該處之輸砂量為零；由於本模式受限於輸砂量為 零時，在求解特徵值時之矩陣運算會出現溢位問題，故本模式只能模 擬接近完全護甲情形。護甲之形成，其過程為一相當複雜之物理現 象，與許多影響因子有關，包括河床質之粒徑分佈、遮蔽效應、水力 強度與輸砂量大小、床形之形成和移動及水力篩選作用等。本文將使

用臧氏(1995)根據 IALLUVIAL 模式護甲理論之修改方法，以河床面

上之不動粒徑所覆蓋之比例來判定護甲程度，並引入一護甲係數 AF，用以代表河床護甲程度之大小，其定義為任一時間，河床表面 無法移動之河床質的比例。可推得護甲係數AF：

( ) ∑

=

−

= ^m

l

i i

cumi

a ds

p V C

t

AF( ) 1 （2.48）

上式，

Ca ：為率定係數；

m_、l ：河床質中共有m_{個粒徑，粒徑}l 至m_{無法被水所傳輸；}

V_cumi ：單位時間顆粒i累積於上層之體積 ds_i ：顆粒i之粒徑。

1. 輸砂之序率特性

由 Gessler (1967)提出之床質起動機率修正C_a值，底床剪應力與

臨界剪應力相同時，床質起動之機率是 50%，則各代表粒徑停留於河

床之機率S_i，可由下式求得：

上式中，θ_c及θ分別為底床處之無因次臨界剪應力及水流之無因次剪 應力。

3. 護甲對輸砂量之修正

河床表面之護甲形成，將降低輸砂能力，防止底床持續沖刷，為 簡化起見，假設成如下之線性關係：

(

^{C AF}

)

q

q_sa = _s 1− ₁⋅ （2.53）

上式中，q_sa ：護甲後之單位寬度輸砂量；

q_s ：單位寬度之輸砂量；

C₁ ：修正係數，介於 0與1之間。

第三章 數 值 方 法

3.1 非均勻質多方式特性法

在天然河道或實驗渠道中，非均勻河床質粒徑之傳遞波、水面波 及底床波之波速差異甚大，造成泥沙粒徑與水流的可蘭數相差頗巨，

為 解 決 此 問 題 ， 本 模 式 採 用 固 定 時 間 間 格(specified time-interval scheme)多方式法，包含隱式法、傳統顯式法、時間延後法及空間延 後法，俾使每一條特性線依其在格網點上的可蘭數 Cr 之大小及設定 之最大延後數M 來決定適宜之方法。可蘭數定義為

t x t x

dt

Cr dx ⁱ

Δ

= Δ Δ

=Δ

/ /

/ λ

( 3.1 ) 其中

dt dx

i =

λ ：為一特徵值，即擾動波傳遞速度；

Δx ：模式演算之空間間距；

Δt ：模式演算之時間間距 。

計算點之格網與可蘭數 Cr 及最大延後數 M 之關係如圖 3.1 所 示。就欲演算的時刻( t=k+1 時)而言，每一格網點有 N+3 條特性線通 過以求解N+3 各未知數（P_bi , h, u, z），特性線除了通過演算時刻的未 知格點外，其另一端交於時間軸或空間軸上的參數值利用線性內差求 解之。今以Φ值表示特性線與格網之交點未知量，上述多方式法經由 線性內差如圖3.2 及圖 3.3 所示，其內差公式可表示如下：

1.隱式法(implicit scheme) (Cr >1) 3.時間延後法(temporal reachback scheme)（ 1 < Cr <1

M ）

上列各式中，

在模式初始演算時，若最初幾個演算時刻數 k(time steps)小於 最大延後數 M，令最大延後數 M 等於 k，其餘時刻之最大延後數

p ：代表待求點；

(

^{, ,} 1¹

)

⁰

1 1

1 ⁺ ₊⁺ =

−+ k

j k j k j

Fi φ φ φ i=1, ..., N+3 ( 3.13 ) 根據上式即可求得格網點上於 k+1 時段之 N+3 條相應方程式。

3.3 特徵值及特徵向量之修正法

由於特徵值及特徵向量的計算結果，會直接影響特性線在格網的 內插點位置及求解方程式[（3.6）及（3.13）式]之係數值，採用較精 確的求法對於模擬之結果有相當程度的改善。原方法[李氏（1992）、

陳氏（1994）、臧氏（1995）]是由單一格網點推求特徵值及特徵向量，

無法反應該點上游或下游流況對此斷面的影響，尤其在斷面變化或上 下游流況差距較大時容易造成質量不守恆。

為避免上述之缺點，陳氏(1996)將特徵值與特徵向量之求法改特 性線傳遞方向之相鄰兩格網點來求解。又特徵值及特徵向量是由控制 方程式[( 2.8 )至( 2.11)式]偏微分項前之係數推導而得，如 u, h, A, B, a 等，由（2.21）式，將特徵值表示如下函數：

^λ_i = ^f

( )

^ϕ i = 1, ..., N + 3 ( 3.14 ) φ= u, h, A, B, a, ...

以亞臨界流為例：介於 j-1 與 j 之正特性曲線λ有 N+2 條，介於 j 與 j+1 之負特性曲線λ有1 條；所以在上游：

_⎟⎟

⎠

⎜⎜ ⎞

⎝

⎛ +

= ⁻

2

1 j

j

i f ϕ ϕ

λ i = 1, ..., N + 3 ( 3.15 )

可得特徵值λi及相對應之特徵向量

[ ]

ai , i =1 ~ N + 3；取其中落於 j-1 與 j 間之 N + 2 條特徵值λ及相對應之特徵向量。

同理，於下游：

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛ +

= ⁺

2

1 j j

i f ϕ ϕ

λ i = 1 , ..., N + 3 ( 3.16 )

可得特徵值λi及相對應之特徵向量

[ ]

ai , i =1 ~ N + 3；取其中落於 j 與 j+1 間之 1 條特徵值λ及相對應之特徵向量。

3.4 邊界條件

由於本模式為一結合演算模式，需在同一時刻求解水流與輸砂方 程式，而在邊界上之求解未知數為流速、水深、底床高程及床質粒徑 組成，必需依據邊界之已知水理及輸砂資料建立關係式，例如流量歷 線、水位歷線及輸砂歷線等，利用這些邊界條件，再配合原特性方程 式（3.13），如此方能得到係數矩陣；若邊界所提供的水理或輸砂資 料並非求解未知數時，可根據邊界上控制體積內質量守恆之原理建立 彼此的關係式，並經由牛頓疊代法使邊界上之未知數在求解過程中，

以同時滿足邊界條件式與特性方程式。

當流況為亞臨界流(subcritical flow)時，河床波、粒徑篩分波之傳 遞方向皆由上游往下游傳，水面波之傳遞方向一與河床波、粒徑篩分 波同向，一則為逆向，即由下游往上游傳；故需一個與水流有關、一

個與河床高程及 N 個粒徑組成相關之上游邊界條件，以及一個與水 流有關之下游邊界條件。當流況為超臨界流(supercritical flow)時，兩 水面波由上游往下游傳，河床波與粒徑篩分波則往上游傳，故需兩個 與水流有關之上游邊界條件及 N+1 個與底床高程及粒徑組成有關之 下游邊界條件。

除了考慮外部的邊界條件，對於合流及分流的交匯區，必須給於 內部邊界條件，邊界條件式之處理與推導過程詳述於後。

1. 上游邊界條件

(1) 流量歷線

一般模式常採用在上游邊界所觀測到之流量歷線Q

( )

t ，作為數值 演算之一邊界條件式。因流量為流速及水深之函數，以牛頓疊代法求 解，則Q= Au可寫為

F^m

( )

u,h =A^mu^m−Q

( )

t_k+1 ( 3.17 )

( )

u h

F , 可是為誤差函數(error function)，其定義為估計值

(

^Amum

)

與 真值^[Q

( )

^{tk+1 ]}之絕對誤差。另定義前後兩次疊代誤差如下：

ΔF =F^m+1

( )

u,h −F^m

( )

u,h ( 3.18 ) (3.17)與(3.18)式中，k 為演算時間，m 為疊代次數。

當求解之未知數收斂時，F^m+1

( )

u,h =0。ΔF 之泰勒展開式如下(忽 略二次以上的項)：

^{u Fm}

( )

^{u h Q}

( )

^{tk Amum}

本模式引進緩衝段法( buffer reach method)加以改進，此法由 Cunge et al.（1980）提出，原應用於有限差分法之分離演算模式，主 要在消除不同 DX 對演算結果之影響。本文依其原理為將上游端河 段引進一緩衝段 ΔXB(如圖 3.4 所示)，並將在 sec0 之計算點移至

本模式引進緩衝段法( buffer reach method)加以改進，此法由 Cunge et al.（1980）提出，原應用於有限差分法之分離演算模式，主 要在消除不同 DX 對演算結果之影響。本文依其原理為將上游端河 段引進一緩衝段 ΔXB(如圖 3.4 所示)，並將在 sec0 之計算點移至

在文檔中 多方式特性法模式應用於動床主支流之研究 (頁 27-0)