模式理論基礎

3-1 顯式有限解析法模式理論

本研究之顯式有限解析法(explicit finite analytic, EFA)數值模式 可分成水理及輸砂兩大部分，於水理計算方面，在滿足de Saint Venant 之基本假設下，應用顯式有限解析法離散動量方程式，並利用沿特性 線積分概念，求解動量方程式，再配合適當斷面處理與差分形式，求 解水流連續方程式。顯式有限解析法有推導容易且程式撰寫較隱式法 簡易之優點，而此法在計算流力與水力計算領域之應用已證實有相當 不錯之成果。

至於輸砂計算方面，因天然河床係由多種不同粒徑之泥砂所組成，

且上下游河段之底床粒徑大小有所差異，因此所採用之模式須具有模 擬非均勻泥砂之特點，以反映河道泥砂部份之不同特性。此外，為考 慮懸浮載與河床載不同之運移機制，故將兩者予以分開計算，並考慮 泥砂在河道底床附近發生沉淤與再懸浮之情形；因此本模式引用許氏 (2002)對懸浮載與河床載間交換速率估算之研究，藉以推估水體中各 懸浮泥砂之濃度變化，以及河床上床質粒徑之組成。本研究於處理輸 砂控制方程式時，採取適用於雙曲線型方程式之特性法求解懸浮載質 量守恆方程式，並與河床載質量守恆方程式及整體河床輸砂之質量守

恆方程式進行結合演算，再利用Newton-Raphson 法疊代聯立求解。

3-1-1 水理控制方程式

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

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

2. 靜水壓分佈：假設河道中水流之垂向流線曲率很小而且忽略 其垂直加速度，因此水深方向速度梯度為零，可忽略垂向加 速度，則假設成立。

3. 河道定量流摩擦損失估計：渠底摩擦與紊流效應對水流所造 成的損失，可以定量流摩擦律估算。

4. 底床坡度甚小：當假設成立時，重力沿河道所造成的分力將 會很小，甚至可忽略不計，亦即水深可以垂向水面與渠底高 程差表示。

5. 忽略柯氏力及風力的影響。

對於不可壓縮水流之控制方程式，包括水流連續方程式與水流動 量方程式，為如下形式。

水流連續方程式

 0

部份沈滓顆粒即向上運移超越參考高度，形成懸浮載(suspended load) 傳輸型態；同樣的，懸浮泥砂在受重力的影響下亦會逐漸沈降至底床，

進而轉變為河床載傳輸型態，圖3- 1 為泥砂傳輸之示意圖。透過以上

敘述可知，泥砂會因所受作用力的不同而改變其傳輸機制，因此可將

3-2 顯式有限解析法模式數值方法

為連續方程式與動量方程式之時間加權因子(time weighting factor) ，

其範圍在［0,1］之間。

q

_li為支流流量，合流時

q

_li為負，分流時

q

_li為 向及流況採用不同的差分方式。當流況為亞臨界流(-1<Fr<1)時，則

r

=i+1，

l

=i-1，n_d=2，代表中央差分；當流況為正（向下游）之超臨界 流(Fr＞1)，則

r

=i，

l

=i-1，n_d=1，代表後項差分；當流況為負（往 上游）之超臨界流時(Fr＜-1)，

r

=i+1，

l

=i，n_d=1，代表前項差分，i 為計算點位置。上標者為時間點，下標者為空間位置，

 t

為時間間 距，

 x

為二斷面之間距，下標符號

^

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

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

當超亞臨界混合流發生時，超臨界流下游端水躍發生處，會因上、

下游水深之共軛特性自然形成一限制條件，對整個流場而言，此限制 條件即為一內部邊界條件，其存在使超臨界流區域得到多於之邊界條

件，導致數值計算發生over-determined 現象。鍾(2009)利用內部相鄰 計算點之水位高程，透過預測-修正數值方法先求取超亞臨界混合流 區域內各斷面通水面積及流量，以減少地形驟變對模式演算之影響與 降低數值震盪幅度，並增加之後代入顯式有限解析法運算之穩定性。

