橋墩水深平均二維沖刷模式之發展

(1)

國

立

交

通

大

學

土木工程學系

碩士論文

橋墩水深平均二維沖刷模式之發展

Development of a Depth-Averaged Two-Dimensional

Pier Scouring Model

研究生：陳彥酉

指導教授

：

楊錦釧博士

指導教授

：

謝德勇博士

中華民國一 ○ ○ 年七月

(2)

橋墩水深平均二維沖刷模式之發展

Development of a Depth-Averaged

Two-Dimensional Pier Scouring Model

研究生：陳彥酉 Student：Yan-You Chen 指導教授：楊錦釧 Advisor：Jinn-Chuang Yang 指導教授：謝德勇 Advisor： Te-Yung Hsieh

國立交通大學土木工程研究所碩士論文 A Thesis

Submitted to Department of Civil Engineering College of Engineering

National Chiao Tung University in partial Fulfillment of the Requirements

for the Degree of Master in

Civil Engineering July 2011

Hsinchu, Taiwan, Republic of China

中華民國一 ○ ○年七月

(3)

I

橋墩水深平均二維沖刷模式之發展

學生：陳彥酉指導教授：楊錦釧指導教授：謝德勇國立交通大學土木工程學系碩士班

摘要

本研究旨在發展一橋墩水深平均二維沖刷模式，使水平二維模式除了可反映橋墩兩側的束縮沖刷外，亦可以合理模擬垂向水流在橋墩周圍底床造成的局部沖刷現象。本研究以謝氏(2003)的水深平均二維動床模式為主體，採用浸沒式邊界法處理計算區域中因水流無法穿越橋墩導致的內部邊界問題，並在模式中加入垂向流速之估算式與相對應的局部沖刷輸砂公式。在垂向流速推估部份，首先針對重要影響參數進行無因次分析，再根據參數的合理範圍設計多組模擬案例，以三維 CFX 模式進行三維水理分析，將模擬結果進行迴歸分析後獲得橋墩周圍的最大垂向流速分布之經驗公式。局部沖刷輸砂公式部分，使用多組橋墩沖刷實驗案例進行沖刷參數的檢定分析，並利用檢定結果與對應案例之水理條件進行迴歸得到沖刷參數的經驗式。最後以其他的橋墩沖刷實驗案例進行驗證分析，結果顯示模式能合理模擬橋墩周圍局部沖刷深度變化的歷程。關鍵字:橋墩沖刷、水深平均二維模式、局部沖刷、CFX

(4)

II

Development of a Depth-Averaged Two-Dimensional

Pier Scouring Model

Student：Yan-You Chen Advisor：Jinn-Chuang Yang

-Hsiang Hung Advisor：Te-Yung Hsiehaaa

Department of Civil Engineering National Chiao-Tung University

Abstract

This study develops a depth-averaged two-dimensional pier scouring model to simulate both the contraction scour and the local scour properties due to vertical flow. The framework of the model developed is based on the depth-averaged two-dimensional mobile bed model by Hsieh (2003), in which employed the immersed boundary method is to deal with the internal pier boundary. The maximum vertical flow velocity and the corresponding sediment transport formula are taken into account to simulate the local scour phenomenon. The maximum vertical velocity regression formula is established from the velocity field simulated by a commercial 3-D model, CFX, through several hypothetical cases designed. Several bridge pier scouring experiment results are used to conduct the model calibration and parameter regression analysis to develop the empirical formula of local scour. Furthermore, the model is verified by using experiment results other than those for calibration. The verification studies show that the variation of the scour depth with respect to time around the bridge piers could be well simulated.

Keywords: piers couring, depth-averaged two-dimensional model, local scour, CFX

(5)

III

誌謝

在結束研究所階段的學習，即將踏入人生的下一階段的現在，心中有太多的感謝想要表達。在這段期間，承蒙恩師楊錦釧教授與謝德勇博士的悉心指導，使我得以增長學識並順利完成學位論文，兩位在學術及生活上的積極態度，更是我的學習典範。感謝口試委員盧昭堯教授、許銘熙教授、賴進松博士在論文審查期間給與諸多寶貴的意見，使本論文能夠更加完善。學生在此向各位老師致上萬分的感謝與敬意。在學期間，特別感謝和政學長、胤隆學長、世偉學長、浩榮學長、建華學長、弘恩學長、昇學學長、仲達學長、仁凱學長、俊宏學長、昀軒學長、仁猷學長、明儒學長、鏡如學姐、筱萍學姐、群玲學姐及綺雯學姐在學業、研究以及生活上，給予我的協助與照顧。感謝國網中心的李明龍博士在模式使用上的排疑解惑，使我的研究得以進行。感謝同窗好友聖翔、家偉、緁玲、唯泰、東洲、紹唐、彥瑜，學弟妹建翔、昀直、舒勤與芳綺，豐富了我的研究所生活並與我互相鼓勵一同學習。感謝高中好友信文、弘毅及大學好友博誠、柏傑、旻璋、宇勝，在我迷惘或遭遇困難時聆聽並給予我中懇的建議。感謝家人在這段時間的體諒並給予我鼓勵，感謝爸爸和老哥照顧媽媽，使我可以將心思放在完成我的論文上；感謝爸媽給我的關懷，使我即使遠在新竹仍備感溫暖，感受到家人對我的愛；感謝阿淳，聰慧如妳才能如此包容任性妄為的我，在我最低潮的時候陪伴我並給予我無限的支持。最後，謹將此論文獻給你們，與你們一同分享我的喜悅。

(6)

IV 目錄

目錄

摘要 ... I Abstract ... II 誌謝 ... III 目錄 ... IV 表目錄 ... VII 圖目錄 ... VIII 符號表 ... X 第一章緒論 ... 1 1.1 研究動機與方向 ... 1 1.2 文獻回顧 ... 1 1.2.1 局部流場 ... 2 1.2.2 局部沖刷深度 ... 4 1.3 研究目的與方法 ... 5 1.4 章節簡介 ... 6 第二章理論基礎 ... 9 2.1 水理部分 ... 9 2.1.1 控制方程式 ... 9 2.1.2 輔助關係式 ... 10 2.1.3 水理邊界條件 ... 11 2.2 沉滓運移部分 ... 12 2.2.1 控制方程式 ... 12 2.2.2 輔助關係式 ... 13 2.2.3 沉滓運移邊界條件 ... 17 第三章數值方法 ... 20

(7)

V 3.1 水理部分 ... 20 3.1.1 分割操作趨近法 ... 20 3.1.2 數值差分式 ... 21 3.2 沉滓運移部分 ... 23 第四章橋墩周圍垂向流速推估 ... 26 4.1 前言 ... 26 4.2 ANSYS-CFX 三維模式簡介 ... 26 4.2.1 模式原理簡介 ... 26 4.2.2 模式測試 ... 27 4.3 因次分析 ... 27 4.4 案例設計 ... 28 4.5 模擬結果分析 ... 28 4.5.1 垂向流速估算式 ... 29 4.5.2 適用性分析 ... 29 4.5.3 估算式之驗證 ... 30 第五章模式檢定與驗證 ... 40 5.1 模擬參數敏感度分析 ... 40 5.2 模擬參數之決定與分析 ... 40 5.2.1 懸浮載源項時間係數 γ ... 40 5.2.2 沖刷參數 α ... 41 5.3 模擬參數之驗證分析 ... 42 5.3.1 Yanmaz 試驗案例 ... 42 5.3.2 Dargahi 試驗案例 ... 43 第六章結論與建議 ... 55 6.1 結論 ... 55 6.2 建議 ... 56

