局部極值消減法應用於二維水深平均模式之研究

(1)

國

立

交

通

大

學

土木工程學系

碩

士

論

文

局部極值消減法應用於二維水深平均模式之研究

Study on Local Extremum Diminishing Scheme

Applied to 2D Depth-averaged Model

研究生：羅琦雯

指導教授：楊錦釧博士

連和政博士

中華民國九十九年七月

(2)

局部極值消減法應用於二維水深平均模式之研究

Study on Local Extremum Diminishing Scheme

Applied to 2D Depth-averaged Model

研究生：羅琦雯 Student: Chi-Wen Lo

指導教授：楊錦釧 Advisor: Jinn-Chuang Yang

連和政 Ho-Cheng Lien

國立交通大學

土木工程學系

碩士論文

A Thesis Submitted to Civil Engineering

College of Engineering

National Chiao Tung University

in Partial Fulfillment of the Requirements

for the Degree of Master

in

Civil Engineering

July 2010

Hsinchu, Taiwan, Republic of China

中華民國九十九年七月

(3)

謝誌

承蒙恩師楊教授錦釧與連博士和政於論文研究期間之悉心指導，同時感謝吳博士祥禎於研究上之協助，使得本論文得以順利完成，在此謹致由衷的敬意與感謝。在論文審定期間，感謝口試委員許教授少華與賴教授進松之細心匡正與寶貴建議，使本論文更臻完善。感謝就學期間研究室胤隆學長、世偉學長、浩榮學長、建華學長、弘恩學長、振家學長、歆淳學姐、俊宏學長、全謚學長於課業與生活上的提攜與照顧，亦感謝同學群玲、昀軒、仁猷、明儒、俊騰、勁頤、新詠在研究生活中的陪伴。同時，由衷感謝交通大學土木系97級同學們的友情扶持與砥礪。最後謹以本論文獻給我親愛的父母與弟弟，感謝你們長久以來對我無條件的付出與支持，給我最大的力量與溫暖。

(4)

II

局部極值消減法應用於二維水深平均模式之研究

學生：羅琦雯指導教授: 楊錦釧連和政

國立交通大學土木工程學系

摘要

本研究旨在以局部極值消減法及配合有限體積法與 Runge-Kutta 多時階法求解二維水深平均方程式。本研究採用局部極值消減法，其具有最大值不會再增加與最小值不會再減少的特性，可以有效抑制數值震盪。同時，為了提升數值的穩定度，在時間離散化上採用 4 階 Runge-Kutta 時階法。本研究以潰壩模擬、穿臨界流流經底床突起之模擬、束縮段產生交波之模擬等案例，探討控制方程式之慣性項、重力項及摩擦項對模擬精度之影響，最後本模式也實際應用於現場案例東埔蚋溪之模擬。另外，針對消散黏滯係數之指數 β 值的探討，本研究採用交波模擬案例作敏感度分析，並以模擬值與實驗值之效率係數、均方根誤差、最大值之相對差異百分比三種指標作比較。經過計算值與解析解及實驗數據的比對，本研究所發展的局部極值消減法具有更佳的求解精度，而且就東埔蚋溪現場應用來看，局部極值消減法亦能處理穿臨界流等複雜流況。

關鍵詞:二維水深平均模式、局部極值消減法、4 階 Runge-Kutta 時階法、

穿臨界流場

(5)

Study on Local Extremum Diminishing Scheme

Applied to 2D Depth-averaged Model

Student：Chi-Wen Lo Advisors：Dr. Jinn-Chuang Yang

Dr. Ho-Cheng Lien

Department of Civil Engineering

National Chiao-Tung University

ABSTRACT

The purpose of this study is to use the Local Extremum Diminishing (LED) scheme integrated with finite volume method and Runge-Kutta

time stepping method for solving 2D depth-averaged equations.The LED

scheme preserves the properties of that the local maxima would not increase and the local minima would not decrease to restrain the numerical oscillation efficiently. Meanwhile, the fourth order Runge-Kutta method is also adopted in time difference term to improve the numerical stability.

In this study, the dam-break flow, the transcritical flow with hump and the cross-wave flow with contraction are simulated to discuss the accuracy affected from inertia force term, gravity force term and bed friction term in the governing equations receptively. Furthermore, the proposed model is also applied to the practical field simulation of Tung-Pu-Ruei creek. The discussions for the exponent β of coefficient of artificial viscosity with cross-wave simulation are presented. The sensitivity analysis of β is analyzed by using Coefficient of Efficiency

(6)

IV

(CE), Root-Mean-Square Error (RMSE), and Error for the Peak Value (EVP), which are analyzed on the basis of comparing the simulation and experimental data.

By comparing with analytic solution and experimental data, the computed results from proposed model with LED scheme can provide adequate accuracy. In addition, for the practical application of Tung-Pu-Ruei creek, the LED scheme can also handle the complex flow.

Keywords: 2D depth-averaged model, local extremun diminishing scheme, fourth order Runge-Kutta method , transcritical flow.

(7)

目錄

誌謝………I 摘要………II Abstract………III 目錄………V 表目錄………VII 圖目錄………VIII 符號表………XII 第一章緒論………1 1.1 研究動機與目的………1 1.2 文獻回顧………1 1.3 研究方法與流程………4 1.4 本文結構………6 第二章理論基礎………7 2.1 基本控制方程式………7 2.2 水深二維平均模式……… 8 第三章數值方法………12 3.1 有限體積法………12 3.2 局部極值消減法………12

(8)

VI 3.3 Runge-Kutta 多時階法……… 16 3.4 邊界條件設定………16 第四章模式驗證分析與東埔蚋溪河道模擬……… 20 4.1 水流經潰壩的情………20 4.2 水流在穿臨界流場中遇到突起的情………22 4.3 超臨界水流遇到河道束縮產生交波的情………25 4.4 β 值敏感度分析……… 28 4.5 現地案例:東埔蚋溪河道模擬………31 第五章結論 … … … 3 3 參考文獻……… 35

(9)

表目錄

表4.1 消散項敏感度分析各指標的結果……… 38

(10)

VIII

圖目錄

圖 3-1 邊界條件示意圖….………39 圖 4-1 潰壩案例初始條件示意圖………39 圖 4-2 β=0.2 時潰壩第 20 秒的水位比較圖……… 40 圖 4-3 β=0.4 時潰壩第 20 秒的水位比較圖……… 40 圖 4-4 β=0.6 時潰壩第 20 秒的水位比較圖……… 41 圖 4-5 潰壩第 20 秒模擬示意圖………41 圖 4-6 潰壩第 40 秒模擬示意圖……… 42 圖 4-7 β=0.2 時潰壩第 40 秒的水位比較圖………42 圖 4-8 β=0.4 時潰壩第 40 秒的水位比較圖………43 圖 4-9 β=0.6 時潰壩第 40 秒的水位比較圖………43 圖 4-10 潰壩第 20 秒流速比較圖……… 44 圖 4-11 潰壩第 40 秒流速比較圖……… 44 圖 4-12 穿臨界流場底床突起示意圖………45 圖 4-13 β=0.2 時穿臨界流未發生震波之水位比較圖………45 圖 4-14 β=0.4 時穿臨界流未發生震波之水位比較圖………46 圖 4-15 β=0.6 時穿臨界流未發生震波之水位比較圖………46 圖 4-16 穿臨界流未發生震波第 1 秒的水位瞬時變化…………47 圖 4-17 穿臨界流未發生震波第 4 秒的瞬時水位變化…………47

(11)

圖 4-18 穿臨界流未發生震波第 10 秒的瞬時水位變化…………48 圖 4-19 穿臨界流未發生震波第 16 秒的瞬時水位變化…………48 圖 4-20 穿臨界流未發生震波第 32 秒的瞬時水位變化…………49 圖 4-21 穿臨界流未發生震波 3000 次迴圈後的模擬結果…… 49 圖 4-22 穿臨界流未發生震波之福祿數變化………50 圖 4-23 穿臨界流未發生震波的剩餘值變化………50 圖 4-24 β =0.2 時穿臨界流發生震波之水位比較………51 圖 4-25 β =0.4 時穿臨界流發生震波之水位比較圖………51 圖 4-26 β =0.6 時穿臨界流發生震波之水位比較圖………52 圖 4-27 穿臨界流發生震波第 10 秒的水位瞬時變化………52 圖 4-28 穿臨界流發生震波第 35 秒的水位瞬時變化………53 圖 4-29 穿臨界流發生震波第 60 秒的水位瞬時變化………53 圖 4-30 穿臨界流發生震波 4000 次迴圈的模擬結果……… 54 圖 4-31 穿臨界流發生震波之福祿數變化………54 圖 4-32 穿臨界流發生震波的剩餘值變化………55 圖 4-33 束縮段示意圖………55 圖 4-34 β =0.2 時束縮渠道之邊牆水位比較圖……… 56 圖 4-35 β =0.4 時束縮渠道之邊牆水位比較圖………56 圖 4-36 β =0.6 時束縮渠道之邊牆水位比較圖………57

