• 沒有找到結果。

水庫沉滓運移模擬邊界條件之影響分析

N/A
N/A
Protected

Academic year: 2021

Share "水庫沉滓運移模擬邊界條件之影響分析"

Copied!
93
0
0

加載中.... (立即查看全文)

全文

(1)

土木工程學系

碩士論文

水庫沉滓運移模擬邊界條件之影響分析

Effect of Boundary Conditions on Sediment Transport

Modeling for Reservoirs

研 究 生:黃建翔

指導教授

楊 錦 釧

博士

指導教授

謝 德 勇

博士

(2)

Effect of Boundary Conditions on Sediment Transport

Modeling for Reservoirs

研 究 生:黃建翔 Student:Chien-Hsiang Huang 指導教授:楊錦釧 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 August 2012

Hsinchu, Taiwan, Republic of China

(3)

致謝

在為期兩年的這段研究所求學期間,承蒙恩師楊錦釧教授與謝德勇 博士的悉心指導,使我得以增長學識並順利完成學位論文。兩位在學術 上的專業指導及為人處事態度的榜樣,更是令我獲益良多。感謝口試委 員賴進松教授、連和政博士在論文審查期間給予諸多寶貴的意見,使本 論文能夠更加完善。學生建翔在此向各位老師致上萬分的感謝與敬意。 交大在學期間,特別感謝文祿學長、胤隆學長、世偉學長、浩榮學 長、建華學長、弘恩學長、昇學學長、仲達學長、仁凱學長、彥酉學長、 聖翔學長、家偉學長、緁玲學姐、唯泰學長、紹唐學長、彥瑜學姐、東 洲學長在學業、研究以及生活上,給予我的協助與照顧。感謝同窗好友 昀直、舒勤、芳綺、信富、韋豪、家榮、岱玲,學弟妹醇國、瑋廷、蓉 瑩、亞雯、于軒、健賓,豐富了我的研究所生活並與我互相鼓勵一同學 習。感謝好友可忠、國高中同窗、大學同窗、交大撞球社好友們,在我 迷惘或遭遇困難時聆聽並給予我中懇的建議。 最後感謝家人在這段時間的體諒並給予我鼓勵,感謝親愛的爸媽和 哥哥,你們是我最安心的避風港,不但讓我生活安逸,所有難關更因你 們而化解。感謝爺爺奶奶、外婆、姑姑姑丈、叔叔嬸嬸給我的愛與關懷, 使我即使遠在新竹仍備感溫暖。亦十分感謝中正八八期的各位長輩們一 路以來對建翔的照顧。最後,謹將此論文獻給你們,與你們一同分享我 的喜悅。

(4)

就水庫泥砂課題而言,泥砂入庫後之運移現象主要包含上游段粗顆粒 造成之三角洲淤積行為,以及細顆粒入庫區後潛入形成異重流之運移行 為。以數值模式的角度來看,沉滓運移現象模擬結果之準確度常受上游河 床載邊界條件之設定,以及底床邊界再懸浮機制之影響。本研究之目的, 即 擬 引 用 謝(2003) 及 鍾 (2012) 所 發 展 之 二 維 (RESED2D) 以 及 三 維 (RESED3D)沉滓運移模式,探討上游及底床邊界條件對沉滓運移之影響。 為探討水庫沉滓運移模擬邊界條件之影響,分別先以三維模式模擬實 驗案例,探討底床邊界條件對濃度傳遞的影響;再以二維模式模擬實驗以 及現地案例,探討入流邊界條件對底床變遷的影響。最後,在簡單垂向入 流濃度的假設條件下,將三維模式應用於石門水庫,以檢視模式之實用 性。 關鍵字:三維模式、沉滓運移、水庫、懸浮載、底床載  

(5)

In respect of sediment issue in reservoirs, the transport of sediment mainly involves with coarse particles in upstream that forms delta and fine articles in reservoirs that result in transport of density current. From numerical modeling view point, the accuracy of simulated results of sediment transport is often influenced by the bed-load boundary condition in upstream and resuspension mechanism of the riverbed boundary. This research adopts Hsieh’s (2003) two-dimensional (RESED2D) and Zhong’s (2012) three-dimensional (RESED3D) mobile-bed models to discuss the effect of upstream and riverbed boundary conditions on sediment transport. To examine the effect of boundary conditions of sediment transport modeling for reservoirs, this study first simulated experimental and on-site case with 2D model for discussion of inflow boundary condition’s effect on the change of riverbed. Then, the experimental cases were simulated with use of 3D model for discussion of the effect of riverbed boundary on delivery of concentration. Finally, with a hypothetical inflow concentration profile, the 3D model was applied to the Shi-men Reservoir to show its applicability.

Keyword: 3D model, sediment transport, reservoir, suspended load, bed load

(6)

致謝 ... I  摘要 ... II  ABSTRACT ... III  目錄 ... IV 圖目錄 ... VII 符號表 ... X 第一章 緒論 ... 1 1.1 研究動機與方向 ... 1 1.2 文獻回顧 ... 2 1.3 研究目的與方法 ... 3 1.4 章節介紹... 4 第二章 理論基礎 ... 6 2.1 水理控制方程式 ... 6 2.1.1 三維部分 ... 6 2.1.2 水平二維部分 ... 7 2.1.3 垂直部分 ... 8  2.2 動床控制方程式 ... 10  2.3 輔助關係式 ... 11 2.3.1 層流與紊流剪應力... 11

(7)

2.3.3 河床載通量 ... 12 2.3.4 懸浮載源 ... 13 2.3.5 作用層源 ... 14 2.3.6 作用層厚度 ... 15 2.4 邊界條件... 15 2.4.1 水理部分 ... 15 2.4.2 動床部分 ... 16  第三章 數值架構 ... 19 3.1 水理部分... 19 3.1.1 求解架構 ... 19 3.1.2 數值方法 ... 22 3.2 沉滓運移部分 ... 24 3.2.1 求解架構 ... 24 3.2.2 數值差分 ... 26 第四章 河床載入砂邊界條件的影響分析... 32 4.1 沖刷案例 ... 32  4.2 淤積案例 ... 33 4.3 沖淤平衡案例 ... 34

(8)

5.1 淤積案例 ... 54 5.2 沖刷案例 ... 55 第六章 模式在石門水庫的應用 ... 59 6.1 石門水庫上游段二維模擬 ... 59 6.2 石門水庫三維模擬 ... 61 第七章 結論 ... 70 7.1 結論 ... 70 7.2 建議 ... 70 參考文獻 ... 72  

(9)

