三維大渦模式在明渠水流之應用

全文

(1)國立交通大學 . 土木工程學系碩士論文 . 三維大渦模式在明渠水流之應用 Use of a Three-Dimensional LES Model for Open-Channel Flow . 研究生：陳蓉瑩指導教授：楊錦釧博士指導教授：謝德勇博士. 中華民國一 ○ 二年七月.

(2) 三維大渦模式在明渠水流之應用 Use of a Three-Dimensional LES Model for Open-Channel Flow. 研究生：陳蓉瑩指導教授：楊錦釧指導教授：謝德勇. Student：Jung-Ying 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 2013 Hsinchu, Taiwan, Republic of China. . 中華民國一 ○ 二年七月 . .

(3) 誌謝 . 承蒙恩師楊教授錦釧及謝博士德勇之悉心指導，使我得以順利完成學業，於. 此致上萬分感謝。感謝口試委員傅教授武雄、黃教授良雄及賴教授進松給予寶貴意見與指教，使本研究更臻完善。感謝王教授克漢多次給予許多指正與建議，幫助完善本論文成果。 . 於兩年研究生涯中，特別感謝和政學長、胤隆學長、世偉學長、浩榮學長、. 文祿學長、建華學長、弘恩學長、仲達學長、仁凱學長、聖翔學長、家偉學長、昀直學長、建翔學長、舒勤學姐與芳綺學姐給予的教導與照顧。感謝同窗醇國、亞雯、瑋廷、健賓、于軒、威辰與莉玲在研究所生活互相扶持與關懷，及感謝學弟妹怡君、瑋宸、佩衡、承儒與伯華的協助，在此表示誠摯謝意。感謝摯友麗嘉、唄楒、京穎、婕恩以及一路走來支持與陪伴著我的朋友們，讓我能正向積極從容的面對生活中的挑戰，最後將此論文獻給親愛的家人朋友們，一同分享這份喜悅。 I. .

(4) 三維大渦模式在明渠水流之應用學生：陳蓉瑩. 指導教授：楊錦釧. 學生：洪聖翔. 指導教授：謝德勇國立交通大學土木工程研究所. 摘要本研究基於垂直水平分離演算 (VHS)之概念，利用大渦模擬 (large-eddy simulation, LES)發展一三維大渦模式，其中紊流黏滯係數採用 Smagorinsky 次網格模式計算，並以雷諾平均假設下之兩種零方程紊流模式相互比較，探討明渠流之紊流流場現象。在座標系統上，水平方向與垂直方向分別採用正交曲線座標系統與  座標系統，使模式能便利處理不規則的渠道與底床邊界。水理模式之架構係先由深度平均之控制方程式求解水平二維每個計算格點之水面高程以及水深平均流速，並將之代入流速差異量方程式，求得三維每個計算格點相對於水深平均流速之流速差異量，即可求解出三維空間上每個計算格點之流速分佈。最後，本研究分別採用突擴渠道流速場案例以及直渠道渦度場案例來模擬紊流流場，探討不同計算網格大小對不同型態紊流模式模擬結果之影響。. 關鍵字：三維模式、正交曲線座標系統、  座標系統、大渦模擬、明渠紊流 . . II.

(5) Use of a Three-Dimensional LES Model for Open-Channel Flow Student：Jung-Ying Chen. Advisor：Jinn-Chuang Yang. Student：Sheng-Hsiang Hung. Advisor：Te-Yung Hsiehaaa. Department of Civil Engineering National Chiao-Tung University. ABSTRACT To. investigate. the. turbulence. of. open-channel. flow,. a. hydrostatic. three-dimensional model based on a vertical-horizontal splitting (VHS) method with large-eddy simulation was used in this study. The eddy viscosity was treated by Smagorinsky sub-grid model. The results were compared with two commonly used Reynolds-averaged zero-equation models. By following the vertical and horizontal splitting concept, the shallow water flow governing equations were split into two parts including the depth-averaged 2D equations and velocity defect equation in vertical direction. The depth-averaged 2D equations were transformed into orthogonal curvilinear coordinate system, to solve the elevation of every water column and depth-averaged velocities. The velocity defect equations can be derived by subtracting the 2D depth-averaged equations from the 3D Navier-Stokes equations. The sigma coordinate was used to fit the complex geometry in channel bed. Incorporating with the continuity equation, the three-dimensional velocity field can therefore be solved. Two experimental cases including sudden expansion and straight channel were simulated to investigate the effect of computational grid size on the model’s accuracy. Keywords：3D model, orthogonal curvilinear coordinate system, sigma coordinate system, large-eddy simulation, open-channel turbulent flows . III.

(6) 目錄誌謝 .............................................................................................................................................I 摘要 ........................................................................................................................................... II ABSTRACT ............................................................................................................................. III 目錄 .......................................................................................................................................... IV 表目錄 ...................................................................................................................................... VI 圖目錄 .................................................................................................................................... VII 符號表 ...................................................................................................................................... IX 第一章緒論 .............................................................................................................................. 1 1.1 研究動機與方向 ......................................................................................................... 1 1.2 文獻回顧 ..................................................................................................................... 1 1.3 研究目的與方法 ......................................................................................................... 3 1.4 章節介紹 ..................................................................................................................... 4 第二章理論基礎 ...................................................................................................................... 5 2.1 控制方程式 ................................................................................................................. 5 2.1.1 三維方程式 ...................................................................................................... 5 2.1.2 二維水深平均方程式 ...................................................................................... 7 2.1.3 流速差異量方程式 .......................................................................................... 8 2.2 輔助關係式 ............................................................................................................... 10 2.2.1 底床剪應力 .................................................................................................... 10 2.2.2 層流與紊流剪應力 ........................................................................................ 10 2.3 紊流黏滯係數 ........................................................................................................... 11 2.3.1 次網格模式 .................................................................................................... 11 2.3.2 零方程式 I ...................................................................................................... 12 2.3.3 零方程式 II .................................................................................................... 12 2.4 邊界條件 ................................................................................................................... 12 2.4.1 二維水深平均部分 ........................................................................................ 12 2.4.2 流速差異量部分 ............................................................................................ 12 第三章數值架構 .................................................................................................................... 15 3.1 求解架構 ................................................................................................................... 15 3.1.1 二維水深平均部分 ........................................................................................ 15 3.1.2 流速差異量部分 ............................................................................................ 17 3.2 數值差分 ................................................................................................................... 18 IV.

(7) 3.2.1 二維水深平均部分 ........................................................................................ 18 3.2.2 流速差異量部分 ............................................................................................ 19 第四章結果與討論 ................................................................................................................ 23 4.1 突擴渠道流速場分析案例 ....................................................................................... 23 4.2 直渠道渦度場分析案例 ........................................................................................... 25 第五章結論與建議 ................................................................................................................ 66 5.1 結論 ........................................................................................................................... 66 5.2 建議 ........................................................................................................................... 67 參考文獻 .................................................................................................................................. 68 附錄 A ...................................................................................................................................... 74 附錄 B ...................................................................................................................................... 76. V.

(8) 表目錄表 4.1 突擴渠道案例之模擬物理條件 ................................................................................... 29 表 4.2 直渠道案例之模擬物理條件 ....................................................................................... 29 表 4.3 直渠道案例不同寬深比(B/H)條件之模擬耗費時間 .................................................. 30. . . . . VI.

(9) 圖目錄圖 2.1 正交曲線座標轉換示意圖 ........................................................................................... 13 圖 2.2 σ 座標轉換示意圖 ........................................................................................................ 13 圖 2.3 水深方向流速剖面示意圖 ........................................................................................... 14 圖 3.1 模式計算流程圖 ........................................................................................................... 20 圖 3.2 水平二維模式控制體積法示意圖 ............................................................................... 21 圖 3.3 交錯格網(staggered grid)示意圖 .................................................................................. 22 圖 3.4 垂直模式控制體積法示意圖(計算區域) .................................................................... 22 圖 4.1 Babarutsi(1989)實驗水槽設置圖 ................................................................................. 31 圖 4.2 Babarutsi(1989)環流案例實驗量測流速分佈圖 ......................................................... 31 圖 4.3 Babarutsi(1989)實驗量測數據校正示意圖 ................................................................. 31 圖 4.4 Babarutsi(1989)計算網格 361  62 佈置圖 .................................................................. 32 圖 4.5 Babarutsi(1989)計算網格 721  123 佈置圖 ................................................................ 32 圖 4.6 Babarutsi(1989)計算網格 721  184 佈置圖 ................................................................ 33 圖 4.7 Babarutsi(1989)計算網格 721x  367 佈置圖 .............................................................. 33 圖 4.8 突擴案例 U / U 0 於各量測斷面之比較圖(網格數為 361 62 ) ................................... 34 圖 4.9 突擴案例 U / U 0 於各量測斷面之比較圖(網格數為 721 123 ) .................................. 35 圖 4.10 突擴案例 U / U 0 於各量測斷面之比較圖(網格數為 721 184 ) ............................... 36 圖 4.11 突擴案例 U / U 0 於各量測斷面之比較圖(網格數為 721 367 )............................... 37 圖 4.12 突擴案例 U / U 0 於各量測斷面不同計算網格之次網格模式結果比較圖 ............. 38 圖 4.13 突擴案例 U / U 0 於各量測斷面不同計算網格之零方程式 I 結果比較圖 .............. 38 圖 4.14 突擴案例 U / U 0 於各量測斷面不同計算網格之零方程式 II 結果比較圖 ............ 39 圖 4.15 Nezu and Rodi(1986)實驗案例設置圖 ...................................................................... 40 圖 4.16 明渠紊流之垂向流速分佈示意圖 ............................................................................ 41 圖 4.17 邊界層之各層定義範圍示意圖 ................................................................................ 41 圖 4.18 邊界層內格網點分佈對模擬結果之影響示意圖 .................................................... 42 圖 4.19 Nezu and Rodi(1986)計算網格 74  5 12 佈置圖 ...................................................... 43 圖 4.20 Nezu and Rodi(1986)計算網格 363  13  23 佈置圖 .................................................. 43. VII.