(8)

VI

參考文獻 ... 57 附錄一 CFX 模式之數值方法簡介 ... 59

(9)

VII

表目錄

表 4.1 局部流場模擬案例設計一覽表 ... 32 表 4.2 各參數標準化後之迴歸關係表(忽略坡度S₀) ... 32 表 4.3 各參數取自然對數後之迴歸關係表 ... 32 表 5.1 沖刷公式迴歸分析所採用之各實驗案例設置 ... 45 表 5.2 沖刷參數檢定結果 ... 45 表 5.3 沖刷參數檢定結果(續) ... 46 表 5.4 各無因次參數對沖刷參數影響程度之迴歸關係表 ... 46

(10)

VIII

圖目錄

圖 1.1 橋墩周圍流場示意圖(Melville 2000) ... 7 圖 1.2 墩前控制體積示意圖(Shen et al. 1969) ... 7 圖 1.3 清水沖刷與含滓沖刷之沖刷歷程(Raudkivi 1986) ... 8 圖 2.1 浸沒式邊界法示意圖 ... 18 圖 2.2 狄拉克脈衝函數 4 點法 ... 18 圖 2.3 作用層示意圖(Spasojevic 1988) ... 19 圖 3.1 控制體積法示意圖 (a)實際區域；(b)計算區域 ... 25 圖 4.1 墩前垂向流速沿水深方向之分布情形(Raudkivi 1986) ... 33 圖 4.2 VOF 法體積分率函數示意圖 ... 33 圖 4.3 CFX 控制體積示意圖 ... 34 圖 4.4 定床實驗案例格網配置圖(Ahmed 案例) ... 34 圖 4.5 垂向流速實驗值與 CFX 模擬結果之比較 ... 35 圖 4.6 距離因子 (xr r) 與角度因子(

  

 ) / 示意圖 ... 36 圖 4.7 無因次垂向流速w U 估算值與模擬值之比較_m ₀ ... 36 圖 4.8 平床時刻垂向流速w 之修正示意圖m₀ ... 37

圖 4.9 墩前渦流最大切線速度v_b_max(Nakagawa & Suzuki 1974) ... 37

圖 4.10 無因次垂向流速w U 修正結果_m ₀ ... 38 圖 4.11 無因次垂向流速w U 估算結果與實驗值之比較_m ₀ ... 39 圖 5.1 垂直射流沖刷示意圖(Clarke 1962) ... 47 圖 5.2 沖刷參數 與時間係數 對沖刷歷程影響之比較 ... 47 圖 5.3 射流沖刷歷程與射流強度之關係 ... 48 圖 5.4 時間係數 值與垂向流速強度之關係 ... 48 圖 5.5 沖刷參數之檢定(Yanmaz 案例) ... 49 圖 5.6 沖刷參數之檢定(Chang 案例) ... 49 圖 5.7 沖刷參數之檢定(Nakagawa 案例) ... 50

(11)

IX 圖 5.8 沖刷參數與各無因次參數之關係 ... 50 圖 5.9 水深對計算值之影響 ... 51 圖 5.10 檢定值與計算值之比較 ... 51 圖 5.11 動床模擬格網配置圖(Yanmaz 案例) ... 52 圖 5.12 墩前沖刷歷程模擬結果與 Yanmaz 實驗值之比較(Run3) ... 52 圖 5.13 墩前沖刷歷程模擬結果與 Yanmaz 實驗值之比較(Run18) ... 52 圖 5.14 動床模擬格網配置圖(Dargahi 案例) ... 53 圖 5.15 墩前沖刷歷程模擬結果與 Dargahi 實驗值之比較 ... 53 圖 5.16 墩前沖刷坑剖面與 Dargahi 實驗值之比較 ... 54

(12)

X 符號表 A渦流截面積； B 渠道寬度； C  濃度； B 摩擦係數； c Chezy 係數； 1 c  顆粒蔡司係數； D橋墩直徑； k D 顆粒 k 之粒徑； m D  不產生移動之最小顆粒粒徑； k D_  無因次顆粒粒徑； d  水深； s d  底床沖刷深度； se d 平衡沖刷深度； 0 d  上游來流水深； m E  作用層厚度； rms E 均方根誤差； p Fr  橋墩福祿數； g 重力加速度； H 平均水深； 1 h 、h2 



、 方向轉換係數；



 Von Karman’s 係數； n 曼寧糙度係數； p孔隙率 i b q  i 方向某一粒徑之河床載通量； R 水力半徑；

(13)

XI p Re  橋墩雷諾數 S  懸浮載源； s 砂比重； f S  作用層源； v S  垂向流場導致之懸浮載源； 0 S  渠道坡度； k T  輸送參數； 11 T 、T 、12 T22 有效剪應力項； t 時間； f t  實驗或模擬結束時間； U  方向平均速度；



u  方向速度；



s U  遠方之表面流速； 0 U  上游平均來流速度； u_  剪力速度； u_c 臨界剪力速度； V  方向平均速度；



v  方向速度；  max b v  墩前最大切線速度； lk w  顆粒 k 之躍起速度； fk w  顆粒 k 之沉降速度； m w  最大垂向流速； 0 m w  平床時刻之最大垂向流速； b z  底床高程； s z  水面高程；   沖刷參數；

(14)

XII



 粒徑百分比；



懸浮載源項時間係數； 1  _、₂  、



方向之亂流傳輸係數； l   層流黏滯係數； t   紊流黏滯係數；   流體密度； s   泥砂密度； 1 b  _、 2 b   底床剪應力在



、方向之分量；



、  平面上兩正交座標方向；   隱藏因子； t   時間間距；



 、 



、方向之格網間距。上標 n  時刻之已知變數； n t 1 ( 1) n   n t時刻之未知變數； 1 2 ( 1) n  n t與n t 間之未之變數；

 

時間平均；

 

水深平均；

 

 時間平均瞬時擾動量。

(15)

1 1.

第一章

緒論

1.1 研究動機與方向 橋梁連結了河岸的兩端，是公路運輸之重要命脈，而河川底床沖淤本是一動態平衡之現象，由於橋梁等結構物所造成的斷面縮減使流場發生變化，應力集中的結果對周圍底床造成過度沖刷。橋墩附近可能產生的沖刷形態，除了局部應力集中在橋墩兩側產生的束縮沖刷外，墩前壅水、下向流、馬蹄形渦流等複雜的紊流流場，亦會對底床的沖淤平衡產生影響，甚至是在橋墩周圍產生較大沖刷量的主因。因此如何有效地計算局部流場所造成之最大沖刷深度，作為建造橋梁的設計參考，或在颱洪來臨之際提供即時的防災預警之用，是諸多水利學者與相關領域人員努力的目標。以往在橋墩沖刷的研究，多以物理實驗來探究沖刷深度與各項水理參數之間的關係，衍生了許多的沖刷深度經驗公式。近年來電腦演算效率的提升，以水理數值模式搭配泥砂運移模式進行橋墩周圍的底床沖刷模擬之研究得以發展蓬勃，並進一步應用於實務工程參考。在數值模式應用方面，三維水理模式求解完整的動量方程式及連續方程式，結果相對準確但求解時間冗長；而二維模式藉由特定假設條件，簡化了控制方程式，求解效率相對較快且維持了一定的精確程度。本研究為了上述之實用目的，在既有之二維動床模式架構中，加入水流在垂直方向的流速估算及垂向流場對底床的沖刷效應，使模式具備有模擬局部沖刷現象之能力，期望能更合理且迅速地表現出局部水流變化對橋墩周圍底床之影響。 1.2 文獻回顧 橋墩沖刷之研究由來已久，國內外學者通常針對不同的橋墩因素(尺寸、形狀、數目、形狀均勻與否)、水流因素(流速、水深)及底床性質(泥砂粒徑、顆粒級配)造成的流場改變以及底床沖淤變化進行相關探討。根據成因的不同，橋墩周圍可能產生的沖刷種類分為: (1) 一般沖刷(general scour):不論結構物存在與否，單純以水流之驅動使得泥砂被沖刷而移動至下游，並且淤積在下游某處，沖淤交替使得整體河床一直處於動態之變化。此一型式之沖刷可由一維以上維度之動床模式等進行模擬； (2) 束縮(contraction)沖刷:建造在渠道中的結構物使通水斷面縮減而發生迴水效