圖2.1 正交曲線座標轉換示意圖 ... 17  圖2.2 Σ座標轉換示意圖 ... 17 圖2.3 水深方向流速剖面示意圖。 ... 18 圖2.4 作用層示意圖 ... 18 圖3.1 模式計算流程圖 ... 29 圖3.2 水平二維模式控制體積法示意圖。 ... 30 圖3.3 交錯格網(staggered grid)示意圖 ... 31  圖3.4 垂直模式控制體積法示意圖(計算區域)。 ... 31  圖4.1 Suryanarayana 試驗之沉滓粒徑分佈曲線 ... 37  圖4.2(a) 非均勻沉滓沖刷案例之底床沖刷時變圖(經驗法) ... 38  圖4.2(b) 非均勻沉滓沖刷案例之底床沖刷時變圖(實測法) ... 39  圖4.3(a) 淤積案例之底床沖刷時變圖(1)(假設上游邊界河床載=0) 40  圖4.3(a) 淤積案例之底床沖刷時變圖(2)(假設上游邊界河床載=0) 41  圖4.3(b) 淤積案例之底床沖刷時變圖(1)(率定法) ... 42  圖4.3(b) 淤積案例之底床沖刷時變圖(2)(率定法) ... 43  圖4.3(c) 淤積案例之底床沖刷時變圖(1)(率定法,並同時讀取上游邊 界懸浮載濃度為0) ... 44  圖4.3(c) 淤積案例之底床沖刷時變圖(2)(率定法,並同時讀取上游邊

(10)

圖4.3(d) 淤積案例之底床沖刷時變圖(1)(經驗法) ... 46  圖4.3(d) 淤積案例之底床沖刷時變圖(2)(經驗法) ... 47  圖4.4(a) 沖淤平衡案例之底床沖刷時變圖(假設為清水入流) ... 48  圖4.4(b) 沖淤平衡案例之底床沖刷時變圖(經驗法) ... 49 圖4.4(c) 沖淤平衡案例之底床沖刷時變圖(實測法) ... 50 圖4.5(a) 沖淤平衡案例之底床沖刷時變圖(假設為清水入流) ... 51  圖4.5(b) 沖淤平衡案例之底床沖刷時變圖(經驗法) ... 52 圖4.5(c) 沖淤平衡案例之底床沖刷時變圖(實測法) ... 53 圖5.1 淤積案例之濃度剖面圖 ... 56 圖5.1 淤積案例之濃度剖面圖(續) ... 57 圖5.2 沖刷案例之濃度剖面圖 ... 58 圖6-1 石門庫區斷面暨模式適用範圍示意圖 ... 62 圖6-2 入庫輸砂量與入庫流量關係圖 ... 63 圖6-3 庫區斷面4- 28底泥取樣粒徑分析結果 ... 63 圖6-4(a) 民國93年之底床沖刷時變圖(假設上游邊界河床載為0) ... 64  圖6-4(b) 民國93年之底床沖刷時變圖(經驗法) ... 64 圖6-4(c) 民國93年之底床沖刷時變圖(率定法) ... 65 圖6-5 民國95年之底床沖刷時變圖(驗證) ... 66

(11)

圖6-7 韋帕颱風垂向濃度分布圖 ... 69  圖6-7 韋帕颱風垂向濃度分布圖(續) ... 69            

(12)

B 渠道寬度; C 懸浮沉滓濃度; 懸浮沉滓水深平均濃度; 近底床沉滓平衡濃度; 近底床沉滓濃度; Chezy 糙度係數; 摩擦係數; 懸浮沉滓中值粒徑; 水深; 河床沉滓代表粒徑; 編號第m 組沉滓之粒徑; 懸浮載向下通量; 懸浮沉滓代表粒徑; 無因次顆粒粒徑; 作用層厚度; 均方根誤差; 懸浮載向上通量; 福祿數; 、 主要與側方向受密度差異量影響之移流作用力; 重力加速度; 正常水深: 高含砂效應下,最大差異發生處之水深; C  , C Cc c  2 / f c Cg c  50 dDb Dm Ds Dms D*m Dp Erms Es Er Ff f   gHm H

(13)

、 、 方向轉換係數; 多項式流變關係之水深平均層流阻力係數; 粗糙高度; 卡門係數(Von Kármán’s) 沉滓代表粒徑總數目; 曼寧糙度係數; 動量方程式之壓力項: 孔隙率; 總流量; 單寬流量: 方向某一粒徑m 之河床載通量: 彎道任一點之曲率半徑: 彎道左岸之曲率半徑: 水力半徑; 彎道中心線曲率半徑; 渠道坡降; 某一粒徑m 之懸浮載源; 砂比重; 作用層源; 水面坡降; 某一粒徑m 水流挾砂能力; 輸送參數; 、 、 有效剪應力項; 1 h h2   l Ks k    MnP  pQqm bi qi RL Rh Rc r  0 Sm s Sg sE Sw S  *m Sk T  11 T T12 T22

(14)

主流方向深度平均速度; 主流方向速度; 斷面平均速度; 剪力速度; 高含砂效應下,最大差異發生處之水深平均速度; 清水流下,最大差異發生處之水深平均速度; 臨界剪力速度; 側方向深度平均速度; 側方向速度; z 方向速度; 沉滓沉降速度; 濃度影響下沉滓沉降速度; 卡式座標系統之垂直方向座標; 底床高程; 水面高程; 最接近底床之垂直格網厚度; 懸浮沉滓傳輸係數; 、 流變關係參數常數項係數; 第m 組代表粒徑在作用層中所佔百分比; 、 流變關係參數指數項係數; 、 懸浮載連續方程式中水平、垂直方向之紊流擴散係數; 層流黏滯係數; 紊流黏滯係數; uuU  * Um Uw U  *c Uvvwf wfh wzb zs zl zu Z  1  2m   1  2sH  sVl   t  

(15)

垂直方向渦流黏滯係數; 動力黏滯係數; 賓漢黏滯係數; 第m 組代表粒徑承受之無因次剪力; 第m 組代表粒徑之無因次臨界剪力; 含砂水流流體密度; 乾砂密度; 賓漢降伏應力; 高含砂效應下,最大差異發生處之底床剪力; 清水流下,最大差異發生處之底床剪力; 、 底床剪應力在ξ、η 方向之分量; 泥砂顆粒啟動之臨界剪力; 含砂水流比重量; 、 平面上兩正交座標方向; 水深分層之垂直方向座標,無因次水深; 座標方向之垂向速度; 時間間距; 、 、η 方向兩相鄰格網點之格網間距。 上標 時刻之已知變數; 時刻之未知變數; 與 間之未知變數; 時間平均; V   t   B   m   cm     s   B   m   w   1 b  2 b   c           t        n n t  1 ( 1) n  n t 1 2 ( 1) n  n t n t

 

(16)

 

差異量(流速為與深度平均流速之差,密度為與清水密度差) 紊流擾動量。 下標 方向格點編號; 方向格點編號; 方向格點編號; 1、2 = 、 方向變數; p、e、w、n、s = 通量所在格點位置; 近底床高度5%水深處變數; 近底床變數。

 

 

i j k   ab

(17)

第一章 緒論

1.1 研究動機與方向 沖淤是河川中常見的重要現象。天然河道的沖刷受水文、地文、河 床質及水理特性等自然因素的影響。在河道中,通常因坡陡湍急,使得 河床載對沖淤行為影響甚小。然而水庫的沉滓運移,大量的入砂以及迴 水效應使流速減慢,導致落淤嚴重時,則河床載勢必不可忽略。每年颱 風帶來的大量降雨對台灣地區而言雖為必要且珍貴的水資源,但降雨挾 帶大量泥砂入庫減少水庫庫容以及於壩址附近形成渾水潭影響取水口 與淨水廠的正常運作,亦為急需解決的問題。 然而面對一欲解決之工程問題,一般常用物理模型試驗或數值模擬 來協助瞭解這些物理現象。設置物理模型試驗為科學家常用之方法,其 為參考原型並依比例縮尺建造一模型並進行物理試驗。試驗結果可信度 雖高,然而其設置所需經費龐大、試驗時間長以及大量消費人力與空間 的問題,再加上模型試驗存在尺度效應(scale effect)、邊界與量測儀器產 生之誤差,使得物理模型在應用上常有相當的限制和困難。數值模擬乃 利用數學模式獲得欲知的物理現象,在某些物理假設條件下,仍可獲得 相當合理之精確度與可靠性。相較於物理模型試驗,數值模擬具有較經 濟及高效率之優點,或具有再利用性,因此許多科學家不斷致力於數值 模擬之發展。 對於河川、湖泊及海洋等大範圍之流場分析,一維與二維模式已廣 泛應用在實務上,其相對於三維模式具有計算快速及容易收歛的優點。 但是當流場在深度方向的分佈為欲探討之問題時,垂向的資訊則相對重 要而不可忽略。相較於全三維模式,擬似三維模式則具有計算快速和容 易收歛的優點,且亦能獲得物理量深度方向分佈的資訊,適合水利工程 實務上的應用。 而粗顆粒形成的河床載一直以來量測不易而缺乏現地數據。故本研

(18)

究嘗試以模式模擬水庫沉滓時,率定出河床載與懸淨載之間存在之比例 關係,做為邊界條件,使之能適當反應出水庫底床淤積行為導致床形的 變化。更進一步再以率定結果做為三維模式模擬沉滓運移時的邊界條 件,期以擬似三維模式計算垂向濃度分布,反應出水庫泥砂落淤的機 制,提供日後相關研究作一參考。 1.2 文獻回顧 河道沉滓之運移行為影響著底床沖淤變化,河口與海岸型態的變遷 等,進而影響水工結構物設計與操作策略研究等相關的課題。因此,如 何正確與合宜的分析河道沉滓的運移歷程,為水利工程一重要的研究課 題。誠如 Dawdy and Vanoni (1986)所言,沒有任何一個沉滓運移模式能 合理解決所有輸砂相關的課題,所以每個模式均有其特定的發展方向, 以解決其特定的輸砂問題。

一般而言,大型沖積河流河床載輸運率(bed-load transport rate)約佔

總沉滓載輸運率之5%~25%,底床質由粗顆粒(coarse materials)組成者所

佔比例更高。若以沉滓傳輸型態的方式來加以分類,河床載(bed material load)模式如 TABS2(1985)及 GSTARS(2001)等,為最普遍的模式型態, 在泥顆粒較粗的情況下,此類模式仍具有相當的實用價值,缺點是類模

式無法有效區分懸浮載與河床載之運移。懸浮載模式如 Celik and

rodi(1988)、Van Rijn(1990)等及 Ziegler and Nisbet(1995)等,較適用於河 床載所佔比例不高的情況,或可用於探討凝聚性沉滓的相關案例。Bell and Sutherland(1983)發展之河床載模式,較適用於粗顆粒沉滓模擬及懸 浮相對不重要的情況。CHARIMA(1990)、Spasojevic and Holly(1990)、 MIKE-11(1992)、NETSTARS(1996)等模式以河床載及懸浮載分開計算的 方式,可模擬懸浮載濃度分佈及底床沖淤的情況外,並可合理模擬非平 衡輸砂的問題。

在數值模式的發展上,水深平均二維模式已普遍應用在天然河川的 模擬,且有許多研究均已驗證其模擬結 果的合理性。二維模式如

(19)

TABS2(1985)、Celik and Rodi (1988)、Spasojevic and Holly (1990)、Ziegler and Nisbet (1995)等,可以用來描述物理變量在側方向(水平二維)或深度 方向(垂直二維)的變化情形,在河川寬深比一般均很大的情況下,河川 沉滓運移模擬應以深度平均之水平二維模式為較佳。 近來許多三維模式開始應用在大型流場如河川、湖泊及海洋的模 擬,相較於二維模式,其能直接提供水深方向的資訊。然而三維模式仍 需要花費不少的時間來模擬,因此有許多擬似三維模式的研究提出 (Lardner and Cekirge, 1988;Jin and Kranenburg, 1993;Wang, 1994; Blanckaert and de Vriend 2003;Hung et al., 2008;Lin and Huang, 2008; Herzfeld et al., 2010;Zhang et al., 2011),除了能降低計算成本,在一些 假設條件下,亦能提供合理的三維流速分布資訊。

擬似三維模式中,Lardner and Cekirge (1988)提出的垂直水平分離演

算法 Vertical Horizontal Splitting (VHS)受到許多學者引用(如 Wang, 1994;Lin and Huang, 2008 等),其將水平與垂直流場分開求解,首先 利用水深平均二維模式計算水位分佈與水深平均流速分量,再透過一子 模式獲得流速在水深方向之分佈。VHS 的概念已經應用在海岸、河口 及湖泊等大型水體的流場分析(Jin and Kranenburg, 1993;Wang, 1994; Lin and Huang, 2008;Herzfeld et al., 2010)。Zhang et al. (2011)延用並擴 展 Muneta and Shimizu (1994)發展之擬似三維模式,將原模式之卡氏座 標系統轉換為非正交曲線座標系統,並考慮原模式忽略的延散項進行運 算。

本研究在三維模擬上主要關心底床邊界之設定,文獻中學者藉由解 析解之濃度垂直剖面直接積分(Rouse, 1937; Zhang et al.,2001),量化懸浮 載運移量,並且以經驗公式推估近底床濃度,作為濃度剖面在近底床處 之邊界條件。而關於泥砂再懸浮機制之探討,錢(2011)等發現在近岸河

口之再懸浮機制發生時,粒徑16μm 以下之細懸浮顆粒所佔比例增加。

近來更有學者應用一種新的沉積物再懸浮裝置,模擬中國太湖水域不同 風浪下懸浮物的垂向分布,如尤(2007)等。利用此法求得水柱總懸浮量

(20)

與擾動頻率的定量關係,是目前較適用於淺水水體懸浮物再懸浮之模擬 方法。 1.3 研究目的與方法 本研究之目的為引用鍾(2012)所發展之三維沉滓運移模式,探討底 床沖淤機制。基於水深平均二維水理模式謝(2003),應用垂直水平分離 演算的概念,加入垂直水理模式,以動量方程式求解水深方向之濃度分 佈,另外並在模式中加入表面風剪力及科氏力之影響。水平方向座標系 統採用正交曲線座標,其能適當表達不規則渠道形狀,且座標軸主軸方 向(本研究使用之軸方向)即為主流方向;而在水深方向則採用 σ 座標

系統(Blumberg and Mellor, 1983),其能解決自由水面在固定格網上變動 而影響模式無法準確計算水面之壓力邊界條件的問題,也能將因為起伏 底床產生之不規則格網,轉換為便於計算的矩形格網,如此可得到精確 度較高的模擬結果。為驗證模式的正確性,本文分別採用具有解析解或 實驗量測數據的案例比較模擬結果。 1.4 章節介紹 前面已闡述本研究之動機與方向、文獻回顧、研究目的與方法,本 節將扼要說明本文章之內容,其中第二章理論基礎及第三章數值架構大 致上均摘自鍾(2012)之論文。 第一章為緒論,說明本研究之背景與目的,並回顧相關模式發展的 文獻,提出本研究之方法與研究之重點。 第二章為理論基礎,在正交曲線座標系統下,由三維那威爾-史托克 司(Navier-Stokes)方程式導出模式控制方程式、輔助方程式的使用及邊 界條件的設定,均於本章介紹。 第三章為數值架構,水平二維水理模式及垂直水理模式之數值方法 於本章說明,並簡述模式之計算流程。

(21)

第四章針對模式發展部分,採用具有解析解或實驗量測數據的案例, 來展現模式在三維運算上預測泥砂濃度分布的功能,並驗證模式預測之 合理性與正確性。 第五章將藉由簡單的實驗案例,以及具有實測數據之現地案例,來 探討河床載邊界條件對床型變化的影響,並驗證模式模擬之合理性與正 確性。 第六章使用石門水庫現地實測資料,結合第四與第五章之概念,分 別以二維與三維模試模擬石門水庫洪颱期間之底床變化,並做一綜合分 析。 第七章為結論與建議,對研究成果作綜合性之歸納說明,並針對研 究尚未考量、不盡完備或日後可繼續研究之處提出建議。

(22)

第二章 理論基礎

模式的座標系統方面,為了能適當表達天然河道不規則的幾何形

狀,模式在水平方向採用正交曲線座標系統,如圖2.1 所示;而在垂直

方向則採用σ 座標系統(Blumberg and Mellor, 1983),如圖 2.2 所示,如

此能將不規則的計算區域轉換至矩形計算區域求解,且能便利地處理渠 道的側壁、自由液面及底床邊界。 控制方程式方面,透過座標系統轉換,將三維控制方程式轉換為正 交曲線座標系統方程式,再將此控制方程式作時間平均及水深平均後, 即可推得水深平均二維控制方程式(謝,1994)。將水平二維水理控制方 程式從三維控制方程式中扣除(VHS),即可求得垂直控制方程式。 2.1 水理控制方程式 2.1.1 三維部分 基於不可壓縮流之假設下,對那威爾-史托克司(Navier-Stokes)方程 式取時間平均後,得控制方程式如下: 連續方程式 2 1 1 2 (h u) (h v) (h h w) 0 z       (2.1) 動量方程式 ξ 方向: 2 2 2 1 2 2 1 2 1 2 1 2 1 2 1 1 2( ) ( ) ( ) ( )u ( u ) ( uv) ( uw) uv h u h v h t h h z hh hh hh                        2 2 2 1 2 2 1 2 1 2 1 2 1 2 1 1 2( ) ( ) ( ) ( u ) ( u v) ( u w) u v h u h v h h h z hh hh hh                                   12 1 22 2 2 11 1 12 1 2 13 1 1 2 1 2 1 2 1 1 ( ) ( ) ( ) c h h p f v h h h h h h h z h h h h                                 (2.2) η 方向:

(23)

2 2 2 2 1 1 2 1 1 2 1 2 1 2 1 1 2( ) ( ) ( ) ( v) ( v ) ( uv) ( vw) uv h v h u h t h h z h h h h h h                        2 2 2 2 1 1 2 1 1 2 1 2 1 2 1 1 2( ) ( ) ( ) ( v ) ( u v) ( v w) u v h v h u h h h z h h h h h h                                   12 2 11 1 1 22 2 12 1 2 23 2 1 2 1 2 1 2 1 1 ( ) ( ) ( ) c h h p f u h h h h h h h z h h h h                                  (2.3) z 方向: 在一般淺水的天然水道,垂直方向之動量方程式可用靜水壓分佈來 簡化, 0 p g z    (2.4) 而垂直方向之流速可透過連續方程式求得,以減少模式之計算量。 以上諸式中,、 、z為三維正交曲線座標方向,其中、 為水 平方向,z為水深方向;下標1、2、3 分別代表物理量在、 、z方向 代號;h1h2分別為、 方向之轉換係數;uvw分別為、、z 方向流速;g為重力加速度;t 為時間; 為層流剪應力;( ) 表時間平 均;( ) 表時間平均瞬時擾動量;fc(=2 sin )為科氏力係數; 為地球 自轉角速度; 為緯度;  為密度;p為壓力。 2.1.2 水平二維部分 將 2.1.1 節之三維水理控制方程式(2.1)、(2.2)、(2.3)利用萊布尼茲 法則對深度方向積分,加上運動邊界條件及動力邊界條件,並取深度平 均值,可得水平二維水理控制方程式。 連續方程式 1 2 ( 2 ) ( 1 ) 0 d hh h ud hvd t

   (2.5) 動量方程式 ξ 方向:

(24)

2 1 2 1 2 1 2 1 2 1 s c z h h u u u v u uv v g f v t h h h h h h h                         1 1 2 1 2 12 11 11 12 22 1 2 2 1 ( ) 1 1 1 2 s b h h h T T T T T dh h dh dh d                             (2.6) η 方向: 2 2 1 2 1 1 2 1 2 2 s c z h h v v v u v uv u g f u t h h h h h h h                 2 2 1 2 1 12 22 22 12 11 1 2 1 2 ( ) 1 1 1 2 s b h h h T T T T T dh h dh dh d                             (2.7) 式中, 2 2 11 ( 11 ) s b z z T

 u u dz (2.8) 2 2 22 ( 22 ) s b z z T

 v vdz (2.9) 12 ( 12 ) s b z z T

 uvu v dz  (2.10) 以上諸式中,d 為水深;zs為水面高程;zb為底床高程; 為 i 方si 向水面剪應力; 為 i 方向底床剪應力;bi ( )表水深平均;( ) 表物理量 之空間微變量(例:u u u  );下標 sb分別為水面及底床代號;T

有效剪應力項(effective stress term),其包含層流剪應力、延散剪應力與 紊流剪應力。 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),即可得到垂直水理控制方程式。 ξ 方向: 1 1 1 2 2 2 u u u u u u u v u v u v u u t hhhhhh                       