(10) 圖 4.21 Nezu and Rodi(1986)計算網格 725  25  25 佈置圖 ................................................. 44 圖 4.22 直渠案例於量測斷面中間水柱縱向流速分佈之比較圖(網格數為 74  5 12 ) ..... 45 圖 4.23 直渠案例於量測斷面中間水柱縱向流速分佈之比較圖(網格數為 363  13  23 ) . 45 圖 4.24 直渠案例於量測斷面中間水柱縱向流速分佈之比較圖(網格數為 725  25  25 ) 46 圖 4.25 直渠案例各斷面中間水柱縱向流速之次網格模式模擬結果 ................................ 47 圖 4.26 直渠案例各斷面中間水柱縱向流速之零方程式 I 模擬結果 ................................. 48 圖 4.27 直渠案例各斷面中間水柱縱向流速之零方程式 II 模擬結果 ............................... 49 圖 4.28 直渠案例模擬結果之渦度比較圖(網格數為 74  5 12 ) ......................................... 50 圖 4.29 直渠案例模擬結果之渦度比較圖(網格數為 363  13  23 ) ..................................... 51 圖 4.30 直渠案例模擬結果之渦度比較圖(網格數為 725  25  25 ) .................................... 52 圖 4.31 直渠案例次網格模式模擬結果之向量及渦度圖(網格數為 74  5 12 ) ................. 53 圖 4.32 直渠案例次網格模式模擬結果之向量及渦度圖(網格數為 363  13  23 ) ............. 54 圖 4.33 直渠案例次網格模式模擬結果之向量及渦度圖(網格數為 725  25  25 ) ............ 55 圖 4.34 直渠案例零方程式 I 模擬結果之向量及渦度圖(網格數為 74  5 12 ) ................. 56 圖 4.35 直渠案例零方程式 I 模擬結果之向量及渦度圖(網格數為 363  13  23 ) ............. 57 圖 4.36 直渠案例零方程式 I 模擬結果之向量及渦度圖(網格數為 725  25  25 ) ............. 58 圖 4.37 直渠案例零方程式 II 模擬結果之向量及渦度圖(網格數為 74  5 12 ) ................ 59 圖 4.38 直渠案例零方程式 II 模擬結果之向量及渦度圖(網格數為 363  13  23 ) ............ 60 圖 4.39 直渠案例零方程式 II 模擬結果之向量及渦度圖(網格數為 725  25  25 ) ........... 61 圖 4.40 直渠案例次網格模式在不同計算網格下於量測斷面中間水柱之 t 分佈圖 ........ 62 圖 4.41 直渠案例零方程式 I 在不同計算網格下於量測斷面中間水柱之 t 分佈圖 .......... 62 圖 4.42 直渠案例零方程式 II 在不同計算網格下於量測斷面中間水柱之 t 分佈圖 ........ 63 圖 4.43 直渠案例次網格模式紊流強度. uˆ 模擬結果之垂向分佈圖 ................................... 64 U*. 圖 4.44 直渠案例次網格模式紊流強度. vˆ 模擬結果之垂向分佈圖 ................................... 64 U*. 圖 4.45 直渠案例次網格模式於近壁區之紊流強度. VIII. uˆ 模擬結果 ....................................... 65 U*.

(11) 符號表 B  渠道寬度. C f  摩擦係數 ( g / c 2 ) Cs  大渦模擬係數. c  Chezy 係數 d  水深 d e  突擴長度 f c  科氏力係數 Fr  福祿數 g  重力加速度 H  下游水深. h  正交曲線座標轉換係數 k s  粗糙高度.   Von Karman’s 係數 Lr  環流長度 p  壓力. Q  流量 Re  雷諾數. S  渠道縱向坡度. T  有效剪應力項 t  時間. U 0  入流均勻流速 U *  剪力速度. u   方向流速 v   方向流速 w   方向流速 IX.

(12) z  卡式座標系統之垂直方向座標. zb  底床高程 zs  水面高程 zl  近底床格網之厚度.   動力黏滯係數.  l  層流黏滯係數  t  紊流黏滯係數  H  水平方向黏滯係數  V  垂直方向黏滯係數.   流體密度.   層流剪應力   緯度   地球自轉角速度   網格幾何平均值. t  時間間距 t LES  大渦模擬時間.  、  、 z   、 、  方向之格網間距。. 上標 n  nt 時刻之已知變數. n  1  (n  1)t 時刻之已知變數 n  1/ 2  (n  1)t 與 nt 間之未知變數 b  底床代號. s  水面代號.   空間平均.   水深平均 X.

(13)    差異量(流速為與深度平均流速之差).    次網格擾動量  ˆ  紊流擾動量下標 i 、 j   、 方向格點編號. 1  物理量在  方向的值 2  物理量在 方向的值. 3  物理量在  方向的值 p 、 e 、 w 、 n 、 s  通量所在格點位置. XI.

(14) 第一章緒論 1.1 研究動機與方向紊流(turbulent flow)是自然界與工程中十分普遍之流場現象，譬如河川中之水流及大氣中之風場。在大多數之情況下，起因為流場之雷諾數過大而造成不穩定之紊流流場，在紊流發展之過程中，黏滯力會破壞來自於因動能消耗而增加之內能，故紊流需一直提供能量來彌補因黏滯力所造成之損失；而在高雷諾數流場之狀況下，各種尺度之渦流(eddy)結構均被激發出來，於時間與空間上呈不規則性，屬三維維度之流場現象。隨著電腦之計算速度及容量已大幅度提高，至今已有許多研究著重於利用紊流數值模擬作為工程實務上之應用。計算紊流之方法可分為三大類：直接數值模擬 (DNS, direct numerical simulation) 、雷諾平均數值模擬 (RANS, Reynolds-averaged Navier-Stokes simulation)、大渦數值模擬(LES, large eddy simulation)，其中 DNS 所需之格點與計算時間龐大，僅限於模擬低雷諾數及簡單幾何形狀之流場；而 RANS 係將流場性質分為平均值及擾動值，代入原本之控制方程式中，衍生出不同之假設式，但由於其為時均化之結果，無法有效表現出紊流流場中流體性質隨時間變動之特性；有鑒於以上兩種方法之限制，LES 將渦流尺度分離，分別使用不同方法求解。由於紊流動能並非均勻地分佈於不同尺度之渦流，而是大多集中於低頻率、大尺度之渦流上，故只需將主導流場之大尺度渦流以 DNS 之方式求解，小尺度渦流則以適用性較廣之次網格模式計算其能量耗散，以減少計算格網點之使用。. 1.2 文獻回顧至今已有許多研究利用物理模型以及數值模式針對渠道紊流之流場現象進行分析探討，在物理模型方面有 Nezu and Nakagawa(1981)、Komori et al.(1982)、 Nezu and Rodi(1986)、Bashidi and Banerjee(1988)、Tominaga et al.(1989)、Komori et al.(1990)等人；而發展紊流之數值模式更是不勝枚舉，以下將文獻中較常見之幾種紊流模式作一概略介紹： Reynolds(1895)提出雷諾平均下之 Navier-Stokes 方程式，後人根據不同之 1.