(16)

2 應，局部流速增加及水流應力集中亦使附近底床產生較大之沖刷；此類型沖刷受到側向水流之影響，需要擬似二維或二維以上之動床模式方可處理。 (3) 局部沖刷(local scour):結構物的出現阻礙了水流前進，水流在其附近形成紊流流場，對底床質造成擾動後將其淘起帶往下游，使底床高程發生變化。此類型沖刷包含了三維之垂向流場對底床之影響，因而一般靜水壓假設之水理模式無法對此進行模擬，需要使用三維水理模式或擬似三維水理模式搭配動床模式處理之。由於橋墩局部沖刷深受周圍流場變化之影響，因此以下將分別就橋墩處局部流場與局部沖刷深度兩部分進行相關研究文獻的蒐集與歸納。 1.2.1 局部流場 橋墩出現造成了壓力梯度變化及水流分離現象，並在橋墩周圍形成局部流場。如圖 1.1 所示，包含墩前的墩前壅水(bow wave)及下向流(downward flow)、墩側的馬蹄形渦流(horeshoe vortex)以及墩後的尾跡渦流(wake flow)。而根據前人之研究，墩前下向流與馬蹄形渦流乃是造成橋墩周圍底床變化之主要因素(Shen et al 1969；Melville 1975；Baker 1980)。針對橋墩前方形成的下向渦流流場，Shen et al.(1969)提出墩前單一渦流的簡化模型概念。假設橋墩面與上游入流點的距離為無窮遠時，在墩前控制體積 ABCD 區域內(如圖 1.2)，水流將在墩前形成順時針朝向底床之強制渦流(force vortex)。利用斯托克斯定律(Stokes’ Theorem)求得橋墩出現所造成的環流量變化，並推論渦流核心強度與橋墩雷諾數(Re_p U D₀ ) 兩者具有函數關係。 Nakagawa and Suzuki(1974)延續 Shen et al.(1969)的理論，藉由實驗數據得到平床

時刻的墩前渦流最大切線速度，實驗結果顯示切線速度約為 0.5 倍上游平均流速。Raudkivi(1986)之實驗結果顯示墩前垂向流速在平床時刻為來流速度的 40%，在發生最大沖刷深度時可達 80%。Melville(1975)針對沖刷坑發展前中後時刻之速度場及渦度變化進行量測，發現在沖刷坑發展過程中，其坑內之水流速度會隨著沖刷坑之發展而增加，最大可達來流速度的 80%；且沖刷平衡時刻之馬蹄形渦流大小及坑內環流量變化受制於沖刷坑幾何形狀，與平床時刻差異相當

(17)

3

大。Baker(1980)則提出不同之看法，沿用 Shen et al.(1969)之理論模型，他推論渦流半徑將隨著沖刷深度加深而增大，但渦度將隨著沖刷坑面積增加而消散減少，此消彼長則系統中的環流量在沖刷過程中將維持恆定   ，並根據此一₀ 假設得到渦流切線速度與沖刷深度之變化關係。Ahmed(1995)設計數組橋墩沖刷試驗，分別對平滑定床、粗糙定床以及動床進行量測，在動床試驗沖刷過程中，墩前下向流之最大流速可達上游速度的95%。

Graf et al. (1998)使用聲波都卜勒流速儀(Acoustic Doppler Velocity Profiler)對 橋墩周圍定床流場進行水平二維流場 u 、 v 及垂向流速 w之量測，量測點位置與 墩前水流方向之夾角 分別為 0°、45°、90°及 157.5°，由流速剖面觀察到二維流場水流分離點(separation point)大致在夾角 135°與 157.5°之間，在這之後便進入尾流區(wake region)；而垂向流場的水流分離點與橋墩中心的距離則隨著角度

 往下游增加而變大，即渦流半徑越往下游越大。Kent and Francis(2003)使用質點影像測速儀(Particle Image Velocimetry；PIV)對圓形橋墩周圍流場進行觀測並分析量測結果，認定橋墩周圍水位之變化與橋墩福祿數Fr 有關，但是否對橋墩_p 周圍速度場造成影響則無法定論。在數值模擬方面，針對橋墩周圍之流場，二維模式可以透過一些數值手段處理內部固體邊界之問題，如有限元素及有限體積法使用非結構性格網以避開內部不透水之邊界。Peskin(1972)則提出浸沒邊界法，將內部固體邊界視為流體的一部分，以作用力與反作用力之概念對物理量進行重新分配。黃氏(2009)曾以此方法進行橋墩及丁壩附近之流場模擬，此法可以合理地模擬在橋墩周圍的縱向、側向流場。然而針對橋墩周圍的垂向流場，二維模式不具有求解能力而無法模擬；三維模式的數值模擬研究方面，針對橋墩周圍之流場，Olsen and Melaaen(1993) 採用 k 紊流模式，對不可壓縮流穩態流場之質量守恆方程式與動量方程式求 解，進行橋墩周圍流場之模擬；曾氏(1994) 採用大渦紊流模擬 (Large eddy simulation；LES)處理大尺度渦流之動量傳輸，而受流場中流體黏滯性影響之小尺度之渦流則以簡化的紊流模式 Smargorinsky(SGS)模式處理。在實務應用上，現今已有多套流體力學商用模式可供選用，如 FLUENT、FLOW-3D 及 ANSYS-CFX 等。其中本研究擬採用之三維商用模式 ANSYS-ANSYS-CFX 使用了有限體積法，

(18)

4

結合有限元素之概念進行格網點間的離散，求解質量守恆方程式與動量方程式。 CFX 模式中包含了多種紊流模式，且具有穩定的多相流模式，使用容積分率法 (Volume of fluid method；VOF)處理兩相之間的質量、動量以及能量傳遞問題，在國內外之研究及應用上皆有相當程度的成果。

1.2.2 局部沖刷深度