(25)

2 1 1 1 2 2 1 2 1 2 1 2 1 2 1 2 2 h h h h h uv uv uv vv v h hh hh hh hh h                       13 ( 1 1) 1 ( ) s b c f v Horizontal Diffusion in d d               (2.11) η 方向: 2 2 2 1 1 1 v v v v v v v u v u v u v v t hhhhhh                                     2 2 2 2 1 1 1 2 1 2 1 2 1 2 1 2 2 h h h h h uv uv uv uu u h hh hh hh hh h                       23 ( 2 2) 1 ( ) s b c f u Horizontal Diffusion in d d                (2.12) 式中, 1 1 2 2 u u v v w t h h h h z                                (2.13) b z z d    (2.14) d t d t         (2.15) 1 zb d d d              (2.16) 1 zb d d d              (2.17) 1 z d    (2.18) 2 1 2 11 12 22 1 2 1 ( Horizontal Diffusion in ) T h 2T h T h dh h                    2 12 11 2 2 2 2 2 1 1 1 1 2 1 2 2 2 1 T 1 T u H 1 H h H u dh dh h h h h h                                 1 2 1 1 1 2 2 1 2 2 2 1 2 H v h H h v H h u h h h h h h h h h                                2 1 1 2 2 2 2 2 2 2 1 1 2 2 1 1 2 2 1 1 2 H H H h u v u h v h h h h h h h h h                                     (2.19)