(15) 假設，逐漸發展出各種類型之紊流模式，若以使方程式封閉之偏微分方程式數目來區分，則有零方程模式、單方程模式、雙方程模式及多方程模式。零方程模式是最早也最簡單之紊流模型，其代表有 Boussinesq(1877)之渦流黏性模型、 Prandtl(1925) 之混合長度模型、 Von Karman(1931) 之局部相似模型以及 Prandtl(1942)之自由剪力層模型。零方程模式具有形式簡單、實用方便之優點，在傳統水力學有廣泛之應用，但未考慮動量交換及擴散傳輸，故缺乏通用性。單方程模式是由 Kolmogorov(1942)及 Prandtl(1945)所提出，增加了一個以紊流動能平均值 k 為變數之紊流傳輸方程式，但其紊流特徵長度不易估算，故其應用僅限於剪力層型流動。因應紊流模型發展之需要，衍生出雙方程模式來符合工程實際之應用；但在雙方程模式當中，標準 k   模型(Harlow and Nakayama 1967； Launder and Spalding 1974；Rodi 1981；Chen 1983；Lumley 1983； Markatos 1986； Takemitsu 1986)其理論假設較為簡單，得到廣泛之應用與驗證，但不適用在非等向性之流場，Rodi(1980)指出等向性之假設無法準確地模擬紊流，且模式假設動能消散率  只與大尺度渦流有關，反而與小尺度渦流無關，此假設與紊流能量消散理論不符，另外對於壁面附近之流場也須強加處理；另一種較為廣泛應用之雙方程模式為 k   模型(Rodi and Spalding 1970；Saffman 1970； Fisher et al. 2000； Rameshwaran and Naden 2003)，適用於黏性底層積分和預測逆壓梯度方面之流場；為了彌補標準 k   模型之缺陷，(Pope 1975；Nisizima and Yoshizawa 1987； Speziale 1987；Rubinstein and Barton 1990)提出非線性 k   模式，適用於具有分離流動、高低雷諾數之流場。隨著計算能力之躍進，進而發展出雷諾應力多方程模型(Hanjalic and Launder 1972；Launder et al. 1975；Lumley 1975；Reece 1977； Basara and Cokljat 1995；Cokljat and Younis 1995；Hyeongsik and Sung 2006)，雷諾應力模型係繼雙方程模型更為理想之模型，其考慮雷諾應力之傳輸，可用來解決更複雜之流場問題，但其需求解之微分方程過於複雜，以致於應用範圍狹隘。隨著電子計算機之功能快速成長，可考慮由直接數值模擬以及大渦數值模擬進一步發展紊流模式。由於直接數值模擬是基於直接求解控制方程式之變數，所以其計算網格以及時間步階必須非常小，才能計算所有尺度之渦流對流速項以及壓力項造成之變化，故受限於低雷諾數及簡單邊界條件下之流場狀況 (Moser and Moin 1984；Spalart 1986；Kim et al. 1987；Mansour et al. 1988)。 2.

(16) 大渦數值模擬係將流場中之大尺度渦流以直接數值模擬之方法求解，小尺度渦流則以選用之次網格模式計算其效應。在次網格模式方面，最早之次網格模式為 Smagorinsky(1963)提出，基於混合長度假設下應用於大氣現象之二維模式，而後開始廣泛應用於有關大氣邊界層流以及水力機械內部流場； Deardroff(1970)採用 Smagorinsky 之次網格模式應用於不考慮壁面效應之渠道流，並以 6720 之網格數與理論值相互驗證，之後的學者亦開始採用大渦模式於低雷諾數之流場(Schumann 1975；Moin and Kim 1981)，從此奠定次網格模式運用在渠道紊流模擬之地位。此後，更將次網格模式應用於二維及三維自由液面之渠道模擬(Rogallo and Moin 1984；Dommermuth and Novikov 1993；Lesieur and Metais 1996；Shi et al. 1999；Meneveau and Katz 2000；Broglia et al. 2003； Armenio 2005；Yue et al. 2005；Rodi et al. 2008；Chung and Pullin 2009)；以及彎道(Stoesser et al. 2008；Balen et al. 2009)、底床突起(Calhoun and Street 2001； Falcomer and Armenio 2002)、流經結構物(Koken and Constantinescu 2008； Mccoy et al. 2008)之應用。由前述紊流模式之發展回顧，可知相關研究已相當成熟，但以目前文獻中之三維大渦模式來說，計算量依然龐大也相當耗時，若能採用擬似三維之模式，將水平及垂直流場分開求解，將可降低計算成本，亦能提供合理的三維流速分佈資訊；另外，目前大渦模式尚未看到以正交曲線座標之型式來處理複雜之幾何邊界條件。因此本研究將運用水平垂直分離演算法之概念，並採用正交曲線座標以便利處理不規則幾何邊界的問題。. 1.3 研究目的與方法本研究之目的為應用大渦模式中之次網格紊流模式，探討明渠水流之紊流流場。基於洪(2011)之靜水壓三維水理模式，應用垂直水平分離演算之概念，將次網格模式轉換成正交曲線座標形式，另外以鍾(2012)雷諾平均架構下之兩種零方程紊流模式並列討論。水平方向座標系統採用正交曲線座標，其能適當表達不規則渠道形狀，且座標軸主軸方向(本研究使用之  軸方向)即為主流方向；而在水深方向則採用  座標系統(Blumberg and Mellor,1983)，其能解決自由水面在固定格網點上變動而影響模式無法準確計算水面之壓力邊界條件的問題，也能將因底床起伏而產生之不規則格網，轉換為便於計算之矩形格網，如此可得到 3.

(17) 精確度較高之模擬結果。為驗證模式之正確性，分別將次網格模式以及兩種零方程模式之模擬結果與實驗量測值互做比較。. 1.4 章節介紹前面已闡述本研究之動機與方向、文獻回顧、研究目的與方法，本節將扼要說明本文章之內容。第一章為緒論，說明本研究之背景與目的，並回顧相關模式發展之文獻，提出本研究之方法與研究之重點。第二章為理論基礎，在正交曲線座標系統下，由三維那威爾-史托克司 (Navier-Stokes)方程式導出模式控制方程式、輔助方程式及紊流黏滯係數之使用及邊界條件之設定，均於本章介紹。第三章為數值架構，水平二維水理模式及垂直水理模式之數值方法於本章說明，並簡述模式之計算流程。第四章為模式驗證，針對模式建構部分，採用具實驗量測數據之案例進行模式的檢定工作，並簡述應用之案例。第五章為結論與建議，對研究成果作綜合性之歸納說明，並針對研究尚未考量、不盡完備或日後可繼續研究之處提出建議。. 4.

(18) 第二章理論基礎本研究以靜水壓為主要假設之三維動量方程式為基礎，以及垂直水平分離演算概念，分為水平二維方程式與流速差異量方程式(洪，2011)，採用大渦紊流模擬(large-eddy simulation , LES)，並結合 Smagorinsky(1963)所提出之次網格模式(sub-grid scale model , SGS)，建構一非零方程之紊流模式；另外並以鍾(2012) 發展之雷諾平均架構下之兩種零方程(zero-equation)紊流模式相互比較。為了能適當表達天然河道不規則的幾何形狀，模式在水平方向採用正交曲線座標系統，如圖 2.1 所示；在垂直方向則採用  座標系統(Blumberg and Mellor,1983)，如圖 2.2 所示，如此能將不規則的計算區域轉換至矩形計算區域求解，且能便利地處理渠道的側壁、自由液面及底床邊界。. 2.1 控制方程式在大渦模式之理論基礎當中，紊流場中之大尺度渦流主導著全流場的物理機制，而大尺度渦流與小尺度渦流之差異在於以下幾點： (1). 大尺度渦流之行為為主要影響平均流場之質量、動量及能量變化，而小尺度渦流與大尺度渦流之間為非線性之交互作用，對平均流場影響甚小。. (2). 流場本身之幾何形狀及性質大大影響大尺度渦流之行為，因此具有非等向性(anisotropic)之特性，而小尺度渦流則趨於等向性(isotropic)。. (3). 大尺度渦流之時間尺度與平均流場近似，相較之下小尺度渦流則快速地生成及消減，以能量耗散之形式存在於流場中。. 由以上理論可知，在控制方程式中僅須求得大尺度之流場行為，故以直接求解之方式計算，而小尺度之行為則以適用性廣泛之次網格模式模擬出能量耗散之效應即可。 2.1.1 三維方程式基於不可壓縮流之假設下，對那威爾-史托克司(Navier-Stokes)方程式取空間平均後，得水理控制方程式如下：連續方程式 5.

(19)    (h2u )  (h1v )  (h1h2 w)  0   z. (2.1). 動量方程式 ξ 方向： 1  1  2(  uv ) h1 (  u 2 ) h2 (  v 2 ) h2     (u )  ( u 2 )  (  uv )  (  uw)  t h1  h2  z h1h2  h1h2  h1h2  . 1  1  2(  u v) h1 (  u2 ) h2  v2 h2    (  u 2 )  (  u v)  (  uw)  h1  h2  z h1h2  h1h2  h1h2 . . 2  1  1  (  vu)  (  wu)   (  uu )  (  uv)  (  uw)  h1  h2  h2  z z. . 2  (uv  uv ) h1  (2uu) h2  (2vv) h2   h1h2 h1h2  h1h2  .   fcv . 1 p 1  h1  h1h2.    12 h1  22 h2     (h2 11 )   (h1 12 )  z (h1h2 13 )   h h   h h    1 2 1 2. (2.2). η 方向： 1  1  2(  uv ) h2 (  v 2 ) h1 (  u 2 ) h1     (v )  (v 2 )  (  uv )  (  vw)  t h2  h1  z h1h2  h1h2  h1h2  . 1  1  2(  u v) h2 (  v2 ) h1 (  u 2 ) h1    (  v2 )  (  uv)  (  vw)  h2  h1  z h1h2  h1h2  h1h2 . . 1  2   1   (  uv)  (  vv)  (  wv)  (  vu )  (  vw) h1  h2  h1  z z. . 2  (uv  u v ) h2  (2vv) h1  (2uu) h1    h1h2 h1h2  h1h2 .    fcu .   h2  11 h1 1 p 1     (h1 22 )  (h2 12 )  (h1h2 23 )   12    h2  h1h2    z  h1h2  h1h2 . (2.3). z 方向：在一般淺水的天然河道，垂向流速遠小於水平方向流速時，垂直方向之動量方程式可用靜水壓分佈來簡化， p  g  0 z 6. (2.4).

