二次流效應影響下湧波傳遞之模擬研究

國立交通大學

土木工程學系

碩士論文

二次流效應影響下湧波傳遞之模擬研究

Numerical Simulation of Surge-Propagation with

Secondary-Flow Effect

研 究 生：鍾浩榮

指導教授：楊錦釧、謝德勇 博士

二次流效應下湧波傳遞影響之模擬研究

Numerical Simulation of Surge-Propagation with

Secondary-Flow Effect

研 究 生：鍾浩榮 Student: Hau-Rong Chung

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

謝德勇 Te-Yung Hsieh

國 立 交 通 大 學

土 木 工 程 研 究 所

碩 士 論 文

A Thesis Submitted to Civil Engineering

College of Engineering

Nation Chiao Tung University

in Partial Fulfillment of the Requirements

for the Degree of Master

in

Civil Engineering

June 2006

Hsinchu, Taiwan, Republic of China

I

誌 謝

感激老師楊教授錦釧、謝博士德勇於研究所兩年間之諄諄教誨， 耐心教導，得以使學位順利完成。 感謝實驗室東霖學長、恩添學長、祥禎學長、夢祺學長、胤隆學 長、昇學學長、世偉學長、曉萍學姊於學業和日常之指導，亦感謝同 學欣瑜、仲達、宣汝、力瑋的共同努力。 另外感謝交通大學土木系 93 級的老同學們的互相砥礪，友情扶 持。 最後謹將本論文獻給逐漸康復的媽媽，支持我的懿軒，還有家中 支持我的爸爸、哥哥、大嫂、姊姊。

二次流效應影響下湧波傳遞之模擬研究

學生：鍾浩榮 指導教授：楊錦釧 謝德勇 國立交通大學土木工程研究所

摘 要

本研究旨在分析二次流效應對彎道湧波模擬之影響。模式在發展 過程中，採用 de Vriend (1977) 的速度剖面假設以求解延散剪應力考 量二次流之效應。在數值處理上，採用雙階分割操作趨近法，將水理 控制方程分割成延散步驟和傳播步驟求解，分別採用顯示和隱示差分 式求解。比較混和型上風法、MacCormack 法和總變量削減法（TVD） 三種方法後，移流項數值方法採用較合適的 MacCormack 法。本研究 採用 Bell 等(1992)所作的實驗來作比較，驗證模式的正確性。本研究 採用的二次流速度剖面的假設為非急彎的條件，因此較小二次流強度 因子案例之模擬結果會較接近實驗值。另外，在不同二次流強度因子 案例中，考慮二次流效應之模擬結果均會較不考慮二次流效應之模擬 結果更接近實驗值。

關鍵字：水深平均、

MacCormack

、二次流、彎道、湧浪

III

Numerical Simulation of Surge-Propagation with

Secondary-Flow Effect

Student：Hau-Rong Chung Advisors：Jinn-Chuang Yang Te-Yung Hsieh Department of Civil Engineering

National Chiao-Tung University

ABSTRACT

The purpose of this study is to numerically analyze the effect of the secondary flow on the surge propagation through a bend. The secondary flow effect was embedded in the dispersion stress terms of the flow momentum equations by adopting the vertical velocity profile proposed by de Vriend (1977). Based on the two-step split-operator approach, the momentum equations were decomposed into dispersion and propagation steps and were solved by using both explicit and implicit techniques. Compared with hybrid upwind scheme and TVD (Total Variation Diminishing) scheme, the MacCormack scheme is most suitable for solving the advection term on simulating the positive surge. The experimental data conducted by Bell al. (1992) was compared with simulating results to examine the accuracy of the proposed model. The vertical shape functions embedded in the model is not appropriate for the flow simulation in the sharp bend. Contrarily, results of the case with the smaller secondary intensity show the better agreement. Results indicate that the simulations either in a sharp bend or not, the simulation considering the secondary flow effect shows better agreement.

Key words: depth averaged, MacCormack scheme, secondary flow, surge

目錄

誌謝... I 中文摘要... II 英文摘要...III 目錄... IV 表目錄... VII 圖目錄...VIII 符號表...X 第一章 緒論...1 1.1 研究動機與方向 ...1 1.2 文獻回顧 ...1 1.3 研究目的與方法 ...2 1.4 章節介紹 ...3 第二章 理論基礎...5 2.1 控制方程式 ...5 2.2 輔助關係式 ...7 2.3 邊界條件 ...8 第三章 數值架構...11 3.1 雙階分割操作趨近法 ...11 3.2 數值差分式 ...13 3.2.1 混合型上風法 ...13 3.2.2 MacCormack 方法 ...14 3.3 一維總變量消減法 ...15 第四章 數值方法比較分析...19

VI 4.1 波傳方程式案例 ...19 4.2 Burgers 方程式案例 ...20 4.3 一維潰壩案例 ...22 第五章 彎道湧波模擬及二次流效應分析 ...35 5.1 彎道湧波模擬 ...35 5.2 二次流影響比較分析 ...37 第六章 結論...62 參考文獻...63

表目錄

表 4.1 數值方法於不同可蘭數的速度波峰位置(波動方程式) ...23 表 4.2 數值方法於不同可蘭數的波峰速度(波動方程式) ...23 表 4.3 數值方法於不同可蘭數的相位誤差(波動方程式) ...24 表 4.4 數值方法於不同可蘭數的消散程度(波動方程式) ...24 表 4.5 不同時間間距的不連續面位置(Burgers 方程式於第 1.8 秒)..24 表 4.6 不同時間間距的相位誤差(Burgers 方程式於第 1.8 秒) ...25 表 4.7 一維潰壩案例結果比較表(第 30 秒)...25 表 5.1 模擬結果與實驗值相關係數比較 ...42 表 5.2 模擬結果與實驗值均方根差比較 ...42 表 5.3 模擬結果與實驗值信賴指標比較 ...42 表 5.4 二次流影響比較湧波深度(Test 21n) ...43 表 5.5 二次流影響相關係數比較(Test 21n) ...44 表 5.6 二次流影響均方根差比較(Test 21n) ...44 表 5.7 二次流影響信賴指標比較(Test 21n) ...44 表 5.8 二次流影響比較湧波深度(Test 21) ...45 表 5.9 二次流影響相關係數比較(Test 21) ...46 表 5.10 二次流影響均方根差比較(Test 21) ...46 表 5.11 二次流影響信賴指標比較(Test 21) ...46

VIII

圖目錄

圖 2.1 彎道二次流示意圖...10 圖 3.1 控制體積法示意圖 (a)實際區域；(b)計算區域 ...18 圖 4.1 初始速度高斯分佈(波動方程式)...26 圖 4.2 解析解(波動方程式)...26 圖 4.3 上風法各可蘭數分佈圖(波動方程式) ...26 圖 4.4 MacCormack 法各可蘭數分佈圖(波動方程式) ...27 圖 4.5 TVD 法各可蘭數分佈圖(波動方程式) ...28 圖 4.6 初始速度的不連續面分佈(Burgers 方程式) ...29 圖 4.7 上風法速度不連續面分佈(Burgers 方程式) ...30 圖 4.8 MacCormacn 法速度不連續面分佈(Burgers 方程式) ...32 圖 4.9 TVD 法速度不連續面分佈(Burgers 方程式) ...33 圖 4.10 一維潰壩案例初始水深狀態...34 圖 4.11 一維潰壩案例第 15 秒水深狀態...34 圖 4.12 一維潰壩案例第 30 秒水深狀態 ...34 圖 5.1 Bell 等(1992)實驗水槽...47 圖 5.2 測站 1 水深量測值...47 圖 5.3 不同θb之 * MaxU 與SI關係圖(謝，2002) ...48 圖 5.4 測站 2 水深比較圖(Test 21n、Test 21) ...49 圖 5.5 測站 4 水深比較圖(Test 21n、Test 21) ...50 圖 5.6 測站 6 水深比較圖(Test 21n、Test 21) ...51 圖 5.7 測站 8 水深比較圖(Test 21n、Test 21)...52 圖 5.8 參考時間解釋圖...53 圖 5.9 測站 2 二次流影響(Test 21n) ...54

圖 5.10 測站 4 二次流影響(Test 21n)...55 圖 5.11 測站 6 二次流影響(Test 21n)...56 圖 5.12 測站 8 二次流影響(Test 21n)...57 圖 5.13 測站 2 二次流影響(Test 21)...58 圖 5.14 測站 4 二次流影響(Test 21)...59 圖 5.15 測站 6 二次流影響(Test 21)...60 圖 5.16 測站 8 二次流影響(Test 21)...61

X

符號表

a =波速； B = 渠道寬度； Cr = Courant number； f C = 摩擦係數； 0 C = 常數； c = Chezy 係數； d = 水深； 0 d = 平均水深； g = 重力加速度； H = 平均水深； 1 h 、h₂ = ξ、η 方向轉換係數； k = von Karman’s 係數； g k =幾何信賴指標 s k =統計信賴指標 L = 渠道長度； N = 模擬案例總計算格點數目； n = 曼寧糙度係數； R = 水力半徑； c R = 渠道中心線平均曲率半徑； r = 曲率半徑； c r = 渠道中心線曲率半徑； SI = 二次流強度因子； 0 S = 渠道坡度；