(26)

1 2 1 22 12 11 1 2 1 ( Horizontal Diffusion in ) T h 2T h T h dh h                    2 22 12 1 2 2 2 2 1 2 2 1 2 2 2 2 2 1 T 1 T v H 1 H h H v dh dh h h h h h                              2 2 1 2 1 2 1 1 2 1 2 1 2 H u h H h v H h u h h h h h h h h h                                2 2 2 1 2 2 2 2 2 1 2 1 2 2 1 1 2 2 1 1 2 H H H h v v u h u h h h h h h h h h                                      (2.20) ; 為水平黏滯係數H   l t; 為層流黏滯係數;l  為紊流黏滯係數t * / 6 u d

 (Falconer 1980);u*為剪力速度; 為von Karman’s 係數(約等

於0.41);d為水深。

2.2 動床控制方程式

沉滓運移部分在懸浮載、河床載分開求解概念下,求解懸浮載之質 量傳輸(mass transport)方程式、作用層連續方程式(active-layer continuity equation)、底床連續方程式(bed-layer continuity equation) (謝德勇, 2002),其中懸浮載之質量傳輸為三維控制方程式。由於本研究模式各 項變數之間的空間關係較為複雜,包含懸浮載(suspended load) 濃度、 懸浮載底床濃度、懸浮載向上與向下通量、河床載(bed load)、作用層 (active layer)等,將各變數於深度方向之關係表示於圖 2.4。本研究將沉 滓運移計算分為懸浮載與河床載兩部分,以底床高程 為分界。懸浮載 運移在 之上,河床載運移在 之下,且河床載運移發生於作用層內。 作用層主要是一理論上的假設物理量,模式假設河床載運移與粒徑變化 均僅發生於作用層當中,頂面位置為底床 所在高程。懸浮載之底床濃 度則是位於 ,交換通量也發生於 。以下逐一介紹圖2.4 中的各變數 與計算方式。 沉滓運移控制方程式將沉滓粒徑分為數種代表粒徑,分開計算懸浮 載與河床載,計算各粒徑在懸浮載中的濃度以及在作用層當中之比例, 並計算整體沉滓在交換後所造成之底床沖淤量,將方程式說明如下 : b z b z zb b z b z zb