(20) 而垂直方向之流速可透過連續方程式求得，以減少模式之計算量。以上諸式中，  、 、 z 為三維正交曲線座標方向，其中  、 為水平方向， z 為水深方向；下標 1、2、3 分別代表物理量在  、 、 z 方向代號； h1 、 h2 分別為  、  方向之轉換係數； u 、 v 、 w 分別為  、 、 z 方向流速； g 為重力加速度； t 為. 時間；  為層流剪應力； ( ) 表空間平均； ( ) 表次網格擾動量； fc (= 2 sin  )為科氏力係數；  為地球自轉角速度；  為緯度；  為密度； p 為壓力。 2.1.2 二維水深平均方程式將 2.1.1 節之三維水理控制方程式(2.1)、(2.2)、(2.3)利用萊布尼茲法則 (Leibuniz rule)對深度方向積分，加上運動邊界條件及動力邊界條件，並取深度平均，可得水平二維水理控制方程式。連續方程式 h1h2. d    (h2ud )  (h1vd )  0 t  . (2.5). 動量方程式 ξ 方向：  u u u v u uv h1 v 2 h2  1 ps  g zs         fc v   t h  h  h h  h h  h  h1       1 2 1 2 1 2 1  .  . h h  1 T12 1  h2 1 T11 ( s1   b1 )  2T12 1  T22 2      R  T11 dh1h2   d    dh2  dh1 . (2.6). η 方向：  v v v u v uv h2 u 2 h1  1 ps  g zs          fcu  h2  h2   t h2  h1  h1h2  h1h2  .  . h h  1 T12 1  h1 1 T22 ( s 2   b 2 )  2T12 2  T11 1      R (2.7)  T22 dh1h2   d    dh1  dh2 . 式中， zs. T11   [ 11   u 2   (2uu   u 2 )] dz zb. 7. (2.8).

