國
立
交
通
大
學
土木工程學系
博
士
論
文
擬似三維分層積分淺水波模式之發展
Development of a Semi-3D Layer-Integrated
Shallow Water Wave Model
研 究 生：洪夢祺
指導教授：楊錦釧 博士
擬似三維分層積分淺水波模式之發展
Development of a Semi-3D Layer-Integrated
Shallow Water Wave Model
研 究 生：洪夢祺 Student：Meng-Chi Hung 指導教授：楊錦釧 Advisor：Jinn-Chuang Yang
國 立 交 通 大 學
土 木 工 程 學 系
博 士 論 文
A DissertationSubmitted to Department of Civil Engineering College of Engineering
National Chiao Tung University in Partial Fulfillment of the Requirements
for the Degree of Doctor of Philosophy
in
Civil Engineering June 2007
Hsinchu, Taiwan, Republic of China
擬似三維分層積分淺水波模式之發展
學生：洪夢祺 指導教授：楊錦釧 國立交通大學土木工程學系博士班摘
要
本研究發展一交錯網格有限體積擬似三維分層積分淺水波模式，提出二 次形狀函數近似層內速度分布以確保層介面之流速與剪力連續，毋需二維水 深平均子模式決定水面高程，可有效降低模式計算量。模式之驗證首先透過 封閉水池風吹之垂直環流假設案例與解析解比較，除觀察模式對層數之敏感 度外，亦檢視模式對環流現象水深方向流速剖面發展過程並反映水面坡降變 化之機制，及垂直水深方向黏滯性分布改變對流速剖面形狀之影響；此外， 由於底床坡降與粗糙度密切影響明渠流況，經由明渠流假設案例之功能測試 檢驗，本研究假設均勻流條件推導渦流滯性與水深、粗糙度、剪力速度之關 係，模擬出水深方向均勻紊流強度之拋物線流速剖面，並應用實驗觀測水深 方向紊流強度拋物線分布，模擬出對數型流速剖面。 關鍵字：三維模式、分層積分、淺水波模式、風剪環流
Development of a Semi-3D Layer-Integrated Model
for Shallow Water Flow
Student: Meng-Chi Hung Advisor: Jinn-Chuang Yang
Department of Civil Engineering National Chiao Tung University
ABSTRACT
A novel semi-three-dimensional (semi-3D) layer-integrated approach was proposed in this study for the shallow water wave computation. A quadratic shape function proposed to approximate the velocity distribution in layer ensures the continuity of velocities and shear stresses at interfaces. To discretize the governing equations, the finite volume formulation with staggered grid is used. As the
two-dimensional (2D) sub-model for locating the free surface is not needed in this approach, the computational time consumption has been dramatically reduced. The model was verified through the wind-induced circulation in a closed basin by
comparing the velocity profiles to the analytical solution. The formation of the velocity profiles due to change of viscosity distribution and the evolution process of water surface slope were also investigated. Further, the eddy viscosity is estimated by a function related to discharge, water depth, and Manning’s n under the assumption of parabolic distribution of longitudinal velocity along the vertical direction for open channel flows. Two designed hypothetically cases were conducted to examine the proposed eddy viscosity relation and to demonstrate the capabilities of the proposed model.
Keywords: three-dimensional model, layer-integrated, shallow water wave model, wind-induced circulation
誌 謝
承蒙恩師 楊教授錦釧的悉心指導與諄諄教誨，於治學與為人處世多所 啟發，使本文得以順利完成，謹致由衷謝忱與敬意。感謝口試委員台大農工 系許教授銘熙、台大土木系黃教授良雄、成大水利系蔡教授長泰、成大水利 系賴教授泉基、交大應數系賴教授明治細心匡正、惠賜卓見，使本文更臻完 善。 感謝大學長張教授哲豪與同窗羅正工程司克信不間斷地鼓勵，謝博士德 勇與蔡博士東霖兩位學長於論文的研討與寶貴意見，以及研究室伙伴們與求 學過程所有同學伴我一起學習、一起成長，並陪我度過低潮而能持續論文研 究，各位的扶持與砥礪永銘我心。 爸媽的無怨無悔的照顧與關懷，哥哥、姊姊的啟發與鼓勵，是求學與研 究過程中最大的動力與精神支柱，謹以本文與父母兄姊分享這份喜悅，並表 無限感恩之心。目 錄
摘要 ... i Abstract ... ii 誌謝 ... iii 目錄 ... iv 表目錄 ... vi 圖目錄 ... vii 符號表 ... ix 第一章 緒 論...1 1.1 研究背景與動機 ...1 1.2 文獻回顧 ...3 1.3 研究方法 ...7 1.4 章節介紹 ...8 第二章 分層積分淺水波模式理論架構 ...10 2.1 雷諾平均控制方程式 ...10 2.2 分層積分控制方程式 ...12 2.3 邊界條件 ...17 2.4 明渠流之垂直紊流黏滯係數推估式 ...17 第三章 分層模式(一)－不考慮層介面連續 ...23 3.1 有限體積法與數值差分式 ...23 3.2 數值演算方法 ...26 3.3 模式檢討 ...27 第四章 分層模式(二)－考慮層介面連續 ...32 4.1 層介面剪力連續條件限制式 ...32 4.2 數值架構 ...384.2.1 有限體積法 ...38 4.2.2 數值差分式 ...40 4.2.3 水柱代數方程組 ...40 4.2.4 求解步驟 ...42 4.2.5 數值穩定限制條件 ...44 4.3 小結 ...44 第五章 風剪垂直環流場驗證 ...47 5.1 層數對流速剖面之影響 ...47 5.2 模式對流場物理現象之反應機制 ...49 5.3 黏滯係數分佈型態對流速剖面之影響 ...50 第六章 明渠流之計算與比較 ...61 6.1 紊流黏滯係數關係式之驗證 ...61 6.2 紊流黏滯係數分佈型態對流速剖面之影響 ...62 6.3 底床層厚度對流速剖面之影響 ...64 6.4 底床邊壁效應對流速剖面之影響 ...65 6.5 風剪對流速與水位剖面之影響 ...66 第七章 結論與建議...83 7.1 結論 ...83 7.2 建議 ...84 參考文獻 ...85 附錄 修正湯瑪斯法...89
表 目 錄
圖 目 錄
圖 2.1 分層模式層切割方式與變數配置示意圖 ... 22 圖 3.1 交錯網格有限體積法變數之控制體積於 x-y 於平面之配置圖... 30 圖 3.2 分層模式(一)求解步驟流程圖... 31 圖 4.1 分層模式(二)求解步驟流程圖... 46 圖 5.1 不同層數對風剪環流場流速剖面之影響 ... 52 圖 5.2 不同層數對風剪環流場水面坡降之影響 ... 53 圖 5.3 風剪環流場流速剖面之發展過程(縱向中間長度處) ... 54 圖 5.4 風剪環流場流速剖面之發展過程(靠近背風邊界處) ... 55 圖 5.5 風剪環流場流速剖面之發展過程(靠近迎風邊界處) ... 56 圖 5.6 風剪環流場水面坡降之發展過程 ... 57 圖 5.7 風剪環流場不同水深處之流速發展過程 ... 58 圖 5.8 不同型態之黏滯係數分佈對流速剖面之影響(均勻 vs 拋物線型)... 59 圖 5.9 不同型態之黏滯係數分佈對流速剖面之影響(均勻 vs 梯型)... 60 圖 6.1 不同曼寧 n 值對直線渠道之水位剖面之影響 ... 69 圖 6.2 不同曼寧 n 值對直線渠道之流速剖面之影響(渠道中間長度處) .... 70 圖 6.3 三種不同黏滯係數分佈型態 ... 71 圖 6.4 黏滯係數分佈型態對流速剖面之影響 ... 72 圖 6.5 黏滯係數分佈型態對無因次流速剖面之影響 ... 73 圖 6.6 使用層數對流速剖面之影響(黏滯係數均勻分佈) ... 74 圖 6.7 使用層數對流速剖面之影響(黏滯係數拋物線型分佈) ... 75 圖 6.8 使用層數對無因次流速剖面之影響(黏滯係數拋物線型分佈) ... 76 圖 6.9 底床邊界速度對流速剖面之影響 ... 77 圖 6.10 底床邊界速度對無因次流速剖面之影響 ... 78圖 6.12 風剪對流速剖面之影響(黏滯係數拋物線型分佈) ... 80 圖 6.13 風剪對水位剖面之影響(黏滯係數均勻分佈) ... 81 圖 6.14 風剪對水位剖面之影響(黏滯係數拋物線型分佈) ... 82
符 號 表
一、變數 h= 層厚度 n= 曼寧粗糙度係數 p= 靜水壓力 u= x-方向流速 u* = 剪力速度 v= y-方向流速 A= 水柱動量離散代數方程組係數矩陣 C= 任意常數 F= 層介面通量 H= 水深 K= 總層數 S= 坡降(底床或能量) Su, Sv= 動量方程式之源項 U= 水深平均流速 α= Elder 黏滯係數關係式之係數 γ= 水體單位體積重量 ηB = 底床高程 ηk = 層介面高程 ηs = 水面高程 μ= 動力黏滯係數 ν= 運動黏滯係數 ρ= 水體密度τs = 風剪應力 二、下標 e,s,w,n= 控制體積之各方向控制面 h= 水平方向 k= 位於第 k 層與第 (k-1) 層間之層介面 k+1/2= 第 k 層之層平均物理量 m= 啞變數，以連續方程式計算垂向速度時使用 v= 垂直方向 E,S,W,N= 與該計算控制體相鄰之各控制體中心點 P= 該計算控制體中心點
第一章 緒 論
1.1 研究背景與目的
自然界流動之水體表面與大氣相通，液面壓力保持與大氣相等，液體內 部壓力不平衡時可反映在液面之升降，為管流所沒有之自由度，故名為自由 液面流，自然界中的水體流動，舉凡漫地流、河川、湖泊、水庫、河口、海 岸、洋流等，多為自由液面流；相對於深海洋流，河川、湖泊與海岸水體之 運動受限於底床邊界，其運動行為近乎水平方向，即所謂淺水波運動，淺水 波運動為重力所帶動而加速，當重力與底床阻力達到平衡時，則流態不再改 變而形成均勻流，然而天然河道、湖泊或水庫由於幾何邊界不規則之影響， 水體因重力與邊界阻力不平衡、或匯入不同密度之水體而形成各種複雜流 況，如河川彎道水面超高造成之二次流、水庫挾砂水流因密度不同而造成之 異重流或密度層變流、河口因海水鹽度而形成之垂直環流等，均直接或間接 影響我們日常生活。此外，在水利工程的領域裡還有許多淺水波運動之流場 為工程師欲瞭解或解決之問題，例如潰壩波的傳遞、高程差所引起的水躍、 流動水體遭遇水工結構物所激起的碎波、地震時因搖晃所引起之水體震動、 海洋水工物結構因地形變化所產生之淺化碎波作用、明渠變量流之水面急劇 變化、水流經過水工結構物如橋墩、潛堰、丁壩…等之流場。 為解決水利工程相關問題而尋求改善方案時，一般都先經由物理模型試 驗或數值模擬來檢討改善方案的功效，物理模型一直是科學家研究流體相關較易被認同，但因其需要較多的經費、較大的試驗場所、較長的試驗時間及 龐大人力的配合，使其在試驗進行時缺乏變通的彈性，一般均僅能針對一、 兩個既定的方案進行模擬試驗，因此，常用之物理模型試驗係以縮小尺度之 模型，觀察並量測流體的物理性質及運動特性，然而受限於模型相似律問 題，對於輸砂或水躍等複雜流況，量測結果往往易受尺度效應(scale effect) 與邊界影響而產生誤差，且量測儀器亦不可避免地造成局部流場、溫度場或 濃度場的改變，而影響觀測精度。 數值模擬乃利用數學模式預測物理現象，在有限差分法、有限元素法、 有限體積法與有限解析法等數值方法提出後，使計算流體力學蓬勃發展，與 物理模型試驗相較之下，數值模式計算具有經濟、效率、資料完整與功能強 大等優點，近年來由於電子計算機的高速化與普遍化，電腦記憶體容量與計 算量的限制已大幅降低，利用電腦演算的數值模式，將可克服物理模型試驗 的限制，並可充份發揮數值模擬的預測功能，設計不同的可行方案進行模擬 分析，以提供改善方案參考。 就淺水波運動現象，水體運動受底床邊界限制，其運動方向幾為水平， 垂直方向之運動相對於水平甚小，垂直動量方程式通常簡化為靜水壓分佈， 就所欲分析問題之需要，又進一步將三維之控制方程式依斷面或水深積分簡 化為一維或二維模式。對於河川與海岸等大範圍之流場模擬，由於一維與二 維模式相對於三維模式具有計算快速、模式容易收斂等優點，較適合作長時 間之模擬，目前已廣泛使用於工程實務上；然而，一維模式僅能看出洪水波 沿河道方向之運動現象，平面二維水深積分模式可提供流場與泥砂濃度之平 面分布，但缺乏深度方向資訊，相對地，垂直二維雖可提供垂直面上之流場 與濃度分布，但缺乏側向資訊。許多水利工程常見問題，其運動方向雖近似 水平方向，但其流場及濃度等物理量之垂直分佈卻為影響傳輸行為之重要因
子，例如水庫泥沙運動、湖泊風剪環流場與河口垂直環流場及其鹽度分佈 等，無論是一維或平面二維模式均無法描述其物理現象，擬似三維淺水波運 動模式除能描述其自然界具自由液面之水體運動行為外，並能提供深度方向 之分佈資訊，相對於全三維流體動力模式(fully three-dimensional hydro- dynamic model)，具有計算經濟與收斂快速，適合作為水利工程相關問題之 分析工具。
1.2 文獻回顧
二維水深平均淺水波數值模式因忽略水深流速、濃度等物理量在水深方 向之相關資訊而計算快速，已經被廣泛應用在河川、水庫與海岸中之非定常 自由液面流模擬，但是，當物理量之水深方向變化為欲探討問題或現象之關 鍵資訊時，如河川彎道二次流、水庫泥砂異重流、與河口海水入侵等問題， 水深平均模式並無提供相關資訊以瞭解其物理現象或提供決策之依據。例如 河道中水深方向流速非均勻分佈以及側向二次流造成之延散剪應力影響則 相對重要，考慮與否以及如何處理將影響模擬結果， Molls 與 Chaudhry [25] 結合層流、紊流與延散剪應力的等效應力觀念，忽略二次流垂直方向之速 度。但彎道中水流實為三維流場，Flokstra [8]指出延散剪應力對彎道水流之 影響，Lien 等人[21]進一步指出未處理延散剪應力，在模擬彎道水流時，流 場呈現如勢流(potential flow)中之自由渦流(free vortex)，外岸流速小，內岸 流速大。Lardner 與 Cekirge [17]將水平與垂直流場分開求解，首先利用二維 深度平均模式(depth-averaged model)計算水位分佈與深度平均速度分量，再 透過求解二水平方向之動量方程式，以獲得水平速度分量在垂直方向之分 佈，並用以計算深度平均模式中因垂直方向不均勻速度分佈所形成之延散擬似三維海岸水動力模式之發展，探討河口與海灣之環流場與鹽度分佈 (Wang [34];吳等人[36])。 隨著電腦計算能力快速發展，許多求解雷諾平均(Reynolds averaged) Navier- Stokes 方程式之三維流體動力模式開始應用到天然河川之模擬，對 於河川或水庫之自由液面流動壓模式，大多對自由液面邊界條件進行簡化， 以降低複雜度、提高穩定度。Demuren 與 Rodi(1983)使用鋼性蓋板(rigid-lid) 概念，自由液面可視做對稱平面，液面高程變化變化可由對稱平面上的壓力 呈現，並將此概念應用至蜿蜒河道污染質擴散之三維模擬。Miyata 等人 [24] 在移動網格座標系統上採交錯式網格(staggered grid)分佈，以有限差分法模 擬船舶前進時之周圍流場，並假設自由液面之黏滯性很小可忽略，將自由液 面上之正向剪力邊界條件簡化為大氣壓力自由液面邊界條件，並刪除切線剪 力邊界條件，求解層流狀態下之自由液面流模式，為解決自由液面在計算過 程中產生不合物理現象之擾動，接引用自由液面平滑化技巧，使計算獲得收 斂解。Neary [27]發展完全三維分流模式，假設水面變化不大採固定水面網 格，考慮紊流流場特性，探討取水工流場分布特性。Gessler 等人[9]發展三 維動床模式探討彎道流場特性與底床變動現象，考慮紊流效應與懸浮載分 布，並探討底床變動與河床沖蝕、懸浮載沉降、河床載傳輸等機制。Wu 等 人[35]以 FAST3D 為基礎，將有限體積法三維動壓模式應用至河道水流與河 床沖淤模擬上。
Hervouet 與 Van Haren [13]提出移動邊界之三維 Navier- Stokes 模式，內 建二維靜水壓假設之淺水波方程，其計算方式係先由二維淺水波方程求解水 面高程，依據所得之水面高程利用移動邊界法建立三維網格來求解三維流 場，並依據三維計算所得水面邊界之壓力分布，再次修正水面高程。Meselhe 與 Sotiropoulos [23]同樣以雷諾平均 Navier- Stokes 模式，應用同一概念但容
許自由液面層之網格變形，讓模擬結果網格變形調整下達到收斂解，因此毋 需每一迭代均重建一次網格，較 Hervouet 與 Van Haren [13]有效降低計算 量，但對於非定常流之模擬，對每一時間步驟Δt，其網格仍需針對水深變化 而重建。Ahsan 與 Blumberg [1]發展三維水理模式探討湖泊內水質問題與溫 度分布，假設水庫內靜水壓分布，由水平方向動量方程式求解水平方向流 場，代入連續方程式得到垂直方向速度，並考慮紊流擴散之效應。此外，由 於目前大部份靜壓三維模式對於渦流滯性係數常引用經驗剖面假設(Davis & Rodi [4]; Heap & Jone [12])，且因這些經驗剖面假設並不具有通用性，故選 用不同的經驗剖面假設，將產生不同之計算結果。 Osher等人[29]提出的等位函數(Level set)法通過設置一個距離函數隱式 跟蹤自由液面，並且將介面條件融入控制方程式，可在整個計算域內用統一 的差分式來解，因此整個計算的過程變的簡單許多，而且變的具有較大的通 用性；由於該方法不涉及座標轉換，可使用的差分形式較多，因此容易使用 較高精度的差分式，已被成功地用於解決一些二相流的問題，其基本概念 是，將兩種界質用統一的N-S方程求解，在界面兩側分別採用兩相流各自的 密度、粘性，並在自由液面上給與一適當的函數，定義φ為空間中某點到分 界面有方向的垂直距離，目前該方法已應用於鋒面運動(Osher與Sethian [29])、渦漩運動模擬(Harabetian et al. [10])以及不可壓縮二相流計算(Sussman
et al. [32])等。Hirt及Nichols [14]提出流體體積法(Volume of Fluid, VOF)，將 液體在一個網格中佔有之體積比率當作指標，以訂定流場中氣、液分離相之 介面，計算體積比率之改變即可追蹤該分離界面之運動，即利用一fractional volume來標記某一計算網格為水(φ＝1)或空氣(φ＝0)或為自由液面(0＜φ＜ 1)，由網格內流體體積函數值，可決定網格內之液面斜率，進而可決定液面 形狀，流體體積法可用於不同相位接觸界面流場，且對碎波亦具處理能力， 因此廣泛應用在不同領等位函數法模擬三維自由液面流域。已廣泛應用於化
模組化之概念並預留模組供使用者延伸定義對實務應用保有相當彈性， Faure等人[7]將入流、出流等邊界處理嵌入CFX-4用以模擬河川之自由液面 流流場，提升全三維流體動力模式之實用性。 二維模式在計算時間上相對於三維模式節省許多，較適合作長時間之模 擬，但由於二維模式須引入二次流速度分布，然引入之速度分部是否符合實 際無狀況並無法得知，三維模式能計算河道任意點之流場，但其所花費之計 算時間（CPU time）甚鉅，且固定網格無法反映水面之變化，分層(multi-level or multi-layer)模式則是對擬似三維模式分層積分推導而得，若層數夠多則可 得與擬似三維模式相近之精度，且亦可減少層數適度降低擬似三維模式之計 算量並避免尺度效應，正好可達成二維模式與三維模式無法完成之部分，而 且，由於分層模式在水深方向使用層厚度(或水深分割)之概念取代三維網 格，容許自由液面隨時間或迭代變動，因此沒有每次迭代均須重新產生網格 的問題，有效降低模式計算量，已經成功應用在大範圍之海岸與河口問題。 Choi [3]採消散 Galerkin 有限元素法發展分層積分泥石流模式，採用移動與 固定格網並行方式探討泥石流進入水庫之前進速度與淤積現象；Li 與 Zhan [20]提出分層有限元素模式並用於東京灣灣流之模擬，Kim 與 Lee [16]則應 用基於多層概念之有限差分模式 BACHOM-3 模擬韓國 Suyoung 灣，Li 與 Gu [19]使用分層積分之概念探討港灣內之潮汐沖刷與底床汙染物交換之問 題，Shankar 等人[31]使用相同概念發展河道之有限差分水理模式應用於河道 突縮、突擴與水門下游之模擬。 分層模式在河口或海岸問題在河口或海岸問題之模擬已經相當成熟，但 因河口或海岸問題中，水位變動固定於海平面附近，其變動範圍遠小於水 深，並無流量－水位是否合理之困擾，若應用至應用至河川或水庫問題時， 自由液面位置之決定仍未解決。Lin 與 Falconer [22]則將分層概念應用至英
格蘭東北海岸 Hummer 河口感潮帶中乾濕交替區域之模擬，但因模式缺少流 量－水深平衡機制，仍需使用二維水深平均子模式來決定自由液面之位置。 顏[38]發展一假設靜水壓分佈之擬似三維曲線座標分層紊流模式，能模擬湧 浪、潰壩等自由液面流，並應用於丁壩周圍流場之模擬。
1.3 研究方法
應用分層概念求解三維模式係於水深方向使用層厚度(或水深分割)之 概念取代三維網格，容許自由液面隨時間或迭代變動，因此沒有每次迭代均 須重新產生網格的問題，有效降低模式計算量，此概念應用於河口或海岸等 淺水波問題之模擬已經相當成熟，但在河川與水庫等問題存在水位－流量之 關係，水深之決定仍須仰賴二維子模式等，本研究擬發展一擬似三維分層積 分淺水波模式，毋須藉助於二維子模式決定水深，能應用於河川或水庫之模 擬，降低完全三維模式對自由液面處理之困難及耗費大量計算時間，並改善 深度平均二維模式無法得知深度方向之流場分布，期能對河川與水庫內流場 作較精確之分析。 蔡[37]將分層積分概念應用於區域地下水與地層下陷模擬，並引入二次 形狀函數確保層介面間之地下水水頭及層介面間地下水通量連續，本研究延 伸此概念應用至自由液面流之計算，透過二次形狀函數與層介面速度、剪力 連續之概念，發展經濟、效率之分層積分數值模式，模式中流體在水深方向 將劃分為若干層，將每一層之控制方程式沿垂直方向積分，可得二維層厚度 平均控制方程式；然後應用二次形狀函數於層內部之速度分佈，基於層介面力分布，進一步藉由層介面剪力連續之條件，推得層介面速度限制式，將各 層積分控制方程式與各層介面速度限制式連結可得水柱方程組。模式建立在 卡氏座標系統上之結構網格(structure grid)上，於交錯網格(staggered grid)上 布置模式之變數，使用具守衡特性、已廣泛應用於計算流體力學之有限體積 法(Finite Volume Method)離散控制方程式。首先在垂直方向利用隱式法求解 各水柱方程組可得各水柱垂直方向之速度剖面，然後逐次求解各水柱方程 組；接著，求解層連續方程式，透過層連續方程式可計算得各層界面間之通 量及層界面上之垂直速度，由底層開始逐次求解至自由液面層，可得水位變 動量；並透過迭代方式以求得每一時間步驟之收斂解。 本研究透過與 Heaps [12]針對定常層流底床不滑移狀況，推導出封閉邊 界水體如湖泊或水庫表面受風吹引起之垂直環流場沿水深方之流速剖面解 析解之比較，並檢視流速剖面之發展過程是否能反映，水面受到風剪應力帶 動，表面水體朝與風同方加速，進而帶動更內部之水體，並在表面水體會產 生堆升，進而造成底部水平方向壓力不平衡而造成底部水體朝與風反向運動 之物理現象，來驗證模式是否具模擬自由液面流場之機制；此外，進一步透 過假設之直線渠道案例驗證文中推導之紊流黏滯係數估算式，並計算比較不 同黏滯係數分佈型態、液面與底床邊界等條件等對流速剖面與水深之影響， 展現模式模擬明渠流之能力。
1.4 章節介紹
前述以闡明本研究之背景與目的、文獻回顧、研究方法，茲扼要說明本 文章節之內容。第一章 緒論：闡述本研究之背景與目的，並回顧三維自由液面流之發展過 程、計算方法與應用範圍，及其遭遇問題之相關文獻，提出本文之 研究方法與解決方式。 第二章 分層積分淺水波模式理論架構：由三維 Navier-Sokes 方程式基於靜 水壓假設，簡化為擬似三維模式，經分層積分得層積分控制方程 式，並簡述邊界條件之處理與推導紊流黏滯係數關係式。 第三章 分層模式(一)－不考慮層介面連續：為本研究所嘗試但未成功之分 層模式，簡述其數值架構與計算流程，並檢討該分層模式所遭遇之 問題。 第四章 分層模式(二)－考慮層介面連續：經由分層模式(一)之經驗，改進 分層模式數值架構，為本研究目前使用分層模式，簡述其數值架 構、引進二次形狀函數以推導層介面剪力連續條件限制式與模式計 算流程及收斂條件。 第五章 風剪垂直環流場模式驗證：透過風剪環流場之解析解比較，驗證模 式之機制與精度，並透過不同紊流黏滯係數分佈測試模式對不同流 速剖面之計算能力。 第六章 明渠流之計算與比較：透過假設渠道案例之計算與比較，驗證本研 究提出之紊流黏滯係數推估式，並針對不同底床粗糙度、黏滯係數 垂直分佈型態、底床滑移邊界及液面風剪作用等，測試模式之特性。 第七章 結論與建議：簡述本研究突破現有三維自由液面流發展所遭遇之難 題，與後續研究所需改進之建議及研究方向。
第二章 分層積分淺水波模式理論架構
自由液面為水與空氣兩相之交界，自然界之水流如河川、水庫、海岸等 均為自由液面流，與壓力流最大差異在於壓力可反應於水深，本章由三維 Navier- Sokes 方程式基於靜水壓假設，簡化為擬似三維淺水波模式，經分層 積分得層積分控制方程式，並簡述邊界條件之處理與推導紊流黏滯係數推估 式。2.1 雷諾平均控制方程式
三維黏性流體之 Navier-Stokes 方程如后： 連續方程 0 ) u ( x t i i = ρ ∂ ∂ + ∂ ρ ∂ (2.1) 動量方程 i j i j i i j j i g x u x x p ) u u ( x ) u ( t ⎟⎟⎠+ρ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ μ ∂ ∂ + ∂ ∂ − = ρ ∂ ∂ + ρ ∂ ∂ (2.2) 其中，ui為流體之速度，μ為流體之動力黏滯係數，g 為重力加速度，p 為壓 力。 自然界中之自由液面流，多為紊流，可由雷諾平均方程得到雷諾剪應力，若”－”代表時間平均，則瞬時流速 ui為時間平均流速ui與擾動量 ui’之 組合，可表示如下 ' u u ui = i + i ，p=p+p' 其中，時間平均之流速ui與壓力p定義為
∫
Δ + Δ = t t t i i udt t 1 u ，∫
Δ + Δ = t t t pdt t 1 p 則物理量具有下列時間平均關係 0 ' ui = ，uj'ui =0，ui =ui 自然界流動之水體可視為不可壓縮流，因此不考慮密度之擾動量，基於 不可壓縮流之假設，Navier-Stokes 方程經由雷諾時間平均，可得雷諾方程式 如下 連續方程 0 x u i i = ∂ ∂ (2.3) 動量方程 i j i j i j i i j j i g ' u ' u x u x x p 1 ) u u ( x t u + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − ∂ ∂ ν ∂ ∂ + ∂ ∂ ρ − = ∂ ∂ + ∂ ∂ (2.4) 其中，ui'uj'為時間平均過程之流速擾動相關量，乘上密度則為擾動引起之 動量傳輸，即為雷諾剪應力(Reynolds stresses)。分佈來簡化，而垂直方向之速度可透過連續方程式計算以減少三維模式之計 算量，即擬似三維模式(semi- three dimensional model)，因此基於靜水壓假設 之不可壓縮三維流體運動控制方程式可改寫如下： 0 z w y v x u = ∂ ∂ + ∂ ∂ + ∂ ∂ (2.5) z y x x p 1 z ) wu ( y ) vu ( x ) uu ( t u xx yx zx ∂ τ ∂ + ∂ τ ∂ + ∂ σ ∂ + ∂ ∂ ρ − = ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ (2.6) z y x y p 1 z ) v w ( y ) v v ( x ) v u ( t v xy yy zy ∂ τ ∂ + ∂ σ ∂ + ∂ τ ∂ + ∂ ∂ ρ − = ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ (2.7) ) z H ( g p=ρ +ηB− (2.8) 其中，t為時間，x、y、z為卡氏(Cartesian)座標，u、v、w分別為x、y、z方 向上之速度分量，p為壓力，ρ為水體密度，g為重力加速度， σxx、σyy、τyx、 τxy、τzx、τzy、分別為應力張量於x-y、y-z、z-x平面上之x、y方向之分量，H 為水深，ηB為底床高程。
2.2 分層積分控制方程式
分層積分模式之概念建立在擬似三維模式上，其網格建立方式在水平 方向與水深平均二維模式相同將模擬範圍內之水體分為單一水柱，可用卡氏 座標、正交曲線座標或非正交曲線座標，而單一水柱之水深方向上進一步簡 化分為數層，將每一層內之擬似三維模式控制式沿層厚度積分，得到數組二 維層積分控制方程式，分層方式之概念與變數配置示如圖2.1，其中包含液 面層(top layer)、內部層(internal layer)、底床層(bottom layer)，基於不可壓縮 流與靜水壓之假設，將擬似三維流體運動控制方程式(2.5)~(2.7)對層厚度積 分，可得分層積分控制方程式。首先定義層平均物理量如下：
∫
+ η η + = φ φ ) y , x ( ) y , x ( 2 / 1 k 1 k k dz ) t , z , y , x ( h 1 (2.9) 其中，φ代表控制方程式之物理量如流速 u、v 與壓力 p，φk+1/2代表第 k 層 之積分平均物理量，ηk+1為介於第(k+1)層與第 k 相鄰之層介面高程，ηk為 介於第(k-1)層與第 k 相鄰之層介面高程，K為總層數，h為層厚度。 動量方程式之對流量為非線性，其積分關係如下：∫
∫
∫
∫
∫
∫
∫
∫
∫
+ + + + + + + + + η η + + η η + + η η + + + + η η + + η η + + η η + + η η + + η η + + + + η η φ φ φ φ + φ φ φ + φ φ φ + φ φ = φ φ φ φ + φ φ φ + φ φ φ + φ φ = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ φ − φ + φ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ φ − φ + φ = φ φ ) y , x ( ) y , x ( 2 1 k , j j 2 1 k , i i ) y , x ( ) y , x ( 2 1 k , i i 2 1 k , j ) y , x ( ) y , x ( 2 1 k , j j 2 1 k , i 2 1 k , j 2 1 k , i ) y , x ( ) y , x ( 2 1 k , j j 2 1 k , i i ) y , x ( ) y , x ( 2 1 k , i i 2 1 k , j ) y , x ( ) y , x ( 2 1 k , j j 2 1 k , i ) y , x ( ) y , x ( 2 1 k , j 2 1 k , i ) y , x ( ) y , x ( 2 1 k , j j 2 1 k , j 2 1 k , i i 2 1 k , i ) y , x ( ) y , x ( j i 1 k k 1 k k 1 k k 1 k k 1 k k 1 k k 1 k k 1 k k 1 k k dz ) -)( -( h 1 dz ) -( h dz ) -( h dz ) -)( -( h 1 dz ) -( h 1 dz ) -( h 1 dz h 1 dz ) ( ) ( h 1 dz h 1 (2.10) 由於 ( - )dz 0 ) y , x ( ) y , x ( 2 1 k , i i 1 k k = φ φ∫
+ η η + ，(2.10)式重新整理如後∫
∫
+ + η η + + + + η η φ φ φ φ + φ φ = φ φ ) y , x ( ) y , x ( 2 1 k , j j 2 1 k , i i 2 1 k , j 2 1 k , i ) y , x ( ) y , x ( j i 1 k k 1 k k dz ) -)( -( h 1 dz h 1 (2.11) 第 k 層動量方程式沿層厚度積分可得層積分動量方程式k zx 1 k zx yx xx 2 / 1 k B k 1 k 2 / 1 k 2 / 1 k 2 / 1 k 2 / 1 k 2 / 1 k ) ( ) ( ) T h ( y ) T h ( x ] x ) z H ( gh [ ) Fu ( ) Fu ( ) u hv ( y ) u hu ( x ) hu ( t ρ τ − ρ τ + ρ ∂ ∂ + ρ ∂ ∂ + ∂ − η + ∂ − = − + ∂ ∂ + ∂ ∂ + ∂ ∂ + + + + + + + + (2.12) k zy 1 k zy yy xy 2 / 1 k B k 1 k 2 / 1 k 2 / 1 k 2 / 1 k / 1 k 2 / 1 k ) ( ) ( ) T h ( y ) T h ( x ] y ) z H ( gh [ ) Fv ( ) Fv ( ) v hv ( y ) v hu ( x ) hv ( t ρ τ − ρ τ + ρ ∂ ∂ + ρ ∂ ∂ + ∂ − η + ∂ − = − + ∂ ∂ + ∂ ∂ + ∂ ∂ + + + + + + + + (2.13) 其中，下標(k+1/2)表示層平均物理量，下標(k+1)與 k 分別表示層之上、下 層介面物理量，h 為層厚度，H 為水深，η B為底床高程，F 為通過層介面之 通量，則第(k-1)層與第 k 層間之層介面之通量 Fk 可表示為 y v x u w F k k k k k k ∂ η ∂ − ∂ η ∂ − = (2.14) 其中，wk為介於第(k-1)層與第 k 相鄰之層介面垂直速度。 同理，對第 k 層連續方層式沿層厚度積分可得層積分連續方程式 0 F F ) hv ( y ) hu ( x k 1/2 ∂ k 1/2 + k 1 − k = ∂ + ∂ ∂ + + + (2.15) (2.14)式滿足底床不透水邊界條件，可得底床速度 y v x u w B 1 B 1 1 ∂ η ∂ + ∂ η ∂ = (2.16) 其中，zB為底床高程。則第 k 層與第(k+1)層間之層介面垂直速度可由層積 分連續方程式(2.15)由底床層逐次向上計算，表示如下：
∑
= + + + + + + + ∂ ∂ + ∂ ∂ − ∂ η ∂ + ∂ η ∂ = k 1 m 2 / 1 m 2 / 1 m 1 k 1 k 1 k 1 k 1 k (hv )] y ) hu ( x [ y v x u w (2.17) 其中，m 為啞指標。而自由液面之水位變動量可由第 K 層（液面層）之頂 部垂直速度對時間積分而得，即ΔH = wK+1Δt，或表示如下0 )] hv ( y ) hu ( x [ t H H K 1 m 2 / 1 m 2 / 1 m t t t = ∂ ∂ + ∂ ∂ + Δ −
∑
= + + Δ + (2.18) 其中，Δt為時間間距，K為總層數。 方程式(2.12)與(2.13)中，Txx、Tyx、Txy、Tyy為層流剪應力、紊流剪應力 與延散剪應力和張量之各分量，可表示如下：∫
+ + − ρ − ρ − ∂ ∂ μ = k1 k z z 2 1/2 k 2 xx u' (u u ) ]dz x u [ h 1 T (2.19)∫
+ + + − − ρ − ρ − ∂ ∂ μ = k1 k z z k 1/2 k 1/2 yx u'v' (u u )(v v )]dz y u [ h 1 T (2.20)∫
+ + + − − ρ − ρ − ∂ ∂ μ = k1 k z z k 1/2 k 1/2 xy u'v' (u u )(v v )]dz x v [ h 1 T (2.21)∫
+ + − ρ − ρ − ∂ ∂ μ = k1 k z z 2 1/2 k 2 yy v' (v v ) ]dz y v [ h 1 T (2.22) 根據Boussinesq紊流渦動黏性觀念 t l +τ τ = τ j i l l x u ∂ ∂ μ = τ j i t j i t x u u u ∂ ∂ μ = ρ − = τ 其中，τ為總剪應力，τl為層流剪應力，τt為紊流剪應力；μl為層流動力黏滯 係數，μt為紊流動力黏滯係數。若忽略延散剪應力，則應力張量各分量可簡 化如下： x u T k 1/2 h xx ∂ ∂ ν = ρ + y u T k 1/2 h yx ∂ ∂ ν = ρ +x v T k 1/2 h xy ∂ ∂ ν = ρ + y v T k 1/2 h yy ∂ ∂ ν = ρ + (2.23) 其中， t h l h =ν +ν ν 為水平層流與紊流運動黏滯係數之和。 將(2.23)式代入方程式(2.12)與(2.13)，則層積分動量方程式可簡化為 k zx 1 k zx 2 / 1 k h 2 / 1 k h 2 / 1 k B k 1 k 2 / 1 k 2 / 1 k 2 / 1 k 2 / 1 k 2 / 1 k ) ( ) ( ) y u h ( y ) x u h ( x ] x ) z H ( gh [ ) Fu ( ) Fu ( ) u hv ( y ) u hu ( x ) hu ( t ρ τ − ρ τ + ∂ ∂ ν ∂ ∂ + ∂ ∂ ν ∂ ∂ + ∂ − η + ∂ − = − + ∂ ∂ + ∂ ∂ + ∂ ∂ + + + + + + + + + + (2.24) k zy 1 k zy 2 / 1 k h 2 / 1 k h 2 / 1 k B k 1 k 2 / 1 k 2 / 1 k 2 / 1 k / 1 k 2 / 1 k ) ( ) ( ) y v h ( y ) x v h ( x ] y ) z H ( gh [ ) Fv ( ) Fv ( ) v hv ( y ) v hu ( x ) hv ( t ρ τ − ρ τ + ∂ ∂ ν ∂ ∂ + ∂ ∂ ν ∂ ∂ + ∂ − η + ∂ − = − + ∂ ∂ + ∂ ∂ + ∂ ∂ + + + + + + + + + + (2.25) 茲將第 k 層之積分控制方程式重新整理如後 層積分動量方程式 k zx 1 k zx 2 / 1 k h 2 / 1 k h 2 / 1 k B k 1 k 2 / 1 k 2 / 1 k 2 / 1 k 2 / 1 k 2 / 1 k ) ( ) ( ) y u h ( y ) x u h ( x ] x ) z H ( gh [ ) Fu ( ) Fu ( ) u hv ( y ) u hu ( x ) hu ( t ρ τ − ρ τ + ∂ ∂ ν ∂ ∂ + ∂ ∂ ν ∂ ∂ + ∂ − η + ∂ − = − + ∂ ∂ + ∂ ∂ + ∂ ∂ + + + + + + + + + + (2.24) k zy 1 k zy 2 / 1 k h 2 / 1 k h 2 / 1 k B k 1 k 2 / 1 k 2 / 1 k 2 / 1 k / 1 k 2 / 1 k ) ( ) ( ) y v h ( y ) x v h ( x ] y ) z H ( gh [ ) Fv ( ) Fv ( ) v hv ( y ) v hu ( x ) hv ( t ρ τ − ρ τ + ∂ ∂ ν ∂ ∂ + ∂ ∂ ν ∂ ∂ + ∂ − η + ∂ − = − + ∂ ∂ + ∂ ∂ + ∂ ∂ + + + + + + + + + + (2.25) 層積分連續方程式 0 F F ) hu ( y ) hu ( x k 1/2 ∂ k 1/2 + k 1− k = ∂ + ∂ ∂ + + + (2.15)
水位變動方程式 0 )] hv ( y ) hu ( x [ t H H K 1 m 2 / 1 m 2 / 1 m t t t = ∂ ∂ + ∂ ∂ + Δ −
∑
= + + Δ + (2.18)2.3 邊界條件
茲將模式所需之邊界條件及其處理方法敘述如後： z 固體邊界條件：底床邊界條件採不滑移邊界條件，即 u 0 B z=η = 、 0 v B z=η = ，與不透水邊界條件wz=ηB =0；側壁邊界條件採不透水邊界條 件vy=0 =0、vy=B=0，與滑移邊界條件∂u/∂yy=0 =0、∂u/∂yy=B =0；其中， B為渠道寬度。 z 自由液面邊界條件：自由液面運動邊界條件與風剪力邊界條件。 z 上游邊界條件：採單寬流量邊界條件，縱向流速於水深方向均勻分布， 即 t q/H u= ，q為單寬流量，Ht為任一計算時刻任一迭代之即時水深，無 側向與垂向流速，即ux=0 =0、vx=0 =0。 z 下游邊界條件：假設均勻流條件，即∂u/∂x=0、∂v/∂x=0。2.4 明渠流之垂直紊流黏滯係數推估式
在河川一維或二維水理模式中，水深之決定可將水柱視為自由體(free三者間相互作用平衡，透過計算底床剪應力之摩擦公式即可估計水深。但在 三維流場模擬中，底床剪應力僅作用於緊鄰底床邊界之最底層，而水體內部 層與層之間剪應力須由紊流黏滯係數與流速梯度決定，然紊流黏滯係數係流 場之函數而非流體之性質，其分佈型態係由流場決定，因而流速剖面與黏滯 係數分佈型態係交互作用影響；又，任意流量下基於連續方程式，流速剖面 又與水深相互影響，由此可知水深則由流量、流速剖面與紊流黏滯係數分佈 三者之間相互作用平衡決定，流量或物理邊界幾何改變都將影響紊流黏滯係 數之分佈，因此，在實際模擬中紊流黏滯係數並由無法透過案例率定而得。 在紊流模式尚未建立前，必需建立紊流黏滯係數與流量、水深、底床粗 糙度之關係，方可模擬明渠流之流場；Elder [6]透過實驗觀察提出之黏滯係 數關係式ν=αHu*，在實務上已廣泛被使用在污染質模擬，透過Elder紊流黏 滯係數關係式可建立紊流黏滯係數與水深、流量、底床粗糙度間相互影響之 機制。 本文發展之模式尚未考慮紊流模式，因此如前所述需建立紊流黏滯係 數關係式，於此處將推導過程敘述如後。首先，考慮三維流體運動之Navier- Stokes方程式 x 2 2 2 2 2 2 g ) z u y u x u ( 1 x P 1 z u w y u v x u u t u + ∂ ∂ + ∂ ∂ + ∂ ∂ μ + ∂ ∂ ρ − = ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ (2.26) 假設靜水壓分佈，且水流為定常均勻流，即∂/∂t=0、∂/∂x=0，假設水流在側 向上均勻分佈，即∂/∂y=0，且底床坡度不大，方程式可簡化為 o 2 2 v gS z u − = ∂ ∂ ν (2.27) 其中So為底床坡度。假設紊流運動黏滯係數沿水深方向均勻分佈，方程式 (2.27)對水深積分可得 1 v o C z gS z u + ν − = ∂ ∂ (2.28)
方程式滿足底床剪應力邊界條件 v o 0 z z u ρν τ = ∂ ∂ = ，可得係數C1 = τo/ρ，則水平速 度u於垂直方向之梯度可表示為 v o v o z gS z u ρν τ + ν − = ∂ ∂ (2.29) 若水面受風剪應力τs作用，則水面流速梯度為 v s H z z u ρν τ = ∂ ∂ = 將水面梯度代入方程式可得下列關係 v o v o v s H gS ρν τ + ν − = ρν τ 或表示成 ρ τ = ρ τ + − o s oH gS 由於水流為均勻流況，因此底床剪應力可由水深、底床坡度、水面風剪應力 三者決定 o s o =τ +γHS τ (2.30) 再對方程式(2.29)沿水深方向積分，可得 2 v o 2 v o z z C 2 gS u + μ τ + ν − = 滿足底床不滑移邊界條件，則C2 = 0，因此，上式簡化為 z z 2 gS u v o 2 v o μ τ + ν − = (2.31) 方程式(2.31)即為流速剖面分佈函數。 基於質量守衡，流速剖面方程式(2.31)沿水深方向積分即為單寬流量q， 即
將方程式(2.31)代入方程式(2.32)，可得下列關係式 q H 2 H 3 gS 2 v s 3 v o = μ τ + ν (2.33) 由於水深H、紊流運動黏滯係數νv均為未知數，因此無法單由方程式(2.33) 決定紊流運動黏滯係數νv之值。若忽略液面風剪應力之作用，則紊流運動黏 滯係數νv可表示為 3 o v H q 3 gS = ν (2.34) 或Elder紊流黏滯係數關係式之形式νv =αHu*，則係數α為 U 3 u* = α 假設渠道寬深比夠大，水力半徑可以水深取代，渠道水深平均流速可由曼寧 公式計算如下 * 6 / 1 2 / 1 o 3 / 2 u g n H S H n 1 U = = (2.35) 則Elder紊流黏滯係數關係式之係數α可改寫為 6 / 1 * H 3 g n U 3u = = α (2.36) 於是，Elder紊流黏滯係數關係式之係數α可由曼寧粗糙度係數n、水深H、底 床剪力速度u*決定。回到方程式(2.36)，若自由液面受到風剪應力作用於上， 則底床剪應力會隨之增加。另，對於非均勻流況，方程式(2.35)中之底床坡 降 So 可由摩擦坡降 Sf 取代。 因此，在紊流模式尚未建立前，紊流黏滯係數可藉由Elder關係式推得， 其中係數 α 則可透過方程式(2.36)由水深、流量、曼寧粗糙係數 n 計算得 到；現場常見之流況，如當曼寧 n 介於 0.02 ~ 0.05，水深介於 0.5m ~ 10m，
由方程式(2.36)計算之係數 α 約介於 0.014 ~ 0.058，該範圍略小於 Nezu與 Rodi [26] 實驗量測範圍(水深平均值約0.06左右)，研判應為假設紊流黏滯係 數均勻分佈，致流速梯度於底床與水面附近過小而半水深處流速梯度過大之 影響。
2 / 1 k u+ uk uk+1 k zx τ 1 k zx + τ wk+1 wk Layer K Layer k Layer 1 H 圖 2.1 分層模式層切割方式與變數配置示意圖
第三章 分層模式(一)－不考慮層介面連續
3.1 有限體積法與數值差分式
在不可壓縮流與靜水壓假設下，非定常流擬似三維分層積分淺水波模式 之控制方程式已詳述於第二章，為方便閱讀，茲將本節相關之動量方程式、 連續方程式、水位變動方程式重述如下： 層積分動量方程式 k zx 1 k zx 2 / 1 k h 2 / 1 k h 2 / 1 k B k 1 k 2 / 1 k 2 / 1 k 2 / 1 k 2 / 1 k 2 / 1 k ) ( ) ( ) y u h ( y ) x u h ( x ] x ) z H ( gh [ ) Fu ( ) Fu ( ) u hv ( y ) u hu ( x ) hu ( t ρ τ − ρ τ + ∂ ∂ ν ∂ ∂ + ∂ ∂ ν ∂ ∂ + ∂ − η + ∂ − = − + ∂ ∂ + ∂ ∂ + ∂ ∂ + + + + + + + + + + (2.24) k zy 1 k zy 2 / 1 k h 2 / 1 k h 2 / 1 k B k 1 k 2 / 1 k 2 / 1 k 2 / 1 k / 1 k 2 / 1 k ) ( ) ( ) y v h ( y ) x v h ( x ] y ) z H ( gh [ ) Fv ( ) Fv ( ) v hv ( y ) v hu ( x ) hv ( t ρ τ − ρ τ + ∂ ∂ ν ∂ ∂ + ∂ ∂ ν ∂ ∂ + ∂ − η + ∂ − = − + ∂ ∂ + ∂ ∂ + ∂ ∂ + + + + + + + + + + (2.25) 層積分連續方程式 0 F F ) hv ( y ) hu ( x k 1/2 ∂ k 1/2 + k 1 − k = ∂ + ∂ ∂ + + + (2.15) 水位變動方程式 0 )] hv ( y ) hu ( x [ t H H K 1 m 2 / 1 m 2 / 1 m t t t = ∂ ∂ + ∂ ∂ + Δ −∑
= + + Δ + (2.18)對於非定常流不壓縮紊流，分層模式之動量方程式，可以下列通式表示 φ + + − + + + + + + + ∂ φ ∂ ν + ∂ φ ∂ ν ∂ ∂ + ∂ φ ∂ ν ∂ ∂ = φ ρ + φ ρ ∂ ∂ + φ ρ ∂ ∂ + φ ρ ∂ ∂ S ) z ( ) y h ( y ) x h ( x ) F ( ) v h ( y ) u h ( x ) h ( t b -t v 1/2 k h 1/2 k h b t 1/2 k 1/2 k 1/2 k 1/2 k 1/2 k 1/2 k (3.1) 其中，φ為 u 或 v，Sφ為水壓項。由分層動量方程通式對對控制體進行體積積 分，可得分層模式有限體積方程式，即 φ − − + − + − + − + + Δ + + + Δ Δ ∂ φ ∂ ν + Δ ∂ φ ∂ ν + Δ ∂ φ ∂ ν = Δ Δ ρ + Δ φ ρ + Δ φ ρ + Δ Δ Δ φ ρ − Δ Δ φ ρ S ) y x z ( ) x y h ( ) y x h ( ) y x F ( ) x F h ( ) y F h ( t / ] ) y x h ( ) y x h [( b t v s n 1/2 k h w e 1/2 k h b -t s n 1/2 k y w e 1/2 k x t 1/2 k t t 1/2 k (3.2) 其中，φ為 u 或 v，Fx、Fy為控制面上的通量，下標 (e,w) 與 (n,s) 表示控制 體積於 (東,西) 與 (南,北) 控制面上之動量通量，e、w、n、s 控制面位置 如圖 3.1 所示。相對應於 u、v 之水壓項 Sφ分別表示如下 w -e 2 / 1 k B u ] ) z H ( y gh -[ S = Δ +η − + (3.3) s -n 2 / 1 k B v ] ) z H ( x gh -[ S = Δ +η − + (3.4) 同理，層積分連續方程式為 0 ) y x F ( ) x v h ( ) y u h ( t ) y x h ( ) y x h ( b -t R s -n 1/2 k w -e 1/2 k 0 n = Δ Δ ρ + Δ ρ + Δ ρ + Δ Δ Δ ρ − Δ Δ ρ + + (3.5) 而水位方程式則可表為
0 ] ) x v h ( ) y u h [( t ) y x H ( ) y x H ( m -s n 1/2 k K 1 m w -e 1/2 k t t t = Δ ρ + Δ ρ + Δ Δ Δ ρ − Δ Δ ρ + = + Δ +
∑
(3.6) 將(3.2)式內擴散項中位於控制面上之梯度值，利用中央差分法近似，以 e 控制面上為例: ) ( x 1 ) y h ( ) y x h ( P , 2 1 k E , 2 1 k e h e 1/2 k h + + + φ −φ Δ Δ ν ≅ Δ ∂ φ ∂ ν (3.7) ) ( y 4 1 ) x h ( ) x y h ( S , 2 1 k N , 2 1 k SE , 2 1 k NE , 2 1 k e h e 1/2 k h + + + + + φ −φ +φ −φ Δ Δ ν ≅ Δ ∂ φ ∂ ν (3.8) 並將對流項中位於控制面上之物理變數，以 QUICK(Quadratic Upstream Intepolation for Convective Kinetics, [18])法推求，QUICK 法可表示如下當Gx,e >0 ] 2 3 [ 8 1 W , 2 1 k P , 2 1 k E , 2 1 k P , 2 1 k e , 2 1 k+ =φ + + φ + − φ + −φ + φ (3.9) 當Gx,e <0 ] 2 3 [ 8 1 EE , 2 1 k E , 2 1 k P , 2 1 k P , 2 1 k e , 2 1 k+ =φ + + φ + − φ + −φ + φ (3.10) 其中 Gx為 x 方向上控制面上之通量，其他控制面與層界面可類比(3.7)式至 (3.10)式，將之代入(3.2)式後，再減去 n P , 2 1 k 4 3 + φ 乘連續方程式，經整理即可化 為有限體積差分通式 φ + + + + + + + + φ + φ + φ + φ + φ + φ = φ S A A A A A A A B , 2 1 k B T , 2 1 k T S , 2 1 k S N , 2 1 k N W , 2 1 k W E , 2 1 k E P , 2 1 k P (3.11)
3.2 數值演算方法
3.1 節將動量方程式轉為有限體積積分通式，φ以 u 或 v 代入，即可得到 x、y 方向動量方程式。∑
= + + =nb E,W,N,S,T,B + u nb , 2 1 k u nb P , 2 1 k u Pu A u S A (3.12)∑
= + + =nb E,W,N,S,T,B + v nb , 2 1 k v nb P , 2 1 k v Pv A v S A (3.13) 由(3.5)、(3.6)、(3.12)與(3.13)式形成一代數方程組，用以求解流場。本研究採 Patankar 與 Spalding [30]提出壓力連結擬似隱式法(Semi- Implicit Method for Pressure-Linked Equations)求解，此法已被廣泛應用至各 類不可壓縮流場與明渠流之計算，其包含一個預測步驟與一個校正步驟。因 此對於正確流場 u、v、h、H 之分佈，SIMPLE 引入修正量(u’、v’、h’、H’)， 用以修正猜值(u*、v*、h*、H*)， u=u*+u' (3.14) v =v*+v' (3.15) h=h*+h' (3.16) H=H*+H' (3.17) 首先，將方程式(3.3)、(3.4)代入(3.12)、(3.13)，並以起始猜值(u*、v*、h*、 H*)取代(u、v、h、H)，則可得到對應 H*之動量方程式
[
B k 1/2 e-w]
u P B , T , S , N , W , E nb 2 ,nb 1 k u nb u P P , 2 1 k A h*[(H* z) ] y g -* u A A 1 * u + = + + +η − Δ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ =∑
(3.18)[
B k 1/2 n-s]
v P B , T , S , N , W , E nb 2 ,nb 1 k v nb v P P , 2 1 k ] ) z * H ( [ * h A x g -* v A A 1 * v + = + + +η − Δ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ =∑
(3.19) 將(3.14)至(3.15)式代入(3.12)至(3.13)式並減掉(3.18)至(3.19)式忽略高次項可 得速度修正量[
B k 1/2 e-w]
u P B , T , S , N , W , E nb 2 ,nb 1 k u nb u P P , 2 1 k A h*[(H' z) ] y g -u' A A 1 u' + = + + +η − Δ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ =∑
(3.20)[
B k 1/2 n-s]
v P B , T , S , N , W , E nb 2 ,nb 1 k v nb v P P , 2 1 k ] ) z ' H ( [ * h A x g -v' A A 1 v' + = + + +η − Δ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ =∑
(3.21) 將(3.20)至(3.21)式代入自由液面層連續方程式(3.6)，可得 H' P S S N N W W E E P PH' A H' A H' A H' A H' M A = + + + + (3.22) 其中，質量殘餘量 H' P M 表示如下 m -s n 1/2 k K 1 m w -e 1/2 k 0 n H' P ] ) x * v * h ( ) y * u * h [( t ) y x H ( ) y x * H ( M Δ ρ + Δ ρ + Δ Δ Δ ρ − Δ Δ ρ = + = +∑
(3.23) H’可由(3.22)式求解，然後代入(3.20)與(3.21)計算 u’與 v’，並重複所有步驟 迭代至質量殘餘量 H' P M 消失，即流場收斂。其演算流程如圖 3.2 所示。3.3 模式檢討
擬似三維為能合理決定自由液面位置，目前仍需仰賴水深積分二維淺水 波子模式，一方面，就實際物理現象而言，流場、水位與水深方向流速剖面 三者是相互影變化的，對於存在垂直環流之流場，或流場垂向分佈之不均勻 性高，二維子模式計算所得之水深是否合理，對於水利工程領域所關心的問 題，如模擬結果之水深與流速剖面是否合理，水位、流況是否反映底床粗糙 度及流量變化等問題，鮮少被討論，另一方面，二維子模式亦增加模式計算 量。 本研究試圖以擬似三維分層積分之方式，改善上述三維模式不易應用於 水利工程領域之缺點，嘗試過程中許多關鍵問題無法解決，因而分層模式(一) 之發展終未成功建立，茲將幾個關鍵原因整理如後，供後續相關研究參考。 1. 自由液面流流速垂向分佈(剖面)與流場內部之紊流強度相互影響，本研 究分層模式(一)缺乏反應此現象之機制，且模式收斂受流場起始值設定 之影響非常大。 2. 對於淺水波而言，由於水平格點與層厚度之尺度差異相當大，對於任意 點之流速而言，層與層之間流速無論在物理上或數值上之關聯均比水平 相鄰格點密切，分層模式(一)在層與層之間垂向上的連結係由層介面上 之剪應力項τt、τb，該剪應力之離散化由上下兩層速度中央差分而得， 模式並無法在垂向上反應實際物理現象。 3. 傳統擬似三維之模擬方法，水位已先由二維子模式決定，三維模式在計 算邊界固定之情況下，僅計算流速之分佈，本研究摒棄二維子模式，使 得每次迭代計算上所得之流場分佈均會反映在水深變化上，因此，數值 計算並不穩定。 4. 三維模式邊界條件設定繁雜，現場觀測資料無法提供模式所需，邊界條 件之設定是否合理將影響模擬結果之正確性，而起始值之設定亦密切影 響數值解之收斂性，分層模式(一)之數值解收斂仍明顯受起始值與上下 游邊界處理方式之影響。
經由而分層模式(一)之發展經驗，本研究提出分層模式(二)－考慮層介 面連續，能有效解決解決分層模式(一)所遭遇之難題，其建立方式詳述於第 四章。
第四章 分層模式(二)－考慮層介面連續
本研究最先嘗試分層模式(一)之建立，但於案例測試上收斂困難，究其 原因，一是數學上，控制方程式缺乏對反應對水柱層與層之間的流場密切關 聯之特性，一是數值上採取用於等向性流場使用計算流體力學SIMPLE求解 方法，不容易反應淺水波流場非等向之特性，分層模式(二)係針對分層模式 (一)之發展經驗，修正其物理上與數學上之不足，茲將分層模式(二)之建立 方式詳述如後。4.1 層介面剪力連續條件限制式
在不可壓縮流與靜水壓假設下，非定常流擬似三維分層積分淺水波模式 之控制方程式已詳述於第二章，為方便閱讀，茲將本節相關之動量方程式、 連續方程式、水位變動方程式重述如下： 層積分動量方程式 k zx 1 k zx 2 / 1 k h 2 / 1 k h 2 / 1 k B k 1 k 2 / 1 k 2 / 1 k 2 / 1 k 2 / 1 k 2 / 1 k ) ( ) ( ) y u h ( y ) x u h ( x ] x ) z H ( gh [ ) Fu ( ) Fu ( ) u hv ( y ) u hu ( x ) hu ( t ρ τ − ρ τ + ∂ ∂ ν ∂ ∂ + ∂ ∂ ν ∂ ∂ + ∂ − η + ∂ − = − + ∂ ∂ + ∂ ∂ + ∂ ∂ + + + + + + + + + + (2.24)k zy 1 k zy 2 / 1 k h 2 / 1 k h 2 / 1 k B k 1 k 2 / 1 k 2 / 1 k 2 / 1 k / 1 k 2 / 1 k ) ( ) ( ) y v h ( y ) x v h ( x ] y ) z H ( gh [ ) Fv ( ) Fv ( ) v hv ( y ) v hu ( x ) hv ( t ρ τ − ρ τ + ∂ ∂ ν ∂ ∂ + ∂ ∂ ν ∂ ∂ + ∂ − η + ∂ − = − + ∂ ∂ + ∂ ∂ + ∂ ∂ + + + + + + + + + + (2.25) 層積分連續方程式 0 F F ) hv ( y ) hu ( x k 1/2 ∂ k 1/2 + k 1 − k = ∂ + ∂ ∂ + + + (2.15) 水位變動方程式 0 )] hv ( y ) hu ( x [ t H H K 1 m 2 / 1 m 2 / 1 m t t t = ∂ ∂ + ∂ ∂ + Δ −
∑
= + + Δ + (2.18) 然而，方程式(2.24)與方程式(2.25)除了求解變數 u、v、H 外，尚有未 知數(τzx)k、(τzx)k+1、(τzy)k、(τzy)k+1等層上、下邊界之x、y方向剪應力存在， 這些剪應力項自由液面流於水深與流速分佈剖面計算之關鍵項目，未適當處 理將造成模式之收斂問題，經過分層模式(一)之發展經驗，本研究引入蔡 (2001)成功應用於地下水流模擬之二次形狀函數與層介面水頭、水流通量連 續之觀念，來評估層介面剪應力項，其處理方式詳述如後。 方程式(2.23)與(2.24)中之未知數(τzx)k、(τzx)k+1、(τzy)k、(τzy)k+1等層上、 下邊界之x、y方向剪應力存在，這些透過層積分所得之剪應力項自由液面流 於水深與流速分佈剖面計算之關鍵項目，本研究計算三維流場時，水深與流 速剖面係由控制方程式計算收斂，未使用二維淺水波子模式先行決定自由液 面位置，因此層介面剪應力若未合理估算，則將造成水位劇烈震盪。部水頭之內插函數，並將模式應用於濁水溪沖積平原之地層下陷模擬；本研 究引入該二次形狀函數，並將層介面水頭、水流通量連續之觀念，修改為自 由液面流中層介面速度與剪力連續。本研究使用之方法為由層介面上配置控 制方程式中物理量 u 與 v ，並以二次形狀函數描述物理量於 u 與 v 層內 部之分布情形，滿足層介面之速度連續條件可求得層內部二次形狀函數之係 數，因此，層上、下介面之流速垂向梯度便可由形狀函數求得，又同一層介 面之剪力可分別由相鄰上、下兩層之流速垂向梯度與垂向紊流黏滯係數求 得，滿足層介面上、下兩層不滑移可得到層介面剪力連續條件限制式；同理， 將上述步驟應用至水柱各層，得到水柱方程組。詳細推導過程如下： 假設水柱各層內部物理量 φ 之分佈可用二次形狀函數描述，φ 為流速 u 或 v ，層上、下介面位置分別為 ηk 與 ηk+1 ，因此，二次形狀函數可表 示如下 2 cz bz a ) z ( = + + φ (4.1) 二次形狀函數具有三個未知係數，除層積分物理量外，須於層上、下介面額 外配置物理量，因此，形狀函數於上、下層介面須滿足實際物理量，即 k a ) 0 z ( = = =φ φ (4.2) 1 k 2 z c z b a ) z z ( =Δ = + Δ + Δ =φ + φ (4.3) 定義層平均物理量 φk+1/2 dz ) cz bz a ( z 1 z 0 2 2 / 1 k
∫
Δ + = Δ + + φ (4.4) (4.1) 式之係數滿足(4.2)至(4.4)式之條件，可得如下k a=φ (4.5) ) 6 4 2 ( z 1 b − φk+1− φk + φk+1/2 Δ = (4.6) ) 6 3 3 ( z 1 c 2 φk+1 + φk − φk+1/2 Δ = (4.7) 其中，層下部介面位置 ηk 當做原點，Δz=ηk+1−ηk，將(4.5)至(4.7)式代入(4.1) 並取一次導數，則第 k 層上、下層介面流速垂向梯度可表示如下 第 k 層內上介面流速垂向梯度 ) 6 2 4 ( z 1 z k 1 k k 1/2 1 k + + η φ − φ + φ Δ = ∂ φ ∂ − + (4.8) 第 k 層內下介面流速垂向梯度 ) 6 4 2 ( z 1 z k 1 k k 1/2 k + + η φ + φ − φ − Δ = ∂ φ ∂ + (4.9) 將上下層介面流速垂向梯度(4.8)式與(4.9)式代入層積分動量方程式 (2.23)式與(2.24)式，層積分動量方程式可改寫如下 ) u 6 u 12 u 6 ( ) y u h ( y ) x u h ( x ] x ) z H ( gh [ ) Fu ( ) Fu ( ) u hv ( y ) u hu ( x ) hu ( t k 2 / 1 k 1 k k , v 2 / 1 k h 2 / 1 k h 2 / 1 k B k 1 k 2 / 1 k 2 / 1 k 2 / 1 k 2 / 1 k 2 / 1 k 2 1 − + ν + ∂ ∂ ν ∂ ∂ + ∂ ∂ ν ∂ ∂ + ∂ − η + ∂ − = − + ∂ ∂ + ∂ ∂ + ∂ ∂ + + + + + + + + + + + + (4.10)
) v 6 v 12 v 6 ( ) y v h ( y ) x v h ( x ] y ) z H ( gh [ ) Fv ( ) Fv ( ) v hv ( y ) v hu ( x ) hv ( t k 2 / 1 k 1 k k , v 2 / 1 k h 2 / 1 k h 2 / 1 k B k 1 k 2 / 1 k 2 / 1 k 2 / 1 k 2 / 1 k 2 / 1 k 2 1 − + ν + ∂ ∂ ν ∂ ∂ + ∂ ∂ ν ∂ ∂ + ∂ − η + ∂ − = − + ∂ ∂ + ∂ ∂ + ∂ ∂ + + + + + + + + + + + + (4.11) 對單一層而言，x 方向層積分動量方程式(4.10)除了層平均流速 uk+1/2 未知 外，尚有上、下層介面流速 uk+1、uk 未知，因此，需有額外之方程式方能求 解流場。 考慮第 (k-1) 層與第 k 層間介面，由於緊鄰層介面兩層不滑移，則必 須滿足速度連續條件 k u u u -k k ≡ = η η+ (4.12) 與剪力連續條件 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ φ ∂ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ν = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ φ ∂ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ν + η + η k -k z z v,k 21 2 1 -k v, (4.13) 將第 (k-1) 層上介面與第 k 層下介面流速垂向梯度，參考(4.8)式與(4.9) 式，代入剪力連續條件(4.13)式，可得 x 方向之剪力連續條件限制式
(
1 r)
u 6r u 2r u 0 4 u 6 u 2 k 1 k k k k k k k 1 2 1 2 1 + + − + = − − + + − (4.14) 其中， 2 1 2 1 k , v k , v k r − + ν ν = 為上、下相鄰兩層垂直紊流黏滯係數之比。同理，重複上述 步驟， y 方向之剪力連續條件限制式(
1 r)
v 6r v 2r v 0 4 v 6 v 2 k 1 k k k k k k k 1 2 1 2 1 + + − + = − − + + − (4.15) 另外，底床之不滑移邊界即為 u1= 0、v1= 0，而底床之滑移邊界條件則可以 u 和 v 取代(4.9)式中之 φ 改寫為 2 3 v, x B, 1 2 / 3 2 z ) u 4 u 6 u 2 ( z 1 z u B ρν τ = − + − Δ = ∂ ∂ η = (4.16) 2 3 v, y B, 1 2 / 3 2 z ) v 4 v 6 v 2 ( z 1 z v B ρν τ = − + − Δ = ∂ ∂ η = (4.17) 同理，可將自由液面剪力條件 2 1 K v, x S, z S z u + η = ρν τ = ∂ ∂ 與 2 1 K v, y S, z S z v + η = ρν τ = ∂ ∂ 亦可以u 和 v 取代(4.8)式中之 φ改寫為 2 1 K v, x S, K 2 / 1 K 1 K 6u 2u ) u 4 ( z 1 + + + ρν τ = + − Δ (4.18) 2 1 K v, y S, K 2 / 1 K 1 K 6v 2v ) v 4 ( z 1 + + + ρν τ = + − Δ (4.19) 將(4.10)式、(4.14)式應用至各層，配合底床滑移(4.16)式或不滑移邊界 條件 u1= 0 及自由液面剪力邊界條件(4.18)式，重新整理可得到 u 之水柱層 積分控制方程組；同理，將(4.11)式、(4.15)式應用至各層，配合底床滑移(4.17) 式或不滑移邊界條件 v1= 0 及自由液面剪力邊界條件(4.19)式重新整理可分 得到 v 之水柱層積分控制方程組。4.2 數值架構
4.2.1 有限體積法
有限體積法(Finite Volume Method)在方程式數值離散時採用控制體積 之觀念，可適用於各種形狀之邊界，且介面相鄰兩元素所採用之通量相等， 守恆性良好，已廣泛應用於計算流體力學。本研究在垂直方向上採用分層積 分之方式離散控制方程式，在水平方向則採用有限體積法離散控制方程式， 以得控制方程式最佳之守恆性。在離散控制方程式時，因控制方程式具有三 個變數 u、v、w(或H)，這三個變數既是待求變數，亦用以推求元素邊界通 量，一般常用之方法有交錯格點(staggered grid)與非交錯格點(non-staggered grid)，非交錯格點將所有待求變數均放置於相同點上，因此僅使用一套格 點，在有限差分法中計算方便，但已有許多研究指出非交錯格點求解之流場 容易出現數值震盪之現象(checker board)，使用交錯格點在有限體積法中方 便計算元素邊界通量，亦可避免數值震盪現象，交錯格點物理配置方式及其 相對之控制體積元素如圖3.1所示。 考慮已經由積分之層積分動量方程式(4.10)與(4.11)式，將變數分別配置 於圖3.1不同之格點上，使用有限體積法將動量方程式離散化，可表示成以 下之形式
(
)
y ] ) z H ( gh [ y x ) u 6 u 12 u 6 ( x ) y u ( h y ) x u ( h y x ) Fu ( ) Fu ( ) x u hv ( ) y u hu ( t / ] ) y x hu ( ) y x hu [( w e 2 / 1 k B k 2 / 1 k 1 k v s n 2 / 1 k h w e 2 / 1 k h k 1 k s n 2 / 1 k 2 / 1 k w e 2 / 1 k 2 / 1 k t 2 / 1 k t t 2 / 1 k Δ − η + − Δ Δ + − ν + ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ Δ ∂ ∂ ν + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ Δ ∂ ∂ ν = Δ Δ − + Δ + Δ + Δ Δ Δ − Δ Δ − + + + − + − + + − + + − + + + Δ + + (4.20)(
)
x ] ) z H ( gh [ y x ) v 6 v 12 v 6 ( x ) y v ( h y ) x v ( h y x ) Fv ( ) Fv ( ) x v hv ( ) y v hu ( t / ] ) y x hv ( ) y x hv [( s n 2 / 1 k B k 2 / 1 k 1 k v s n 2 / 1 k h w e 2 / 1 k h k 1 k s n 2 / 1 k 2 / 1 k w e 2 / 1 k 2 / 1 k t 2 / 1 k t t 2 / 1 k Δ − η + − Δ Δ + − ν + ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ Δ ∂ ∂ ν + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ Δ ∂ ∂ ν = Δ Δ − + Δ + Δ + Δ Δ Δ − Δ Δ − + + + − + − + + − + + − + + + Δ + + (4.21) 其中，上標 t 與 t+Δt 表示計算之時刻，下標 (e,w) 與 (n,s) 表示控制體積 於 (東,西) 與 (南,北) 控制面上之動量通量，e、w、n、s 控制面位置如圖 3.1所示。或可表示成下列之傳輸通式(
)
y ] ) z H ( gh [ y x ) 6 12 6 ( ) y ( D ) x ( D ) G ( ) G ( ) G ( ) G ( t / ] ) ( ) [( w e 2 / 1 k B k 2 / 1 k 1 k v s n 2 / 1 k x w e 2 / 1 k x k z 1 k z s n 2 / 1 k y w e 2 / 1 k x t 2 / 1 k t t 2 / 1 k Δ − η + − Δ Δ φ + φ − φ ν + ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ φ ∂ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ φ ∂ = φ − φ + φ + φ + Δ φ ∀ − φ ∀ − + + + − + − + + − + − + + Δ + + (4.22) 其中， φ 為物理量 u 或 v ，∀為控制體體積，Gx、Gy、Gz 分別為 x、y、 z 方向上之對流通量，表示如下 y ) 2 u u h( Gx,e E,k 1/2 P,k 1/2 Δ + = + + ) y 2 u u h( Gx, w W,k 1/2 P,k 1/2 Δ + = + + x ) 2 v v h( Gy,n N,k 1/2 P,k 1/2 Δ + = + + ) x 2 v v h( Gy,s S,k 1/2 P,k 1/2 Δ + = + + y x F Gz,k+1= k+1Δ Δ Gz,k =FkΔxΔy Fk、Fk+1定義如(2.13)式，而擴散通量Dx、Dy則可表示為4.2.2 數值差分式 對 方 程 式 (4.20) 而 言 ， 可 將 對 流 項 (Gxφ)e,w、 (Gyφ)n,s 與 擴 散 項 (Dx∂φ/∂x)e,w、(Dyφ)n,s 均位於控制體邊界面上，而非控制體中心上，須由內 插而得，茲將本研究採用之內插方式，敘述如後。 一階上風法(upwind scheme)雖然精度不高，但因其穩定性佳，常用在工 程實務之模擬，由於三維模式數值處理複雜，因此對流通量(Gxφ)e,w、(Gyφ)n,s 採用一階上風法 ⎪ ⎩ ⎪ ⎨ ⎧ < ⋅ φ > ⋅ φ = φ 0 ) n v ( if 0 ) n v ( if e E e P e (4.23)
擴散通量(Dx∂φ/∂x)e,w、(Dy∂φ/∂y)n,s 採用中央差分法(central difference scheme)
) ( x 1 ) y h ( ) y x h ( h e h e φE−φP Δ Δ ν ≅ Δ ∂ φ ∂ ν (4.24) 其中，物理量 φ 為 u 或 v, 上標 P 表示控制體積 CV 中心點 P，上標 E 表示相鄰控制體積 CV 中心點 E，其相關位置請參考圖 3.1。 4.2.3 水柱代數方程組
將(4.23)式與(4.24)式代入流速 u、v 之傳輸通量通式(4.22)式，可得到 控制體之動量離散方程式 u t 2 / 1 k , 1 j , i S t 2 / 1 k , 1 j , i N t 2 / 1 k , 1 i E t 2 / 1 k , j , 1 i W t t k , j , i b t t 2 / 1 k , j , i P t t 1 k , j , i t S u A u A u A u A u A u A u A + + + + = + + + − + + + + + − Δ + Δ + + Δ + + (4.25) v t 2 / 1 k , 1 j , i S t 2 / 1 k , 1 j , i N t 2 / 1 k , 1 i E t 2 / 1 k , j , 1 i W t t k , j , i b t t 2 / 1 k , j , i P t t 1 k , j , i t S v A v A v A v A v A v A v A + + + + = + + + − + + + + + − Δ + Δ + + Δ + + (4.26) 其中，下標 i、j 與 k 代表變數沿 x、y、與 z 方向之離散指標，上標 u 與 v 用以區分 u 與 v 離散動量方程式之通量係數，下標 P 代表該控制體中心 點，下標 N、E、W、S 則代表相鄰控制體之中心點，參考圖 3.1，S 為其 餘各項總和。通量係數AP、AE、… 等，可由(4.23)式與(4.24)式代入(4.22) 計算而得，並表示如後 e , x e e x, E (1 )G x D A − −α Δ = w x,w w x, W G x D A +α Δ = n , x n n y, N (1 )G y D A − −α Δ = s x,s s y, S G y D A +α Δ = y x 6 G At = z,k+1 − νvΔ Δ Ab =Gz,k −6νvΔxΔy y ) D D ( G ) 1 ( G x ) D D ( G ) 1 ( G y x 12 t A s , y n , y s , y s n , y n w , x e , x w , x w e , x e v P Δ + + α − − α + Δ + + α − − α + Δ Δ ν + Δ ∀ =
x ] ) z H ( gh [ v t S B k 1/2 n s t 2 / 1 k v − +η − Δ Δ ∀ = + + − (4.27) 其中， ⎪⎩ ⎪ ⎨ ⎧ < > = α 0 G if 0 0 G if 1 e , x e , x e ⎪⎩ ⎪ ⎨ ⎧ < > = α 0 G if 0 0 G if 1 w , x w , x w ⎪⎩ ⎪ ⎨ ⎧ < > = α 0 G if 0 0 G if 1 n , y n , y n ⎪⎩ ⎪ ⎨ ⎧ < > = α 0 G if 0 0 G if 1 s , y s , y s 以方程式(4.25)取代方程式(4.10)，並應用至各分層，則(4.25)式、(4.14) 式應用至各層，配合底床滑移(4.16)式或不滑移邊界條件 u1= 0 及自由液面 剪力邊界條件(4.18)式，可得到 u 之水柱層積分動量離散代數方程組。因 此，對 K 層之水柱而言，共有 K 個層積分方程式、K-1 個剪力連續條件限 制式、兩個邊界條件，與未知變數同樣共有(2K+1)個，則形成單一水柱流速 u 之完整代數方程組。同理，以方程式(4.26)取代方程式(4.11)，並應用至各 分層，(4.26)式、(4.15)式應用至各層，配合底床滑移(4.17)式或不滑移邊界 條件 v1= 0 及自由液面剪力邊界條件(4.19)式，可得到(2K+1)個代數式形成 之 v 之水柱層積分動量離散代數方程組。 4.2.4 求解步驟 對水柱之 u 或 v 動量離散代數方程組，水深方向採全隱式法(fully
implicit scheme)，水平方向採顯式法(explicit scheme)，因為剪力連續條件限 制式(4.14)式與(4.15)式均具有五個未知變數，因而水柱層積分動量離散代數 方程組係數矩陣形成具五條對角線之帶狀矩陣，可用附錄 A 之修正湯瑪斯 法(modified Thomas algorithm)求解。由於流場中同一水柱之水深、流速剖 面、紊流黏滯係數相互密切影響，因此求解程序上先由垂直方向將各個水柱 離散方程組求解，然後代入連續方程式求解垂向速度與水深變化，不僅符合 物理機制，數值上也相對簡單與穩定，求解步驟如圖 4.1 所示並說明如後： a. 設定起始流速與水深。 b. 使用Elder紊流黏滯係數關係式計算紊流黏滯係數，其係數 α 以(2.35) 式估算。 c. 以方程式(4.25)、(4.14)式並應用至各層，配合底床滑移(4.16)式或不滑 移邊界條件 u1= 0 及自由液面剪力邊界條件(4.18)式重新整理可得到 (2K+1)個代數式形成之 u 之水柱層積分動量離散代數方程組求解 u。 d. 以方程式(4.26)、(4.15)式應用至各層，配合底床滑移(4.17)式或不滑移 邊界條件 v1= 0 及自由液面剪力邊界條件(4.19)式重新整理可得到 (2K+1)個代數式形成之 v 之水柱層積分動量離散代數方程組求解 v。 e. 將 (c)、(d) 求得之流速 u、v 代入連續方程式求解 垂向速度 w 與 水深 H。 f. 重複步驟 (b)-(e) 迭代修正非線性對流項。 g. 判對 (f) 步驟是否達到要求之收斂條件，若收歛則跳至下一時刻。
4.2.5 數值穩定限制條件 本研究之數值方法採採垂向隱式與水平顯式之半隱式法求解，對擬似 三維之淺水波方程組必須滿足 CFL 收斂條件(Courant- Friedrichs-Levy condition) 2 2 y 1 x 1 gH t 2 CFL Δ + Δ Δ = (4.28) 且透過連續方程式求解垂向速度亦採顯式法，故 w z t≤ Δ Δ (4.29)
4.3 小結
分層模式(二)之發展方式係針對分層模式(一)所遭遇之困難作改進，其 改進方式包含數學上及數值上之改進，使得分層模式之數值解具有穩定、快 速收斂之特性，茲將改進方式之依據與模式收斂特性之關係敘述如後： a. 考慮淺水波之物理非等向性流場之物理機制，淺水波流場因水平與垂 直尺度之差異，流場垂向變化大，因此，引入層內部二次函數，利用 剪力連續條件導出層介面剪力連續條件限制式，對一水柱建立水柱水 柱方程式，改進其數學特性。 b. 底床邊界條件使用不滑移邊界條件，即Dirichlet邊界條件，使解在數 學上具唯一性，較底床、液面都使用Neumann邊界條件穩定。 c. 數值差分使用垂向全隱式、水平顯式之半隱式法求解，確保垂向上之 數值穩定，並簡化水平上之計算方式，使模式容易建立卻兼顧數值解之收斂性。 d. 求解步驟放棄水深二維平均模式常用之水平ADI數值方法，先逐一對 單一水柱方程式求解，建立水柱上符合物理意義之暫態流速剖面，進 而利用迭代方法讓數值解收斂。 e. 使用有限體積法離散控制方程式，上游流速邊界設定為單寬流量均勻 非配於任意時刻任意迭代之水深上，使上游流量邊界之處理正好符合 有限體積法之守恆特性，巧妙避免流速邊界引起之數值震盪。