局部沖刷可進一步分為清水沖刷(clear-water scour) 及含滓沖刷 (lived-bed scour)兩種型態(Chabert and Engeldinger 1956)。Laursen(1956)以單位時間內沖刷坑的體積變化率d V dt 來定義底床沖刷的型態: d V 1 2 s s q q dt   (1.1) (1)水流在單位時間內將泥砂帶出沖刷坑的運移量q 等於由上游補充進入坑內的s₁ 運移量q ，代表底床無沖刷；s₂ (2) qs1  且0 0qs2 qs1 為清水沖刷；(3) 1 2 0 s s q q  則為含滓沖刷。根據前人研究，清水沖刷相對於含滓沖刷，具有兩項顯著特性 (1)清水沖刷之平衡沖刷深度(equilibrium scour depth)較含滓沖刷大，Shen et al.(1969)認為含滓沖刷的最大沖刷深度約比清水沖刷小 10%；(2)清水沖刷較含滓沖刷需要較多的時間來達到平衡沖刷深度(如圖 1.3)(Shen 1969；Raudkivi 1983；Melville and Chiew 1999)。綜合以上所述，並考量沖刷程度對橋墩安全之影響，國內外橋墩沖刷研究大多對清水沖刷之情況做探討。橋墩局部沖刷深度會受到周圍水流強度的影響，當水流速度大於顆粒臨界起動速度，底床質被水流帶起運移至下游，沖刷坑深度將隨著時間漸漸加深。由於沖刷坑截面積擴大，Baker(1980)相信這將造成渦流渦度的擴散而導致坑內水流速度的下降，進而使坑內環流量維持定值；Melville(1975)則認為馬蹄型渦流的環流量在沖刷過程中會隨著沖刷坑發展而增加，其沖刷深度變化率及墩前下向流速度則會隨著時間漸漸變小。橋墩周圍沖刷行為與時間之演進有密切關係，根據 Melville and Chiew(1999)的定義，若 24 小時內產生的沖刷深度變化未達橋墩直

徑的 5%，則認定此時刻沖刷達到平衡；其實驗結果顯示大約在所需平衡時間的

(19)

5 驗中，局部沖刷在實驗開始 12 小時之後達至平衡，且在實驗開始後 2 小時內即達成最大沖刷深度的60%。由於泥砂運移之機制太過於複雜，前人在動床模擬之研究多以泥砂傳輸經驗式搭配質量守恆方程式求解底床變化，而不同的經驗式有其適用限制，採用不同的經驗公式會對模擬結果造成很大的差異。Shen et al.(1969)根據前人之實驗數據，得到平衡沖刷深度與橋墩雷諾數之關係式。Baker(1980)延續 Shen(1969)之理論及假設，藉由實驗數據推導出其沖刷公式為雙曲函數之型式。Kothyari et al. (1992)探討含滓沖刷情況下墩前最大沖刷深度之歷程，由實驗得到墩前渦流的尺寸與橋墩尺寸之關係，藉由底床剪力得到渦流尺寸隨時間之變化。張氏(1998)曾改寫 Clarke(1962)垂直射流沖刷公式，將墩前下向水流對底床之沖刷視同垂直射流之效應，提出墩前下向流沖刷公式。Clarke(1962)之沖刷公式乃是基於深水假設，忽略水深對底床沖刷程度之影響，且射流強度在作用時間中維持定值，不同於橋墩周圍流場隨時間以及床形變化之情況。張氏(1998)曾以作用時間與沖刷深度之關係，對時間項指數係數進行修正。以上兩位學者皆未在公式中考慮水深之影響，然而根據 Raudkivi and Ettema(1983)之研究，水深 d 與橋墩直徑D之比值

d D 在小於 3~4 之間時，墩前下向流及墩側馬蹄形渦流會與反向之壅升水流產

生交互作用，使得底床沖刷之程度受到影響；Melville and Sutherland(1988)則認 為，無因次水深 d D 小於1.43 時，會影響底床沖刷。由於局部流場對於底床之沖淤變化影響很大，尤其以垂向流速分量對底床造成之影響最為直接，因此垂向流場資訊對於局部沖刷尤其重要。回顧前人之文獻時，發現前人之實驗量測對於此方面資訊較為缺乏，無法找到足夠的數據用來決定垂向流場在空間上的變化趨勢。因此，本研究決定使用三維模式進行一系列之流場模擬，藉此獲得合理的垂向流場推估經驗式之後，在二維模式中加入垂向流場所造成之局部沖刷效應，適切地表現水流對底床所造成的沉滓交換機制，使其將底床變動之資訊反覆疊代的同時仍保有二維模式在計算速率上的優勢。 1.3 研究目的與方法 以謝氏(2003)所發展的二維水深平均 RESED2D 動床模式為主體，內部邊界以浸沒式邊界法處理(黃氏，2009)；垂向流場資訊方面，本研究擬以三維 ANSYS-CFX 模式為工具，進行多組案例之流場模擬並對結果進行統計迴歸得到垂向流場在空間上分布之經驗公式，提供二維模式做為推估垂向流場之理論依

(20)

6 據。動床模擬部份，在輸砂公式加入垂向流場對底床產生之局部沖刷效應，並針對 Clarke(1962)之沖刷公式進行修改，增加水深對底床沖刷之影響，最後以實驗案例對相關沖刷參數進行檢定及驗證。 1.4 章節簡介 本章前幾節已就本研究之動機與方向、文獻回顧、研究目的與方法進行說明，本節將簡單描述各章節之內容。第一章為緒論，即對本研究之動機與方向做說明，並回顧前人相關之研究文獻，進而提出本研究之目的與方法，最後介紹各章節之內容架構。第二章為理論基礎，將分別就正交曲線座標下之基本控制方程式、輔助關係式以及邊界條件，分為水理及沉滓運移兩大部份進行說明。第三章為數值方法，有關模式中水理及沉滓運移兩部份之控制方程式所採用的數值處理方式，皆在本章做說明。第四章為橋墩周圍垂向流速推估，首先對所使用之三維模擬工具做簡介，爾後進行數組案例之流場模擬，接著說明對模擬結果進行迴歸分析後得到垂向流速推估公式之過程，並以定床實驗案例結果與計算值進行比較。第五章為模式檢定與驗證，本章首先將對模式中待定係數做敏感度分析與檢定，然後以實驗數據來決定沖刷公式中的各項參係數，最後就整體模式進行驗證分析。第六章為結論與建議，除了對本研究之整體成果做一歸納說明與討論外，並對不盡完備之處提出建議以供日後改進與後續研究之參考。

(21)

7 圖 1.1 橋墩周圍流場示意圖(Melville 2000) 圖 1.2 墩前控制體積示意圖(Shen et al. 1969)

(22)

8 圖 1.3 清水沖刷與含滓沖刷之沖刷歷程(Raudkivi 1986) s d se d live-bed scour clear-water scour

(23)

9 2.

第二章

理論基礎

模式透過座標轉換將基本控制方程式轉換為正交座標系統方程式，並對控制方程式取時間平均與水深平均後，即二維正交曲線座標模式所需之控制方程式。其中水理控制方程式包括水流連續方程式及水流動量方程式；而沉滓運移部份之控制方程式，則將輸砂通量區分為非均勻質之懸浮載與底床載，控制方程式包含了某一粒徑懸浮載之質量守恆方程式、某一粒徑於作用層(active layer)之質量守恆方程式以及整體河床沉滓運移之質量守恆方程式。其中，為反映底床之局部沖刷現象，本研究將對控制方程式中懸浮載源項進行修正。以下將就水理及沉滓運移部份之理論基礎進行說明: 2.1 水理部分 2.1.1 控制方程式 為求控制方程式之簡化，需對數值模式做出若干假設，分別為:(1)流體為不可壓縮之牛頓流體；(2)水壓為靜水壓分布；(3)忽略風剪力之效應；(4)忽略地球自轉之科氏力效應。根據以上假設進行座標轉換可得正交曲線座標下之水理控制方程式: (1) 水流連續方程式 1 2 ( 2 ) ( 1 ) 0 d h h h ud h vd t    _  _  _    (2.1) (2) 水流動量方程式 ξ方向： 1 2 1 2 1 2 1 2 1 2 1 1 2 1 11 1 12 12 22 1 2 1 2 1 2 1 2 2 11 2 11 1 12 1 12 1 2 1 1 ( ) ( ) 1 1 1 1 ( ) ( ) 1 ( ) ( ) ( ) ( ) b b s b s s b s u u u v u h h g uv v f z d t h h h h h h h h h h T h T T T h h d h h d h h d h h d z z z h h h h d h h d                          _  _  _  _  _ _  _                               b b z      _    (2.2) η方向：