(12)

X 圖 4-37 束縮渠道第 0.17 秒的水位瞬時變化……… 57 圖 4-38 束縮渠道第 0.5 秒的水位瞬時變化………58 圖 4-39 束縮渠道第 0.7 秒的水位瞬時變化………58 圖 4-40 束縮渠道第 1.3 秒的水位瞬時變化………59 圖 4-41 束縮渠道經過 2000 次迴圈之後的模擬結果……… 59 圖 4-42 渠道束縮模擬之水位等高線圖………60 圖 4-43 渠道束縮實驗之水位等高線圖………60

圖 4-44 1995 年 Molls and Chaudhry 所模擬之水位等高線圖 60 圖 4-45 與其他模式之邊牆水位比較模擬…….……… 61 圖 4-46 三個模式與實驗值的誤差百分比……… 61 圖 4-47 束縮渠道模擬結果之福祿數變化……… 62 圖 4-48 束縮渠道模擬結果之剩餘值變化……… 62 圖 4-49 束縮渠道之中心線水位比較圖………….……….………...63 圖 4-50 三個模式與實驗值的中心線水位誤差比較……..………...63 圖 4-51 束縮渠道模擬結果之效率係數變化………….……… 64 圖 4-52 束縮渠道模擬結果之最大相對誤差變化………….………64 圖 4-53 束縮渠道模擬結果之均方根誤差變化………….…………65 圖 4-54 β與CPU 耗時變化……….………65 圖 4-55 β與CPU 效能變化………..………….……… ..66

(13)

圖 4-56 東埔蚋溪流域概況圖………67

圖 4-57 東埔蚋溪河道示意圖………68

圖 4-58 東埔蚋溪河道模擬之水位變化………68

圖 4-59 東埔蚋溪河道模擬之福祿數變化………69

(14)

XII

符號表

A 面積 Fs 水面運動邊界 Fb 底床運動邊界 n N 法線方向 n 曼寧係數 g 重力加速度 g = 9.81 ( m/sec2 ) H 水面高程 H = 水深 + 底床高程 = h + Zb h 水深=水面至渠道底床之距離 P 壓力 R 殘餘值 Sox

、

Soy x

、

y 方向的底床坡降 Sfx

、

Sfy x

、

y 方向摩擦坡降 t 時間 V 流速 u

、

v

、

w x

、

y

、

z 方向之平均流速 uj ,uk 為格網點j、k 的 x 方向的流速 x

、

y

、

z 正交卡式座標 x Δ

、

Δy、Δz x.y.z 方向之格網間距

(15)

Zb 底床高程 β 消散黏滯係數之指數 ε 人工黏滯係數 ρ 流體密度 τ 應力張量 υ 消散黏滯係數 Γs 固體邊界 in Γ 入流邊界 out Γ 出流邊界

(16)

1

第一章緒論

1.1 研究動機與目的

當河川特性為長度短、坡度大，而且河床坡度隨山區地形變化而有陡坡與緩坡互相接連的現象，常造成亞臨界與超臨界流同時存在於同一流場中，如發生水躍、跌水、乾床等不連續面之流況，這是使得一般河川數值模式應用於天然河道中常發生不易收斂的最大原因。不易收斂的情況可能導致模擬結果的誤差擴大，甚至是程式發散。為解決不連續面以及數值震盪所產生的影響，連和政等(2007-2009)將局部極值消減法應用於一維淺水波方程式之求解，本研究採用局部極值消減法結合有限體積法及 Runge-kutta 多時階法，進一步求解二維淺水波方程式。本研究模擬兩種不連續面情況以及複雜流場的交波案例，分別是潰壩模擬、穿臨界流流經底床突起的案例以及超臨界水流流經束縮段產生交波的案例。討論三個案例的模擬結果，藉此驗證模式的精確度及穩定度。同時也將模式初步應用於東埔蚋溪之河道模擬，提供未來模擬現地案例時，解決不連續面流場問題的較佳方法。針對消散項之指數做敏感度分析，以超臨界水流流經束縮段產生交波為例，對控制消散項係數做敏感度分析，找出適當的係數範圍。

1.2 文獻回顧

隨著數值方法的演進，針對雙曲線偏微分方程式所面臨的不連續問題，已經有很多數值方法可以克服如Toro（2001）發展的總變量消

減法（total variation diminishing，簡稱 TVD）及 Jameson（1981）發展的局部極值消減（local extremum diminishing，簡稱 LED)。已經有

(17)

許多學者採用包含完整源項之二維水深平均模式，求解穿臨界流場問題，而得到不錯的成果 Toro(2001)以總變量消減法求解淺水波方程式， TVD 法廣泛被應用到航太與水利工程上的數值求解上。優點為其方法可以捕捉複雜震波問題，求解不連續解有不錯的結果，且可抑制非物理性的震盪，近年來在計算水力學上也有很大的發展。古孟晃(2003)以水深平均為基本假設，採用總變量消減法(TVD)之有限差分數值模式，模擬一維、二維淺水波方程式之變量流問題，探討乾濕交界之處理、能量坡降與地形變化之處理。林宗賢(2007)利用高精度TVD格式初步模擬台北三重地區淹水範圍，求得淹水水深資料及範圍，利用GIS 系統加以套疊處理，提供都市防洪資訊，再配合淹水潛勢資料可作為防災計畫與應變措施訂定之依據，來測試模式對三重地區的適用性。吳銘順(2006) 採用混合型的總變量消減法，搭配限制式模擬各種潰壩的情況，討論乾濕交界的處理，並實際應用於台北新店溪流域的模擬。 Lin et al.等(2003)比較四個有限體積混和型智慧 TVD 模式，分

別為 the Liou–Steffen splitting (LSS) ， van Leer splitting ，

Steger–Warming splitting and local Lax–Friedrichs splitting schemes，以

各種潰壩及水躍之流況來模擬驗證，結果顯示以 LSS 模式的模擬結

果最佳。Lin et al.(2005)以高解析二階 TVD 模式進行理想之潰壩模擬，

其中flux vector splitting (FVS)和 direct MUSCL-Hancock(DMH) slope

limiter approach 結合後之模式的模擬最為相當理想。

總變量消減法在處理不連續面上有不錯的成效，但是在計算上的過程較為繁複，應用時必須針對所採用的控制方程式推導特徵向量矩陣（eigenvector matrices），如此一來，求解的計算量很大，模擬的過程將會很耗時，而且未來模式功能需要擴充時，如增加處理濃度傳輸

(18)

3

方程式或動床泥砂傳輸方程式等，便會相當不容易。

Jameson(1984)發展局部極值消減法求解機翼上穿音速流場（transonic flow）的震波問題（shock wave problem），同時也使用4 階Runge-Kutta時階法改善收斂速度。局部極值消減法主要精神就是加入人工黏滯項來抑制不連續面所產生的震盪。李隆正(1988)也曾採用Jameson的局部極值消減法求解多孔性翼剖面在穿音速流場的問題。局部極值消減法後來被引用在水利工程方面，Gharangik and Chaudhry(1989) 在MacCormack模式中加入了消散項來模擬水躍的形成。Fennema and Chaudhry (1990) 加入人工黏滯項(artificial viscosity) 來減少MacCormack 的震盪現象，並且探討在二維矩形渠道的潰壩湧浪現象。Younus and Chaudhry(1994)在ADI模式中加入人工黏滯項，分別模擬束縮渠道的水躍以及圓形水躍。Chaudhry and Khan(1995)將人工黏滯項加入ADI 隱式法並和 MacCormack 顯式法一同模擬在 spur-dike附近的水流情況。Thomas molls & chaudhry (1995 )利用ADI 法並加上人工黏滯項，模擬超亞臨界流、束縮段流況、180度彎道流況和潰壩模擬。Meselhe et al(1997)發展自調性消散項，用以模擬穿臨

界流場的情況。Rahman and Chaudhry(1998)加入人工黏滯項模擬潰壩，

並且調整格網討論適應性。連和政等(2007)採用局部極值消減法求解具源項之一維淺水波方程式，以不規則底床坡度變化之流場為例驗證模式的準確度。連和政等(2008) 以Jameson et al. (1981)所提出的自調式人工黏滯法求解急變流場不連續面之問題，以及採用4階Runge-Kutta改善程式之收斂速度。模擬假設的一維複雜流場案例作為模式測試。連和政等(2009) 採用 Qy模式直接求解流量與水位的變數，數值方法採用局部極值消減法，