顯式有限解析法利用已預測修正後之通水面積、流量等數值繼續運算，

得到真正所需超亞臨界混合流區域之通水面積、流量及超臨界流區域 之下游邊界水深。鍾(2009)以此方式改良修正顯式有限解析法於超亞 臨界混合流流場計算之限制，詳細理論可參考鍾(2009)之碩士論文。

3-2-2 輸砂方程式之數值方法

由於輸砂方程式中各物理量具有高度之相關性，如河床質與懸浮 質藉著懸浮載源連結，河床高程及作用層內之粒徑百分組成又因懸浮 載源而變動等，故有賴結合演算方式求解此三條方程式，所採用之方 法為半隱式法差分聯立求解。因懸浮載之質量守恆方程式依物理特性 可分割為移流及反應項(advection-reaction term)與擴散項(diffusion term)二部份，故首先將懸浮載質量守恆之移流及反應項與作用層質量 守恆方程式、整體河床質載守恆方程式，利用 Newton-Raphson 法疊 代聯立求解；然後，根據所獲得之變數值，再與懸浮載質量守恆方程 式之擴散項反覆疊代至收斂為止。各控制方程式離散後之形式如下：

作用層質量守恆差分式：

理 以 求 得 較 正 確 之 移 流 軌 跡。假 設 移 流 軌 跡 從 離 開 點 D 至 到

為 一 源 項 ， 為 影 響 沈 滓 濃 度 分 佈 的 重 要 物 理 量 。

以上三式為非線性代數式，可線性化後利用 Newton-Raphson 法疊代

求解： 單位河川流功率。根據前述之文獻回顧，以Annandale (1995)之水力 沖蝕理論與Sklar and Dietrich (2004)之泥砂磨蝕理論較能定量描述河 道中軟岩侵蝕問題。而軟岩受到水流作用之塊狀崩落或風化沖洗等其 餘侵蝕現象，目前仍未有較具明確之經驗公式或數學方程式可描述，

多為實驗或物理過程之觀測研究，數模應用上較為困難。有鑑於此，

本研究挑選 Annandale (2006)之水力沖蝕機制與 Sklar and Dietrich (2004)之泥砂磨蝕機制作為軟岩侵蝕模式之理論基礎，本研究定義水 力作用於軟岩表面之沖刷為沖蝕，泥砂作用於軟岩表面之沖刷為磨蝕，

以下則針對相關理論加以說明。

3-3-1 水力沖蝕模式理論

一般來說，水流作用力常對軟岩節理處產生瞬間或逐漸沖蝕破壞，

或對單獨軟岩塊體表面產生沖擊性之沖蝕。此二種機制均屬水力沖蝕 之範疇。此外在紊流條件下，瞬間剪應力與流速將隨著河川平均流功 而改變，與穩流條件下之穩定作用力相比，此瞬間作用於岩床之剪應 力或壓力差為造成岩床破壞之主要原因(Annandale, 2006)。

Annandale (1995)曾以水流沖蝕能量為概念，探討河川流功與沖 蝕指數之相對關係。沖蝕指數之概念最早係由Kirsten (1982)所提出，

沖蝕指數可以以下形式所表示：

s d b s

h

M K K J

K

 (3-21)

式中，

M s

為材料強度；

K b

為顆粒/塊體尺寸；K

d

為弱面/或顆粒間之抗 剪強度；J

s

為地盤構造條件。各指數之物理特性及計算方式如下所 述：

材 料 強 度 參 數

M s

主 要 代 表 岩 石 密 度 以 及 無 圍 單 壓 強 度 (unconfined compressive strength, UCS)等特性，無圍單壓強度可由

ASTM D-7012 標準試驗獲得。Kirsten 對於計算 M

s

之方式表式如下：

若該材料為非黏性粒狀材料(non-cohesive granular material)，則可簡化 為

K b

=1000d

³