(24)

10 2 2 2 1 1 2 1 2 1 2 2 1 2 2 12 1 22 11 12 1 2 1 2 1 2 1 2 2 12 2 12 1 22 1 22 1 2 1 1 ( ) ( ) 1 1 1 1 ( ) ( ) 1 ( ) ( ) ( ) ( ) b b _s _b _s _b s b s b v u v v v h h g uv u f z d t h h h h h h h h h h T h T T T h h d h h d h h d h h d z z z z h h h h d 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  底床高程；z_s  水面高程； i b   底床剪應力在



與 方向之分量；

 

 時間平均；

 

 水深平均；

 

  時間平均瞬時擾動量； 下標 s 、 b 分別代表變數在水面與底床的數值；T 、₁₁ T 、₁₂ T₂₂  有效剪應力項

(effective stress term)，包含層流剪應力、亂流剪應力與延散剪應力(dispersion stresses)，由於本研究乃針對直線渠道，在此暫不考慮延散剪應力項； f( )



及

( )

f



分別為



方向與方向的虛擬反饋力。 2.1.2 輔助關係式

(1) 底床剪應力

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

1 2 2 1/2 ( ) b Cf u u v     (2.7) 1 2 2 1/2 ( ) b Cf v u v     (2.8) 其中C_f g c/ 2為摩擦係數；C chezy係數 (2) 層流與亂流剪應力模式採用Boussineq 之渦流黏性理論，將層流與亂流剪應力合併表示為

(25)

11 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)；u  剪力速度； k  von Karman’s 係數(約等於 0.4)

2.1.3 水理邊界條件 本研究考量渠道入流、渠道出流以及固體邊界三種邊界型式，其中渠道入流邊界設定為單位寬度入流量，渠道出流邊界則為水位高程。固體邊界則分為外部邊界及內部邊界兩種型式: (1) 外部固體邊界外部邊界採用主流方向流速為零之滑移條件(slip condition)或速度梯度為零之非滑移條件(non-slip condition)處理。 (2) 內部固體邊界對於渠道內部固體邊界如橋墩，本研究使用浸沒式邊界法(黃，2008)來模擬內部邊界所造成的局部流場集中現象。如圖 2.1 所示，內部邊界格網點為 Lagrangian 格點，結構格網部分的則為 Eulerian 格網點。實際上，內部邊界並不具有實體，浸沒式邊界法的原理，是使用虛擬力(virtual force)的概念，來近似流體撞擊固體邊界時，固體反饋給周圍流場的反作用力。計算虛擬力 f( )



及 f( )



的方式是假設內部固體邊界為非滑移，也就是邊界上流速為零，利用邊界與周圍流場的流速差，計算兩者間的反作用力。虛擬力與反作用力的關係可表示為 1 ( ) ( ) ( ) m m M i i h i i k f x F X  x X s  



  (2.12) ( ) ( ) ( ) h xi Xi d xh  X d y h  Y      (2.13)

(26)

12 2 2 1 (3 2 1 4 4( ) ) 8 1 (3 2 7 12 4( ) ) 2 8 0 i i i i i i h i i i i i r r r r h h h h h r r r d r h h h h h otherwise          _           (2.14) 上方各式中， f x 為( )_i Eulerian force； ( ) m i F X 為Lagrangian force；i 下標為_m



方向與方向的變數；M 為標記點總數； s 為標記點之間間距；_h為狄拉克脈衝

函數(Lai and Peskin 2000)；h 為尤拉正交格網點間距。_i (2.14)式為狄拉克脈衝函

數 4 點法，計算標記點與格網之間權重時，在



方向與方向往標記點上下左右方向四個格網點的距離做計算如圖 2.2。 Eulerian 流速與 Lagrangian 流速之關係如下 1 2 ( ) ( ) ( ) m m i i h i i i U X u x  x X  h h _(2.15)

本模式以下列步驟進行Lagrangian force 之求解(參考 Su et al. 2007)

1 2 ( ) ( ) ( ) ( ) m m m i i i h i i x U X U X f x x X h h t   _     



(2.16) 將(2.13)式帶入後移項整理 1 2 1 ( ) ( ) ( ) ( ) ( ) m m m m m M i i i i h i i i k x U X U X x X x X s h h F X t     _ _  _      _      (2.17)

求解上式即可得Lagrangian force，並將其代回(2.12)式即可得到 Eulerian force。

2.2 沉滓運移部分 2.2.1 控制方程式 依照運移方式的不同，輸砂通量可分為非均勻質之懸浮載與底床載。因此模式中的控制方程式包含了某一粒徑懸浮載質量守恆方程式、某一粒徑於作用層之質量守恆方程式以及整體河床中泥砂運移之質量守恆方程式，茲說明如下 (1) 某一粒徑懸浮載之質量守恆方程式

(27)

13 1 1 1 2 1 2 1 2 2 1 2 2 1 1 dC u dC v dC d h dC d h dC S dt h d h d h h d d  h d h h d d  h d d        _ _ _ _     (2.18) (2) 某一粒徑於作用層之質量守恆方程式 1 2 2 1 1 2 ( ) ( ) ( ) (1 ) m b b 0 s f h q h q E p h h S S t                 (2.19) (3) 整體河床泥砂運移之質量守恆方程式 1 2 2 1 1 2 ( ) ( ) (1 ) _b b b 0 s h q h q z p h h S t           _   _ 



_   _ (2.20) 以上式子中，_s  泥砂密度；



 粒徑百分比；p 孔隙率；E_m  作用層厚度； i b q  i 方向某一粒徑之河床載通量； S  懸浮載源；S_f 

作用層源(active-layer floor source)。 2.2.2 輔助關係式 (1) 亂流傳輸係數₁、₂ 亂流傳輸係數採用Elder(1959)之經驗公式表示 1 5.93u d 0.232 u d   _ ，  _ (2.21) (2) 河床載通量q_b 凝聚性泥砂之通量可假定為零，而非凝聚性沉滓之河床載通量，模式中採用 van Rijn(1984a)之輸砂經驗式計算(以



方向為例): 1 1 2.1 0.3 ( ) 0.053 ( 1) k k b b k s k k T q q D s gD D D      _(2.22) 式中， 1 3 2 ( 1) k k s g D D       _ _    無因次顆粒粒徑；Dk  顆粒粒徑； 2 2 2 ( ) ( ) c c k k k u u T u       輸砂參數； 1 c u g u c    臨界剪力速度； 1 90 12 18log( ) 3 d c D   顆粒蔡司參數； s s    砂比重。

(28)

14 另外，泥砂運移過程中，較細之顆粒除了以懸浮載型式移動外，由於泥砂為非均質，細顆粒可能隱藏於較粗顆粒之間而不易被水流帶走。為體現此一機制，模式進一步將非均勻粒徑所產生的遮蔽效應納入考慮並進行修正如下 ( ) i k i b h k b k q   q D (2.23) 由於水流在渠道中運行時，底床載源之變化將受到:(1)縱向與橫向之底床坡降S₀ 、S0；(2)縱向及橫向之流速 u 、 v 的影響；模式採用 Struiksma(1985)之公式，對某一粒徑之底床載通量進行修正(張氏，2005)