(21) zs. T22   [ 22   v 2   (2vv  v2 )] dz zb. zs.     (uv  vu   u v)]dz T12   [ 12   uv zb. (2.9) (2.10). 以上諸式中，d 為水深； zs 為水面高程； zb 為底床高程；  s1 、 s 2 為主流與側方向之水面剪應力； b1 、 b 2 為主流與側方向之底床剪應力； ( ) 表水深平均； (  ) 表物理量之空間微變量(例：u  u  u )；下標 s 、b 分別為水面及底床代號；T 為有效剪應力項(effective stress term)，其包含層流剪應力、延散剪應力與紊流剪應力；另外，本研究以 R 、 R 表示為水深平均後的殘餘項，如附錄 B 所示，假設其數值計算結果可忽略，在模式計算時並未考慮。 2.1.3 流速差異量方程式將 2.1.1 節之三維水理控制方程式(2.2)、(2.3)以靜水壓假設代入，並令 u  u~  u 、 v  v~  v (如圖 2.3)，得到之方程式減去 2.1.2 節之水平二維水理方程. 式(2.6)、(2.7)，即可得到流速差異量方程式。 ξ 方向： u u u u u u u v u v u v u u        t h1  h1  h1  h2  h2  h2   .  h1 uv   h1 2vv h2 v 2 h2 uv h1 uv     h1h2  h1h2  h1h2  h1h2  h1h2 .  f c v . 1  13 ( s1   b1 )   ( Horizontal Diffusion in  )  d  d. (2.11). η 方向： v v v v v v v u v u v u v v        t h2  h2  h2  h1  h1  h1    h2 uv   h2 2uu h1 u 2 h1 uv h2 uv      h1h2  h1h2  h1h2  h1h2  h1h2 .   f c u . 1  23 ( s 2   b 2 )   ( Horizontal Diffusion in  )  d  d. 式中， 8. (2.12).

(22) .  u  u  v  v       w t h1  h1  h2  h2  z. (2.13). . z  zb d. (2.14).   d  d t t. (2.15). 1 zb  d     d  d . (2.16). 1 zb  d     d  d . (2.17).  1  z d. (2.18). 其中，  為水深分層之垂直方向座標，無因次水深。 ( Horizontal Diffusion in  )  . h h  1  h2  2T12 1  T22 2   T11  dh1h2     . 1 T12 1 T11 u  2 H   1  2 H h2  2 H  2u         dh2   dh1    h1   h1  h12 h2   h12  2. . 2 H   v  h1  H   h2        h1   h1h2   h2   h1  .  H h1  2  u   2    2 H h2  2  h1 .  v   H   h1    u         h2  h2   h2    h1 .  1   v  1   u   h1 2 H v h2  2  2   2    h1   h2  h2   h1    h1h2  . ( Horizontal Diffusion in  )  . (2.19). h h  1  h1  2T12 2  T11 1   T22  dh1h2     . . 1 T22 1 T12 v  2 H   1  2 H h1  2 H  2 v        dh2   dh1    h2   h2  h1h22   h22  2. . 2 H   u  h2  H   h2        h2   h1h2   h1   h1  . .  H h2  2  v     2 H h12  2  h2 .  v   H   h1    u         h2  h1   h2    h1 .  1   v  1   u   h2 2 H u h1  2  2   2    h1   h2  h2   h1    h1 h2  . (2.20). 其中， H 為水平方向黏滯係數   l  t ； l 為層流黏滯係數，約為 1.12 106 m 2 / s ；. 9.

(23)  t 為紊流黏滯係數。. 2.2 輔助關係式. 2.2.1 底床剪應力底床剪應力採用 French (1986)之經驗式：   30 z1    b1  ub ub  2.5ln      2.72ks  . 2.   30 z1    b2  vb vb  2.5ln     2.72k s   . 2. (2.21). (2.22). 式中， ub 、 vb 分別為  、 方向之近底床流速； z1 為近底床流速之格網與底床間之垂直距離； ks 為粗糙高度。 2.2.2 層流與紊流剪應力採用 Boussinesq(1877)之渦流黏性理論，層流與紊流剪應力可合併表示為  1 u  11 v h1   (2uu   u2 )  2 H      h1  h1h2  . (2.23).  1 v  22 u h2   (2vv  v2 )  2 H      h2  h1h2  . (2.24).  12  (uv  vu   u v)  2 H .  h2   v  h1   u         h1   h2  h2   h1  . (2.25). 垂向剪應力以流速梯度表示時，若考慮正交曲線座標以及  座標時，可表示為 13 1 u 1  w  V (  ) d  h1  . (2.26).  23 1 v 1  w  V (  ) d  h2  . (2.27). 10.

(24) 式中 V 為垂直方向黏滯係數   l  t 。. 2.3 紊流黏滯係數紊流黏滯係數理論上也應歸類於輔助方程式中，但因此部分為本研究探討之重點，因此特別將此部分的處理方式以專節之方式彙整說明。如前所述，本研究大渦模擬在紊流係數之處理上係使用次網格模式，另外雷諾平均模擬採用兩種零方程模式求解。茲將上述三種類型紊流黏滯係數之處理方式說明如下： 2.3.1 次網格模式本研究採用 Smagorinsky(1963)之次網格紊流模式，此模式理論之假設有二，其一為採用 Boussinesq 提出之渦流黏性理論，將紊流剪應力表示為紊流黏滯係數與剪應變率之關係式，  ij  2 t Sij. (2.28) 1 ui u j  )。 2 x j xi. 上式中 t 為 H   l  t 與 V   l  t 中之 t ； Sij 為剪應變率  (. 其二之假設為基於混和長度假說(mixing length hypothesis)，可將 t 表示為  t  (l ) 2 S. (2.29). 其中， l  Cs  S  (2Sij Sij )1/2. 上式中 Cs 為大渦模擬係數，本研究取 0.1 (Deardoff 1970)；   (  x  y  z )1/3 為能直接計算之最小渦流尺寸。綜合以上假設之次網格紊流黏滯係數可合併表示為  t  (Cs  ) 2 (2 Sij Sij )1/ 2. (2.30). 由於 Smagorinsky(1963) 次網格模式之假設， t 為純量無方向性，因此本研究採  H 與 V 為相同值。. 11.

(25) 2.3.2 零方程式 I 採用 Elder(1959)提出之紊流黏滯係數經驗式  t  U * d / 6. (2.31). 其中，  為 Von Karman’s 係數(約等於 0.41)； U * 為剪力速度； d 為水深。 2.3.3 零方程式 II 另一種採用 Jobson and Sayre(1970)提出之紊流黏滯係數經驗式  t  U *d (1   ). (2.32). 其中，(2.29)式為(2.28)式之積分結果。. 2.4邊界條件. 2.4.1 二維水深平均部分水平二維水理模式考量三種邊界條件設定，分別為渠道入流、渠道出流與固體邊界。一般而言，渠道入流邊界條件設定為單位寬度入流量，渠道出流邊界條件則採用水位高程設定。在固體邊界處，沿固體邊界法線方向採不透水邊界條件；而沿固體邊界切線方向可分為滑移與非滑移邊界條件。 2.4.2 流速差異量部分垂直水理模式考量渠道入流、渠道出流、自由液面及底床邊界條件。在渠道入流及渠道出流處假設均勻流邊界條件。自由液面採風剪力邊界條件： u d  s1 v d  s2 、     V    V . 而底床則採用底床剪力邊界條件： u d  b1 v d  b2 、     V    V . 12.

(26) 圖 2.1 正交曲線座標轉換示意圖. 圖 2.2 σ 座標轉換示意圖. 13.

(27) u. u. u  u  u. z. u 圖 2.3 水深方向流速剖面示意圖 ( u 為時間平均流速； u 為水深平均流速； u 為流速差異量). 14.

(28) 第三章數值架構水理部分之整體架構基於垂直水平分離演算概念，首先由水平二維模式求解各水柱之水面高程以及水深平均流速，其中之有效剪應力項須依靠三維流場計算而得到，將求得的水面高程以及水深平均流速代入流速差異量方程式求解水深方向流速剖面，數值計算流程如圖 3.1 所示。三維流場經流速差異量方程式求解後，與二維模式疊代收斂。. 3.1 求解架構. 3.1.1 二維水深平均部分本研究基於隱式雙階分割操作之觀念，將水深平均動量方程式分割成二個步驟(延散步驟及傳播步驟)，並利用隱式數值方法求解。延散步驟求解移流項和擴散項，傳播步驟求解壓力項、底床剪應力項和連續方程式。水理控制方程式先就時間部分離散如下：延散步驟  V     t . n. 1 2.  (V n )V. n. 1 2. . 1. .  T. n. 1 2. (3.1). 傳播步驟  V     t . n 1.  V     t . n. 1 2.   g ( zb  d )n 1 .   V n 1  0. b d. (3.2) (3.3). 式中，V 表示速度向量；T 表示擴散及延散項；n +1 表示(n +1)Δt 時刻之未知變數；Δt = tn+1 −tn；n 表示 nΔt 時刻之已知變數；n + 1/2 表示在(n +1)Δt 與 nΔt 間之未知變數。 . 15.

(29) (3.1)~(3.3)的一般式可表示成延散步驟 u u v u  uv h1 v 2 h2  u      h1  h2   h1h2  h1h2   t . h h  1  h2 1 T12 1 T11  2T12 1  T22 2     T11     dh2   dh1   dh1h2  . . 1  dh1h2.  zs zb zs zb   (h2 11 ) s   (h2 11 )b   (h1 12 ) s   (h1 12 )b    . (3.4). v u v v v  uv h2 u 2 h1       h1  h2   h1h2  h1h2   t . h h  1  h1 1 T12 1 T22  2T12 2  T11 1     T22     dh1   dh2   dh1h2  . . 1  dh1h2.  zs zb zs zb   (h2 12 ) s   (h2 12 )b   (h1 22 ) s   (h1 22 )b    . (3.5). 傳播步驟 u g  ( zb  d )  b1   h1 d t . (3.6). v g  ( zb  d )  b 2   h2 d t . (3.7). 和 h1h2. d    (h2ud )  (h1vd )  0 t  . (3.8). 針對 n +1 時刻的水深值(d n+1)做線性化處理，且僅保留一階項，(3.8)式可改寫成 h1h2.     ( Δd )  d    ( Δd )   1 Δd   1     2 Δd   2   0  1 2 t        . (3.9). h2 n  12 h2 gΔt  zbn 1 d n  h2 gΔt n h1 gΔt n n u   式中，1   d ；1  d ；  ； 1  1d ； 2   C h1 C h2 C C h1     C f (u h1 gΔt  z d  n    ；  2   2 d ； C  1  Δt C h2     Δd  d n 1  d n ； C f  g / c 2 為摩擦係數；c 為 Chezy 係數。. 2 . h1 v C. 1 n 2. . n 1 b. n. 16. n. 1 2 2. )  (v. d. n. n. 1 2 2. ). ；.

(30) 3.1.2 流速差異量部分三維流場之求解係以水柱(water column)為單位發展數值結構，在求解流速差異量時，可由水深判斷有必要求解三維流場之範圍，在需求解範圍內針對個別水柱求解流速剖面。在此概念之下，流速差異量 u 、 v 在離散時，可將方程式區分為垂直方向(沿水深)與水平方向兩方向之物理量，分別以隱示法與疊代法求解。在垂直水理模式方面，每一水柱以隱示法分別求解，將(2.11)與(2.12)式改寫為 ξ 方向：  h1 uv   h1  V  u  u  u u u uv       t  h1   h1h2  h1h2  d 2   . n 1.  M un. (3.10).  M vn. (3.11). η 方向：   h2  V  v  v  v v v uv h2 uv       t  h2  z h1h2  h1h2  d 2   . n 1. 式中，  (   ) u u u u v u v u v u uv h1 M un   f c v  s1 b1       d h1  h1  h2  h2  h2  h1h2   2vv h2 v 2 h2  V  w     ( Horizontal Diffusion in  h1h2  h1h2  dh1  .  ) . n. (3.12).   h2 (   ) v v v v u v u v u v uv M vn    f c u  s 2 b 2       d h2  h2  h1  h1  h1  h1h2   2uu h1 u 2 h1  V  w     ( Horizontal Diffusion in  h1h2  h1h2  dh2  .  ) . n. (3.13). 其中 Horizontal Diffusion in  與 Horizontal Diffusion in  為水平方向之剪應力，如第二章式(2.19)、式(2.20)所述。. 17.

(31) 3.2 數值差分. 3.2.1 二維水深平均部分本模式採用控制體積(control volume)法的觀念來離散控制方程式，控制體積法的基本概念如圖 3.2 所示，其中(a)圖為實際區域，(b)圖為計算區域，E、W、 N、S 表相鄰格點，e、w、n、s 表控制面。模式計算之變數則放置在交錯網格 (staggered grid)中，如圖 3.3。在控制方程式中，除了移流項採用一階精度混合型上風法(hybrid scheme) (Spalding1972)差分外，所有空間差分均採用二階精度的中央差分法。另外，時間項則採用簡單的前向差分方法。中央差分法可表示成  Ψ  Ψe  Ψw    Δ    p. (3.14).  Ψ  Ψn  Ψs    Δ    p. (3.15). 式中， Ψe  0.5 (ΨE  ΨP )  0.5 (Ψi 1, j  Ψi , j ) ； Ψw  0.5 (ΨP  ΨW )  0.5 (Ψi , j  Ψi 1, j ) ； Ψn  0.5 (ΨN  ΨP )  0.5 (Ψi , j 1  Ψi , j ) ； Ψs  0.5 (ΨP  ΨS )  0.5 (Ψi , j  Ψi , j 1 ) ； Ψ 可表為. u 、 v 、 h1 、 h2 、 d 、 zs 和 zb ； i 、 j 分別代表水平格網上任一點之縱向及橫向. 位置。混合型上風法為上風法(upwind scheme)與中央差分法組合而成，當移流效應重要時，採用上風法；移流效應不重要時，則採用中央差分法。至於移流效應重要性的判斷，則採用格網雷諾數(mesh Reynolds number) Rx、Ry 作為判斷的因子，當|Rx|或|Ry|大於 2 時，代表移流效應重要，差分方法採用能反映方向性的上風法；|Rx|或|Ry|小於等於 2 時，移流效應可視為不重要，差分方法採用中央差分法。混合型上風法應用於本研究移流項的處理可表成 uin, j u n  Φ n 1     0.5 h1    h1i , j.   Φin1,1 j  Φin, j 1   Φin, j 1  Φin1,1 j  (1   x )    (1   x )  Δ Δ    . 18.     . (3.16).

(32) vin, j v n  Φ n 1   0.5   h2    h2i , j.   Φin, j 11  Φin, j 1   Φin, j 1  Φin, j 11     (1  ) (1   y )     y     Δ Δ      . (3.17). 其中 0  x   1 1 . 上列諸式中， Rx . uin, j h1i , j Δ. /. Rx  2 0  Rx  2 ；  y   1 1 Rx  2 . ； Ry . vin, j h2i , j Δ. /. Ry  2 Ry  2 Ry  2. (3.18). ；  為流體動力黏滯係數(dynamic. viscosity)； Φ 可表成 u 或 v 。 3.2.2 流速差異量部分垂直水理模式亦採用控制體積法的觀念離散控制方程式，如圖 3.4 所示，E、 W、N、S、 P k 1 、 P k 1 為相鄰格點，e、w、n、s、t、b 為控制面。模式計算之變數亦放置在交錯網格(staggered grid)上，如圖 3.3。(3.10)及(3.11)等式左邊採用 Crank-Nicolson method，其時間差分為二階經度；而等式右邊的 M u 和 M v 之空間差分則採用中央差分法。. 19.

(33) 輸入幾何資料與水理條件. 1.二維水深平均方程式. 出入流邊界條件. 否. 二維模式收斂是. 2.流速差異量方程式. 底床邊界條件紊流黏滯係數. 3.垂直模式計算垂向流速和剪應力項. 4.求解水深方向流速剖面. 否. 否. 垂直模式收斂是. 5.水平二維(延散、傳播). 垂直及二維模式收斂是. 輸出模擬結果. 圖 3.1 模式計算流程圖. 20.

(34) Ni , j +1 n e. P. w Wi -1 , j. E. i +1 , j. i,j. s. y S i , j -1. x. (a) N i , j +1. n. w. W. P i,j. i -1 , j. e. Ei +1 , j. s η. ξ. S i , j -1. (b) 圖 3.2 水平二維模式控制體積法示意圖 i 、 j 分別代表水平格網上任一點之縱向及橫向位置；(a)實際區域；(b)計算區域 (謝德勇，2002). 21.

(35) u & u v & v w & zs. u w. (洪聖翔,2011). 圖 3.3 交錯格網(staggered grid)示意圖. z. T i , j , k +1. η ξ t N W. i -1 , j , k. P. w. k. n i,j,k. s. e. i , j +1 , k. E i +1 , j , k. S i , j -1 , k b. Bi , j , k -1. (洪聖翔,2011). i 、 j 分別代表水平格網上任一點之縱向及橫向位置； k 代表在垂向格網上的位置. 圖 3.4 垂直模式控制體積法示意圖(計算區域). 22.

(36) 第四章結果與討論本研究係以大渦紊流模擬(LES)結合 Smagorinsky(1963)之次網格模式建構一非零方程之紊流模式，在本章所要探討的是計算網格大小對次網格模式模擬結果之影響，除了與實驗量測數據驗證，並以另兩種零方程紊流模式相互比較。由於大渦紊流模式屬時間相依，在迭代過程須確認結果已達統計不變性，也就是經過數個大渦模擬時間 t LES (large eddy turnover time)之後，前後取樣值已在極小的誤差範圍內，其中大渦模擬時間可定義為(Y.Cheng et al. 2003) t LES .  U*. (4.1). 其中  為渠道一半之長度， U * 為剪力速度。另外，在數值已達收斂值後之模擬結果當中，取數個 t (time step)之數據將其作時間平均，所得之結果即為大渦紊流模式最終之結果，其數值收斂所需之時間長短須經由案例測試而定(案例部分會再次詳述)。若依大渦模擬之理論，大尺度渦流係以直接求解之方式，網格尺寸係依據流場本身之雷諾數(Reynolds number)來給定(E.M. Matos and M.Mori 2012)，其定義為   Re 3/ 4. (4.2). 其中  為最小網格尺寸。由此可見，模擬所需的網格數須為 Re9/4 之多，但以目前之電腦配備而言，若以如此的網格大小來計算恐耗費太多時間。於此，本研究在選擇網格尺寸上先以較疏之網格作模擬測試，而後再作不同程度之加密，待模擬結果與理論值相接近時，便可確立此案例所適用之網格大小。 (Lindsey,Jasim and Hanif 2013). 4.1 突擴渠道流速場分析案例本研究採用 Babarutsi et al.(1989)之環流流場實驗案例，模擬案例為實驗 3，其試驗水槽平面佈置如圖 4.1 所示，均勻入流從左側上方流入試驗渠道，在流經渠道突擴處時即會產生環流流場，其中 Lr 表環流流場之環流長度(recirculating length)。上游均勻入流流速 U 0 為 0.145m/s，下游水深 H 為 0.0819m，渠道的突擴長度 d e 為 0.305m，突擴後渠道的寬度 B 為 0.61m，渠道縱向坡度 S 為 0.0098， 23.

(37) 摩擦係數 C f 為 0.00527，福祿數 Fr 及雷諾數 Re 分別為 0.162、10 4 。圖 4.2 為沿水流方向各量測斷面之無因次水平流速分佈圖，圖中虛線部分內為環流產生區域，文獻當中指出此實驗之目測環流長度約為 2.4m。圖 4.3 為實驗量測數據之校正示意圖，由於此實驗量測儀器之限制而無法量測出反方向之流速值，須將反曲之量測值修正為正確之流速值。本案例模擬之側壁採用非滑移邊界條件，上游邊界條件為固定入流量，下游邊界條件為給定水位資料。並採用四種不同大小之計算網格，分別為 361 62  19 均勻網格( 、 方向均為 0.01m)、721 123  19 均勻網格(  、 方向均加密一倍至 0.005m)、721 184  19 非均勻網格(  方向維持在 0.005m， 方向在環流產生區域加密至 0.0025m)、 721 367 19 非均勻網格(  方向維持在 0.005m， 方向在環流產生區域加密至 0.001m)，其中三數字依序為三方向之網格數，分別為主流  方向，側方向 ，水深  方向，網格佈置圖如圖 4.4 至圖 4.7 所示。此案例模擬之物理條件整理如表 4.1。比照前述取樣之理論，此案例之 t LES 約為 11.7 秒，在模式計算 5400 秒(約 460 個 t LES )之後進行取樣，取樣方式為每隔 5000 秒計算時間取樣一次，共取十次作時間平均，當作最終之結果。圖 4.8 至圖 4.11 分別為四種不同大小計算網格下之無因次水平流速模擬結果與實驗值比較圖，模擬結果具次網格模式及零方程式 I、II 三種，圖中縱軸為渠寬，橫軸為主流方向之水深平均流速與入流流速之比值。圖 4.8 為 361 62 之均勻網格，主流方向與側方向之網格大小相同，圖 4.8(a)到(e)分別為沿主流方向之各量測斷面之比較結果，可看出次網格模式之環流產生區域較實驗值為後，在圖 4.8(d)才出現較明顯的負流速現象，整體而言其結果與零方程式 II 差別不大；而零方程式 I 在環流區域所反應出之負流速值較符合實驗量測值。圖 4.9 為 721123 之均勻網格，  、 方向網格大小加密至 0.005m，可明顯看出網格大小. 對次網格模式之影響，在環流產生過程之斷面均有不錯的模擬結果。圖 4.10 為 721 184 之非均勻網格，如圖 4.6 所示，在環流產生區域進行加密，將側方向網. 格大小加密至 0.0025m，如此使得次網格模式模擬結果與實驗量測值更為接近；而在加密之後，由圖 4.10(b)、(c)可看出零方程式 I 及 II 在環流區域之預測結果與量測值出現誤差，環流之效應開始減弱。圖 4.11 為 721 367 之非均勻網格，如圖 4.7 所示，將環流區域之側方向網格大小加密至 0.001m，由圖 4.10 與圖 4.11 可看出，次網格模式在環流區域部分之模擬結果已不太有明顯改變；反觀零方 24.

(38) 程模式在網格加密之條件下，模擬結果反而不如先前條件理想，原因可能為垂直方向網格未同時加密所致。圖 4.12 至圖 4.14 分別為次網格模式、零方程式 I 及 II 在  方向網格大小相同之下， 方向加密不同程度之模擬結果比較，可以更清楚地比較出網格加密對結果之影響。由圖 4.13 及圖 4.14 可看出，零方程式 I、II 同樣在距上游 250cm 處之模擬結果與量測值出現較大差距，而在縱方向 150cm 處均可看出網格加密反而使環流區域之負方向流速效應減小；而由圖 4.12，在次網格模式的結果當中，環流產生區域之網格加密程度愈大，愈能適當反應出環流之效應，在在說明網格尺寸對於次網格紊流模式預測結果之重要性。. 4.2 直渠道渦度場分析案例本研究採用 Nezu and Rodi(1986)之實驗案例，模擬案例為實驗 8 (寬深比為 16.9)，其試驗水槽佈置如圖 4.15 所示，渠道長度 20m，渠寬 B 為 0.6m，下游水深 H 為 0.036m，平均流速為 0.417m/s，渠道坡度 S 為 0.00142，福祿數 Fr 及雷諾數 Re 分別為 0.707、4.8  104 ，其實驗量測斷面為距上游入口 18m 處之中間水柱，探討垂向分佈資訊。本案例模擬之側壁採用非滑移邊界條件，上游邊界條件為固定入流量，下游邊界條件為給定水位資料。圖 4.16 為明渠紊流之垂向流速分佈示意圖，一般可將其速度剖面分為內層(inner layer)及外層(outer layer)，圖中內層垂直壁面向上又可依序分為黏性次層(viscous sublayer)、緩衝層(buffer layer)及紊流對數層(log-law region)，邊界層之各層定義範圍示意圖如圖 4.17 所示；外層為紊流分佈區，在此展現出紊流之特性。在高雷諾數的流場中，靠近固體邊界之邊界層較薄(也就是圖 4.16 中之內層範圍)，若要精確地計算邊界層內之流況需要非常細小之網格，在三維流場計算中，若忽略該邊界層之影響，則會造成模擬結果不符合實際流場之狀況。 Anderson(1995)提出在模擬紊流流場時，須將大量且密集之格網點分佈於速度梯度較大之區域內，如此不但能夠減小截斷誤差(truncation error)，亦可獲得流場中真正之物理特性。文獻中以一平板邊界層流之例子說明邊界層內格網點分佈之重要性，如圖 4.18 所示。為了能明確分辨出三方向網格大小相互之影響，本案例採用三種不同大小之計算網格，分別為 74  5 12 非均勻網格、 363  13  23 及 725  25  25 之非均勻網. 25.

(39) 格，網格佈置圖如圖 4.19 至圖 4.21 所示。此案例模擬之物理條件整理如表 4.2。比照先前取樣之理論，此案例之 t LES 約為 20.1 秒，在模式計算 18000 秒(約 900 個 t LES )之後進行取樣，取樣方式為每隔 6000 秒計算時間取樣一次，共取十次作時間平均，當作最終之結果。圖 4.22 至圖 4.24 分別為三種不同大小計算網格下於量測斷面中間水柱之縱向流速在垂向分佈之模擬結果與實驗值比較圖，模擬結果具次網格模式及零方程式 I、II 三種，圖中縱軸為水深，橫軸為縱向流速。圖 4.22 為 74  5 12 之非均勻網格，由於近底床之格網點僅少許落於內層上方，故三種模擬結果之趨勢較和層流之流速分佈類似，未能模擬出紊流之現象。圖 4.23 為 363  13  23 之非均勻網格，主流方向與側方向之網格大小相同，並在近底床部分之網格加密，使加密之格網點增加落於過渡層及紊流對數層之範圍內，少數落於黏滯邊界層內。由圖 4.22 與圖 4.23 可看出，次網格模式以及零方程式 II 之模擬結果和先前有明顯差異，此時次網格模式之結果稍微顯現出紊流特性，但尚未能主宰成紊流流場；而零方程式 II 之模擬結果在近底床部分與實驗量測值略為接近。圖 4.24 為 725  25  25 之非均勻網格，主流方向與側方向之網格大小相同，較先前加密一倍至 0.025m，並在近底床部分加密，使加密之格網點增加落於黏滯邊界層內，可看出次網格模式流速分佈之趨勢已近似於紊流之分佈，在近底床部分也有不錯之結果；由圖 4.23 與圖 4.24 可看出，零方程式 II 之結果在底床部分因為加密網格而更近於實驗量測值，但整體之流速分佈已無明顯改變；而零方程式 I 仍保持在層流之流速分佈。圖 4.25 至圖 4.27 分別為次網格模式、零方程式 I 及 II 在不同計算網格大小之下沿主流方向上各斷面之縱向流速分佈圖。以網格數為 725  25  25 為例，可看出距離上游愈遠，其邊界層內之流速梯度愈大，邊界層成長有愈厚之趨勢，此結果與實際物理現象吻合。圖 4.28 至圖 4.30 為次網格模式、零方程式 I 及 II 在不同計算網格大小之下於量測斷面之渦度圖，圖中縱軸為該比較點之水深除以總水深，橫軸為該比較點距側壁之長度除以總水深。此渦度圖之縱橫軸比例為 16.9，為實驗案例之寬深比，但從圖中較難看出渦度的細部結構，所以將渦度圖另繪製成縱橫軸等比例表示，如圖 4.31 至圖 4.39 所示。圖 4.31 至圖 4.33 為次網格模式在不同計算網格大小下於量測斷面(距入口 18m 處)之 v、w 合向量圖以及渦度圖，圖中縱軸為該比較點之水深除以總水深，橫軸為該比較點距側壁之長度除以渠寬，繪製. 26.

(40) 成一正方形之示意圖。由圖 4.31 可看出在網格數為 74  5 12 之狀況下，合向量值幾乎趨近於零，也無法模擬出具有渦流現象之流場。由圖 4.32 及圖 4.33 可看出網格愈密之情況下，所反應出之速度向量值愈大，渦度峰值也愈大，而渦度圖之渦流形狀也愈趨於對稱。另外，在網格數為 363  13  23 及 725  25  25 條件下之合向量圖中，會出現看似於一正一負之向量，原因為本案例渠道之寬深比為 16.9，當以正方形之示意圖表示時會有渦流形狀被壓縮之可能，故容易出現此一現象；而在壁面附近之流場現象不明顯，有一可能原因為近底床之網格加密過後之尺寸與側方向網格之尺寸比例失衡，小尺寸渦流在過程中易被過濾掉。圖 4.34 至圖 4.36 為零方程式 I 在不同計算網格大小下於量測斷面之 v、w 合向量圖以及渦度圖；圖 4.37 至圖 4.39 為零方程式 II 在不同計算網格大小下於量測斷面之 v、w 合向量圖以及渦度圖。由圖 4.34 至圖 4.39 可看出兩種零方程式均能模擬出流場之垂向斷面具對稱之特性。圖 4.40 為次網格模式在不同計算網格大小下於量測斷面之紊流黏滯係數分佈圖，圖中縱軸為該比較點之水深除以總水深，橫軸之無因次參數為紊流黏滯係數值除以水深以及剪力速度，在次網格模式計算下之紊流黏滯係數可看出較不具規則性，但在網格愈密之情況下，上下數值差距之程度愈小，整體之分佈愈趨均勻。圖 4.41 為零方程式 I 在不同計算網格下於量測斷面之紊流黏滯係數分佈圖，係一沿水深均勻分佈之函數，其結果也較不受網格大小之影響。圖 4.42 為零方程式 II 在不同計算網格下於量測斷面之紊流黏滯係數分佈圖，係一沿水深呈現近似拋物線分佈之函數，大約在水深一半處為最大值。圖 4.43 至圖 4.45 為計算網格為 725  25  25 之下的紊流強度之垂向分佈圖，紊流強度定義為紊流擾動量除以剪力速度，其中紊流擾動量 uˆi 定義為 uˆi  (ui ) 2. (4.3). 其中， ui 為各方向之次網格擾動量， i 代表 1、2、3，分別為  、 、  三方向。由於此網格數之下的流速分佈較近似紊流特性，故以此網格數為例作圖。圖 4.43 為主流方向紊流強度之垂向分佈圖，圖中縱軸為主流方向紊流擾動量與剪力速度之比值，橫軸為該比較點之水深除以總水深。圖 4.44 為側方向之紊流強度之垂向分佈圖，圖中縱軸為側方向之紊流擾動量與剪力速度之比值，橫軸之定義同圖 4.43。由圖 4.43 及圖 4.44 可看出主流方向之紊流強度遠大於側方向之紊流強度，從整體趨勢來看，靠近壁面會有一極值出現而後趨於平緩，此結果與理論相符。圖 4.45 為主流方向之紊流強度在近壁區之垂向分佈圖，縱軸為主流方 27.

(41) 向紊流擾動量與剪力速度之比值，橫軸之無因次參數為該比較點之水深乘上剪力速度除以層流黏滯係數值，由於近壁區之格網點較少，故無法將近壁理論(Law of the wall)之特性完整呈現，僅能展示落於邊界層內之格網點數給作參考。另外，本研究也以 Nezu and Rodi(1986)的直渠道實驗 7 作模擬測試(寬深比為 12.4)，試驗渠道佈置如圖 4.15 所示，下游水深為 0.049m，平均流速為 0.858m/s，渠道坡度 S 為 0.00409，福祿數 Fr 及雷諾數 Re 分別為 1.240 及13.1104。以次網格模式以及零方程式 I、II 分別測試不同寬深比案例下之計算網格尺寸以及所耗費之 CPU 時間均列於表 4.3。實驗 7(寬深比為 12.4)案例之模擬結果在模式總計算時間為 45000 秒時還未達收斂階段，且在相同計算網格大小下時，其所耗費之 CPU 時間已遠大於實驗 8(寬深比為 16.9)，原因可能為本模式理論係基於淺水波、靜水壓之假設，對於寬深比較小之案例會受限於模式假設，除了數值不易穩定之外，模擬之結果也易有較大誤差。. . . 28.

(42) 表 4.1 突擴渠道案例之模擬物理條件 . . .  (m).  (m). z (m). 361. 62. 19. 0.01. 0.01. 4.55  103. Babarutsi. 721. 123. 19. 0.005. 0.005. 4.55  103. et.al(1989). 721. 184. 19. 0.005. 加密至 0.0025. 4.55 103. 721. 367. 19. 0.005. 加密至 0.001. 4.55  103. 實驗案例. 表 4.2 直渠道案例之模擬物理條件 . . .  (m).  (m). z (m). Nezu. 74. 5. 12. 0.25. 0.15. 加密至 1.8 103. and. 363. 13. 23. 0.05. 0.05. 加密至 7.2 106. Rodi(1986). 725. 25. 25. 0.025. 0.025. 加密至 5.4 106. 實驗案例. . 29.

(43) 表 4.3 直渠道案例不同寬深比(B/H)條件之模擬耗費時間 ※模式總計算時間 T=45000s Case B/H. LES. 1. 2. 16.9. 3. 3. TypeI. TypeII. LES. 1. 2. 紊流模式. 12.4. TypeI. TypeII. 網格數. CPU time(s). 74  5 12. 59058.51s. 363 13  23. 225433.37s(約 2.6 天). 725  25  25. 1244967.41s(約 14.4 天). 74  5 12. 29250.31s. 363 13  23. 73557.17s. 725  25  25. 814025.84s(約 9.42 天). 74  5 12. 30127.54s. 363 13  23. 76850.55s. 725  25  25. 814890.82s(約 9.43 天). 74  5 12. 98948.16s(約 1.1 天). 363 13  23. 482379.04s(約 5.6 天). 725  25  25. 3193860.15s(約 36.9 天). 74  5 12. 31035.94s. 363 13  23. 252999.14s(約 2.9 天). 725  25  25. 258914.15s(約 2.99 天). 74  5 12. 36799.71s. 363 13  23. 300217.63s(3.47 天). 725  25  25. 338202.55s(3.9 天). ※註：本研究使用 CPU 為 i7-3930k(3.8G Hz.)/6 cores/Max memory:64GB i7-2600k(3.8G Hz.)/4 cores/ Max memory:32GB. 30.

(44) 圖 4.1 Babarutsi(1989)實驗水槽設置圖. 圖 4.2 Babarutsi(1989)環流案例實驗量測流速分佈圖. 圖 4.3 Babarutsi(1989)實驗量測數據校正示意圖. 31. .

(45) .  ( m). . 圖 4.4 Babarutsi(1989)計算網格 361  62 佈置圖.  ( m). 圖 4.5 Babarutsi(1989)計算網格 721  123 佈置圖. 32.

(46) .  ( m). . 圖 4.6 Babarutsi(1989)計算網格 721  184 佈置圖.  ( m). 圖 4.7 Babarutsi(1989)計算網格 721x  367 佈置圖. 33.

(47) x=100cm. x=150cm 70. 60. 60. 60. 50. 50. 50. 40. 40. 40. 30 20. width(cm). 70. width(cm). width(cm). x=50cm 70. 30 20. 30 20. 10. 10. 10. 0. 0. 0. -0.3 0.0 0.3 0.6 0.9 1.2. -0.4. 0.0. 0.4. 0.8. 1.2. -0.4. 0.0. 0.4. U/U0. U/U0. U/U0. (a). (b). (c). 70. 60. 60. 50. 50. 40. 40. 30 20. 30 20. 10. 10. 0. 0. -0.5. 0.0. 0.5. 1.0. 1.2. x=250cm. 70. width(cm). width(cm). x=200cm. 0.8. 1.5. -1.0 -0.5 0.0 0.5 1.0 1.5. U/U0. U/U0. (d). (e). 圖 4.8 突擴案例 U / U 0 於各量測斷面之比較圖(網格數為 361 62 ) 實驗值(○)；次網格模式模擬結果(─)；零方程式 I 模擬結果(- -)；零方程式 II 模擬結果(─ ─). 34.

(48) x=100cm. x=150cm 70. 60. 60. 60. 50. 50. 50. 40. 40. 40. 30 20. width(cm). 70. width(cm). width(cm). x=50cm 70. 30 20. 30 20. 10. 10. 10. 0. 0. 0. -0.3 0.0 0.3 0.6 0.9 1.2. -0.4. 0.0. 0.4. 0.8. 1.2. -0.4. 0.0. 0.4. U/U0. U/U0. U/U0. (a). (b). (c). 1.2. x=250cm. 70. 70. 60. 60. 50. 50. 40. 40. width(cm). width(cm). x=200cm. 0.8. 30 20. 30 20. 10. 10. 0. 0. -0.8 -0.4 0.0 0.4 0.8 1.2. -0.8 -0.4 0.0 0.4 0.8 1.2. U/U0. U/U0. (d). (e). 圖 4.9 突擴案例 U / U 0 於各量測斷面之比較圖(網格數為 721 123 ) 實驗值(○)；次網格模式模擬結果(─)；零方程式 I 模擬結果(- -)；零方程式 II 模擬結果(─ ─). 35.

(49) x=100cm. x=150cm 70. 60. 60. 60. 50. 50. 50. 40. 40. 40. 30 20. width(cm). 70. width(cm). width(cm). x=50cm 70. 30 20. 30 20. 10. 10. 10. 0. 0. 0. -0.3 0.0 0.3 0.6 0.9 1.2. -0.8 -0.4 0.0 0.4 0.8 1.2. -0.8 -0.4 0.0 0.4 0.8 1.2. U/U0. U/U0. U/U0. (a). (b). (c). x=250cm 70. 60. 60. 50. 50. 40. 40. width(cm). width(cm). x=200cm 70. 30 20. 30 20. 10. 10. 0. 0. -0.8 -0.4 0.0 0.4 0.8 1.2. -0.8. -0.4. 0.0. U/U0. U/U0. (d). (e). 0.4. 0.8. 圖 4.10 突擴案例 U / U 0 於各量測斷面之比較圖(網格數為 721 184 ) 實驗值(○)；次網格模式模擬結果(─)；零方程式 I 模擬結果(- -)；零方程式 II 模擬結果(─ ─). 36.

(50) x=100cm. x=150cm 70. 60. 60. 60. 50. 50. 50. 40. 40. 40. 30 20. width(cm). 70. width(cm). width(cm). x=50cm 70. 30 20. 30 20. 10. 10. 10. 0. 0. 0. -0.3 0.0 0.3 0.6 0.9 1.2. -0.8 -0.4 0.0 0.4 0.8 1.2 1.6. -0.8 -0.4 0.0 0.4 0.8 1.2. U/U0. U/U0. U/U0. (a). (b). (c). x=250cm 70. 60. 60. 50. 50. 40. 40. width(cm). width(cm). x=200cm 70. 30 20. 30 20. 10. 10. 0. 0. -0.6 -0.3 0.0 0.3 0.6 0.9. -0.8 -0.4 0.0 0.4 0.8 1.2. U/U0. U/U0. (d). (e). 圖 4.11 突擴案例 U / U 0 於各量測斷面之比較圖(網格數為 721 367 ) 實驗值(○)；次網格模式模擬結果(─)；零方程式 I 模擬結果(- -) ；零方程式 II 模擬結果(─ ─) 37.

(51) x=50cm. x=100cm. x=150cm. x=200cm. x=250cm. 70 60. width(cm). 50 40 30 20 10 0. -0.4 0.0 0.4 0.8 1.2-0.8 -0.4 0.0 0.4 0.8 1.2-0.8 -0.4 0.0 0.4 0.8 1.2-0.8 -0.4 0.0 0.4 0.8 1.2-0.8 -0.4 0.0 0.4 0.8 1.2. U/U0. U/U0. U/U0. U/U0. U/U0. 圖 4.12 突擴案例 U / U 0 於各量測斷面不同計算網格之次網格模式結果比較圖實驗值(○)；721  123(─)；721  184(─ ─)；721  367(──). x=50cm. x=100cm. x=150cm. x=200cm. x=250cm. 70 60. width(cm). 50 40 30 20 10 0. -0.5 0.0 0.5 1.0 1.5-0.5 0.0 0.5 1.0 1.5-0.5 0.0 0.5 1.0 1.5-1.0 -0.5 0.0 0.5 1.0 1.5-1.0 -0.5 0.0 0.5 1.0 1.5. U/U0. U/U0. U/U0. U/U0. U/U0. 圖 4.13 突擴案例 U / U 0 於各量測斷面不同計算網格之零方程式 I 結果比較圖實驗值(○)；721  123(─)；721  184(─ ─)；721  367(──). 38.

(52) x=50cm. x=100cm. x=150cm. x=200cm. x=250cm. 70 60. width(cm). 50 40 30 20 10 0. -0.5 0.0 0.5 1.0 1.5-0.5 0.0 0.5 1.0 1.5-0.4 0.0 0.4 0.8 1.2-1.0 -0.5 0.0 0.5 1.0-0.8 -0.4 0.0 0.4 0.8 1.2. U/U0. U/U0. U/U0. U/U0. U/U0. 圖 4.14 突擴案例 U / U 0 於各量測斷面不同計算網格之零方程式 II 結果比較圖實驗值(○)；721  123(─)；721  184(─ ─)；721  367(──). 39.

(53) (a). (b). 圖 4.15 Nezu and Rodi(1986)實驗案例設置圖 (摘錄自 Nezu and Rodi 1986). 40.

(54) 圖 4.16 明渠紊流之垂向流速分佈示意圖. z/H . . Z (. zU *. . ). 圖 4.17 邊界層之各層定義範圍示意圖 (摘錄自 Stephen B.Pope, Turbulent Flows). 41.

(55) (a). (b) 圖 4.18 邊界層內格網點分佈對模擬結果之影響示意圖 (a) 邊界層內無格網點 (b)邊界層內設有格網點 (摘錄自 Anderson 1995). 42.

(56) 圖 4.19 Nezu and Rodi(1986)計算網格 74  5 12 佈置圖. Log-scale. 圖 4.20 Nezu and Rodi(1986)計算網格 363 13  23 佈置圖. 43.

(57) Log-scale. 圖 4.21 Nezu and Rodi(1986)計算網格 725  25  25 佈置圖. 44.

(58) 0.04. z(m). 0.03. 0.02. 0.01 experiment LES results TypeI results TypeII results 0.00 0.0. 0.2. 0.4. 0.6. 0.8. U(m/s). 圖 4.22 直渠案例於量測斷面中間水柱縱向流速分佈之比較圖(網格數為 74  5 12 ). 0.04. z(m). 0.03. 0.02. 0.01 experiment LES results TypeI results TypeII results. 0.00 0.0. 0.2. 0.4. 0.6. 0.8. U(m/s). 圖 4.23 直渠案例於量測斷面中間水柱縱向流速分佈之比較圖(網格數為 363  13  23 ). 45.

(59) 0.04. z(m). 0.03. 0.02. 0.01 experiment LES results TypeI results TypeII results. 0.00 0.0. 0.2. 0.4. 0.6. 0.8. U(m/s). 圖 4.24 直渠案例於量測斷面中間水柱縱向流速分佈之比較圖(網格數為 725  25  25 ). 46.