(19)

使得在不規則斷面之數值模擬，仍然可以得到良好的數值守恆性。局部極值消減法相較於總變量消減法具有計算量較小的優勢，配合4階 Runge-Kutta時階法可提升數值穩定度，本研究延續連和政博士等人之研究，採用局部極值消減法配合有限體積法與4階Runge-Kutta時階法求解二維淺水波方程式，模擬不連續面流場之情況。

1.3 研究方法與流程

首先著手於二維水深平均水理模式之應用的文獻蒐集與整理。理論基礎主要利用水深平均之定義及邊界條件配合萊布尼茲積分推導出二維水深平均之控制方程式。本研究之數值方法採用有限體積法、具有抑制震盪作用的局部極值消減法以及提高數值穩定度的四階 Runge-Kutta 時階法。應用案例則選擇模擬潰壩案例、穿臨界流案例及超臨界流經過束縮段產生交坡之案例，分別是為了測試控制方程式中的慣性項、重力項及摩擦項對模擬精度之影響，並且挑選超臨界流經過束縮段產生交波之案例，作為消散項係數之敏感度分析，也針對不同消散項係數所耗費的 CPU 作討論。同時將模式初步應用於東埔蚋溪河道模擬，驗證模式之可行性。結論的部分，針對模擬結果作分析以及誤差原因的推測，並且提出需要改善的事項或未來的發展方向作為建議。

(20)

5 文獻回顧資料蒐集水理模式資料蒐集、整合理論基礎二維水深平均控制方程式公式推導數值方法 1. 有限體積法 2. 局部極值消減法 3. 四階 Runge-Kutta 時階法案例模擬及參數分析 1. 潰壩案例(慣性項) 2. 穿臨界流場案例(重力項) 3. 超臨界流在束縮河段產生交波案例(摩擦項) 4. 針對消散項係數作敏感度分析 5. 東埔蚋溪河道模擬結論 1. 模式在模擬中的精確度與穩定度 2. 敏感度分析的結果 3. 未來建議

(21)

1.4 本文結構

第一章導論: 說明研究的動機、目的、研究流程與回顧前人的貢獻。第二章理論基礎: 推導二維水深平均控制方程式，並將其守恆形式（Conservation form）寫成矩陣型式以利求解。第三章數值方法: 數值模式採用有限體積法做差分;採用局部極值消減法，藉此抑制不連續面上所產生的震盪。採用 Runge-Kutta 時階法改善收斂速度並提高穩定度。第四章數值計算結果及討論: 根據所模擬的三個模擬案例以及現地案例作探討與比較。針對消散項係數之指數，以交波案例作敏感度分析，討論出一個適當的係數範圍。第五章結論與建議: 將本文各章節所得結論作一綜合說明，提出更進一步之研究建議，以利後續研究發展。

(22)

7

第二章理論基礎

2.1 基本控制方程式

對於不可壓縮流，三維的連續方程式和動量方程式如下 (1) 連續方程式 0 u v w x y z ∂ ₊∂ ₊∂ ₌ ∂ ∂ ∂ _(2.1) (2) 動量方程式 ( ) xx yx zx x u u u u p u v w g t x y z x x y z τ τ τ ρ ∂ + ∂ + ∂ + ∂ =ρ −∂ +∂ +∂ +∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ _(2.2) ( v u v v v w v) g_y p xy yy zy t x y z y x y z τ τ τ ρ ∂ + ∂ + ∂ + ∂ =ρ −∂ +∂ +∂ +∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ _(2.3) ( ) xz yz zz z w w w w p u v w g t x y z y x y z τ τ τ ρ ∂ + ∂ + ∂ + ∂ =ρ −∂ +∂ +∂ +∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ _(2.4) x

、

y

、

z =卡式座標系統之座標軸。t =時間。u

、

v

、

w=水流速度於 x

，

y

，

z 方向之分量。ρg_x，ρg_y，ρg_z=單位質量重在 x

、

y

、

z 方向之分 量。P=壓力。ρ=密度。 xx τ 、τ_yx、τ_zx、τ_xy、τ_yy、τ_zy、τ_xz、τ_yz、τ_zz=應力張量，其定義如下所示： 2 ( ) xx xx u x τ = ρε ∂ ∂ (2.5) 2 ( ) yy yy u y τ = ρε ∂ ∂ _(2.6)

(23)

2 ( ) zz zz w z τ = ρε ∂ ∂ (2.7) ( ) xy yx xy u v y x τ =τ =ρε ∂ +∂ ∂ ∂ _(2.8) ( ) xz zx xz u w z x τ =τ =ρε ∂ +∂ ∂ _∂_(2.9) ( ) yz zy yz v w z y τ =τ =ρε ∂ +∂ ∂ ∂ _(2.10) 其中: ε_xx、ε_yy、ε_zz.εxy、εxz、εyz、εyx、εzx、εzy為各方向之紊流黏性擴散係數。

2.2 水深二維平均模式

(1)基本假設對一般天然河川模式而言，通常其寬度遠大於水深，平面的地形地貌改變大，則平面上水理特性的改變將遠大於深度上的改變，所以可將深度方向的水理特性忽略。根據以上特性，假設其水深方向變化遠小於水平方向，並忽略風力、科氏力之影響，利用水深方向之積分，可將三維之控制方程式簡化成二維水深平均方程式(2D-depth- averaged model)。基本假設如下: 1.流場水平尺度遠大於垂直水深尺度： 由於水平方向各種物理變化量均遠大於水深方向的各種物理變化量，因此可忽略水深方向之加速度及對流加速度與水深方向速度分量所產生之剪應力。

(24)

9 2.忽略流場之科氏力( Coriolis force )及風力的影響： 因地球各部的線性速率不同，地球自轉造成氣流偏向的力稱科氏力。科氏力影響不大，僅於長時間和長距離運動有關聯。 3.壓力為靜水壓分佈： 假設渠道中水流之垂向流線曲率很小而忽略其垂直加速度，因此水深方向速度梯度為零，以及水面坡度甚小，可忽略垂向加速度，因此靜水壓分佈的假設成立。 4.假設渠道無側向流量流入及底床無沖刷之影響: 假設渠道無任何側向流量包括：支流、降雨、入滲、蒸發等等。假設底床為一固定且無法刷深的河道。對任意物理量V 取其水深平均之定義： ( , , ) ( , ) 1 ( , , ) ( , , , ) b H x y t Z x y f x y t f x y z t dz h =

_∫

(2.11) 式中， f 為V 之水深平均值，Z_b為底床高度，H 為水面高程，h 為水深，h = H - Z_b。 (2)邊界條件: 進行水深平均時，需引入底床及水面運動邊界條件，水面運動邊界條件: DFs H _u H _v H _w 0 Dt t x y ∂ ∂ ∂ = + + − = ∂ ∂ ∂ _(2.12) s F 為水面運動方程 ( , , , ) ( , , ) 0 s F x y z t = −z H x y t = _(2.13) 底床運動邊界條件:

(25)

0 b b b b DF Z Z Z u v w Dt t x y ∂ ∂ ∂ = + + − = ∂ ∂ ∂ ， 0 b Z t ∂ ₌ ∂ _(2.14) b F 為底床運動方程式， F x y z t_b( , , , )= −z Z x y_b( , ) 0= _(2.15) 將基本假設帶入上述的邊界條件，改寫之後: 連續方程式： h ( )hu ( )hv 0 t x y ∂ ₊∂ ₊∂ ₌ ∂ ∂ ∂ _(2.16) x 方向動量方程式： 2 2 0 ( ) ( ) ₂ ( ) ( _x _fx) gh hu hu huv gh S S t x y ∂ + ∂ ₊ ₊∂ ₌ ₋ ∂ ∂ ∂ _(2.17) y 方向動量方程式： 2 2 0 ( ) ( ) ( ) ₂ ( _y _fy) gh hv hv huv gh S S t x y ∂ + ∂ ₊∂ ₊ ₌ ₋ ∂ ∂ ∂ _(2.18) 其中：h 為水深；u、v 為 x

、

y 方向平均速度；S_ox與S_oy為x

、

y 方 向的底床坡降；S_fx與S_fy為x

、