方向: qb  hk kq Db( ) cosk  (2.24) 方向: qb  hk kq Db( ) sink  (2.25) 在上面兩式中 1 1 sin tan 1 cos b s b s z f z f           _   _     _  _   _    (2.26) 1 tan v u  _        (2.27) 以上，  沉滓運移角度；  底床剪應力之方向； f_s  沉滓之形狀因子， 1 f_s  ；2 2 2 2 50 ( ) ( 1) u v c s D      Shields parameter。 (3) 懸浮載源 S 針對非凝聚性之泥砂而言，懸浮載源可定義為懸浮載向下通量與底床質向上 通量之間交互作用的結果。對某一特定粒徑 k 而言，受到重力作用而沉降至底 床，其向下之通量可表示為下式 k k k d f d S  w C _(2.28) 其中負號代表通量向下，反之為向上之通量； [3.25 0.55ln( k )] k f d k w C C u  _   (Lin 1984)；C_k  顆粒 k 的深度平均濃度； k f w  顆粒 k 的沉降速度，其計算方式依照 粒徑大小採用不同關係式(van Rijn 1984b)如下

(29)

15





2 0.5 ( 1) 1 ( 100 ) 18 0.01( 1) 10 1 1 (100 1000 ) 1.1 1 ( 1000 ) k k k k f k k k k s gD D m s gD w m D m D s gD D m    _ _     _         _ __  _ _         _ _ _ _     (2.29) 在水中未設有結構物的渠道中，底床質向上之通量主要是由於水流剪力對底床之作用，使底床質揚起成為懸浮質，即一般沖刷及束縮沖刷之現象。其向上通量可表示為 k k k e l k e S  w C _(2.30) 上式 1.5 0.3 0.015 k k k k e D T C a D  (van Rijn 1984b)； a砂丘高度之一半； k l w  顆粒 k 躍起 之速度。顆粒躍起速度定義為河床質發生跳躍(saltation)離開底床瞬間時的垂直速度，採用Hu and Hui(1996)提出之經驗公式 3.2 4.5log 1.2 3.1 1.2 k l w u_        _{ }  (2.31) 式中 ( ) b s gDk       (4) 橋墩局部沖刷懸浮載源S_v 受到橋墩周圍局部流場如墩前下向流與馬蹄型渦流之影響，底床質將被淘起成為懸浮載之型式往下游移動，此一沉滓交換機制即橋墩附近發生局部沖刷之緣由，為使模式能夠合理地表現此一現象，本研究加入垂向流場之沖刷效應。若將此沖刷機制視為近似垂直射流之作用，則可運用 Clarke(1962)垂直射流沖刷公式進行沖刷深度之計算 0.43 0.05 0.05 5.5 =(0.21 0.003)D jet jet c u u f f s c U U D gt D gD w w d  _ _{ } _{ } _  _ _ _{ } _{ } _        _ _ _ _{ } _    _  (2.32)

(30)

16 上式中，D_c=沖刷坑之寬度；D_u=射流出口直徑；U_jet=射流出口速度；d_s  沖刷深度。張氏(1998)進一步將(2.32)式改寫為橋墩局部沖刷之經驗式並移項整理如下 0.48 0 0 s m d w U t D U D         _ _   (2.33) 在上式中，；為待定係數，此乃公式推導過程中所產生之參數，用以代表垂向流場之沖刷能力，其大小是由上游平均流速U 、重力加速度₀ g、沉滓落淤速度w 及橋墩直徑_f D決定； 則為時間項之指數係數。將(2.33)式對時間微分(張文鎰 2002)，在此視垂向流場導致之底床沖刷變化為一附加之懸浮載源項， 0.48 ₁ 0 0 0 k k m v k k k w U t S U U D       _{ } _  _ _{ } _     (2.34) 式中S k v  某一粒徑區間 k 在垂向流場作用下形成的沉滓交換率，正值表示通量 向上；U₀  平均來流速度； D為橋墩之直徑；w_m  最大垂向流速；k與k為 相對應不同粒徑 k 之待定係數。 結合(2.28)、(2.30)及(2.34)式，相對應某一粒徑 k 之懸浮載源可以表示為 k k k k e v d S S S S _(2.35) (5) 作用層源S _f 作用層示意圖如圖 2.3。作用層源之產生乃肇因於母層(active stratum)頂面之升降，當其下降時， (1 ) [( ) ( )] f s s k b m S p z E t         (2.36) 式中( )s k  母層中某一顆粒 k 之粒徑百分比；若母層厚度增加，即其頂面上升時，則將(2.36)式中 ( )s k以k取代。 (6) 作用層厚度E_m

(31)

17 1 ( n n) m em b b E  C z  z (2.37) 式中C 為一數值參數，模式暫取為_em 20。當河床表面接近護甲條件時，作用層厚度接近零，在這種情況下，可用 Borah 等(1982)所提出護甲層之厚度，予以修正: 1 1 ( )+ 1 n n m m b b K k m k D E C z z p       



(2.38) 式中，D_m  不產生移動的最小顆粒粒徑。另外，作用層在淤積期間可定義為: 1 ₍ 1 ₎ n n n n m m b b E  E  z  z (2.39) 2.2.3 沉滓運移邊界條件 沉滓運移求解的為之變數為深度平均懸浮載濃度 C 、粒徑百分比



及河床高程z 。渠道入流邊界設定為入流濃度分布，渠道出流邊界條件為_b    。在C  0 固體邊界處，由於控制方程式並不具有動量方程式，因此動床子模式沿用謝氏 (2003)之概念，邊界條件設定為   。C  0



及z 之入流邊界條件則以已知值b 設定，由於下游邊界坡度通常較為平緩，邊界條件可以設為   與  0 0 b z     ，在固體邊界部分同樣設定為   與  0    。 zb  0

(32)

18 圖 2.1 浸沒式邊界法示意圖圖 2.2 狄拉克脈衝函數 4 點法

(33)

19

(34)

20 3.

第三章

數值方法

3.1 水理部分 3.1.1 分割操作趨近法 本研究使用分割操作之概念，將求解動量方程式之過程分成延散、傳播以及浸沒邊界求解三個步驟。其中延散步驟求解移流項以及擴散項，傳播步驟則求解壓力項、底床剪應力項及連續方程式，浸沒邊界求解步驟則進行內部邊界虛擬力之求解。根據以上所述將水理控制方程式改寫為延散步驟 1 1 1 2 2 1 2 ( ) n n n n V V V T t        _{ } _ _{  }  _    (3.1) 傳播步驟 1 2 ( ) n n n d b V V g z d t t d          _  _{  } _ _  _   _      (3.2) 0 n V     (3.3) 浸沒邊界求解步驟 1 n n n V V f t t       _  _  _   _      (3.4) 上方各方程式中，V 代表流速在各方向之分量；T 代表擴散項及延散項；上標 1 n 表示在(n 1) t時刻之未知變數； n 1 n t t  t    ；上標 n 代表在 n t 時刻之已知變數；上標n1 2及n 分別代表求解傳播與浸沒邊界步驟時，在 n t 與 (n 1) t時刻之間的未知變數。將(3.1)至(3.3)之一般式表示如下延散步驟

(35)