(27)

質量傳輸方程式 (2.21) 水平擴散項由於 sigma 座標轉換後產生多項二次微分之乘積,本研 究中假設其值可予以忽略。 作用層連續方程式 (2.22) 底床連續方程式 (2.23) 上列諸式中, 某一代表粒徑的懸浮載體積濃度; 、 垂直、 水平方向紊流擴散係數; 泥砂落淤速度; 含砂水流密度; 乾砂密度; 第m 組代表粒徑在作用層中所佔百分比; 孔隙率; 作用層厚度,模式中假設為 0.2~0.5 (m); 作用層源(source of active layer); 底床高程; 、 水流方向、水流側方向方向第 m 組代表粒徑之河床載通量; 為第m 組代表粒徑懸浮載源(source of suspended load)。 2.3 輔助關係式 2.3.1 層流與紊流剪應力 採用Boussinesq 之渦流黏性理論,層流與紊流剪應力可合併表示為 2 2 1 2 2 2 1 2 1 1 sV sV wfh wfh C C C C C C C C uh vh t h h D D D D                           2 2 2 1 2 2 1 1 2 2 1 2 1 1 1 2 2 2 1 ( ) ( ) sH C h h C sH h sH C sH C h h C sH h sH C hh h h h h h h                                               2 1 1 2 1 2 1 2 ( ) (1 ) ( ) ( ) 0 m m m m p s b b s E E p h q h q S S t h h h h                 2 1 1 2 1 1 2 1 2 (1 ) ( ) ( ) 0 m m m M b s b b s m z p h q h q S t h h h h               

  C sVsHfh w   sm   pp ESEb zqb1m 2 m b qm s S