y 方向摩擦坡降；g 為重力加速度。 (3)二維水深平均模式之控制方程式守恆形式寫成矩陣格式如下 S y G x F t U = ∂ ∂ + ∂ ∂ + ∂ ∂ (2.19) ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = vh uh h U ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ + = uvh gh h u uh F 2 2 2 1 ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ + = 2 2 2 1 gh h v uvh vh G ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ − − = ) ( ) ( 0 fy oy fx ox S S gh S S gh S

(26)

11 式中，h= 水深； u

＝

x 方向之流速； v

＝

y 方向之流速； x z gh S b ox _∂ ∂ − = ＝x 方向之底床坡降； y z gh S b oy _∂ ∂ − = ＝y 方向之底床坡降； ₄_/₃ 2 2 2 h n v u u Sfx + = ＝x 方向之摩擦坡降（friction slope）； 3 / 4 2 2 2 h n v u v S_fy = + ＝y 方向之摩擦坡降（friction slope）；n=曼寧 係數。

(27)

第三章數值方法

3.1 有限體積法

為了使差分方程式能盡量滿足物理的守恆律，本研究採用有限體

積法(finite volume method) 。

由第二章式(2.19)得出的二維水深平均模式之控制方程式如下: S y G x F t U = ∂ ∂ + ∂ ∂ + ∂ ∂ (2.19) 將式(2.19)作控制體積積分，並對含 F、G 的面積分應用於散度定理 (Divevgence Theory)，化為沿其邊界的線積分： ( , ) A B A U dA F G ndB SdA t ∂ ₊ _• ₌ ∂

∫

r

∫

(3.1) 式中，

nr

=各邊的單位法向量。展開之後寫成: 4 1 ( )_k k U A F y G x S A t = ∂ ₊ _{Δ − Δ} _{= ⋅} ∂

∑

(3.2) 對空間之離散化可表為: 4 , , , 1 ( ) i j k i j i j k U A F y G x S A t = ∂ + Δ − Δ = ∂

∑

_(3.3) 最後，對時間之離散化可表為: 4 1 , , , , , 1 ( n n ) / ( n n ) n i j i j i j k i j i j k A U + U t F y G x S A = − Δ +

∑

Δ − Δ = (3.4)

3.2 局部極值消減法

Jameson（1981）的局部極值消減法，其具有最大值不會再增加以及最小值不會再減少的特性，藉此抑制數值震盪。Jameson 應用局部

(28)

13 極值削減法於機翼上穿音速流場的震波問題，我們採用局部極值消減法求解二維水理模式。由於實際地型崎嶇複雜，時常有不連續面流場發生，不連續面的流場容易產生數值震盪，震盪會引起非線性的不穩定，可能導致程式發散，因此本研究採用局部極值消減法抑制數值的震盪。其原理可從偏微分方程式看起，一般偏微分方程式的形式可表為 : ( ) ( ) 0 u f u g u t x y ∂ ₊ ∂ ₊ ∂ ₌ ∂ ∂ ∂ _(3.5) 令uj ，u_k為格網點 j

、

k 的 x 方向的流速，j

、

k 其為鄰近的兩個格網點 將(3.4)使用泰勒展開式將各項展開並整理成(3.5)式 j jk k k du c u dt =

∑

_(3.6) 再將(3.5)式寫成差分型式: j ( ) jk k j k j du c u u dt ≠ =

∑

− (3.7) (3.7)中，係數c_jk的總合需為 0，即 _jk 0 k c =

∑

且(3.7)式中的係數c_jk需為正值 0, jk c ≥ k≠ j 在(3.7)中，若 uj為最大值: 則uk-uj≦0 ， j 0 du dt ≤ ，可知最大值不會再增加。在(3.7)中，若 uj為最小值:

(29)

則uk-uj≧0， duj 0 dt ≥ ，可知最小值不會再減少。經由以上的原理得知，加入人工消散項可以去除過大或是過小的值，藉此來抑制數值的震盪。因此我們在二維水深平均模式中，針對x 和 y 方向各加入人工消散項: 4 1 , , , , , 1 ( n n ) / ( n n ) n n n i j i j i j k i j i j x y k A U + U t F y G x S A D D = − Δ +

∑

Δ − Δ = + + (3.8) (3.8)式中 x D 為 x 方向的人工消散項， Dy 為 y 方向的人工消散項。 x D =(di+1/ 2,j −di−1/ 2,j)_， Dy= (di j, 1/ 2+ −di j, 1/ 2− )_(3.9) U t x d_i _j [ _i _j _i _j _i(4) _j( _i ₃_/₂_,_j 2 _i₁_/₂_,_j _i₁_/₂_,_j)] , 2 / 1 , 2 / 1 ) 2 ( , 2 / 1 * , 2 / 1 + + + + + − + _Δ Δ − Δ − Δ −Δ Δ = ε ε (3.10) (2) (4) 1/ 2, *[ 1/ 2, 1/ 2, 1/ 2, ( 1/ 2, 2 1/ 2, 3/ 2, )] i j i j i j i j i j i j i j x d U t ε ε − − − − + − − Δ = Δ − Δ − Δ − Δ Δ _(3.11) U t x d_i _j [ _i_j _i_j _i(4_j) ( _i_,_j ₃_/₂ 2 _i_,_j₁_/₂ _i_,_j₁_/₂)] 2 / 1 , 2 / 1 , ) 2 ( 2 / 1 , * 2 / 1 , + _Δ + Δ + − + Δ + − Δ + −Δ − Δ = ε ε (3.12) (2) (4) , 1/ 2 *[ , 1/ 2 , 1/ 2 , 1/ 2( , 1/ 2 2 , 1/ 2 , 3/ 2)] i j i j i j i j i j i j i j x d U t ε ε − − − − + − − Δ = Δ − Δ − Δ − Δ Δ _(3.13) j i j i j i1/2,U =U 1, −U, Δ₊ ₊ (3.14) 1/ 2, 1, 2, i− jU Ui− j Ui− j Δ = − (3.15) , 1/ 2 , 1 , i j+ U Ui j+ Ui j Δ = − (3.16) , 1/ 2 , 1 , 2 i j− U Ui j− Ui j− Δ = − (3.17) (2) (2) 1/ 2, max( 1, , , ) i j k i j i j β β ε₊ = υ ₊ υ (3.18)

(30)

15 (2) (2) 1/ 2, max( 1, , 2, ) i j k i j i j β β ε− = υ − υ − _(3.19) (2) (2) , 1/ 2 max( , 1, , ) i j k i j i j β β ε + = υ + υ _(3.20) (2) (2) , 1/ 2 max( , 1, , 2) i j k i j i j β β ε − = υ − υ − _(3.21) ) , 0 max( (2) , 2 / 1 ) 4 ( ) 4 ( , 2 / 1 j i j i+ = k ε+ ε _(3.22) (4) (4) (2) 1/ 2, max(0, 1/ 2, ) i j k i j ε− = ε− (3.23) (4) (4) (2) , 1/ 2 max(0, , 1/ 2) i j k i j ε + = ε + (3.24) (4) (4) (2) , 1/ 2 max(0, , 1/ 2) i j k i j ε − = ε − (3.25) j i j i j i j i j i j i j i h h h h h h , 1 , , 1 , 1 , , 1 , 2 2 − + − + + + + − = υ (3.26) 2, 1, , 1, 2, 1, , 2 2 i j i j i j i j i j i j i j h h h h h h υ + + + + + − + = + + (3.27) , 1, 2, 1, , 1, 2, 2 2 i j i j i j i j i j i j i j h h h h h h υ − − − − − − + = + + (3.28) , 2 , 1 , , 1 , 2 , 1 , 2 2 i j i j i j i j i j i j i j h h h h h h υ + + + + + − + = + + (3.29) , , 1 , 2 , 1 , , 1 , 2 2 2 i j i j i j i j i j i j i j h h h h h h υ − − − − − − + = + + (3.30) * t Δ ＝CFL＝ x t gh V Δ Δ + ) ( ＝1 之時階。本研究採用_k(2)＝1.0，_k(4)＝0.001。

(31)

3.3 Runge-Kutta 多時階法

本研究在空間上的差分採用顯示法，顯示法雖有程式容易撰寫的優點，但此法是屬於線性有條件的穩定，因此必須滿足 CFL（courant number）≦1，但這樣一來就會影響程式的收斂速度。為了改善顯示法的收斂速度，本研究也採用 Jameson （1984）的 4 階 Runge-Kutta 時階法，則式（3.7）可改寫成： 1 , , ,

, (

n n

) /

( ) 0

i j i j i j

Ai j U

+

−

U

Δ +

t R U

=

(3.8) 式中 4 , , , 1

( )

(

n n

)

n n n i j k x y i j i j k

R U

F y G x

D

S A

=

∑

Δ −

Δ

−

式（3.8）可寫成： ) ( , ) 0 ( , n j i j i U U = ( (0)) , , 1 ) 0 ( , ) 1 ( ,j ij ij ij i U tR U U = −α Δ ( (1)) , , 2 ) 0 ( , ) 2 ( ,j ij ij ij i U tR U U = −α Δ ( ) ) 2 ( , , 3 ) 0 ( , ) 3 ( ,j ij ij ij i U tR U U = −

α

Δ ( ) ) 3 ( , , 4 ) 0 ( , ) 4 ( .j i j i j i j i U tR U U = −

α

Δ 本研究採用α₁=1/4，α₂ =1/3，α₃ =1/2 ，α₄ =1

3.4 邊界條件設定

在自然界的河川中，邊界可區分為固體邊界 (陸地邊界) 與開放邊界 (水域邊界) 兩類。

(32)

17 (1) 固體邊界: 固體邊界 Γs ，可分成滑動 (slip) 與非滑動 (no-slip) 邊界條件。 若流體黏滯性可忽略時(ν = 0 ) ，則可設定為滑動條件 (slip condition)；若有考慮流體黏滯性的影響時必須設定為非滑動條件 (non-slip condition)。定義如下所示: 0, x 0, 0 y n n q h q N N ∂ ∂ ₌ ₌ ₌ ∂ ∂ ， on Γs (3.31) 0, _x 0, _y 0 n h q q N ∂ ₌ ₌ ₌ ∂ ， on Γs (3.32) 其中：Nn為法線( normal )方向。如圖 3-1 顯示，定義 M 為 x-方向網格編號，差分展開後格式如下所示 slip: ( 1) ( ) ( 1) 0 ( 1) ( ) x y y h M h M q M q M q M + = + = + = on Γs (3.33) nonslip: ( 1) ( ) ( 1) 0 ( 1) 0 x y h M h M q M q M + = + = + = on Γs (3.34) 本研究採取二階精度，邊界需再延伸一點，其格式如下:

(33)

滑動邊界: ( 2) ( 1) ( 2) 0 ( 2) ( 1) h M h M qx M qy M qy M + = − + = + = − on Γs (3.35) 非滑動邊界: ( 2) ( 1) ( 2) 0 ( 2) 0 h M h M qx M qy M + = − + = + = on Γs (3.36) 邊界為固體邊界時，在邊界點外需設”虛擬元素(ghost)”，作為鏡射 (mirror)元素之用。當水位在邊界處內、外有著相同值，其流量則會差一負號，屬於反射(reflective)邊界。對於不規則邊界條件，由實際地形模式中水位高程函數反應出來，因為屬於全域計算，所以其不規則邊界即可由固體邊界反應進計算域中，不需額外再給定。 (2)開放邊界 開放邊界條件，分為入流與出流條件。在入流邊界條件，即h=hin， Q=qin， on Γin。上游入流與下游出流邊界屬於開放邊界，可分成穿

透邊界(Transmissive boundary)與反射邊界(Reflective boundary) ，定義如下所示：穿透邊界: 0 0 0 n x n y n h N q N q N ∂ ₌ ∂ ∂ = ∂ ∂ = ∂ (3.37)

(34)

19 反射邊界: 0 0 0 n x y h N q q ∂ ₌ ∂ = = _(3.38) 本研究採取二階精度，邊界需再延伸一點，其格式如下: 穿透邊界: ( 1) ( ) ( 2) ( 1) ( 1) ( ) ( 2) ( 1) ( 1) 0 ( 2) 0 x x x x y y h M h M h M h M q M q M q M q M q M q M + = + = − + = + = − + = + = _(3.39) 反射邊界: ( 1) ( ) ( 2) ( 1) ( 1) ( ) ( 2) ( 1) ( 1) 0 ( 2) 0 x x x x y y h M h M h M h M q M q M q M q M q M q M + = + = − + = − + = − − + = + = _(3.40) 當M=0 為上游入流邊界，其內、外點，如圖 3-2 所示。當水位在邊界處內、外有著相同值，其流量則會差一負號，屬於反射(Reflective) 邊界。

(35)

第四章模式驗證分析與東埔蚋溪河道模擬

為了驗證模式在各種流場的精確度及穩定性，本研究選擇三種不同的流況來模擬分析，如下表所示: 案例流況特性驗證項潰壩模擬水平底床不考慮摩擦慣性項 2 u h、_{v h}2 穿臨界流流經底床突起考慮底床坡降不考慮摩擦重力項 0 x

S

、

S

_{0 y} 河道束縮段產生交波水平底床考慮摩擦情況摩擦項 fx

S

_、

S

_fy 在模擬的過程中，抑制數值震盪的消散項是關鍵，幅度過大的震盪可能導致程式發散，因此消散項的強弱對數值的穩定性有重要的影響。控制消散項之係數是局部極值消減法公式中人工黏滯係數

υ

的指數，以

β

值代稱。在模擬的過程中，每一個格網點的

β

值都是固定的。在同樣

β

值的條件下，將模擬結果分別與解析解或實驗值作比較。

4.1 潰壩的案例

在潰壩發生時，水流會形成一個極端的不連續面，在不考慮底床坡降與摩擦的情形下，觀察水流流經不連續面時，模式的精確度與穩

(36)

21 定度。本案例目的為測試模式中的慣性項，模擬結果與 Stoker(1957) 之解析解比較。實驗條件為一水平無摩擦之理想潰壩，渠道長度為3000 m，渠道寬度為1250 m，壩址位置假設在渠道中央（即距上游 1500 m 處），上游水位為20 m，下游水位為 1 m，如圖 4-1 所示。當 t > 0 時，瞬時移開閘門模式即進行潰壩演算。觀察模擬結果在潰壩發生後不同時刻的水位與流速。由 Stoker 的解析解可知臨界流況會發生在下游與上游水深比為 0.138 時(即H H2/ 1=0.138)，當H H2/ 1<0.138會發生超臨界流潰壩湧浪，而H H₂/ ₁ >0.138則會發生亞臨界流潰壩湧浪。本案例下游與上游水深比為H H₂/ ₁=0.05，故會發生超臨界流之潰壩湧浪。

β

值設定中，分別以

β

=0.2、

β

=0.4、

β

=0.6 來做測試，取其結果較佳的情況來做分析。採用301x26 個計算格網點，一共 79 次迴圈，

實際模擬時間為40 秒，CPU(Intel Daul Core 2.26GHz) 耗時為 2.16

秒，模擬時間內上下游邊界均未受到影響，因此在模擬過程中上游邊界入流量保持為零且下游水位保持為1 m。模擬結果討論不同

β

值對數值震盪的抑制能力有所差異，以潰壩發生第 20 秒的水位比較案例來討論，圖4-2 為

β

=0.2 時，由於消散項效應強，抑制震盪的能力也強，因此未出現數值振盪的情況，在負湧浪的上游出現最大誤差，約 5.5%；圖 4-3 為

β

=0.4 的情況，數值未出現震盪，而其最大誤差發生在負湧浪的上游，約4%；圖 4-4 為β =0.6 的情況，在正湧浪部分出現了明顯的數值振盪。根據三種

β

值的模擬結果，

β

=0.4 的情況最符合解析解，因此以

β

=0.4 進行模擬。接著以潰壩發生第 40 秒的水位比較案例來討論，圖 4-7 為

β

=0.2 的情況，最大誤

(37)

差出現在負湧浪，誤差約為5%;圖 4-8 為

β

=0.4 情況，最大誤差出現在負湧浪，約為4%;圖 4-9 為

β

=0.6 情況，在正湧浪的部分出現了明顯的震盪，因此

β

=0.6 的時候，模擬結果較不穩定，綜合以上結果比較，

β

=0.4 的情況最符合解析解。圖 4-1 是潰壩的初始狀態，壩址位於渠道中央，觀察潰壩開始之後，針對渠道中心線的水位及流速做討論。當閘門打開後，由於上下游水位的差異而形成超臨界流湧浪，並且產生兩個震波，分別是一個向下游傳遞的正波(positive wave)及一個向上游傳遞的負波(negative wave)。模擬結果在正負湧浪的上下游產生誤差。圖 4-3 為第 20 秒的渠道中心線水位與解析解比較，負湧浪部分所產生的最大誤差約4%; 正湧浪部分的最大誤差為3%。圖 4-8 第 40 秒的渠道中心線水位與解析解比較中，負湧浪部分產生的最大誤差為4%;正湧浪部分出現的最大誤差為1%。模擬結果的趨勢比解析解來的和緩。流速的部分，產生誤差的地方亦是正負湧浪的區域。圖 4-10 第 20 秒渠道中心線流速模擬結果與解析解比較中，最大誤差約為 5.2%; 圖4-11 第 40 秒渠道中心線流速模擬結果與解析解比較中，最大誤差約為 4%。在震波尚未到達的區域，其值都與解析解重合，維持在初始狀態的部分。

4.2 穿臨界流遇到底床突起的案例

在現地模擬中，由於地形及流況多變，時常遇到亞臨界流場與超臨界流場相接而產生的穿臨界流場，所以模式將模擬穿臨界流場的情況。由於本案例考慮到底床的突起，因此案例目的為測試模式源項中跟底床坡降有關的重力項(S0 x.S0 y)。底床突起地形尺寸引用 Zhou et al. (2001)的實驗，突起地形尺寸如下

(38)

23 2 0 if < 8 m ( ) 0.2 0.05( 10) if 8 x 12 m 0 if x > 12 m b x Z x x ⎧ ⎪ =_⎨ − − ≤ ≤ ⎪ ⎩ 圖4-12 為底床突起示意圖，實驗條件為無摩擦的矩形渠道，長 25m、寬1m，分別模擬未發生震波及發生震波兩種情況。 (1) 模擬穩定水流流經突起且未發生震波的情形 單寬入流量上游初始水位下游水位 1.53m2/s 0.6m 0.4m 格網點為 201x3，一共 3000 次迴圈，總模擬時間為 62.69 秒， CPU 耗時為 5.56 秒。觀察模式在穿臨界流場中的穩定度與精確度，模擬結果與解析解比較水位、福祿數及剩餘值變化。模擬結果討論不同β 值對數值震盪的抑制能力有所差異，在穩定水流流經突起且未發生震波的案例，分別以β =0.2， β =0.4， β =0.6 三種情況來討論，圖 4-13 為當β =0.2 的情況，在不連續面靠近上下游的部分，皆可看出數值震盪。圖4-14 為β =0.4 的情況，數值未出現明顯震盪，大致符合解析解。圖 4-15 為β =0.6 的情況，不連續面的上下游部分可看出些微數值震盪。根據以上三種模擬結果，以β =0.4 的情況最符合解析解。因此本案例以β =0.4 進行模擬。圖 4-16 至圖 4-21 是穿臨界流在各時刻的瞬時水位變化。圖 4-14 是模擬水位與解析解的比較，大致與解析解吻合。圖 4-23 是福祿數變化的趨勢與比較，上游為亞臨界流，經過底床突起之後，福祿數急速上升，在底床突起的部分發生穿臨界流場，水流流過突起的底床後，

(39)

則變成超臨界流。模擬數值跟解析解大致吻合。圖 4-24 是 3000 次迴圈的剩餘值變化，殘餘值公式計算可表為: 1 , ( ) Re N i j i R U t s N = Δ =

∑

(3.41) 在 0 至 1000 次迴圈的區域，剩餘值逐漸下降;超過 1000 次之後，剩餘值些微上升。第 3000 次迴圈的剩餘值控制在 1.00E-4 以下，是合理的範圍。 (2) 模擬穩定水流流經突起且發生震波的情形 單寬入流量上游初始水位下游水位 0.18m2/s 0.25m 0.33m 格網點為 201，一共 4000 次迴圈，總模擬時間為 154.3 秒，CPU 耗時為7.5 秒。觀察模式在發生震波的穿臨界流場中的情況，與解析解比較水深、福祿數及剩餘值變化。模擬結果討論不同β 值對數值震盪的抑制能力有所差異，在穩定水流流經突起且未發生震波的案例，分別以β =0.2，β =0.4，β =0.6 三種情況來討論，圖 4-25 為β =0.2 的情況，明顯的看出在不連續面在 x=7m 的部份開始出現明顯的震盪以及 x=12m 至 x=13m 的部分也有明顯震盪；其模擬出的最低水位約為0.22m，遠高於解析解的 0.15m。圖 4-26 為 β =0.4 的情況，在不連續面部分的震盪較之前微小許多，誤差主要出現在最低水位的部分，誤差約 17%，其餘部份大致符合解析解。圖 4-27 為β =0.6 的情況，不連續面的水位在上游的部分，約是 x=6m 的地方開始出現微幅震盪，誤差大約在 0.2%左右，最低水位的誤差約為18%。以上三種情況相比，以β =0.4 的情況最符合解析解。因此本案例以β =0.4 進行模擬。

(40)

25 圖 4-26 是模擬水位與解析解的比較，除了最低水位的部分出現誤差，其餘部分都跟解析解大致吻合。模擬結果的最低水位是0.181m，而且解析解的最低水位是 0.15m，誤差約在 17%。圖 4-28 至圖 4-31 是穿臨界流在各時刻的瞬時水位變化。圖 4-32 是模擬結果的福祿數與解析解的福祿數比較，本案例是從亞臨界流經過底床抬升之後變成超臨界流，又從超臨界流回到亞臨界流的狀態，所以經歷了兩次的穿臨界流場。模擬結果不連續面的部分發生誤差，模擬結果在x=11.4m 出現最大的福祿數，而解析解的最大福祿數則是發生在x=11.6m，誤差約為 1.7%；而解析解的最大福祿數為2.63，模擬結果的最大福祿數是 2.30，誤差約為 12.5%。圖 4-33 是 4000 次迴圈的剩餘值變化，在第 1500 次之後，剩餘值都控制在 1.00E-4 以下，是合理的範圍。

4.3 超臨界水流遇到河道束縮產生交波的案例

當超臨界流流經束縮段渠道時，震波受邊牆的反射而交會，稱之為交波(cross wave)。在工程實務上，若能將邊牆的走向做適度的布置，便可以利用波的交互干擾特性來消除不需要的水面震波。本案例引用 Coles and Shintaku(1943)的實驗案例，一個由寬變窄的束縮渠道，渠道斷面為矩形，上游渠寬為0.629 m，下游渠寬為0.314 m，上游水深固定為0.0314m，上游單寬流量0.0717m2/s， Chzey C=84.1m1/2/s，模擬採用41x21個計算格點，一共2000次迴圈，實際模擬時間為30.67 秒， CPU耗時6.67秒。渠道示意圖如圖4-34所示。由於考慮到底床摩擦，因此案例目的為測試模式中源項的摩擦項(S_fx、S_fy)。觀察模式在超臨界水流因河

(41)

道束縮而產生交波時，模擬邊牆水位與實驗值的比較並觀察剩餘值變化。單寬入流量上游初始水位摩擦係數 0.0717m2/s 0.0314m Chzey C=84.1m1/2/s 模擬結果討論不同β 值對數值震盪的抑制能力有所差異，分別以β =0.2， β =0.4， β =0.6三種情況來討論，圖4-35為β =0.2的結果，是三種情況中跟實驗值誤差最大的，由於消散項效應較大，將較大幅度的震盪去除，最大誤差出現在第一個交波，約為12%；圖4-36為β =0.4的模擬結果，最大誤差出現在第二個交波，約為9%；圖4-37為β =0.6的結果，最大誤差出現在第二個交波後期，約為11%。根據以上三種模擬情況相比，以β =0.4的情況最符合實驗值。因此本案例以β =0.4進行模擬。圖 4-38 到圖 4-42 為束縮段產生交波的模擬情況。圖 4-36 是束縮段邊牆水位模擬結果與實驗值的比較，模擬結果的趨勢大致與實驗值符。模擬結果在第一個交波之後，水位上升至 0.0429m 的地方然後呈現水平狀；實驗值在第一個交波的地方，水位急速上升至 0.047m 形成一個高峰，之後的幾個實驗值微幅滑落至0.0435m，這部分模擬值的最大誤差約為 9%。另一個誤差部分是在第二個交波時，模擬結果在 x=1.783m 的地方有一個高峰的水位為 0.0709m，實驗值的第二次抬升不明顯。高峰之後，水位先是些微下降然後上升至 0.0729m，接著再往下降。實驗值在第二交波後期，最後的兩個實驗值都大於模擬結果，最大的誤差約在6%。本案例中，一共發生三次水位的抬升。第一個交波發生後，由於

(42)

27 水流進入束縮段，導致流量變化而使得水位發生第一次抬升，第二個交波發生在束縮段結束進入矩形渠道的部分，估計第二次水位抬升是由於渠道變窄所致，而後水位些微下降又開始了第三次抬升，可能是受到了前面交波的影響所致。實驗值的第二次抬升較不明顯，可能是實際操作時量測不易所導致。

本模式結果分別與 1995 年的 Molls and Chaudhry 以及 1992 年的 Bhallamudi and Chaudhry 的模擬結果做比較。圖 4-46 是三個模式的模

擬水位比較，圖 4-47 則是三個模式與實驗值的誤差比較，，

Bhallamudi and Chaudhry 的結果在第一個交波的模擬值跟實驗值接近，而第二個交波後期的結果起伏震盪大，誤差較明顯；Molls and Chaudhry 則是在第二個交波前期的模擬值與實驗值偏離較多，第一個交波的模擬結果與實驗值接近。交波的案例中，發生了三次水位的抬升。第二次抬升的情況，在本模式以及另外兩個模式的模擬情況中都有出現，而實驗值的第二次抬升卻不明顯，估計可能是實際操作中有現實因素未被考慮進模式中，因此與實驗結果有出入；亦可能是實驗中量測水位不易，量測時間與水位很難準確掌握。圖 4-43 是模擬的水位等高線圖，圖 4-44 是實驗量測出來的等高線圖，由模擬值看來，呈現出兩個明顯的交波型態。兩圖比較中，在 X=1.3m 以前，模擬結果跟實驗值接近。在束縮段的後半段，第二個交波之後，實驗值靠近邊牆的位置水位較模擬值來的高。圖 4-45 為

1995 年的 Molls and Chaudhry 所模擬之水位等高線圖，其第二個交波的型態不甚明顯，第二交波後期的邊牆水位值較接近實驗值。圖 4-48

為福祿數的變化情形，在約X=0.3m 時，渠道開始束縮，流速受到影

(43)

渠道在X=1.79m 的位置停止束縮，福祿數約在 1.5m處就開始往下降，估計是由於接近束縮口的關係，流速受到影響。渠道停止束縮之後，福祿數降為 2.5 左右。圖 4-49 所示，約在第 1170 次之後，剩餘值相當趨近於零，是理想的狀態。本案例同時也模擬了渠道中心線水位變化，實驗值的水位發生一次抬升，中心線水位發生了一次交波。本模式模擬的結果，模擬水位在交波尚未發生時，與實驗值相當接近，發生交波之後，模擬水位上升至 0.065m 至 0.07m 左右即停止上升，而實驗值則是上升至 0.08m 才往下降。

本模式結果分別與 1995 年的 Molls and Chaudhry 以及 1992 年的 Bhallamudi and Chaudhry 的模擬結果做比較。圖 4-49 為三個模式與實驗值的水位比較圖，1992 年與 1995 年在交波後期的水位明顯比實驗值來的高出許多，1992 年的模擬結果在交波尚未發生時即發生誤差，發生交波之後，水位與實驗值比較則有明顯的誤差；1995 年的模擬結果，亦是在交波尚未發生時即產生誤差，交波發生後誤差變大。圖 4-50 是三個模式與實驗值的誤差比較，由圖可之，1992 年的模式在交波發生時誤差最大；本模式與 1995 年的誤差趨勢接近，本模式在交波發生前期的誤差較小，1995 年模式在交波發生後期較接近實驗值，整理來看，本模式的誤差相對較小。

4.4 β 值敏感度分析

本研究中，抑制數值震盪的消散項具有相當重要的影響。幅度過大的震盪可能導致程式發散，因此消散項會影響數值的穩定度。與抑制震盪能力有關的人工黏滯項υ ，其指數β 值的大小在0~1 之間，β 值越小，消散項效應越強，抑制震盪的能力越強;反之β 值越大，消

(44)

29 散項效應越弱，抑制震盪的能力趨弱。β 值的合理範圍，以敏感度分析來決定。本研究所採用的指標有以下三種。 1. 效率係數 (Coefficient Of Efficiency，CE) 當 CE 值越趨近於 1 時，表示模擬結果與實際觀測之資料越吻合，模擬之精確度越高。其關係表示式如下所示：

[

]

{

}

{

}

2 1 2 1 ( ) ( ) 1 ( ) n i n i EXP H SIM H CE EXP H EXP = = − = − ⎡ − ⎤ ⎣ ⎦

∑

(4.1) ( )

EXP H 為實驗值， SIM H( )為本程式模擬值，EXP實驗值之平均值。

2. 均方根誤差 (Root-Mean-Square error，RMSE) 樣本之各點位推估值與觀測點位值之平均方根誤差，當 RMSE 越趨近於0，表示模式推估值越符合觀測值，其定義如下所示：

[

]

{

}

0.5 2 1 ( ) ( ) n i EXP H SIM H RMSE N = ⎧ ⎫ − ⎨ ⎬ ⎩ ⎭ =

∑

(4.2) ( ) EXP H 為實驗值， SIM H( )為本程式模擬值。

3.最大值之相對差異值(Error for the Peak Value， EVP%)

當模擬最大值的相對差異值越小，即代表模擬結果與實際觀測資料越吻合，其定義如下所示:

EPV=(模擬最大值-實驗最大值)/實驗最大值＊100%

利用上述三種指標針對超臨界流流經束縮段產生交波這個案例做

(45)

來做敏感度分析。由於在β =0~0.15 的情況模擬結果發散，所以從β =0.2 開始討論。表 4-1 為敏感度分析結果，並將結果繪製成圖以便分析。由圖4-49 CE 效率系數比較圖來看，0.25 到 0.4 這個範圍的效率係數最高，大約在 0.948~0.953 之間，皆趨近於 1。由圖 4-51 RMSE 圖來看，0.3 到 0.4 這個範圍的均方根誤差都在 0.003，皆趨近於 0。圖 4-50 為最大值之相對差異值的比較，一開始是逐漸變小的，直到 0.4 後，差異值又開始變大，到 0.5 之後差異值開始逐漸變小。根據三種數據來看，符合效率係數高、均方根誤差小、最大相對誤差小的條件，綜合結果來看，當β =0.3~0.4 這個範圍時，效率係數(CE)最高，最大相對誤差(EVP%)較低，且均方誤差(RMSE)=0.003 皆趨於最小，所以β =0.3~0.4 應可視為最佳範圍。接著針對不同 β 值的模擬來討論 CPU 的耗時以及將將 CPU 耗時除以實際時間，CPU 以 Intel Daul Core 2.26GHz 來作測試。測試 β 由 0~1 的模擬時間，將其結果繪製成表 4-2 及圖 4-54、圖 4-55。由圖得知， CPU 耗時呈現不規則變化。經由敏感度分析所得出的 β 適合範

圍是0.3-0.4，這個區域之 CPU 耗時約莫是 6.57-6.63 之間，在整體比

較中，CPU 耗時是較低的。將 CPU 的耗時除以時間，當 CPU 耗時越

小，而實際模擬時間越大時，代表效能越好，所得出的值越小；反之，當 CPU 耗時越大，而實際模擬時間越小時，代表效能較弱，所得出的值越大。由圖 4-55 可知，圖型呈現不規則變化，依照敏感度分析所得出的 β=0.3~0.4 之間，計算出的值低於 0.21，在整體的曲線看來是較低的值，顯示在這個範圍中效能偏高。，因此 β=0.3~0.4 中，求解階度較高，CPU 耗時偏低且 CPU 的效能偏高。

(46)

31

4.5 現地案例:東埔蚋溪河道模擬

本章節將模式實際應用於東埔蚋溪之河道模擬，以下針對地理位置、模擬數值條件以及模擬結果分析來討論。【地理位置】：東埔蚋溪為濁水溪重要支流，主流發源於阿里山山脈之領頭山、鳳凰山及大鞍山。流域分佈如圖 4-52 主流全長 18.5km，本研究所模擬的河段是從延平橋至濁水溪匯流口，長約 2.56km，河床平均坡降約為1.136%，屬於急變流河川。【模擬數值條件】：圖 4-53 為延平橋至濁水溪匯流口河段之河道示意圖，模擬範圍是延平橋至濁水溪匯流口之河道內部區域。斷面資料參考東埔蚋溪治理規畫檢討報告(民 91)以及高含砂水流整治規畫條件檢討及計算模式研發應用總成果報告(民 97)，為桃芝颱風後第四河川局修復之堤防與疏濬後之量測斷面，糙度係數以曼寧 n=0.042 為常數進行設定，下游水位為濁水溪斷面100 的洪水位 135.05m。流量部分依照東埔蚋溪治理規劃檢討報告，採用 100 年重現期距的洪水量來模擬，給定的上游入流量為1950cms，上游河寬為 36.4m，單位寬度流量為 53.57 m2/s。由定量流進行河道模擬，格網點為189×15，實際模擬時間為 442 秒， CPU 耗時為 35.45 秒，模擬結果以水位、福祿數、殘餘值作分析。【模擬結果分析】：圖 4-53 為東埔蚋溪河道示意圖，模擬河段起點位於右下角，濁水溪匯流口位於左上角，河川方向是由右下往左上流動。圖 4-54 分別是左岸、右岸、渠道中心線之水位模擬圖。圖4-53 與圖 4-54 作比較，河道經歷了三個明顯的彎道，河段起點的河寬較窄，河道逐漸向左岸拓寬，右岸則是維持不變；因此左岸的水位低於右岸水位，接著

(47)

河道方向由左向右彎，由於向心力的緣故，左岸水位順著灣道逐漸升高，在彎道的部分為左岸水位高於右岸水位。過了第一個彎道之後，河道逐漸變窄，在束縮段的部分，水位明顯上升。束縮段過後，河道進入第二個彎道，由右向左彎，由於向心力的緣故，右岸水位順著彎道逐漸升高，左岸水位則是降低。過了第二個彎道之後，由於河道方向未有明顯改變，所以右岸與左岸的水位大致相同。接著進入第三個彎道，河道方向是由右往左彎，因此右岸水位抬升，左岸水位下降，之後河道方向跟寬度都未再改變，河道左岸與右岸的水位呈現大致相同之平衡的狀態。渠道中心線水位則是大致介於左岸與右岸的中間。根據圖 4-55 福祿數變化顯示，從延平橋開始水流就呈現超臨界流的狀態，河道遇到彎道的部分幾乎都呈現超臨界流。由福祿數變化可知穿臨界流頻繁發生，顯示此河段的流況複雜。圖 4-56 殘餘值的部分，在 3000 次迴圈之後，殘餘值大約低於

10

−5，相當趨近於零，數值呈現穩定狀態。根據以上模擬結果，本研究所採用之局部極值消減法亦能處理穿臨界流等複雜流況。

(48)

33

第五章

結論

本研究旨在發展二維水深平均模式進行二維水理模擬。模式發展: 1. 為了使模式盡量滿足物理的守恆律，採用有限體積法作差分。 2. 為了克服不連續面所造成的數值振盪而採用局部極值消減法 (LED)。 3. 使用 ROUNGE-KUUTA 多時階法來提高數值的穩定度。模擬分析與建議: 1. 根據三個驗證案例的模擬結果分析，經過計算值與解析解及實驗數據的比對，在水位、福祿數與殘餘值上都獲得不錯的結果，潰壩案例的最大誤差控制在6%以內，而與實驗的交波案例，其最大誤差控制在 9%以內。本模式之模擬結果與 1995 年的 Molls and

Chaudhry 以及 1992 年的 Bhallamudi and Chaudhry 的模擬結果比較，顯示本研究所發展的局部極值消減法在不連續面流場中具有不錯的求解精度。 2. 將超臨界流流經束縮段產生交波的案例作敏感度分析，跟消散項黏滯係數有關的

β

在0.3~0.4 之間的範圍，其效率係數、均方根誤差、最大值之相對差異值百分比的綜合考量最佳。將CPU 耗時與效能配合β 值作分析，可得出在 β 等於 0.3~0.4 的效能也是最佳範圍。 3. 在東埔蚋溪河道模擬中，由福祿數變化看出穿臨界流頻繁發生，顯示模擬河段流況複雜，根據模擬結果顯示本模式亦能處理穿臨界流之複雜流況。 4. 建議可以在模式中加入消散項的限制式，利用限制式來決定消散項抑制能力之強弱。在不連續面流場中，當數值震盪較劇烈時，

(49)

需要較大的抑制能力；反之，在震盪輕微甚至未發生震盪的區域，過大的抑制能力反而導致數值誤差。配合上限制式之後，模擬結果會更符合實際情況。 5. 在穿臨界流發生震波的案例中，水位在底床突起的位置出現了物理振盪。局部極值消減法會去除幅度較大的震盪，使得本該存在的物理振盪也有被抑制的趨勢。保留物理振盪且抑制數值振盪才是最接近實際情況的模擬，使用限制式或是其他方法使消散項辨別物理振盪與數值振盪，可作為後續研究的方向。

(50)

35

參考文獻

1. 李隆正，(1988)，「多孔性翼剖面在穿音速流場之數值模擬」，國立成功大學航空太空工程研究所，碩士論文。 2. 東埔蚋溪治理規劃檢討報告，(2002)，經濟部水利署水利規劃試驗所。 3. 古孟晃，(2003) ，「以總變量消減(TVD)法求解淺水波方程式」，逢甲大學土木及水利工程研究所碩士論文。 4. 吳銘順，(2006) ，「實際地形潰壩數值模擬」，逢甲大學水利工程研究所碩士論文。 5. 林宗賢，(2007)，「河川洪水對三重地區淹水之影響」，逢甲大學水利工程研究所碩士論文。 6. 高含砂水流整治規畫條件檢討及計算模式研發應用總成果報告， (2008)，經濟部水利署水利規劃試驗所。 7. 連和政、李隆正、張文鎰，(2008)，“自調式人工黏滯法應用於急變流場之研究”，第十七屆水利工程研討會，台中 8. 連和政、張文鎰、吳祥禎、張哲豪，(2009) ，「局部極值消減法應用於不規則斷面之流場模擬」，第十八屆水利工程研討會， p.H154-p.H160。 9. 連和政、張文鎰、吳祥禎、張哲豪，（2009），「以局部極值消減法求解具源項之淺水波方程式」，九十八年電子計算機於土木水利工

(51)

程應用研討會，p.425-430

10. Fennema R ， Chaudhry MH. 1990. “Explicit Methods for 2-D

Transport Free-Surface Flows” J. Hydr.Engrg. 116: 1013-1034.

11. Gharangik A， Chaudhry MH. 1989. “Numerical Simulation of

Hydraulic Jump” J. Hydr. Engrg. 117: 1195-1211.

12. Jameson， A.， Schmidt， W. and Turkel， E. 1981. “Numerical

Solution of the Euler Equation by Finite Volume Method Using Runge-Kutta Time Stepping Schemes”， AIAA 14th Fluid and Plasma Dynamics Conference， 81:1259.

13. Jameson， A. and Baker， T.J. 1984. “Multigrid Solution of the Euler

Equation for Aircraft Configurations” ， AIAA 22nd Aerospace Sciences Meeting，84:93.

14. Molls， T.， and Chaudhry， M. H. 1995. “ Depth-averaged open

channel flow model.” J. Hydr. Engrg. 121:453-464.

15. Molls， T. and Chaudhry ，M. H. and Khan，W.K. (1995) ” Numerical

simulation of two-dimensional flow near a spur-dike”，Advances in Water Resources.18:227-236.

16. Meselhe EA ， Sotiropoulos F ， Holly FM. 1997.”Numerical

Simulation of Transcritical Flow in Open Channels”. J. Hydr. Engrg. 123: 774-783.

(52)

37

Shallow Flows”， John Wiley & Sons， Ltd.

18. Younus and Chaudhry .1994.” A Depth-Averagedk-e Turbulence

Model for the Computation of Free-Surface Flow”. J. Hydr. Res. 326: 414-436.

(53)

β 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 CE 0.944 0.949 0.953 0.951 0.948 0.944 0.935 0.931 0.944 EVP(%) -5.91 -5.182 -4.651 -4.519 -4.439 -4.519 -4.651 -4.784 -4.917 RMSE 0.004 0.004 0.003 0.003 0.003 0.004 0.004 0.004 0.004 β 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 CE 0.926 0.923 0.919 0.916 0.913 0.911 0.908 0.906 EVP(%) -4.638 -4.28 -3.909 -3.605 -3.393 -3.234 -3.014 -2.796 RMSE 0.004 0.004 0.004 0.004 0.004 0.005 0.005 0.005 表 4-1 消散項敏感度分析各指標的結果 β 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 CPU TIME 6.515 6.625 6.578 6.59 6.625 6.69 7.031 6.734 6.67 REAL TIME 31.79 31.77 31.78 31.68 31.6 31.54 31.51 31.49 31.5 β 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 CPU TIME 6.6 6.61 6.75 6.65 6.875 6.625 6.58 6.58 REAL TIME 31.51 31.5 31.42 31.34 31.28 31.23 31.16 31.06 表 4-2 β值與 CPU 耗時、實際模擬時間

(54)

39 圖 3-1 邊界條件示意圖圖 4-1 潰壩案例初始條件示意圖

(55)

圖 4-2 β =0.2 時潰壩第 20 秒的水位比較圖圖 4-3 β =0.4 時潰壩第 20 秒的水位比較圖 distance (m) 0 500 1000 1500 2000 2500 3000 3500 w a te r lev e l (m ) 0 5 10 15 20 25 simulation exaction solution distance (m) 0 500 1000 1500 2000 2500 3000 3500 w a te r lev e l (m ) 0 5 10 15 20 25 simulation exact solution

(56)

41 圖 4-4 β =0.6 時潰壩第 20 秒的水位比較圖圖 4-5 潰壩第 20 秒模擬示意圖 distance (m) 0 500 1000 1500 2000 2500 3000 3500 w a te r lev e l (m ) 0 5 10 15 20 25 simulation exact solution

局部極值消減法應用於二維水深平均模式之研究

國

立

交

通

大

學

土木工程學系

碩

士

論

文