21 2 1 2 1 2 1 2 1 2 1 11 1 12 12 22 1 2 1 2 1 2 1 2 2 11 2 11 1 12 1 12 1 2 1 1 1 1 1 ( ) ( ) 1 ( ) s ( ) b ( ) s ( ) b s b s b u u u v u h h uv v t h h h h h h h T h T T T h h d h h d h h d h h d z z z z h h h h h h d                                 _  _    _  _                    _    _       (3.5) 2 2 1 1 2 1 2 1 2 2 12 1 22 11 12 1 2 1 2 1 2 1 2 2 12 2 12 1 22 1 22 1 2 1 1 1 1 1 ( ) ( ) 1 ( ) s ( ) b ( ) s ( ) b s b s b v u v v v h h uv u t h h h h h h h T h T T T h h d h h d h h d h h d z z z z h h h h h h d                                 _  _    _  _                    _    _       (3.6) 傳播步驟 2 2 1 (z_b d) C u uf v u g t h  d      _{ } _    _  _ (3.7) 2 2 1 (z_b d) C u uf v u g t h  d      _{ } _    _  _ (3.8) 1 2 ( 2 ) ( 1 ) 0 d h h h ud h vd t    _  _  _    (3.9) 對n 時刻之水深1 dn1做線性化處理，剔除高階項並保留一階項，則(3.9)式可改寫為 1 2 1 1 1 2 2 2 ( ) ( ) 0 d d d h h d d t                _    _{  } _    _{  } _       _  _  _  _ (3.10) 式中 2 1 1 n h g t d C h_     ； 1 1 2 2 2 1 1 [ ] n n n b z h h g t d u C_ C h_              ； 1 1 n d   ； 1 2 2 n h g t d C h_     ； 1 1 1 2 1 2 2 [ ] n n n b z h h g t d v C_ C h_              ； 2 2 n d   ； 1 1 2 2 2 2 ( ) ( ) 1 n n f n C u v C t d        ； n 1 n d d  d    。 3.1.2 數值差分式

(36)

22 延散及傳播兩步驟採用隱式數值方法求解。模式採用控制體積(control volume)法之觀念來離散控制方程式，控制體積法的基本概念如圖 3.1 (a) 圖 3.1 (b)及所示，其中圖 3.1 (a)為實際區域，圖 3.1 (b)則為計算區域，E、W、N、S 代表相鄰之格點，e、w、n、s 代表控制面。在水理控制方程式部分，除移流項採用一階精度之混合型上風法(hybrid scheme)(Spalding 1972)外，其餘的空間差分採用二階精度之中央差分法，時間差分則使用簡單的前項差分法。浸沒邊界求解步驟中之虛擬力為純量，則無須進行差分。中央差分法表示如下: 1 ₁ ₁ n _n _n e w p    _ _  _    _  _   (3.11) 1 ₁ ₁ n _n _n n s p    _ _  _    _  _   (3.12) 式中，代表了在e、w、n、s 四個控制面之u、v、h₁、h₂、d、z 和s zb 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                  混合型上風法由上風法(upwind scheme)與中央差分法組合而成，當移流效應重要時，採用上風法；而移流效應不重要時，則採用中央差分法。移流效應之重要程度以格網雷諾數(mesh Reynolds number) Re_x、Re_y做為判斷之因子，當 Rex

或Re_y 大於 2 時，代表移流效應重要，差分方法採用能反應方向性的上風法； Re_x 或Re_y 小於2 時，移流效應視為不重要，差分方法採用中央差分法。混合型上風法應用於本研究移流項的處理可表示為 , 1 1 1 1 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.13) , 1 1 1 1 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.14) 其中

(37)

23 0 2 0 2 1 2 1 2 1 2 1 2 y x x x y y x y R R R R R R          _  _  _ _{ } _ _{ }  _ ； (3.15) 在上列諸式中， , 1i j, n i j x u h R      ； , 2i j, n i j y v h R      ；  流體動黏滯係數(dynamic viscosity)；可代表 u 或 v 。 3.2 沉滓運移部分 沉滓運移控制方程式中，懸浮載之質量守恆方程式為一延散方程式，求解方式與水流運動方程式相同，分為以下兩個步驟: (1) 移流及反應(advection-reaction)步驟 1 2 a dC u dC v dC S dt h d h d d            (3.16) (2) 擴散(diffusion)步驟 1 1 1 2 1 2 2 1 2 2 1 1 a dC dC d h dC d h dC dt dt h h d d  h d h h d d  h d       _ _  _ _ _ _       (3.17) 分解之後的懸浮載質量守恆方程式，移流及反應步驟為一雙曲線偏微分方程式，擴散步驟則為橢圓偏微分方程式。先求解(2.13)、(2.14)及(3.16)三式之聯立方程，接著將求解得到的各變數之中間值與(3.17)進行反覆之疊代至收斂。沉滓運移控制方程式之處理，是採用控制體積法進行方程式的差分，接著以 ADI 法求解。此法與水理延散步驟相似，於此將不再贅述。總結以上，本模式之計算流程概述如下: (1) 輸入初始底床及水理資訊與格網坐標 (2) 以(3.1)式求解 n t 與(n 1) t時刻之間的未知變數 n 1 2 V  、 n 1 2 T  。 (3) 將上一步驟之結果代入(3.2)及(3.3)式求解未知變數 n V 、Tn。 (4) 接著以上式求得之結果，計算出虛擬力 n f 後代入(3.4)式求解(n 1) t時刻的

(38)

24 1 n V  、 n 1 T  ，反覆疊代求解直至收斂。 (5) 將收斂之流場結果及底床條件代入(3.16)及(3.17)式以求得下一時刻之底床資訊。 (6) 將步驟(4)所得到之結果代回步驟(1)反覆進行運算直到到達設定的運算時間為止。 (7) 輸出底床及水理模擬結果。

(39)

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

(40)

26 4.

第四章

橋墩周圍垂向流速推估

4.1 前言 為合理展現橋墩附近之局部沖刷現象，本研究在模式中加入了懸浮載源修正項，描述垂向水流對底床造成的沖刷起懸作用，此修正項之大小與垂向流場強度有著密切關係。沖刷坑形成後，坑內垂向流速之大小隨著沖刷坑幾何形狀變化。張氏(1998)曾根據 Raudkivi(1986)之試驗結果(圖 4.1)，得到垂向流速隨著沖刷深度線性變化之關係式 0(1 0.33 ) s m m d w w D   (4.1) 上式w_m  最大垂向流速； 0 m w  平床時刻之最大垂向流速；d_s  為沖刷深度； D為橋墩直徑。使用(4.1)式可以快速地估算沖刷過程中的最大垂向流速值。然而，要獲得沖刷過程中的最大垂向流速，除了底床變化之資訊外，如何決定平床時刻的垂向流場乃是相當關鍵之課題。由於本模式建構於二維水理模式架構之下，不具有水流垂向動量方程式。因此本研究選擇以三維 ANSYS-CFX 模式進行流場模擬後，對模擬結果進行迴歸分析，以統計分析之方式獲得平床時刻的最大垂向流速在空間上的分布經驗式，供後續動床模擬之使用。 4.2 ANSYS-CFX 三維模式簡介 4.2.1 模式原理簡介 本研究採用之 CFX 為一架構完整的泛用型流體力學軟體，具有發展成熟之多相流模型，可以模擬相與相之間相互作用傳遞之質量、動量以及能量問題。例如在進行兩相流的模擬時，液氣交互作用的自由液面是最難處理的部分，CFX 使用容積分率法(=Volume of Fluid method；VOF)，處理不交相混合的兩流體間的邊界問題。其基本原理是定義在格網單元內的體積分量(fractional volume)為函數F，如圖 4.2 所示，當單元中充滿了特定流體則體積分量函數F 1，相反則