(28)

2 11 1 1 1 2 1 2 H u v h u h h h               (2.24) 2 22 2 2 1 2 1 2 H v u h v h h h               (2.25) 12 2 1 1 2 2 1 2 H h v h u u v h h h h

                 (2.26) 2.3.2 底床剪應力 底床剪應力採用French (1986)之經驗式 2 1 2.5ln 30 1 2.72 b b b s z u u k                (2.27) 2 2 2.5ln 30 1 2.72 b b b s z v v k                (2.28) 式中,ubvb分別為、 方向之近底床流速;z1為近底床流速之 格網與底床間垂直距離;ks為粗糙高度。 2.3.3 河床載通量qb 凝聚性泥砂之通量可假定為零,而非凝聚性沉滓之河床載通量,模式 中採用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.29) 式 中 , 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   顆粒蔡司參數;ss    砂比重。

(29)

另外,泥砂運移過程中,較細之顆粒除了以懸浮載型式移動外,由 於泥砂為非均質,細顆粒可能隱藏於較粗顆粒之間而不易被水流帶走。 為體現此一機制,模式進一步將非均勻粒徑所產生的遮蔽效應納入考慮 並進行修正如下 ( ) i k i b h k b k q

 

q D (2.30) 由於水流在渠道中運行時,河床載源之變化將受到:(1)縱向與橫向之 底 床 坡 降S0S0 ;(2)縱向及橫向之流速u 、 v 的影響;模式採用 Struiksma(1985)之公式,對某一粒徑之河床載通量進行修正(張氏,2005)

方向: qb  hk k bq D( ) cosk  (2.31)  方向: qb  hk k bq D( ) sink  (2.32) 在上面兩式中 1 1 sin tan 1 cos b s b s z f z f                        (2.33) 1 tan v u        (2.34) 以上,

沉滓運移角度; 底床剪應力之方向;fs  沉滓之形狀因子, 1 fs  ;2 2 2 2 50 ( ) ( 1) u v c s D      Shields parameter。 2.3.4 懸浮載源 S 針對非凝聚性之泥砂而言,懸浮載源可定義為懸浮載向下通量與底 床質向上通量之間交互作用的結果。對某一特定粒徑k 而言,受到重力 作用而沉降至底床,其向下之通量可表示為下式:

(30)

k k k d f d S  

w C (2.35) 其中負號代表通量向下,反之為向上之通量; [3.25 0.55ln( k )] k f d k w C C u   (Lin 1984);Ck  顆粒 k 的深度平均濃度;wfk 顆粒k 的沉降速度,其計 算方式依照粒徑大小採用不同關係式(Van Rijn 1984b)如下

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.36) 在水中未設有結構物的渠道中,底床質向上之通量主要是由於水流 剪力對底床之作用,使底床質揚起成為懸浮質,即一般沖刷及束縮沖刷 之現象。其向上通量可表示為: k k k e l k e S

 

w C (2.37) 上式 1.5 0.3 0.015 k k k k e D T C a D  (Van Rijn 1984b);a 砂丘高度之一半;wlk  顆粒k 躍起之速度。 顆粒躍起速度定義為河床質發生跳躍(saltation)離開底床瞬間時的垂 直速度,採用Hu and Hui(1996)提出之經驗公式: 3.2 4.5log 1.2 3.1 1.2 k l w u           (2.38) 式中 ( ) b s gDk       2.3.5 作用層源 Sf 作用層示意圖如圖2.4。作用層源之產生乃肇因於母層(active stratum)

(31)

頂面之升降,當其下降時, (1 ) [( ) ( )] f s s k b m S p z E t         (2.39) 式中( )s k  母層中某一顆粒 k 之粒徑百分比;若母層厚度增加,即 其頂面上升時,則將(2.36)式中 ( )s k以k取代。 2.3.6 作用層厚度 Em

沖刷現象發生時,根據Bennet and Nordin (1977)之研究,可以下式表

示: 1

(

n n

)

m em b b

E

 

C z

z

(2.40) 式中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.41) 式中,Dm  不產生移動的最小顆粒粒徑。 另外,作用層在淤積期間可定義為: 1 ( 1 ) n n n n m m b b E E zz (2.42) 2.4 邊界條件 2.4.1 水理部分 水平二維水理模式考量三種邊界條件設定,分別為渠道入流、渠道 出流與固體邊界。一般而言,渠道入流邊界條件設定為單位寬度入流 量,渠道出流邊界條件則採用水位高程設定。在固體邊界處,沿固體邊 界法線方向採不透水邊界條件;而沿固體邊界切線方向可分為滑移與非

(32)

滑移邊界條件。 垂直水理模式考量渠道入流、渠道出流、自由液面及底床邊界條 件。在渠道入流及渠道出流處假設均勻流邊界條件。自由液面採風剪力 邊界條件 s V u d       、 V s v d       ;而底床則採用底床剪力邊界條件 b V u d       、 V b v d       。 2.4.2 動床部分 本研究底床濃度為 所在高程,如圖 2.4 所示。計算過程中底床 若為濃度捲昇過程,底床平衡濃度 計算結果會大於上層水體濃度, 因此直接採用下式: (2.43) 若為濃度落淤過程, 會小於上層水體濃度,可判斷為落淤情況

(van Rijn,1985),此時則採用 Neumann 邊界條件:

(2.44) = 0.015D Tk1.5 Dk0.3 (Van Rijn 1984b) ;a砂丘高度之一半;      顆粒 躍起之速度。  如此模式便可合理模擬捲昇與落淤兩種情況的濃度剖面。 b z , a e C , a a e CC , a e C 0 a C    , a e C lk wk

(33)

圖2.1 正交曲線座標轉換示意圖

(34)

圖2.3 水深方向流速剖面示意圖。(u 為實際流速剖面;u 為水深平均流 速;u為流速微變量) 圖2.4 作用層示意圖(Spasojevic 1988)

u

u

u u u

 

u

z

(35)

第三章 數值架構

3.1 水理部分 水理部分之整體架構基於水平垂直分離演算概念,首先由水平二維模 式求解水面高程以及水深平均流速,其中之有效剪力項依靠三維流場計 算而得到。三維流場係經流速差異量方程式求解後,與二維模式疊代收 斂。以下分別就水理部分之求解架構與數值差分方法做說明。 3.1.1 求解架構 水深平均控制方程式 本研究基於隱式雙階分割操作之觀念,將深度平均動量方程式分割