11 T 、T₁₂、T₂₂ = 有效剪應力項； t = 時間； U = ξ 方向平均速度； u = ξ 方向速度； w u = 近固體邊界的水深平均速度； * u = 剪力速度； V = η 方向平均速度； v = η 方向速度； i X = 某一時刻實驗值； i Y = Xi相對應之模擬結果； w y = 固體邊界與鄰近固體邊界格點的距離； s z = 水面高程； sm z = 平均水面高程； 1 ε 、ε2 = ξ、η方向之亂流傳輸係數； μ = 流體動力黏滯係數； l υ = 層流黏滯係數； t υ = 亂流黏滯係數； ρ = 流體密度； ρc＝相關係數 1 b τ 、 2 b τ = 底床剪應力在 ξ 與 η 方向之分量； ξ、η = 平面上兩正交曲線座標方向； ς = 距離底床之高度與水深之比值； t Δ = 時間間距； x Δ 、Δy = ξ、η方向之格網間距。

XII 上標 b = BF 模式之模擬結果； m = 疊代次數； n = n tΔ 時刻之已知變數； 1 n+ = (n+ Δ1) t時刻之未知變數； 1 2 n+ = (n+ Δ1) t與n tΔ 間之未知變數； * = MacCormack 預測與修正中間的值 ( )= 時間平均; ( ) =水深平均； ('_{) = 時間平均瞬時擾動量。} 下標 c = 模擬值； M = A、D之中點值； m = 實驗值； s = 變數在水面的值； b = 變數在底床的值。

第一章 緒論

1.1 研究動機與方向 過去二維湧浪模擬鮮少考慮二次流效應之影響，主要受制於水流 行經彎道誠屬三維現象，不易以二維模式模擬得到精確之結果，但於 工程規劃設計實務之需要，二維模式常有其功能及實用性，因此本研 究嘗試以水深平均二維模式模擬彎道中二次流(secondary flow)效應 對湧波傳遞的影響，希冀能提供日後水利防災、水利工程設計所用。 1.2 文獻回顧 預測洪水水位以及到達時間是模擬湧波傳遞的主要原因，而湧波 傳遞主要受到移流項(advection term)的影響；為了預測湧波的前進， 正湧波的模擬是不可忽略的部份，而正湧波的精度主要是受到移流項 採用數值方法的影響。 湧波的現象被探討最廣的是潰壩的模擬， Stokers(1957)以解析解 來分析一維潰壩；除了解析方法之外，以數值方法求解也可以得到很 好的結果。二階精度的 MacCormack 方法從 Garcia and Chaudhry (1986) 開始廣為使用，經由水平二維的淺水波方程式來分析潰壩的模擬，但 是並沒有考慮彎曲渠道的湧波模擬；Bellos(1987)以一維模式配合 MacCormack 方法來解乾床的潰壩；為了解決震盪的現象，Fennema and Chaudhry(1990) 加 入 人 工 黏 滯 項 (artificial viscosity) 來 減 少 MacCormack 的震盪現象，並且探討在二維矩形渠道的潰壩湧浪現象。

同樣是二階精度的方法，Harten(1983)發展一個不需要人工黏滯 項且不會震盪的方法，稱為總變量削減法(簡稱 TVD 法)。Yang and Hsu(1993)發展了二維的 TVD 法，應用在保守型曲線座標淺水波方程

2

式，對高速非穩定流的急湧浪(bore waves)在經過收縮擴展的變化來 做出分析。Tseng and Chu (2000)提出增進 TVD-MacCormack 方法的 預測-修正(predictor-corrector)TVD 法對一維潰壩做出很好的模擬精 度。Wang 等(2000)以混合型 TVD 法為出發點，比較各種限制函數對 潰壩湧波模擬精度的效果，並且對複雜的卡式座標二維流場做出潰壩 模擬，也得到很好的效果。Lin 等(2005)以有限體積法實做出二維 TVD 法，並且模擬震波以及潰壩湧浪，得到和實驗值很接近的結果。 Dammuller 等(1989)以 MacCormack 應用在水平二維正交曲線座 標系統，模擬 Miller and Chaudhry (1989)所做的彎道湧浪實驗；Molls and Chaudhry(1995)發展水深平均二維模式並對 Bell 等(1992)所作之 彎道潰壩實驗做驗證，都沒有考慮彎道的二次流效應；Soulis(1992) 使用 MacCormack 方法應用在水平二維非正交曲線座標系統，對 Bellos 等(1991)收縮-擴張水槽潰壩實驗做出分析，此研究為縱方向的 潰壩模擬，沒有關於彎道的模擬；謝(2004)以上風法應用在水深平均 二維模式，演算迴水、彎道、環流流場以及潰壩等案例，並且考量延 散剪應力在彎道的應用，同時展示了演算潰壩模擬的結果，但並沒有 對彎道湧浪做出分析。

實 驗 方 面 ， Bell 等 (1992) 、 Miller(1988) 以 及 Miller and Chaudhry(1989)針對矩形斷面 180 度彎道做了潰壩湧浪的實驗，並且 記載了水位的變化、湧浪在彎道內外側的高度以及波前到達的時間， 這對後人所要做的二維彎道湧浪模擬是很重要的實驗成果。 1.3 研究目的與方法 本研究之研究重點旨在延續謝(2003)所發展之水深平均二維模 式，利用此模式在彎道水理模擬功能的完備性，模擬湧浪在彎道的運

移歷程。此模式為一水深平均二維非穩態之水理數值模式，座標系統 採用正交曲線座標，水理控制方程採用連續及動量方程式求解，為使 水理二次流理論更為完善，模式在發展過程中不忽略延散剪應力。在 數值方法之選用上，採用雙階分割操作趨近法，將水理控制方程分割 成顯示延散步驟和隱式傳播步驟。 文獻回顧了三種重要的數值方法：MacCormack 方法、總變量削 減法(TVD 法)以及上風法。MacCormack 方法主要的精神是以預測-修正法對保守型方程式做離散，在非保守型方程式也可以應用，本研 究基於延續既有的水深平均二維模式，應用 MacCormack 法、上風法 探討非保守型正交曲線座標方程式，而總變量削減法則由於必須用保 守型方程式求得特徵速度才可應用，因此本研究以一維模式來探討 TVD 法於明渠水流應用的問題。本研究對既有水深二維水理模式的 主要貢獻為選用三種數值方法之一來增進移流項的精度。 由前面的文獻回顧可以知道，目前仍然沒有研究探討二次流效應 對彎道湧浪傳遞之影響，因此本研究將利用二維彎道模式，針對此課 題進行進一步的分析工作。 1.4 章節介紹 前三節已闡述本研究之研究動機和方向、文獻回顧、研究目的與 方法，以下將簡要敘述本論文各章節之內容。 第一章為緒論，針對本研究之緣起和方向作說明後，回顧相關之 研究文獻後，再提出本研究之目的與研究方法，並於章末作本論文架 構說明。 第二章為理論基礎，分別闡述本研究之二維正交曲線座標系統水

4 理控制方程式、輔助方程式及相關之邊界條件。 第三章為數值架構，說明水理控制方程式採用的雙階分割操作法 和介紹所採用的三種數值差分式。 第四章為數值方法比較分析，針對本研究採用的三個數值方法， 進行不同案例下的比較，選擇消散程度和相位誤差較小的方法，進行 模式在彎道湧浪的模擬。 第五章為比較湧浪在彎道受到二次流效應影響下湧浪到達時間 和水深的變化所表現的差異，以及不同二次流強度因子案例探討分 析。 第六章為結論，針對本研究成果作綜合性之歸納說明。

第二章 理論基礎

透過座標系統轉換將控制方程式轉換為正交曲線座標系統方程 式，再將此控制方程式作時間平均及水深平均後，即可推得水深平 均二維正交曲線座標模式所需之控制方程式。其中水理控制方程包 含水理連續及動量方程式。茲將水理之理論基礎敘述如下： 2.1 控制方程式 為適度簡化複雜的控制方程式，需對數學模式作若干假設，分 別為(1)不可壓縮牛頓流體(incompressible Newtonian fluid)；(2)靜水 壓分布；(3)忽略風剪力；(4)忽略科氏力。則水深平均二維正交曲線 座標水理控制方程式可表示為 (1)水流連續方程式 1 2 ( 2 ) ( 1 ) 0 d h h h ud h vd t

ξ

η

∂ ₊ ∂ ₊ ∂ ₌ ∂ ∂ ∂ (2.1) (2)水流動量方程式

ξ

方向： 2 1 2 1 2 1 2 1 2

1

1

u

u

u

v

u

h

h

uv

v

t

h

h

h h

h h

∂

∂

∂

∂

∂

∂

+

∂ξ

+

∂η

+

∂η

−

∂ξ

2 11 1 12 1 1 2 1 2 1 1 ( _b ) ( ) ( ) g z d h T hT h h h d h h d

∂

∂

∂

∂ξ

ρ

∂ξ

ρ

∂η

= − + + + 1 1 2 12 22 1 2 1 2 1 h 1 h b T T h h d h h d d

τ

∂

∂

ρ

∂η

ρ

∂ξ

ρ

+ − − 2 11 2 11 1 12 1 12 1 2 1 ( ) s ( ) b ( ) s ( ) b s b s b z z z z h h h h h h d

τ

τ

τ

τ

ρ

ξ

ξ

η

η

⎡ ∂ ∂ ∂ ∂ ⎤ + _⎢− + − + _⎥ ∂ ∂ ∂ ∂ ⎣ ⎦ (2.2)

6

η

方向： 2 2 1 1 2 1 2 1 2 1 1 v u v v v h h uv u t h h h h h h

∂

∂

∂

∂

∂

∂

+

∂ξ

+

∂η

+

∂ξ

−

∂η

2 12 1 22 2 1 2 1 2 1 1 ( b ) ( ) ( ) g z d h T hT h h h d h h d

∂

∂

∂

∂η

ρ

∂ξ

ρ

∂η

= − + + + 2 1 2 11 12 1 2 1 2 1 h 1 h b T T h h d h h d d

τ

∂

∂

ρ

∂η

ρ

∂ξ

ρ

− + − 2 12 2 12 1 22 1 22 1 2 1 ( ) s ( ) b ( ) s ( ) b s b s b z z z z h h h h h h d

τ

τ

τ

τ

ρ

ξ

ξ

η

η

⎡ ∂ ∂ ∂ ∂ ⎤ + _⎢− + − + _⎥ ∂ ∂ ∂ ∂ ⎣ ⎦ (2.3) 式中， 2 2 11 [ 11 ' ( ) ] s b z z T =

∫

τ

−

ρ

u −

ρ

u −u dz (2.4) 2 2 22 [ 22 ' ( ) ] s b z z T =

∫

τ

−

ρ

v −

ρ

v −v dz (2.5) 12 21 [ 12 ' ' ( )( )] s b z z T =T =

∫

τ

−

ρ

u v −

ρ

u −u v −v dz (2.6) 以上諸式中，

_ξ

、

_η

= 平面上兩正交曲線座標方向；h₁ =

ξ

方向轉 換係數；h₂ =

η

方向轉換係數； u =

ξ

方向速度； v =

η

方向速 度；ρ = 流體密度； d = 水深； g = 重力加速度； t = 時間；z_b = 底 床高程； s z = 水面高程； i b

τ

=底床剪應力在

ξ

與

η

方向之分量; ( )= 時間平均; ( ) =水深平均；(') = 時間平均瞬時擾動量；下標s、b 分 別代表變數在水面與底床的值； 11 T ,T₁₂,T₂₂ = 有效剪應力項(effective stress term)，包含層流剪應力、亂流剪應力與延散剪應力(dispersion stresses)。 2.2 輔助關係式

(1)底床剪應力

底床剪應力採用Rastogi and Rodi (1978)之經驗式

1 2 2 1/ 2 ( ) b Cf u u v

τ

=

ρ

+ (2.7) 2 2 2 1/ 2 ( ) b Cf v u v

τ

=

ρ

+ (2.8) 式中， 2 f C =g c = 摩擦係數；c = Chezy 係數。 (2)層流與亂流剪應力 採用Boussinesq之渦流黏性理論，層流與亂流剪應力可合併表 示為 2 11 1 1 1 2 1 ' 2 u v h u h h h

τ

_υ

ρ

ξ

η

⎡ ∂ ∂ ⎤ − = _⎢ + _⎥ ∂ ∂ ⎣ ⎦ (2.9) 2 22 2 2 1 2 1 ' 2 v u h v h h h

τ

_υ

ρ

η

ξ

⎡ ∂ ∂ ⎤ − = _⎢ + _⎥ ∂ ∂ ⎣ ⎦ (2.10) 12 2 1 1 2 2 1 ' ' 2 h v h u u v h h h h

τ

_υ

ρ

ξ

η

⎡ ∂ ⎛ ⎞ ∂ ⎛ ⎞⎤ − = _{⎢ ⎜ ⎟}+ _{⎜ ⎟⎥} ∂ _{⎝ ⎠} ∂ _{⎝ ⎠} ⎣ ⎦ (2.11) 式中， l t

υ υ υ

= + ； l

υ

= 層流黏滯係數；

υ

_t = 亂流黏滯係數 = ku d_* / 6 (Falcon 1979)； 1/ 2 * ( b/ ) u =

τ ρ

= 剪力速度； k = von Karman’s 係數(約 等於0.4)。 (3)延散剪應力 由於水流在進入彎道後，流場隨水流而彎曲，致使流線因彎曲 而產生徑向慣性力，水面因而形成超高以產生徑向靜水壓差，得以 與徑向慣性力取得平衡。在這兩種力之作用下，水流除了以縱向方 向流動外，在徑向尚產生兩層水流，上層水流之外岸慣性力大於靜

8 水壓差，下層水流則反之，因此造成上層水流流動方向為朝外岸， 下層水流則為朝內岸流動，稱之二次流，如圖2.1所示。 為積分水深平均所產生之延散剪應力項，須對流速剖面作一適 當假設。本模式在延散剪應力的處理方面，則是僅考量二次流的影 響，並採用 de Vriend (1977)對二次流速度剖面之假設： 1 g g ln _m( ) u u uf kc kc

ζ

ζ

⎡ ⎤ = _⎢ + + _⎥= ⎣ ⎦ _(2.12) 1 2 2 ( ) 2 ( ) ( ) 2(1 ) ( ) m m g g ud v vf F F f k r kc kc

ζ

⎡

ζ

ζ

ζ

⎤ = + _⎢ + − − _⎥ ⎣ ⎦ (2.13) 式中， 1 ₀1 ln ( ) 1 F

ζ

ζ

d

ζ

ζ

= −

∫

； 2 1 2 ₀ ln ( ) 1 F

ζ

ζ

d

ζ

ζ

= −

∫

；ζ =(z−zb)/d = 距 離底床之高度與水深之比值；r = 曲率半徑；。 de Vriend (1977)二次流速度剖面之適用範圍為(1)水深遠小於渠 道寬度；(2)渠道寬度遠小於渠道之曲率半徑；(3)單一二次流(single secondary eddy only)；(4)完全發展流況(developed flow)。

2.3 邊界條件

本模式目前考量三種邊界條件設定，分別為渠道入流、渠道出 流與固體邊界。一般而言，渠道入流邊界條件設定為單位寬度入流 量，渠道出流邊界條件則採用水位高程設定。在固體邊界處，應用 側壁理論(Law of the wall)，設定靠近固體邊界的邊界條件為：

* 1 ln( ) w u Ey U k + = (2.14) 式中，u =近固體邊界的水深平均速度； E =糙度因子=9.0 (Lien 等 _w

1999)；y y Uw */

υ

+ _{= ，} w y =固體邊界與鄰近固體邊界格點的距離。 圖 2.1 彎道二次流示意圖 外 岸 內 岸

11

第三章 數值架構

3.1 雙階分割操作趨近法 本研究基於分割操作之觀念，將動量方程式分割成二個步驟(延散步 驟及傳播步驟)，延散步驟使用顯示法，傳播步驟利用隱式數值方法求解。 延散步驟求解移流項和擴散項，傳播步驟求解壓力項、底床剪應力項和連 續方程式。據此，水理控制方程可改寫成： 延散步驟 1 2 1 ( ) n n n n V V V T t

ρ

+ ∂ ⎛ ⎞ _{= − ⋅∇ + ∇ ⋅} ⎜ _∂ ⎟ ⎝ ⎠ _(3.1) 傳播步驟 1 1 2 1 ( ) n n n b b V V g z d t t d

τ

ρ

+ + + ∂ ∂ ⎛ ⎞ ₋⎛ ⎞ _{= − ∇ + −} ⎜ _∂ ⎟ ⎜ _∂ ⎟ ⎝ ⎠ ⎝ ⎠ _(3.2) 1 0 n V + ∇ ⋅ = _(3.3) 式中，V 表示速度向量；T表示擴散及延散項；n+1表示(n+ Δ1) t時刻之未 知變數； n 1 n t t + t Δ = − ；n表示n tΔ 時刻之已知變數； 2 1 + n 表示在(n+ Δ1) t與n tΔ 間之未知變數。 (3.1)~(3.3)的一般式可表示成： 延散步驟 2 1 2 1 2 1 2 1 u u u v u h h uv v t h

ξ

h

η

h h

η

ξ

⎡ ⎤ ∂ _{= −} ∂ ₋ ∂ ₋ ∂ ₋∂ ⎢ ⎥ ∂ ∂ ∂ _⎣∂ ∂ _⎦ 2 11 1 12 1 2 12 22 1 2 1 2 1 2 1 2 1 (h T ) 1 (hT ) 1 h 1 h T T h h d h h d h h d h h d

ρ

ξ

ρ

η

ρ

η

ρ

ξ

∂ ∂ ∂ ∂ + + + − ∂ ∂ ∂ ∂

2 11 2 11 1 12 1 12 1 2 1 ( ) s ( ) b ( ) s ( ) b s b s b z z z z h h h h h h d

τ

τ

τ

τ

ρ

ξ

ξ

η

η

⎡ ∂ ∂ ∂ ∂ ⎤ + _⎢− + − + _⎥ ∂ ∂ ∂ ∂ ⎣ ⎦ (3.4) 2 2 1 1 2 1 2 1 v u v v v h h uv u t h

ξ

h

η

h h

ξ

η

⎡ ⎤ ∂ _{= −} ∂ ₋ ∂ ₋ ∂ ₊ ∂ ⎢ ⎥ ∂ ∂ ∂ _⎣∂ ∂ _⎦ 2 12 1 22 1 2 11 12 1 2 1 2 1 2 1 2 ( ) ( ) 1 h T 1 h T 1 h 1 h T T h h d h h d h h d h h d ρ ξ ρ η ρ η ρ ξ ∂ ∂ ∂ ∂ + + − + ∂ ∂ ∂ ∂ 2 12 2 12 1 22 1 22 1 2 1 ( ) s ( ) b ( ) s ( ) b s b s b z z z z h h h h h h d τ τ τ τ ρ ξ ξ η η ⎡ ∂ ∂ ∂ ∂ ⎤ + _⎢− + − + _⎥ ∂ ∂ ∂ ∂ ⎣ ⎦ (3.5) 傳播步驟 2 2 1 ( b ) C u uf v u g z d t h

ξ

d + ⎛ ⎞ ∂ _{= −} ∂ + ₋ ⎜ ⎟ ∂ _⎝ ∂ _⎠ (3.6) 2 2 2 ( b ) C v uf v v g z d t h

η

d + ⎛ ⎞ ∂ ∂ + = − _{⎜ ⎟}− ∂ _⎝ ∂ _⎠ (3.7) 和 2 1 1 2 ( ) ( ) 0 d h ud h vd h h t

ξ

η

∂ ₊∂ ₊ ∂ ₌ ∂ ∂ ∂ (3.8) 針對_n+1時刻的水深值( n 1 d + )做線性化處理，且僅保留一階項，(3.8) 式可改寫成 1 2 1 1 1 2 2 2 ( ) ( ) 0 d d d h h d d t

ξ

α

ξ

β

γ

η

α

η

β

γ

⎛ ⎞ ⎛ ⎞ ∂ ₊ ∂ ∂ Δ _{+ Δ + +} ∂ ∂ Δ _{+ Δ + =} ⎜ ⎟ ⎜ ⎟ ∂ ∂ _⎝ ∂ _⎠ ∂ _⎝ ∂ _⎠ (3.9) 式中， 2 1 1 n h g t d C h_τ

α

= − Δ ； 12 1 2 2 1 1 [ ] n n n _b h h g t z d u C_τ C h_τ

∂

∂

β

∂ξ

∂ξ

+ + Δ = − + ; _{1 1} n d

γ

=

β

； 1 2 2 n h g t d C h_τ

α

= − Δ ; 12 1 1 1 2 2 [ ] n n n _b h h g t z d v C_τ C h_τ

∂

∂

β

∂η

∂η

+ + Δ = − + ; _{2 2} n d

γ

=

β

；

13 1 1 2 2 2 2 ( ) ( ) 1 n n f n C u v C t d τ + ₊ + = + Δ ； n 1 n d d + d Δ = − 。 3.2 數值差分式 本模式採用控制體積(control volume)法的觀念來離散控制方程式，控 制體積法的基本概念如圖3.1所示，其中(a)圖為實際區域，(b)圖為計算區 域，E、W、N、S 表相鄰格點，e、w、n、s 表控制面。在水理控制方程 式中，除了移流項採用一階精度混合型上風法(hybrid scheme)(Spalding 1972)差分以及MacCormack方法外，所有空間差分均採用二階精度的中央 差分法。另外，時間項則採用簡單的前向差分方法。 中央差分法可表示成 1 _{1 1} n _{n n} e w p

ξ

ξ

+ _{+ +} ⎛∂Ψ⎞ ₌ Ψ − Ψ ⎜ _∂ ⎟ _Δ ⎝ ⎠ (3.10) 1 1 1 n n n n s p

η

η

+ _{+ +} ⎛∂Ψ⎞ ₌ Ψ − Ψ ⎜ _∂ ⎟ _Δ ⎝ ⎠ (3.11) 式中， 1 1 1 1 1 1, , 0.5 ( ) 0.5 ( ) n n n n n e E P i j i j + + + + + + Ψ = ⋅ Ψ + Ψ = ⋅ Ψ + Ψ ； 1 1 1 1 1 , 1, 0.5 ( ) 0.5 ( ) n n n n n w P W i j i j + + + + + − Ψ = ⋅ Ψ + Ψ = ⋅ Ψ + Ψ ； 1 1 1 1 1 , 1 , 0.5 ( ) 0.5 ( ) n n n n n n N P i j i j + + + + + + Ψ = ⋅ Ψ + Ψ = ⋅ Ψ + Ψ ； 1 1 1 1 1 , , 1 0.5 ( ) 0.5 ( ) n n n n n s P S i j i j + + + + + − Ψ = ⋅ Ψ + Ψ = ⋅ Ψ + Ψ ； Ψ 可表為 u, v, h1, h2, d, zs 和 z_b。 3.2.1 混合型上風法

混合型上風法為上風法(upwind scheme)與中央差分法組合而成，當移 流效應重要時，採用上風法；移流效應不重要時，則採用中央差分法。至 於 移 流 效 應 重 要 性 的 判 斷 ， 則 採 用 格 網 雷 諾 數 (mesh Reynolds number)R_x、R_y作為判斷的因子，當 R_x 或 R 大於2時，代表移流效應重y 要，差分方法採用能反映方向性的上風法； x R 或 R 小於等於2時，移流y 效應可視為不重要，差分方法採用中央差分法。 混合型上風法應用於本研究移流項的處理可表成 , 1 , 1, , , 1, 1 1 0.5 (1 ) (1 ) i j n n n n n n n i j i j i j i j i j x x u u h

ξ

h

α

ξ

α

ξ

+ + − ⎡ ⎛Φ − Φ ⎞ ⎛Φ − Φ ⎞⎤ ⎛∂Φ ⎞_{= − + +} ⎢ ⎜ ⎟ ⎜ ⎟⎥ ⎜ _∂ ⎟ _{Δ Δ} ⎝ ⎠ ⎣ ⎝ ⎠ ⎝ ⎠⎦ (3.12) , 1 , , 1 , , , 1 2 2 0.5 (1 ) (1 ) i j n n n n n n n i j i j i j i j i j y y v v h

η

h

α

η

α

η

+ + − ⎡ ⎛Φ − Φ ⎞ ⎛Φ − Φ ⎞⎤ ⎛∂Φ ⎞_{= − + +} ⎢ ⎜ ⎟ ⎜ ⎟⎥ ⎜ _∂ ⎟ _{Δ Δ} ⎝ ⎠ ⎣ ⎝ ⎠ ⎝ ⎠⎦ (3.13) 其中 0 2 1 2 1 2 x x x x R R R

α

⎧ ≤ ⎪⎪ =_⎨ > ⎪ − < − ⎪⎩ ; 0 2 1 2 1 2 y y y y R R R

α

⎧ ≤ ⎪⎪ =_⎨ > ⎪ − < − ⎪⎩ (3.14) 上列諸式中， , 1 , n i j i j x u h R

ξ

μ ρ

Δ = ； , 2 , n i j i j y v h R

η

μ ρ

Δ = ；

μ

= 流體動力黏滯係數 (dynamic viscosity)；Φ 可表成 u 或 v 。 3.2.2 MacCormack 方法 MacCormack 法為預測-修正的方法，應用在本研究的移流項的處理可 以表示為：

15 , , * , , ₁ 1 1 2 ( ) [ ] [ ] i j i j n _{n n n n} n n i j i j _{i i} j j u u u u u v u u t h

ξ

h

η

+ + − ₋ − = − − Δ Δ Δ (3.15) 1 * 2 , , , 1 ( ) 2 n n i j i j i j u + = u +u (3.16) , , 1 1 2 * * * * , , 1 1 1 2 ( ) [ ] [ ] 1 2 i j i j n n n n i j i j i i j j u u _{u u u v} u u h h t

ξ

η

+ + − − − ₋ − = − − Δ Δ Δ (3.17) , , * , , 1 1 1 2 ( ) [ ] [ ] i j i j n _{n n n n} n n i j i j i i j j v v _{u v v v} v v t h

ξ

h

η

+ + − _{= −} − ₋ − Δ Δ Δ (3.18) 1 * 2 , , , 1 ( ) 2 n _n i j i j i j v + = v +v (3.19) , , 1 1 _{2 * *} * * , , ₁ 1 1 2 ( ) [ ] [ ] 1 2 i j i j n n _{n n} i j i j _{i i} j j v v u v v v v v h h t

ξ

η

+ + − − − − − = − − Δ Δ Δ (3.20) 其中 * u 以及 * v 表示為中間時刻的值，Δt為時間間距，(3.15)式以及 (3.18)式為預測步驟，(3.17)式以及(3.20)式為修正步驟。 3.3 一維總變量削減法 本研究考慮一般雙曲型保守型問題，其方程式的形式為：

0

u

E

t

x

∂

₊

∂

₌

∂

∂

(3.20) E 為通式，在Burgers方程式為流速平方之半，在一階波動方程式 就是流速與波速乘積。 若以通量來分解(3.20)式

1 1 1 2 2

1

(

)

2

n n n n i i i i

t

u

u

h

h

x

+ + −

Δ

=

−

−

Δ

(3.21) (3.21)式中 h 為數值通量，參考 Wang 等(2000)所採用的混合型 TVD

法，以上風法和Lax-Wendroff 法用限制函數(limiter function)組合數值通

量， 2 1 1 1 1 1 2 2 2 2 1 {( ) [ (1 )] } 2 i i i i i i t h E E a u x

φ

φ ψ

+ + + + + Δ = + − + − Δ Δ (3.22) 2 1 1 1 1 1 2 2 2 2 1 {( ) [ (1 )] } 2 i i i i i i t h E E a u x

φ

φ ψ

− − − − − Δ = + − + − Δ Δ (3.23) 其中 a為特徵速度， ₁ 2 i+ u Δ 定義為u_i+1−u_i 1 1 1 2 2 1 2 1 2 ( ) if 0 if 0 i i i i i i E E u u a E u u + + + + + − ⎧ _{Δ ≠} ⎫ ⎪ _Δ ⎪ ⎪ ⎪ = ⎨ ⎬ ⎪∂ _{Δ =} ⎪ ⎪_∂ ⎪ ⎩ ⎭ ， 1 1 1 2 2 1 2 1 2 ( ) if 0 if 0 i i i i i i E E u u a E u u − − − − − − ⎧ _{Δ ≠} ⎫ ⎪ _Δ ⎪ ⎪ ⎪ = ⎨ ⎬ ⎪∂ _{Δ =} ⎪ ⎪_∂ ⎪ ⎩ ⎭

φ為限制函數，使用 Roe 的 Superbee 函數(參考Yee 1986)，如下

2 2 1 1 1 3 3 2 1 2 2 2 2 1 1 2 2 2 2 1 1 1 1 1 1 2 2 2 2 , ( )( ) ( )( ) if 0, if 0 (( )( )) (( )( )) i i i i i i i i i i i i i i i i i i t x a a u u a a u u r a r a a a u u a a u u

λ

λ

λ

λ

λ

− + + − − + + + + + + + + + + Δ = Δ − − − − = > = < − − − −

φ

=MAX(0,min(2r,1),min(r,2))

ψ

為消散函數，定義如下

17 1 1 1 1 2 2 2 2 2 2 2 2 1 1 1 1 2 2 2 2 1 1 2 2 ... ... , ... ... 2 2 i i i i i i i i i i a a a a a a a a ε ε ψ _ε ψ _ε ε ε ε ε + + − − + − + − + − ⎧ ⎫ ⎧ ⎫ ≥ ≥ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ =_{⎨ + ⎬} =_{⎨ + ⎬} ⎪ ⎪ ⎪ ⎪ < < ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎭ ⎩ ⎭ ε為極小且大於零的數，約為0.5~0.05之間。

P i,j e E i+1,j w n s W i-1,j i,j+1 N S i,j-1 P i,j e _E i+1,j s n W i-1,j w S i,j-1 N i,j+1 ξ η X Y (a) (b) 圖 3.1 控制體積法示意圖 (a)實際區域； (b)計算區域

19

第四章 數值方法比較分析

經由第二章與第三章在模式理論基礎及數值方法之說明介紹 後，本章對應用的數值方法分析比較，包括上風法、MacCormack 法 以及總變量削減法(TVD)。 4.1 波傳方程式案例 一維波傳方程式是表現一固定波形以波速a 向縱方向傳遞，理論 上正確的數值模擬不應該有消散，也不應該有相位誤差。以波傳方程 式比較數值方法的目的是分析可蘭數(courant number)對消散程度以 及相位誤差的影響。 一維波傳方程式為： 0 u u a t x ∂ ₊ ∂ ₌ ∂ ∂ (4.1) u 為流速。本設計案例總長為 10000 公尺，總模擬時間為 160 分 鐘，a=0.5=波速，可蘭數為Cr a t x Δ = Δ ，上游邊界條件為高斯分佈如(4.2) 式， 2 ( ) 10exp 528 at u= − (4.2) 模擬開始後邊界條件為此高斯分佈另外一半以波速a向正x方向 進入計算區域。 2 ( ) 10exp 528 x u= − (4.3) 初始條件為(4.3)式,0≤ ≤x 10000，如圖 4.1；解析解為高斯分佈的

波形不改變，以波速 a 向右行進，如圖 4.2。 接下來以上風法、MacCormack法以及 TVD法模擬在不同可蘭數 時的相位誤差和消散程度，相位誤差是指某一時刻模擬速度波峰位置 和解析解速度波峰位置的差距，再除以解析解速度波峰位置；消散程 度是某一時刻模擬速度波峰模擬結果速度和解析解速度波峰值的差 距，再除以解析解速度波峰值。三個方法的波峰位置和速度的比較詳 列在表 4.1 和表 4.2，相位誤差和消散程度列表在表 4.3 和表 4.4，三 個方法在可蘭數為1時都沒有相位誤差和消散。上風法如圖4.3a~c， 可以觀察到上風法隨可蘭數越小則波形越消散； MacCormack 模擬 結果如圖4.4a~c，可以發現有減少消散的程度，但是增加相位誤差； 第三項是以 TVD 法模擬，如圖 4.5a~c，相位誤差和消散程度較 MacCormack方法精度低；以表 4.3和表4.4 的第9600秒可以觀察到 當可蘭數不等於1時，三個方法均存在明顯的誤差，上風法的相位誤 差比 MacCormack 好，而 TVD 法更次之，但是消散程度卻是三個方 法中誤差最大的，以MacCormack 方法來說，其受到可蘭數的影響是 三者中比較小的。 4.2 Burgers 方程式案例 本節案例所用的 Burgers 方程式是無黏性的非線性方程式，保守 形式為 2 ( ) 2 ₀ u u t x ∂ ∂ _{+ =} ∂ ∂ (4.4) 以非保守形式表示為 0 u u u t x ∂ ₊ ∂ ₌ ∂ ∂ (4.5)

21 本設計案例如圖 4.6，初始條件為： 0 2 u=1 m/s 2 4 u=0 m/s x x ≤ < ⎧ ⎨ ≤ ≤ ⎩ ，上游 邊界條件為u=1 m/s，下游邊界條件為 0 m/s，總模擬時間為 2 秒，渠 道總長為 4 公尺。本案例以上風法、MacCormack 法以及 TVD 法來 模擬， (4.4)式中的 u2/2 為(3.20)式中的 E，特徵速度就是 u。 為測試時間間距對 Burgers 方程式不連續速度分佈移動的影響， 以上列設計案例改變時間間距，網格間距固定為0.1m，時間間距 0.5 秒以及 0.1 秒適用在上風法、MacCormack 法以及 TVD 法，另外 0.2 秒可以適用在上風法和 TVD 法，本案例以相位誤差探討不連續面移 動的精度，相位誤差為參考Hoffmann(1993)模擬結果來和本研究模擬 結果做比較；模擬結果如表4.5，本研究 MacCormack 方法在Δt=0.2s 時結果會發散，因此以*表示，而 Hoffmann(1993)所做的模擬因為沒 有 Δt=0.2s 的結果，所以也是以*表示；圖 4.7a~c 列出了上風法在 Δt=0.2s、Δt=0.1s 和 Δt=0.05s 的不連續面模擬結果，可以發現不連續 面 在 三 種 時 間 間 距 的 情 況 都 沒 有 相 位 誤 差 ， 圖 4.8a~b 列 出 MacCormack 方法在Δt=0.1 和 Δt=0.05 的模擬結果，Δt=0.1s 的情況沒 有相位誤差，但是在 Δt=0.05 的情況會有震盪，圖 4.9a~c 列出 TVD 方法在Δt=0.2s、Δt=0.1s 和Δt=0.05s 的不連續面模擬結果，其中Δt=0.1s 到Δt=0.05s 會有相位誤差。 以上結果經由整理，列表各個不連續面位置於表4.5，表 4.6 則詳 列 以 Hoffmann(1993) 為 標 準 和 本 研 究 三 種 方 法 的 相 位 誤 差 ， MacCormack 法不連續面傳播的準確率在比較的方法裡面是較好且不 受到時間間距的影響。本研究的目標為湧浪的到達時間方面的預測， 且基於非守恆的控制方程式的實作方便，本研究並沒有把 TVD 方法 應用在二維模式中。

4.3 一維潰壩案例 本案例旨在比較MacCormack 方法和上風法在正湧波、負湧波以 及波前速度模擬的精確度。 圖4.10 是本案例的初始狀態示意圖，壩址位於渠道中間，距離上 下游各500m，t=0 時水壩閘門為關閉狀態，閘門上游(左側)為水深 10m 的靜止水體，下游(右側)為水深 2m 的靜止水體，當 t>0 時瞬時移開 閘門，渠道假設無摩擦效應。案例中採用 201x36 均勻計算網格，數 值參數為Δx=5m，Δy=1m，Δt=0.3s，總模擬時間為 30 秒。模擬時間內 上下游邊界都沒有受到影響，因此上游入流量為零及下游水位保持 2m。 圖4.11、圖 4.12 分別是第 15 秒和第 30 秒時渠道中心線和解析解 之比較圖，比較正湧浪和負湧浪可以發現，上風法以及MacCormack 法在負湧浪的表現極為近似，呈現數值消散的誤差，而正湧浪的模擬 結果發現 MacCormack 的精度比上風法接近解析解，詳細列表於表 4.7。 綜合前述之分析，本研究採用 MacCormack 方法為移流項的數值 方法，進行後續彎道湧波之模擬分析。

23 表4.1 數值方法於不同可蘭數的速度波峰位置(波傳方程式) 單位：公尺 時刻 可蘭數 差分式 2400 秒 4800 秒 7200 秒 9600 秒 上風法 1200 2400 3600 4800 MacCormack 方法 1200 2400 3600 4800 1 TVD 法 1200 2400 3600 4800 上風法 1200 2400 3600 4800 MacCormack 方法 1000 2200 3400 4600 0.5 TVD 法 1000 2200 3200 4200 上風法 1200 2400 3600 4800 MacCormack 方法 1200 2200 3400 4600 0.375 TVD 法 1000 2000 3000 4000 解析解 1200 2400 3600 4800 表 4.2 數值方法於不同可蘭數的速度波峰值 (波傳方程式) 單位：公尺每秒 時刻 可蘭數 差分式 2400 秒 4800 秒 7200 秒 9600 秒 上風法 10 10 10 10 MacCormack 方法 10 10 10 10 1 TVD 法 10 10 10 10 上風法 7.549 5.9 4.99 4.4 MacCormack 方法 9.88 9.59 9.1 8.62 0.5 TVD 法 8.24 7.32 6.72 6.32 上風法 6.99 5.34 4.48 3.94 MacCormack 方法 9.82 9.45 8.91 8.4 0.375 TVD 法 8.22 7.36 6.85 6.56 解析解 10 10 10 10

表 4.3 數值方法於不同可蘭數的相位誤差(波傳方程式) 時刻 可蘭數 差分式 2400 秒 4800 秒 7200 秒 9600 秒 上風法 0% 0% 0% 0% MacCormack 方法 0% 0% 0% 0% 1 TVD 法 0% 0% 0% 0% 上風法 0% 0% 0% 0% MacCormack 方法 16.67% 8.33% 5.56% 4.17% 0.5 TVD 法 16.67% 8.33% 11.11% 12.5% 上風法 0% 0% 0% 0% MacCormack 方法 0% 8.33% 5.56% 4.17% 0.375 TVD 法 16.67% 16.67% 16.67% 16.67% 表 4.4 數值方法於不同可蘭數的消散程度(波傳方程式) 時刻 可蘭數 差分式 2400 秒 4800 秒 7200 秒 9600 秒 上風法 0% 0% 0% 0% MacCormack 方法 0% 0% 0% 0% 1 TVD 法 0% 0% 0% 0% 上風法 24.57% 41.02% 50.08% 55.95% MacCormack 方法 1.18% 4.08% 9.04% 13.78% 0.5 TVD 法 17.63% 26.8% 32.85% 36.77% 上風法 30.11% 46.58% 55.17% 60.62% MacCormack 方法 1.75% 5.46% 10.93% 16.03% 0.375 TVD 法 17.82% 26.38% 31.45% 34.39% 表 4.5 不同時間間距的不連續面位置(Burgers 方程式於第 1.8 秒) Δt(秒) 方法 0.2 0.1 0.05 Hoffmann(1993) * 2.8 2.8 上風法 2.8 2.8 2.8 MacCormack 法 * 2.8 2.7 TVD 法 2.8 2.5 2.3 單位:公尺

25 表4.6 不同時間間距的相位誤差(Burgers 方程式於第 1.8 秒) Δt 方法 0.2 0.1 0.05 上風法 * 0% 0% MacCormack 法 * 0% 3.571% TVD 法 * 10.71% 17.86% 表4.7 一維潰壩案例結果比較表(第 30 秒) 比較項目 解析解 MacCormack 法 上風法 正湧波波前前進速度(m/s) 5.692 5.51 5.418 波前高度(m) 5.079 5.359 5.467

0 1 2 3 4 5 6 7 8 9 10 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 ₁₀₀₀₀ Diatance(m) V e lo c ity(m/ s) 圖4.1 初始速度高斯分佈(波傳方程式) 0 1 2 3 4 5 6 7 8 9 10 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 Distance(m) V elocit y (m /s ) 0 SEC 2400 SEC 4800 SEC 7200 SEC 9600 SEC 圖 4.2 解析解(波傳方程式)

27 0 2 4 6 8 10 0 1000 2000 3000 4000 5000 6000 7000 Distance(m) V e lo ci ty (m /s) 0 SEC 2400 SEC 4800 SEC 7200 SEC 9600 SEC 圖4.3a 上風法可蘭數=1 分佈圖(波傳方程式) 0 2 4 6 8 10 0 1000 2000 3000 4000 5000 6000 7000 Distance(m) V e lo ci ty (m /s) 0 SEC 2400 SEC 4800 SEC 7200 SEC 9600 SEC 圖4.3b 上風法可蘭數=0.5 分佈圖(波傳方程式) 0 2 4 6 8 10 0 1000 2000 3000 4000 5000 6000 7000 Distance(m) V e lo ci ty (m /s ) 0 SEC 2400 SEC 4800 SEC 7200 SEC 9600 SEC 圖4.3c 上風法可蘭數=0.375 分佈圖(波傳方程式)

0 2 4 6 8 10 0 1000 2000 3000 4000 5000 6000 7000 Distance(m) V e lo ci ty (m /s) 0 SEC 2400 SEC 4800 SEC 7200 SEC 9600 SEC 圖 4.4a MacCormack 法可蘭數=1.0 分佈圖(波傳方程式) 0 2 4 6 8 10 0 1000 2000 3000 4000 5000 6000 7000 Distance(m) V e lo ci ty (m /s) 0 SEC 2400 SEC 4800 SEC 7200 SEC 9600 SEC 圖4.4b MacCormack 法可蘭數=0.5 分佈圖(波傳方程式) 0 2 4 6 8 10 0 1000 2000 3000 4000 5000 6000 7000 Distance(m) V e lo ci ty (m /s) 0 SEC 2400 SEC 4800 SEC 7200 SEC 9600 SEC 圖4.4c MacCormack 法可蘭數=0.375 分佈圖(波傳方程式)

29 0 2 4 6 8 10 0 1000 2000 3000 4000 5000 6000 7000 Distance(m) Vel o ci ty (m/s ) 0 SEC 2400 SEC 4800 SEC 7200 SEC 9600 SEC 圖 4.5a TVD 法可蘭數=1.0 分佈圖(波傳方程式) 0 2 4 6 8 10 0 1000 2000 3000 4000 5000 6000 7000 Distance(m) V e lo ci ty (m /s) 0 SEC 2400 SEC 4800 SEC 7200 SEC 9600 SEC 圖4.5b TVD 法可蘭數=0.5 分佈圖(波傳方程式) 0 2 4 6 8 10 0 1000 2000 3000 4000 5000 6000 7000 Distance(m) V e lo ci ty (m /s) 0 SEC 2400 SEC 4800 SEC 7200 SEC 9600 SEC 圖4.5c TVD 法可蘭數=0.375 分佈圖(波傳方程式)

x u 2.0 0 1.0 4.0 圖4.6 初始速度的不連續面分佈(Burgers 方程式案例) 0 0.2 0.4 0.6 0.8 1 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 Distance(m) Ve lo city (m/ s) 0 0.6 SEC 1.2 SEC 1.8 SEC 圖4.7a 上風法Δt=0.2s 速度不連續面分佈(Burgers 方程式) 0 0.2 0.4 0.6 0.8 1 0.1 0.6 1.1 1.6 2.1 2.6 3.1 3.6 Distance(m) Velocity( m /s ) 0 0.6 SEC 1.2 SEC 1.8 SEC 圖4.7b 上風法 Δt=0.1s 速度不連續面分佈(Burgers 方程式)

31 0 0.2 0.4 0.6 0.8 1 0.1 0.6 1.1 1.6 2.1 2.6 3.1 3.6 Distance(m) Velocity( m /s ) 0 0.6 SEC 1.2 SEC 1.8 SEC 圖 4.7c 上風法Δt=0.05s 速度不連續面分佈(Burgers 方程式)

0 0.2 0.4 0.6 0.8 1 0.1 0.6 1.1 1.6 2.1 2.6 3.1 3.6 Distance(m) Velocity( m /s ) 0 0.6 SEC 1.2 SEC 1.8 SEC

圖 4.8a MacCormacn 法Δt=0.1s 速度不連續面分佈(Burgers 方程式)

0 0.2 0.4 0.6 0.8 1 1.2 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 Distance(m) Vel o ci ty (m/s ) 0 0.6 SEC 1.2 SEC 1.8 SEC 圖4.8b MacCormacn 法Δt=0.05s 速度不連續面分佈(Burgers 方程式)

33 0 0.2 0.4 0.6 0.8 1 0 0 1 1 2 2 3 3 Distance(m) Vel ocity(m /s) 0 0.6 SEC 1.2 SEC 1.8 SEC 圖4.9a TVD 法Δt=0.2s 速度不連續面分佈(Burgers 方程式) 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 3.5 Distance(m) Velo city(m /s) 0 0.6 SEC 1.2 SEC 1.8 SEC 圖4.9b TVD 法Δt=0.1s 速度不連續面分佈(Burgers 方程式) 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 3.5 Distance(m) V elocity (m/s) 0 0.6 SEC 1.2 SEC 1.8 SEC 圖4.9c TVD 法Δt=0.05s 速度不連續面分佈(Burgers 方程式)

圖4.10 一維潰壩案例初始水深狀態 0 2 4 6 8 10 0 100 200 300 400 500 600 700 800 900 1 000 Distance (m) Dept h (m) MacCormack scheme Hybrid upwind analytical solution 圖 4.11 一維潰壩案例第 15 秒水深狀態 0 2 4 6 8 10 0 100 200 300 400 500 600 700 800 900 1000 Distance (m) Dept h (m) MacCormack scheme Hybrid upwind analytical solution 圖4.12 一維潰壩案例第 30 秒水深狀態 壩址

35

第五章 彎道湧波模擬及二次流效應分析

5.1 彎道湧波模擬

本研究採 Bell 等(1992)彎道湧浪實驗數據(Test 21n 以及 Test 21) 來驗證本研究在彎道湧浪模擬的正確性。水平矩形實驗水槽如圖 5.1 所示，渠道寬 0.305m，總長 12.26m，測站 1(station 1)到測站 2 之前 是直線段，長 4.323m，測站 2 到測站 6 為 180 度的彎道，內徑 0.914m， 外徑 1.219m，測站 6 之後到尾端都為直線道，長 3.738m，測站 8 距 離末端 1.448m，測站 4 位於彎道中點 90 度之處。而造成湧波的潰壩 裝置位於測站 1，以一個可以瞬間拉起的鐵板蓄水。 一般認為曲率半徑(R)對渠寬(B)比小於 3(R/B<3)的彎道屬於急 彎，Lien 等(1999）建議以滑移邊界作為急彎案例的模擬條件，結果 比非滑移邊界較符合實驗數據。雖然本案例並不符合此急彎的範疇， 但是 Molls and Chaudhry(1992)認為該案例仍然有很強的二次流特 性，應以滑移邊界設定側邊界條件。本研究參考以上研究，以滑移邊 界作為探討二次流影響的假設條件。 模擬參數的設定，本研究使用 249x26 的網格計算，時間間距 0.005s，依照 Bell 等(1992)的實驗參數，上游流量邊界條件如圖 5.2 所示以測站 1 之水深作為率定，渠道初始水深為 0.0762m，下游水深 邊界為 0.0762m，曼寧 n=0.04，下游直線道多延伸了 4m 使模擬時間 內下游邊界不影響計算果結果，因此下游邊界可以初始水深作為水深 邊界為 0.0762m。曼寧 n=0.04 為案例 Test21n，第二個案例 Test21 的 參數唯一不同是曼寧 n=0.0165。 參考謝(2003)之研究，二次流強度因子(SI)計算式可表示為：

f c C r H SI = (5.1) = _/ 2 f

C

g

c

，圖 5.3(謝，2003)為不同θb之 * MaxU 與 SI 關係圖， ( / 2 ) b L rc

θ

=

π

為彎道長度因子。 * MaxU 為ξ方向速度之最大相對差異 = ( / ) b u u Max Δ ，其中 uΔ =縱向的水深平均速度差， b u =彎道中縱向的水 深平均速度差，由圖 5.3 可以發現在固定θb的情況，SI 值越大， * MaxU

也就越大。在 Bell 等(1992)的兩個實驗案例中，Test 21n 為 SI=0.89 以及 Test 21 為 SI=2.13。 為了比較整體模擬結果和實驗值的差距，本研究使用相關係數 c

ρ

、均方根差 RMS 以及信賴指標

k

_g /

k

_s來作為比較的工具，其中： 2 2 ( )( ) ( ) ( ) i i c i i x x y y x x y y ρ = − − − −

∑

∑

∑

(5.2) − =

∑

2 ( ) [ ] i i N

x

y

RMS

N

(5.3) 幾何信賴指標＝ 2 1 2 1 1 1 1 [ ] 1 1 1 1 [ ] 1 i N i i i i g i N i i i i y x y N x k y x y N x = = − + − = − − −

∑

∑

(5.4) 統計信賴指標＝ 2 1 1 exp( [log( )] ) N i s i i y k N = x =

∑

(5.5) 其中，x_i為某一時刻的實驗值，y_i為同一時刻的模擬值，上標

−

表示平均數，N 表示比較值的數量，幾何信賴指標與統計信賴指標的

37

比值(

k

_g /

k

_s)接近 1則表示模擬結果越接近實驗值；RMS 越小表示模

擬結果的整體誤差越小；

ρ

_c越大表示模擬結果越接近實驗值。

圖 5.4(a)(b)～圖5.7(a)(b)分別為測站2 到測站8 的數值模擬結果

和實驗數據的比較圖，(a)表示案例 Test 21n，(b)表示案例 Test 21。

測站 2和實驗數據比較，模擬結果有分內外側不同的數據，由於 測站2 是彎道入口，在湧波完全進入後會受到外側水深超高的影響， 因此會有外側超高的現象是符合預期的結果；測站4 和測站6 的模擬 值有相同的結論：在大SI值的案例，其水位超高會比小SI值案例大 的現象，而模擬結果也反映出同樣的趨勢，水深的模擬大致符合實驗 數據，而到達時間和實驗值比較也符合實驗的結果；測站8 是直線道 接近終點的測站，模擬結果顯示內外側水深已經相差不明顯。 表 5.1 為兩個案例模擬結果與實驗值用相關係數作比較，測站 2 兩個案例模擬結果顯示相關係數很高，都相當接近實驗值；測站4 以 及測站6 相關係數高，同時小SI值的相關係數比大SI值要高，顯示 在彎道中點和彎道出口，小 SI 值比較接近實驗值；直線道終點的測 站8 也是有小SI值的相關係數較高的情形。 表 5.2為模擬結果均的方根差比較，可以發現小SI值的案例在彎 道中間、彎道出口以及直線道終點的模擬結果比大 SI值的誤差要小。 表 5.3說明信賴指標的比值方面，小SI值的案例在彎道中間、彎 道出口和直線道終點都很接近1，也比大SI值案例更接近1。 以模擬結果和三種統計方法比較來看，本研究兩個案例的模擬結 果與實驗結果都相當吻合，驗證本模式彎道湧浪模擬之正確性。 5.2 二次流影響比較分析

根據謝(2003)對彎道水理二次流所作的分析中，水深平均二維模 式可分為兩種型態，一種為傳統模式(Conventional Model；CN模式)， 一種為彎道模式(Bend-Flow Model；BF模式)，這兩種模式最主要之 差異在於動量方程式中對延散剪應力的處理。CN 模式忽略彎道二次 流之效應，即假設垂直方向之流速為均勻分佈，因此忽略延散剪應力 項；BF 模式則為引進深度方向二次流流速剖面來處理延散剪應力 項，以適切反映彎道二次流效應之影響。本研究的 SI 值計算公式為 (5.1)式，透過謝(2003)之分析，得知當二次流強度因子(SI)小於 0.02 時，二次流效應對流場之影響可忽略不計，本研究的SI值為0.89（案 例Test 21n）以及2.13(案例Test 21)，因此都必須考量二次流效應。 為比較實驗值與模擬結果的誤差量，本研究對誤差的定義： 實驗值-模擬值 誤差＝ 100％ 實驗值 × (5.2) 本研究為比較波前的水深，定義一「參考時間」：參考圖 5.8 某 一測站的外側水深實驗數據有跳起現象(也就是實驗值的明顯波前) 的時間，先找到參考時間後，再分別找到所對應到的水深，其中有內 外側實驗水深以及數值內外側的水深，由此才可比較同一時刻內數值 和實驗的水深誤差。

圖5.9(a)(b)~5.12(a)(b)為案例Test 21n(SI=0.89)考慮與不考慮二次

流影響下湧波在彎道傳遞的模擬結果和實驗數據比較圖；圖5.9(a)(b)

為測站2 的比較圖，整體差距不大，但由外側可以明顯的發現考慮二

次流影響下較為接近實驗值；圖 5.10 中，觀察測站 4 波前到達時可

發現考慮二次流影響下比較接近實測值；圖 5.11(a)(b)，測站 6 的模

39 較接近實驗值；觀察圖 5.12(a)(b)也可以發現測站 8 的水深模擬結果 在二次流影響下比較接近實測值。以下列出表5.4 做波前水深模擬數 據做細部的比較。 表 5.4是案例Test 21n 考慮與不考慮二次流影響下模擬深度和實 驗深度比較的結果。參考時間在測站2～測站8 分別是2.81秒、3.79 秒、4.78秒以及6.46秒。 表 5.4 中剛進彎道的測站 2，考慮二次流並未減少誤差，這是由 於實驗值在彎道入口處沒有外側水位超高的紀錄，因此只能以相同的 實驗值和內外側的模擬值作誤差比較。彎道中點測站4 的水深在考慮 二次流影響下比較接近實驗數據。測站 6彎道考慮二次流影響下，其 誤差都比較小，很明顯的沒有考慮二次流的效應的話會有比較大的誤 差，同樣直線道出口也是同樣的情況。 表 5.5 為相關係數的比較，考慮與不考慮二次流影響都很接近 1，比較考慮與沒有考慮二次流影響的相關係數，可以發現有考慮二 次流效應的模擬結果其相關性都比沒有考慮二次流效應的精確度要 高一些。 表 5.6為均方根差的比較；可以發現除了測站2 以外，彎道的測 站 4、測站 6 以及出口的測站 8，都因為考慮二次流而使得平均誤差 比不考慮二次流小。 表 5.7為信賴指標比值比較，由彎道中間、彎道出口的信賴指標 可以發現，在考慮二次流的情況下，模擬結果會比較接近實驗值。 本案例(SI=0.89)的模擬結果顯示，在模擬波前深度時，二次流的 影響在彎道中不可忽略，由三種統計方法的整體比較可以發現，彎道

中點和彎道出口在考慮二次流影響下，其模擬結果比不考慮二次流更 接近實驗值。

圖5.13(a)(b)~5.16(a)(b)為案例Test 21(SI=2.13)考慮與不考慮二次

流影響下，湧波在彎道傳遞的模擬結果和實驗數據比較圖；圖 5.13 為測站2 的比較圖，由圖中波前到達時刻可發現不考慮二次流影響下 彎道入口較為接近實驗值；圖 5.14 中，測站 4 不考慮二次流影響下 波前會比較慢；而測站 6 波前延遲的現象更明顯，可以由圖 5.15 看 出來；測站8 則是大致符合實驗的結果；以下列出表5.8 做波前水深 數據做細部的比較。 表5.8是案例Test 21考慮與不考慮二次流影響下模擬深度和實驗 深度比較的結果。參考時間在測站 2 到測站 8 分別是 2.64 秒、3.57 秒、4.46秒以及6.07秒。 表 5.8中剛進彎道的測站2，和Test21n相同，考慮二次流並不一 定減少誤差。彎道中點的測站4 的水深在考慮二次流影響下比較接近 實驗數據，由內側的數據明顯的看出來。測站6 彎道考慮二次流影響 下，其誤差都比較小，很明顯的沒有考慮二次流的效應的話會有比較 大的誤差。直線道出口測站8 發現實驗值和渠道內側模擬值差距很大 的情況，這是由於渠道末端模擬結果沒有水位超高的現象，只能以相 同的模擬值和 內外側的實驗值作誤差比較。 表 5.9 為相關係數比較，比較有考慮二次流影響與沒有考慮二次 流影響的相關係數，可以發現有考慮二次流效應會比較接近實驗值。 表 5.10 為均方根差的比較；可以發現除了測站 2 以外，彎道中 的測站 4、測站 6 以及出口的測站 8，都因為開啟二次流而使的平均 誤差比關閉二次流小。

41 表 5.11為信賴指標比值比較，由彎道中間、彎道出口的信賴指標 可以發現，在考慮二次流的情況下，模擬結果會比較接近實驗值。 SI=2.13 的案例模擬結果顯示，在模擬波前深度時，二次流的影 響在彎道中不可忽略，由三種統計方法的整體比較可以發現，彎道中 點和彎道出口在考慮二次流影響下，其模擬結果比不考慮二次流較接 近實驗值。 兩個案例比較，在大SI值的案例(Test 21)，考慮二次流不一定在 所有的比較指標中表現很好的精度，例如信賴指標方面就只有測站6 表現比較好，而在小 SI 值的案例明顯看得出二次流影響在彎道的貢 獻。

表5.1 模擬結果與實驗值相關係數比較 SI=2.13 SI=0.89 彎道內側 彎道外側 彎道內側 彎道外側 測站2 0.956 0.924 0.91 0.936 測站4 0.935 0.961 0.967 0.988 測站6 0.93 0.928 0.961 0.91 測站8 0.94 0.825 0.97 0.964 表5.2 模擬結果與實驗值均方根差比較 SI=2.13 SI=0.89 彎道內側 彎道外側 彎道內側 彎道外側 測站2 0.015 0.012 0.019 0.011 測站4 0.016 0.016 0.01 0.009 測站6 0.014 0.021 0.006 0.009 測站8 0.035 0.027 0.011 0.011 表5.3 模擬結果與實驗值信賴指標比較 SI=2.13 SI=0.89 彎道內側 彎道外側 彎道內側 彎道外側 測站2 1.062 1.054 1.069 1.046 測站4 1.081 1.059 1.042 1.032 測站6 1.052 1.08 1.028 1.038 測站8 1.147 1.091 1.053 1.05

43 表5.4二次流影響比較湧波深度(Test 21n) 單位：公尺 數值模擬 測站 內外側 數值與誤差 實驗值 考慮二次流 不考慮二次流 湧波深度 0.173 0.157 0.155 彎道內側 誤差比較 9.121% 10.532% 湧波深度 0.173 0.165 0.167 第2 站 2.81秒 彎道外側 誤差比較 4.532% 3.185% 湧波深度 0.161 0.15 0.144 彎道內側 誤差比較 6.689% 10.072% 湧波深度 0.174 0.164 0.164 第4 站 3.79s 彎道外側 誤差比較 6.126% 6.126% 湧波深度 0.153 0.136 0.123 彎道內側 誤差比較 10.876% 19.386% 湧波高程 0.164 0.152 0.135 第6 站 4.78秒 彎道外側 誤差比較 7.293% 17.841% 湧波高程 0.147 0.139 0.128 彎道內側 誤差比較 5.442% 13.218% 湧波高程 0.151 0.139 0.128 第8 站 6.5秒 彎道外側 誤差比較 7.947% 15.503%

表 5.5二次流影響相關係數比較(Test 21n) 考慮二次流 不考慮二次流 彎道內側 彎道外側 彎道內側 彎道外側 測站2 0.97 0.964 0.941 0.939 測站4 0.967 0.988 0.966 0.982 測站6 0.961 0.91 0.952 0.83 測站8 0.97 0.964 0.969 0.95 表5.6 二次流影響均方根差比較(Test 21n) 考慮二次流 不考慮二次流 彎道內側 彎道外側 彎道內側 彎道外側 測站2 0.019 0.011 0.017 0.01 測站4 0.01 0.009 0.011 0.009 測站6 0.006 0.009 0.01 0.014 測站8 0.011 0.011 0.013 0.017 表 5.7二次流影響信賴指標比較(Test 21n) 考慮二次流 不考慮二次流 彎道內側 彎道外側 彎道內側 彎道外側 測站2 1.069 1.046 1.062 1.041 測站4 1.042 1.032 1.043 1.043 測站6 1.028 1.038 1.041 1.069 測站8 1.053 1.05 1.041 1.069

45 表5.8 二次流影響比較湧波深度(Test 21) 單位：公尺 數值模擬 測站 內外側 數值與誤差 實驗值 考慮二次流 不考慮二次流 湧波深度 0.18 0.157 0.16 彎道內側 誤差比較 12.533% 10.983% 湧波深度 0.18 0.164 0.171 第2 站 2.64秒 彎道外側 誤差比較 8.711% 5.011% 湧波深度 0.176 0.157 0.13 彎道內側 誤差比較 10.678% 25.678% 湧波深度 0.212 0.173 0.127 第4 站 3.57秒 彎道外側 誤差比較 22.229% 48.219% 湧波深度 0.169 0.15 0.12 彎道內側 誤差比較 11.321% 29.023% 湧波高程 0.222 0.168 0.086 第6 站 4.46秒 彎道外側 誤差比較 24.333% 61.329% 湧波高程 0.155 0.19 0.166 彎道內側 誤差比較 -22.281% -7.028% 湧波高程 0.195 0.19 0.166 第8 站 6.07秒 彎道外側 誤差比較 2.8% 14.872%