為0，而 0F  即代表自由表面所在之格網單元。 1

另外，在使用 CFX 進行多相流之模擬時:(1)流體性質需為不可壓縮或微可壓縮；(2)不考慮不同相之間的輻射熱轉換；(3)即便不同相的流體，仍採用相同的

(41)

27 紊流模式。數值處理方面，CFX 採用以單位元(element)為基礎之有限體積法對控制方程式(質量守恆方程式及動量守恆方程式)進行離散後以共軛型式進行求解。以簡化之二維平面示意圖概略說明控制體積之概念，如圖 4.3 所示，在模擬區域中某一格網點與鄰近格網之間連結成為一個單位元，單位元中心點與鄰近單位元之中心點連結則成為一控制體積(虛線所圍之區域)，所傳遞之物理量皆儲存於格網節點上。數值方法採用SIMPLE (Semi-implicit Method for Pressure-Lined Equations)、 SIMPLEC (SIMPLE consistent) 或 PISO (Pressure-Implicit with Splitting of Operators)等方法來處理不可壓縮流之壓力項的求解。詳細之數值方法將於附錄一中介紹。 4.2.2 模式測試 為求謹慎並確定後續模擬之可信程度，本研究在使用 CFX 模擬自行設計之案例前，以 Ahmed(1995)之定床實驗案例做初步的模式測試。實驗渠道長 18.4m，寬 1.22m，深 0.67m，在距離渠道入口 13.1m 處置放直徑為 0.0889m 且表面光滑之圓形橋墩。上游水深為 0.182m，上游平均流速為 0.2927m/s，渠道底床為光滑之鋁製渠道。計算格網布置採用非均勻的結構格網如圖 4.4 所示(x r  35.95 ~ 35.95， 13.71 ~ 13.71 y r   ，z r0 ~ 15.625)，橋墩中心設於座標原點，其格網點總數目為 273564。在模擬設定方面，為了處理橋墩周圍的紊流流場，選用 k 紊流 模式進行模擬，計算時距則以最小格網間距與上游流速來決定，設定計算時距 0.025 t   秒，收斂條件設定為誤差≤10-4。在邊界設定方面，橋墩與底床固定邊界採非滑移(non-slip) 且光滑之條件，模式中定義其等效粗糙高度 kseq <0.03mm(Schlichting 2000)，換算為曼寧糙度係數約 0.012。上游入流及下游出流則分別以流速及靜水壓分布為邊界，上邊界則為開放式邊界。將最大垂向流速在不同距離點上(x r  2.25 ~ 3.5 )於水深方向分布之模擬結果與實驗量測數值做比較如圖 4.5。模擬結果與實驗值大致趨勢吻合，因此後續自行設計之案例將比照此測試案例的細部設定進行模擬。 4.3 因次分析

(42)

28 本研究探討直線渠道中圓形橋墩之沖刷，而橋墩周圍的局部垂向流場乃是構成局部沖刷的主因。參考明渠流實驗及模擬設置，選用可能的影響因子共四大類，分別為水力性質、流體性質、輸砂特性與橋墩性質。水力性質包括上游平均流速U 、上游水深₀ d 、渠寬₀ B、重力加速度g、渠道坡度S ；流體性質則包含₀ 流體密度、流體動力黏滯係數；輸砂性質則以Chezy糙度係數 c 代表底床粗 糙度，橋墩性質則為橋墩直徑D。將上述眾多影響因子與平床時刻之最大垂向流速w 之函數關係表示如下m₀ 0 ( , , , , , , , , )0 0 0 m w  f U d B g S   c D _(4.2) 式中；cChezy糙度係數R1/6 /n；R 水力半徑； n曼寧係數。由於坡度不同所造成的流場差異太大，局部臨界流況之產生嚴重影響了模擬結果進行迴歸分析時的準確程度。因此本研究為求簡化選擇忽略底床坡度的效應，並將上游平均速度U 、橋墩直徑₀ D、重力加速度g列為主要變數，應用柏金漢(Vachy-Buckingham)理論得到 5 個無因次參數如下 0 0 0 2 p p f U D U d D g Re Fr C D B c gD   ，  ，，，  _(4.3)

其中Re_p  橋墩雷諾數(pier Reynolds number)；Fr_p  橋墩福祿數(pier Froude

number)；d₀ / D 水深與橋墩直徑比； /D B 橋墩直徑與渠寬束縮比；C_f  摩擦因子。由上述無因次參數之組合，可將(4.2)式改寫為 0 0 (Re , , 0 , , ) m p p f w U  f Fr d D D B C _(4.4) 4.4 案例設計 為了探討垂向流速與(4.4)式中各項無因次參數之間的關係，本研究設計 31 組模擬案例，每一組的模擬設定僅針對一樣參數做調整而固定其他參數，觀察特定參數對垂向流場之影響，各組模擬案例之設定條件如表 4.1 所示。模擬案例的橋墩形狀皆為均勻圓柱，CFX 模擬之各項細部設定則參考 4.2 節。 4.5 模擬結果分析

(43)

29 4.5.1 垂向流速估算式 為了消弭三維模擬結果的眾多不確定性在統計過程可能造成的不良影響，以及避免各項無因次參數在數量級上的差異可能導致影響程度之誤判，將針對三維模擬得到的原始結果採取標準化(standardize)並刪減離群值的動作。離群值以標準化後的依變數為判定標準，刪除絕對值|_{ |大於 3 之資料。標準化的方式如下}* 所示: *       (4.5) 式中=代表各項無因次參數，上標(*)代表標準化後之參數，統計上稱之為標準分數(standard score)；  算數平均值；



樣本標準偏差。將參數標準化後進行各獨立變數的敏感度分析。採用線性迴歸模型進行迴歸分析，得到各獨立變數與依變數的迴歸關係如表 4.3。從迴歸關係表中可以看出各參數對垂向流速之影響，由迴歸係數大小可以看出橋墩福祿數Fr 相對其參數p 較為敏感，而橋墩雷諾數Re 則是敏感度最低的參數。且橋墩雷諾數及橋墩直徑_p 與渠寬之束縮比D B 的 p-value 皆大於 0.05，其代表意義為參數與應變數最大垂/ 向速度之關係不顯著，若以相同條件進行100 次的預測則其中會有 5 次的結果不一致。接下來，將較為重要之橋墩福祿數Fr 、水深與橋墩直徑比_p d₀ /D 與糙度 係數C_f ，以及必要之距離因子與角度因子取自然對數後進一步做迴歸得到垂向流速迴歸關係式(迴歸係數列於表 4.3)如下 0 0.479 1.07 0.436 1.49 0.133 0 0 0.308 m p f w _d _x _r Fr C U D r             _ _ _ _{ } _       (4.6) 式中 x計算格點與橋墩中心之距離；r 橋墩半徑；  格網點與水流方向之夾角，以上空間變數之意義請參考圖 4.6。 4.5.2 適用性分析 為了解估算式之預測能力與適用範圍，本研究使用(4.6)式進行估算並將結果與 CFX 模擬結果進行比對如圖 4.7(a)與(b)，可以發現在橋墩兩側至墩後的範圍內，公式預測能力明顯下降。推論此一現象之原因，主要是因為水流分離多發生在與墩前水流方向夾角為 90°~120°之間，水流在此區域形成紊流的尾流區，複

橋墩水深平均二維沖刷模式之發展

國

立

交

通

大

學

土木工程學系

碩士論文