成二個時間步驟,延散步驟(advection and diffusion step)及傳播步驟

(propagation-step),分別求得 與 時刻之流場,並利用隱式數值 方 法 求 解 。 延 散 步 驟 求 解 移 流 項(advection terms)和擴散項(diffusion terms),傳播步驟求解壓力項、底床剪應力項和連續方程式。水理控制方 程式先就時間部分離散如下: 延散步驟: (3.1) 傳播步驟: (3.2) (3.3) 式中, 表示速度向量; 表示擴散及延散項; 表示 時 刻之未知變數; ; 表示 時刻之已知變數; 表示在 1/ 2 nn1 1 2 1 1 2 2 0 1 ( ) n n n n n V V V V T t            1 2 1 1 0 ( ) n n n b b V V g z D t D           1 0 n V     V T n1 (n 1) t 1 n n t tt    n n tn1/ 2

(36)

與 間之未知變數。 (3.1)~(3.3)的一般式可表示成: 延散步驟: (3.4) (3.5) 傳播步驟: (3.6) (3.7) 以及連續方程式: (n 1) t n t 1/2 2 1 2 1 2 1 2 1 2 0 n n f h h u u u u v u uv v t h h h h h h                   1/2 2 1 2 11 12 11 12 22 0 1 2 0 1 0 2 1 1 1 2 n h h h T T T T T Dh h Dh Dh                               1/2 2 1 2 1 1 2 0 n n f h h v v v v u v uv t h h h h                   1/2 1 2 1 22 12 22 12 11 0 1 2 0 2 0 1 1 1 1 2 n h h h T T T T T Dh h Dh Dh                               1 1/2 1 1 1 1 0 1 0 ( ) n n n n n b z D u u g g D D d d t h h                           



   1 1/2 2 1/2 2 1 1 2 1/2 1/2 2 2 0 0 0 ( ) ( ) 8 ( ) ( ) ( ) n n n n n f B B l n n n n n c u u v u u K D D D u v              1 1/2 1 1 1 2 0 2 0 ( ) n n n n n b z D v v g g D D d d t h h                           



   1 1/2 2 1/2 2 1 1 2 1/2 1/2 2 2 0 0 0 ( ) ( ) 8 ( ) ( ) ( ) n n n n n f B B l n n n n n c v u v V v K D D D u v                 

(37)

(3.8) 針對 時刻的水深值 做線性化處理,且僅保留一階項,(3.8) 式可改寫成: (3.9) 式中, ; ; ; ; ; ; ; 。 垂直部分 在垂直水理模式方面,考量數值計算的穩定性,一樣以隱示法求 解,如此將(2.11)與(2.12)式改寫為: ξ 方向: 1 1 1 2 1 1 2 1 2 n n V u h h u u u u uv uv u M t h h h h h d                                 (3.10) 1 1 1 1 2 ( 2 ) ( 1 ) 0 n n n n D D h h h uD h vD t               1 n(Dn1) 1 1 2 1 1 1 2 2 2 ( ) ( ) 0 n n D D D D h h D D t                         2 1 1 n h g t D C h     1 2 2 n h g t D C h     1 1 2 2 2 1 1 0 1 0 n n n b z h u tgh D D tgh d d C C h C h                               



  1 1 1 1 1 2 2 0 2 0 n n n b z h v tgh D D tgh d d C C h C h                               



  1 1 n D   2 2Dn 1/2 1/2 2 2 2 1/2 1/2 2 2 0 0 0 ( ) ( ) 1 8 ( ) ( ) ( ) n n f B B l n n n n n t C u v t t K C D D D u v                    1 n n D DD   

(38)

η 方向: 1 2 2 2 2 1 2 1 2 n n V v h h v v v v uv uv v M t h z h h h h d                                      (3.11) 式中, 1 1 1 1 1 2 2 2 1 2 ( ) n s b u c h u u u u v u v u v u uv M f v d h h h h h h h                                          2 2 2 1 2 1 2 1 2 ( ) n V h h vv v w Horizontal Diffusion in h h h h dh                    (3.12) 2 2 2 2 2 1 1 1 1 2 ( ) n s b v c h v v v v u v u v u v uv M f u d h h h h h h h                                           2 1 1 1 2 1 2 2 2 ( ) n V h h uu u w Horizontal Diffusion in h h h h dh                    (3.13) 式中, 為垂向黏滯係數。 V 3.1.2 數值方法 在數值差分方法選用的考量上,利用顯式數值方法求解時,演算時 間間隔將會受到很大的限制,在模擬天然明渠水流問題時將耗費冗長的 演算時間與龐大的電腦計算量,在應用上有其困難存在。為解決這個問 題,本研究採用隱式數值方法求解。 本模式採用控制體積(control volume)法的觀念來離散控制方程 式,控制體積法的基本概念如圖 3.2 所示,其中(a)圖為實際區域,(b) 圖為計算區域,E、W、N、S 表相鄰格點,e、w、n、s 表控制面。模 式計算之變數則放置在交錯網格(staggered grid)上,如圖 3.3。在控制方 程式中,除了移流項採用一階精度混合型上風法(hybrid scheme) (D.B. Spalding 1972)差分外,所有空間差分均採用二階精度的中央差分法。另 外,時間項則採用簡單的前向差分方法。 中央差分法可表示成

(39)

e w p Ψ Ψ Ψ Δ          (3.14) n s p Ψ Ψ Ψ Δ          (3.15) 式中, Ψe0.5 ( ΨEΨP) 0.5 ( Ψi1,jΨi j, ); , 1, 0.5 ( ) 0.5 ( ) w P W i j i j Ψ ΨΨ ΨΨΨn 0.5 ( ΨNΨP) 0.5 ( Ψi j, 1Ψi j, ); , , 1 0.5 ( ) 0.5 ( ) s P S i j i j Ψ ΨΨ ΨΨ Ψ 可表為 uvh1h2dzszb ij分別代表水平格網上任一點之縱向及橫向位置。 混合型上風法為上風法(upwind scheme)與中央差分法組合而成,當 移流效應重要時,採用上風法;移流效應不重要時,則採用中央差分法。 至於移流效應重要性的判斷,則採用格網雷諾數(mesh Reynolds number)