，d 為材料粒徑大小。

地盤構造條件參數

J s

，其代表岩石材料於地面顯露時，抵抗侵蝕

根據 Annandale (2006)對河川流功與沖蝕指數關係之研究，沖蝕 指數

K h

之數值可將河川流功作用之方式區分為兩種類型，分別為顆 公式型式可表示如下(Greimann, 2008)：

 

係數需由試驗結果以及現場資料進行檢定；

U 為河川斷面平均之水流

(roughness height)。

3-3-2 泥砂磨蝕模式理論

一般於河道表面以河床載形式運移之泥砂顆粒，對軟岩河道同時 具有具侵蝕特性的磨蝕效應與覆蓋保護特性的工具效應(tool effect)。

最大軟岩磨蝕量則受到河床載之臨界供砂量(critical sediment supply) 限制；若運移之河床載超越臨界供砂量，則泥砂將落淤於軟岩河道表 面形成覆蓋，並對河道表面產生具保護特性之工具效應，使其不受磨 蝕。

依據 Sklar and Dietrich (2004)之實驗結果，發現軟岩磨蝕速率與 軟岩表面張力強度成反比，但與泥砂顆粒質量成正比，並在達到最大 泥砂傳輸率後，軟岩磨蝕速率將下降。當泥砂顆粒質量較小時，磨蝕 作用將會增加更多顆粒磨蝕軟岩表面；而超過最大值時，泥砂顆粒則

Sklar and Dietrich (2004)修正 Foley (1980)理論後，提出泥砂磨蝕

為單位河寬之輸砂能力(kg/m/s)；w

s

為泥砂顆粒之沖擊速度(m/s)；Y 為軟岩彈性模數(Pa)；k

v

為軟岩強度參數，該參數值介於 10

¹²

~10

¹³

， 為泥砂非平衡狀態與平衡狀態下之輸砂量，前者為輸砂率(sediment

transport rate)，後者為輸砂能力(sediment transport capacity)。若軟岩 表面覆蓋泥砂顆粒，則由顆粒跳動所引起之軟岩磨蝕量將減少；若輸 砂量等於輸砂能力，磨蝕作用將停止。 1 / 項反應懸浮效應，

當水流強度增大時，屬懸浮狀態之泥砂顆粒比例增高，此時泥砂顆粒 於軟岩表面跳動行為減少，磨蝕量則下降。同時，泥砂磨蝕率亦與

⁄ 1 項成反比，反應其物理特性。

3-4 軟岩侵蝕模組數值方法

結合軟岩侵蝕模組之顯式有限解析法於模擬運算前，需先設定具 軟岩侵蝕特性之斷面，大安溪河段現場案例之設定區域以現地觀測及 航拍圖為主要設定依據。其後須給定軟岩侵蝕模組參數之計算方式，

大安溪河段現場案例以現地試驗實測各軟岩參數計算之沖蝕指數

K h

為參數設定依據，該參數於模式中亦可給定其相關計算參數，由模式 於運算過程中求取各時刻之代表數值。至於軟岩參數如軟岩彈性模數

Y、軟岩張力強度 σ T

，則以該區域之觀測資料為設定依據。

模式運算時，水理模組先進行運算，求取河道中各斷面之流量、

水位等水力特性參數，提供軟岩侵蝕模組計算所需之水理條件數值。

由於顯式有限解析法動床模式採非耦合運算，水理模組運算與動床模 組運算間隔不同，當時間增加至啟動動床模組運算間隔時，動床模組

蝕模組計算所需之動床條件數值。軟岩侵蝕模組運算間隔與動床模式 相同，動床模組運算終了後，軟岩侵蝕模組開始運算，利用水理模組

蝕模組計算所需之動床條件數值。軟岩侵蝕模組運算間隔與動床模式 相同，動床模組運算終了後，軟岩侵蝕模組開始運算，利用水理模組

在文檔中 一維軟弱岩盤河道侵蝕數值模式之研發與應用 (頁 31-49)