水力篩選與河床分層

⎜⎜ ⎞

⎝

⎛ −

=

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 之計算點移至 secb，其初始條件可由 sec0 及 sec1 依緩衝段長短內差而得。

此處以特性線之觀念求得緩衝段之長度，以特性線投射到空間軸 上之交點作為緩衝段之長度。通常特性線判斷標準依可蘭數 Cr 及最 大延後數M 決定，假設最大延後數 M=1 之情況下，由可蘭數之定義

t Cr x ⁱ

Δ

=Δ /

λ ，λ_i 為擾動波傳遞速度，則λ_b 底床擾動波傳遞速度。因

此處定義緩衝段為起始時刻之前一時刻之距離，故緩衝段之長度即為 前一時刻之空間間距與可蘭數 Cr之乘積(如圖 3.5 所示)，即特性線與 空間軸交點至初始時刻(即下一時刻)之距離，其數學關係式可表為

t x

C

x_b = _r ⋅Δ = _b⋅Δ

Δ λ 。緩衝段會隨著輸砂公式的選用而有不同的長度；

由（2.3）式滿足輸砂連續條件差分式而得：

( ) ( ) ( )

[

^{Q Q Q Q}

]

^{t X}B

(

^pr

)

^{B z}

k B k s k

B k

s⁺¹− ⁺¹ + 1−θ − Δ =Δ 1− Δ

θ ( 3.24 )

式中，

Qs ：上游邊界入砂量；

QB ：緩衝段下游端之輸砂量；

ΔXB ：緩衝段長度。

取 θ = 1，而ΔQ_s^k+1 =Q_s^k+1−Q_B^k+1，則由( 3.24 )可得

( )

₁

1 +

+ = _{i k}

k

i P t

P i=1,…,N ( 3.28 ) 因沈滓粒徑百分組成不易獲得，所以數值上仍考慮緩衝段法，以 緩衝段下游端之粒徑百分組成為上游粒徑變化之邊界條件，其求法乃 假設在緩衝段之中，滿足沈滓粒徑組成連續方程式[(2.4)式]：

[ ( ) ( ) ]

i^k k

S k B b k Si k Bi bf k

i C Q Q F Q Q P

P ⁺¹= ⁺¹− ⁺¹ − ⋅ ⁺¹− ⁺¹ + ( 3.29 ) 其中，

(

r

)

B

bf a p x

C t

Δ

−

⋅ Δ

= − 1

F_b ：沖刷時，F_b =P_bi；淤積時，F_b =P_boi； Q_Bi ：緩衝段下游端顆粒 i 之輸砂量；

Q_B ：緩衝段下游端之輸砂量；

Qsi ：上游邊界顆粒 i 之入砂量；

Q_s ：上游邊界輸砂量；

a ：交換層厚度；

P_r ：孔隙率。

2. 下游邊條件

本模式可有如下之選擇：

(1) 水位歷線

當流況為亞臨界流時，以下游邊界之實測水位歷線作為模式演算 之邊界條件：

h^k+1+z^k+1= y

( )

t_k+1 ( 3.30 )

(1) 合流:

將合流渠道分為三段，主流上游段、主流下游段及支流段，水流 流向及示意圖如下所示：

i i+1

j

主流

支流

圖 3.6 合流交匯區斷面示意圖

其中，i 表示主流上游段最後一個斷面點，i+1 表示主流下游段第一 斷面點，j 表示支流段最後斷面點。本研究採用兩點水位相等法，此 方法只需將渠道分為兩段，主流段和支流段，在亞臨界流況下，計算 主流段時上游給定流量，下游給定水位，當計算支流段時，支流之下 游水位未知，故假設 i 點水位和 j 點水位相同，如此合流之邊界給定，

即可算出內部斷面之流量及水深。此方法優點計算上較簡易，且要達 到兩點水位相等的收斂精度也較容易，但缺點為 i 點和 j 點水位相同

即可算出內部斷面之流量及水深。此方法優點計算上較簡易，且要達 到兩點水位相等的收斂精度也較容易，但缺點為 i 點和 j 點水位相同

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