Rx、Ry作為判斷的因子,當|Rx|或|Ry|大於 2 時,代表移流效應重要,差 分方法採用能反映方向性的上風法;|Rx|或|Ry|小於等於 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 Φ hhΔ  Δ                                  (3.16) , 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 Φ hhΔ  Δ                              (3.17) 其中 0 2 1 2 1 2 x x x x R R R           ; 2 0 1 2 1 2 y y y y R R R           (3.18) 上列諸式中, , 1, / i j n i j x u h Δ R     ; , 2, / i j n i j y v h Δ R     ; 為流體動力黏滯係數 (dynamic viscosity);Φ可表成uv 。 垂直水理模式亦採用控制體積法的觀念離散控制方程式,如圖 3.4

(40)

所示,E、W、N、S、T、B 為相鄰格點,e、w、n、s、t、b 為控制面。 模式計算之變數亦放置在交錯網格(staggered grid)上,如圖 3.3。(3.12) 及(3.13)等式左邊採用 Crank-Nicolson method,其時間差分為二階精度; 而等式右邊MuMv的空間差分則採用中央差分法。 3.2 沉滓運移部分 沉滓運移首先求解懸浮載濃度之移流傳輸,求得懸浮載濃度的空間 分布後,再行求解作用層連續方程式與底床連續方程式,下列分述其求 解架構與數值方法。 3.2.1 求解架構 質量傳輸方程式 與流速差異量相同,基於尺度之考量,濃度的移流傳輸採用垂直水 平切割的操作方法,在垂直方向的微分項採用隱式法,水平方向的微分 項採用疊代收斂,首先將式(2.21)分為濃度對垂直方向微分項與濃度對 水平微分項( ),並對時間離散,如下所示: (3.19) 其中 包含其他水平方向的移流擴散項,上標 表示疊代收斂, 其包含的物理項如下: (3.20) 作用層與底床連續方程式 作用層連續方程式與底床連續方程式所求解的粒徑比例與底床沖 c M 1 1 1 2 1 1 1 * 2 2 2 1 n n n n n n n f f n sV sV c w w C C C C C C C M t D D D D                          * n c M n* 2 2 2 2 2 1 2 1 1 1 1 sH c sH sH h h h C C C M h h h h h                            2 1 1 1 2 1 2 2 2 2 sH sH sH h h h C C C C C uh vh h h h                                

(41)

淤量具相似的時間尺度,且作用層方程式的守恆性,必須依賴作用層連 續方程式與底床連續方程式在同一時間間距下同時求解。本研究採用式 (2.39)之作用層源維持各代表粒徑之比例在作用層中的和為 100%,以結 合演算法(Spasojevic and Holley,1990)概念求解如下:

在主要格點P(圖 3.1)之離散方程式共有 m+1 個,若以向量表示則 如下式: (3.22) 若整體運算的沉滓可以 m 個代表粒徑做表示,則包含 m 個 需求 解粒徑之間的比例,剩餘一個變數為底床高程。以更簡單的表示方法如: (3.23) M 表示代表粒徑總數目,m 表示代表粒徑編號。將式(3.19)與式 (3.20)分別以向量形式表示則如下式: (3.24) (3.25) 式(2.22) 與 式 (2.3) 均 為 一 非 線 性 代 數 式 , 加 以 線 性 化 後 , 利 用 Newton-Raphson 法疊代求解: (3.26) (3.27) 式(3.23)與式(3.24)中, 為 Jacobian 係數矩陣中之列向量; 為前一次疊代未知向量; 為疊代過程的向量修正值。解得向量 修正值 後,可得新的 向量: 1 1 1 ( , ,..., ,..., ) n n b m M s   z      m  1 1 1 ( , ) n m s   s s m1,M 1 1( ) 0 n F s  1 1( ) 0 n m F s   m1,M 1 1 1 [ F ] s F s(l n ) s         1 1 1 [ m ] (l n ) m F s F s s           m1,M [F/s] 1 l ns  s s  l1sn1

(42)

(3.28) 當達到收斂條件時,疊代得以結束;一般而言,疊代2 至 3 次即可 達到誤差小於 之精度。 3.2.2 數值差分 採用混合型上風法與中央差分法: (3.29) 其中 (3.30) 與式(3.18)相似,其中 為垂向之格網雷諾數,用 以判斷垂向濃度移流所適用的差分方法。 式(3.26)以中央差分法離散: 1 1 1 lsn l ns   s 8 1 10  1 1 1 1 1 (1 ) (1 ) n n n n n n P P T P P B z z T P P B C C C C C C t                1 1 1 1 (1 ) (1 ) n n n n fh T P P B z z T P P B w C C C C D                       1 1 1 1 1 1 2 2 2 ( ) ( ) ( ) ( ) 1 ( )( ) ( ) n n n n n n sV T sV B sV T B T P P B T B T B T P P B P B C C C C C C D D                                            1 1 ( fh T) ( fh B) * n n P c T B w w C M D            0 1 1 z        2 2 2 z z z R R R     , , / ( / 0) n z i j k R  z   * * * 2 2 * * 2 2 1 2 1 1 1 1 1 ( ) ( ) n n n n n c sH e w sH n s e w n s h h h h M C C C C h hh hh h                                   * * * * * 2 2 1 1 ( n n )[( ) ( ) ] ( n 2 n n ) e w sH e sH w sH E P W h h C C C C C h   h       

數據

圖 2.1 正交曲線座標轉換示意圖
圖 2.3 水深方向流速剖面示意圖。( u 為實際流速剖面; u 為水深平均流 速; u 為流速微變量)  圖 2.4 作用層示意圖(Spasojevic 1988) uuu u u uz
圖 3.1 模式計算流程圖
圖 3.2 水平二維模式控制體積法示意圖。 i 、 j 分別代表水平格網上任一 點之縱向及橫向位置;(a)實際區域;(b)計算區域
+7

參考文獻

相關文件

This study first uses the nine indicators of current domestic green architecture to examine those items needed to be considered in the air force base.. Then this study,

This study hopes to confirm the training effect of training courses, raise the understanding and cognition of people on local environment, achieve the result of declaring

Results of this study show: (1) involvement has a positive effect on destination image, and groups with high involvement have a higher sense of identification with the “extent

explore the effect of the different analysis method on the determination of the safety factor for slope stability under stormy

Investigating the effect of learning method and motivation on learning performance in a business simulation system context: An experimental study. Four steps to

Therefore, this study intends to combine the discussion method with the interactive response system of Zuvio IRS for flipped teaching in the course "Introduction to

For meeting the requirement in aggregate usage and maximizing utilization ratio of E&R resources obtained from the site, this study tried to adjust the optimum

In Type I, one performed comparison of 2D&3D model, impact of intrinsic stress of devices, package thermal effect, and transistor location effect.. On the other hand,