浸沒邊界法在二維淺水波流場模擬之應用探討

國

立

交

通

大

學

土木工程學系

碩

士

論

文

浸沒邊界法在二維淺水波流場模擬之應用探討

Use of the Immersed Boundary Method for

2D Shallow Water Flow Simulation

研 究 生：黃振家

指導教授：楊錦釧 博士

指導教授：

謝德勇 博士

浸沒邊界法在二維淺水波流場模擬之應用探討

Use of the Immersed Boundary Method for

2D Shallow Water Flow Simulation

研 究 生 ：黃振家 Student: Chen-Chia Huang

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

謝德勇 Te-Yung Hsieh

國 立 交 通 大 學

土 木 工 程 研 究 所

碩 士 論 文

A Thesis Submitted to Civil Engineering

College of Engineering

National Chiao Tung University

in Partial Fulfillment of the Requirements

for the Degree of Master

in

Civil Engineering

June 2009

Hsinchu, Taiwan, Republic of China

I

浸沒邊界法在二維淺水波流場模擬之應用探討

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

摘 要

本研究在一既已發展之水深平均二維淺水波模式中，植入浸沒邊 界法(Immersed Boundary Method)處理內部邊界問題，用以模擬流場 內部具非浸沒式水工構造物(橋墩、丁壩)的案例。本研究探討浸沒邊 界法中 2 點、4 點狄拉克脈衝函數(Dirac Delta Function)於淺水波模擬 應用之適用性，並分析內部邊界標記點(IB markers)於格點佈置之關係 性。模式之應用測試主要乃選取流場中具橋墩、丁壩的定床及動床實 驗案例進行模擬，除比對流場之正確性外，並就局部沖淤的相對位置 及刷深深度之比較結果，說明模式之限制條件。

II

Use of the Immersed Boundary Method for

2D Shallow Water Flow Simulation

Student：Chen-Chia Huang Advisors：Jinn-Chuang Yang Te-Yung Hsieh

Department of Civil Engineering

National Chiao-Tung University

ABSTRACT

The present study embeds the Immersed Boundary Method (IB method) in the 2D depth-averaged shallow water flow model being developed to solve the interior boundary problem in open channel. Two point and four point Dirac delta functions used in IB method were studied to examine their suitable equation. In addition, the effect caused by various ratios of IB marker’s gap to grid size were also analyzed. Two experimental cases, including cylinder pier and spur dike, were studied herein to demonstrate the practical applicability and limitation of IB method in fixed-bed and movable-bed open channel flow simulation.

Keywords:Immersed Boundary Method, shallow water flow, depth averaged, two-dimensional model

III

目錄

摘 要 ... I ABSTRACT ... II 目錄 ... III 表目錄 ... VI 圖目錄 ... VII 符號表 ... X 第一章 緒論... 1 1.1 研究動機與方向 ... 1 1.2 文獻回顧 ... 1 1.3 研究目的 ... 3 1.4 研究流程 ... 4 1.5 章節介紹 ... 4 第二章 理論基礎 ... 6 2.1 水理部分 ... 6 2.1.1 水理控制方程式 ... 6 2.1.2 輔助關係式 ... 7 2.2 沉滓運移部分 ... 8 2.2.1 控制方程式 ... 8

IV 2.2.2 輔助關係式 ... 9 2.3 外部邊界條件處理 ... 11 2.3.1 水理邊界條件 ... 12 2.3.2 沉滓運移邊界條件 ... 12 2.4 內部邊界條件處理 ... 12 第三章 數值架構 ... 19 3.1 水理部份 ... 19 3.1.1 分割操作趨近法 ... 19 3.1.2 數值差分式 ... 20 3.2 沉滓運移部分 ... 22 第四章 模式參數適用性分析 ... 25 4.1 橋墩案例 ... 25 4.1.1 案例介紹 ... 25 4.1.2 計算時距 ... 26 4.1.3 標記點間距與格網間距比例 ... 27 4.1.4 標記點間距非均勻性的影響 ... 27 4.1.5 質量守恆分析 ... 28 4.2 丁壩案例 ... 28 4.2.1 案例介紹 ... 28

V 4.2.2 計算時距 ... 29 4.2.3 標記點間距與格網間距比例 ... 29 4.2.4 標記點間距非均勻性的影響 ... 30 4.2.5 質量守恆分析 ... 30 4.3 參數分析彙整說明 ... 31 第五章 定床模擬驗證及動床案例測試 ... 50 5.1 定床模擬驗證 ... 50 5.1.1 橋墩案例 ... 50 5.1.2 丁壩案例 ... 51 5.2 動床案例測試 ... 51 5.2.1 橋墩案例 ... 51 5.2.2 丁壩案例 ... 53 第六章 結論與建議 ... 63 6.1 結論 ... 63 6.2 建議 ... 64 參考文獻 ... 65

VI

表目錄

表 4.1 不同Δt對橋墩流場均方根誤差之關係表(2 點法) ... 33  表 4.2 不同Δt對橋墩流場均方根誤差之關係表(4 點法) ... 33  表 4.3 不同Φ對橋墩流場均方根誤差之關係表(2 點法) ... 34  表 4.4 不同Φ對橋墩流場均方根誤差之關係表(4 點法) ... 34  表 4.5 非均勻分布標記點對橋墩流場均方根誤差之關係表 ... 35  表 4.6 流量比較表(橋墩) ... 35  表 4.7 不同Δt對丁壩流場均方根誤差之關係表(2 點法) ... 36  表 4.8 不同Δt對丁壩流場均方根誤差之關係表(4 點法) ... 36  表 4.9 不同Φ對丁壩流場均方根誤差之關係表(2 點法) ... 37  表 4.10 不同Φ對丁壩流場均方根誤差之關係表(4 點法) ... 37  表 4.11 非均勻分布標記點對丁壩流場均方根誤差之影響表 ... 38  表 4.12 流量比較表(丁壩) ... 39  表 4.13 模擬採用之參數一覽表 ... 39  表 4.14 CPU TIME ... 39 

VII

圖目錄

圖 1.1 研究流程圖 ... 5  圖 2.1 作用層示意圖 ... 16  圖 2.2 浸沒邊界示意圖 ... 17  圖 2.3(a) 狄拉克 2 點脈衝函數法示意圖 ... 18  圖 2.3(b) 狄拉克 4 點脈衝函數法示意圖 ... 18  圖 3.1 控制體積法示意圖 (a)實際區域；(b)計算區域 ... 24  圖 4.1 橋墩實驗渠道配置圖 ... 40  圖 4.2 橋墩實驗格網配置圖 ... 40  圖 4.3 最大可蘭數與橋墩案例實驗值之均方根誤差關係圖 ... 41  圖 4.4 Φ與橋墩案例實驗值之均方根誤差關係圖 ... 41  圖 4.5(a) 案例一標記點配置示意圖(橋墩) ... 42  圖 4.5(b) 案例二標記點配置示意圖(橋墩) ... 42  圖 4.5(c) 案例三標記點配置示意圖(橋墩) ... 43  圖 4.5(d) 案例四標記點配置示意圖(橋墩) ... 43  圖 4.6 流量比較斷面示意圖(橋墩) ... 44  圖 4.7 無因次縱斷面流量比較圖(橋墩) ... 44  圖 4.8 丁壩實驗渠道配置圖 ... 45  圖 4.9 丁壩案例流速比較範圍示意圖 ... 45 

VIII 圖 4.10 最大可蘭數與丁壩案例實驗值之均方根誤差關係圖 ... 46  圖 4.11 Φ與丁壩案例實驗值之均方根誤差關係圖 ... 46  圖 4.12(a) 案例一標記點配置示意圖(丁壩) ... 47  圖 4.12(b) 案例二標記點配置示意圖(丁壩) ... 47  圖 4.12(c) 案例三標記點配置示意圖(丁壩) ... 48  圖 4.12(d) 案例四標記點配置示意圖(丁壩) ... 48  圖 4.13 流量比較斷面示意圖(丁壩) ... 49  圖 4.14 無因次縱斷面流量比較圖(丁壩) ... 49  圖 5.1(a) 橋墩案例無因次縱斷面流速比較圖 ... 54  圖 5.1(b) 橋墩案例無因次橫斷面流速比較圖 ... 54  圖 5.2(a) 丁壩案例無因次縱斷面流速比較圖(y/b=1) ... 55  圖 5.2(b) 丁壩案例無因次縱斷面流速比較圖(y/b=1.5) ... 55  圖 5.2(c) 丁壩案例無因次縱斷面流速比較圖(y/b=2) ... 56  圖 5.2(d) 丁壩案例無因次縱斷面流速比較圖(y/b=3) ... 56  圖 5.2(e) 丁壩案例無因次縱斷面流速比較圖(y/b=4) ... 57  圖 5.3 y/b=1 處之流場示意圖 ... 57  圖 5.4 動床案例實驗量測範圍示意圖(橋墩) ... 58  圖 5.5 水流流經橋墩之流況 ... 58  圖 5.6 橋墩底床高程模擬結果圖(第 5 小時) ... 59 

IX 圖 5.7 墩前底床高程比較圖(第 3.83 小時) ... 59  圖 5.8 墩後底床高程比較圖(第 4.17 小時) ... 60  圖 5.9 墩側底床高程比較圖(第 4 小時) ... 60  圖 5.10 動床案例實驗量測範圍示意圖(丁壩) ... 61  圖 5.11 丁壩底床高程模擬結果圖(第 4.5 小時) ... 61  圖 5.12 丁壩沖刷坑深度比較圖(第 4.5 小時) ... 62 

X

符號表

a = 沙丘高度之一半； b = 丁壩長度； C = 濃度； C = 摩擦係數； f Ck = 顆粒 k 之深度平均濃度； c = Chezy 係數； c1 = 顆粒蔡司係數； Dk = 顆粒 k 之粒徑； Dm = 不產生移動的最小顆粒粒徑； D*k = 無因次顆粒粒徑； d = 水深； 5 0 d = 中值粒徑； E _{= 糙度因子；} Em _{= 作用層厚度；} f x( )_i = 尤拉反饋力； ( ) m i F X = 拉格蘭齊虛擬力； g = 重力加速度； h1 = ξ方向轉換係數； h2 = η方向轉換係數； k = von Karman’s 係數； M _{= 標記點總數量；} p = 孔隙率； qbi = i方向某一粒徑之底床載通量；

XI r = 橋墩半徑； r = 懸浮載源； Sf = 作用層源； s = 砂比重； T11,T12,T22 = 有效剪應力項； Tk = 輸送參數； t = 時間； u =

ξ

方向速度； uw = 近固體邊界的水深平均速度； u* _{= 剪力速度；} u*c = 臨界剪力速度； v = η方向速度； wfk = 顆粒 k 之沉降速度； wlk = 顆粒 k 之躍起速度； yw _{= 固體邊界與鄰近固體邊界格點的距離；} zb = 底床高程； zs = 水面高程； Δs = 標記點間距； Δξ 、Δη _{= 格網點間距；} Δh = 局部加密處格網點間距； Δ z = 沖刷深度； ξ、η = 平面上兩正交曲線座標方向； Φ = 標記點間距與格網間距比例； Φ = 狄拉克脈衝函數； δh = 隱藏因子；

XII β = 粒徑百分比； (βs k) = 母層中某一顆粒 k 之粒徑百分比 ρ = 流體密度； ρs = 泥砂密度； τbi=底床剪應力在

ξ

與

η

方向之分量;

υ

l = 層流黏滯係數； t

υ

= 亂流黏滯係數； ε1、ε2 = ξ、η方向之亂流傳輸係數。 上標 n = n tΔ 時刻之已知變數； n+1 = (n+ Δ1) t時刻之未知變數； 1 2 n+ = (n+ Δ1) t與n tΔ 間之未知變數； ( ) = 時間平均; ( ) = 水深平均； ('_{) = 時間平均瞬時擾動量。} 下標 s = 變數在水面的值； b = 變數在底床的值； i = 尤拉格網點 ξ 、η方向之變數； im = 拉格蘭齊標記點 ξ 、η方向之變數。

1

第一章 緒論

1.1 研究動機與方向

台灣為一狹長的島國，並有中央山脈橫亙其中，因此河川特性多 為短小流急，且為防洪、水資源規劃及交通需求，沿河道設置諸多水 工構造物(如橋墩、丁壩等)，使得局部流速加大，造成局部沖刷的現 象，因此改變原有河性，影響河川的穩定發展。 近年來電腦科技蓬勃發展，為探討前述河性及水理之改變，甚多 學者致力於數值模式的發展，數值模式在某些假設條件下，大致上皆 能獲得合理的模擬結果，但於處理跨河設施及水工構造物(內部邊界) 時，模式之發展常遭遇不同之問題，針對流場內部邊界問題，於淺水 波模式之發展中主要不外乎採用處理不規則邊界之非結構性格網，或 貼壁座標轉換之處理。 本研究為擴充原已採用有限差分法之水深平均二維淺水波模 式，使其能夠合宜的模擬河川中具水工構造物的流場問題，於原模式 中植入浸沒邊界法(Immersed Boundary Method)，探討其於明渠流應 用之適用性及限制性。

1.2 文獻回顧

就明渠水流而言，水理模式之發展已漸趨成熟，但在處理流場中 具內部邊界(如橋墩、丁壩)的案例時，仍具相當程度的困難性。 在數值方法處理方面，有限元素法及有限體積法可利用非結構性 格網(unstructured grids)來貼合不規則內部邊界，在遇內部邊界時將格 網繞過不作計算。李明旭(1991)應用有限元素法，以 FESWMS – 2DH 探討橋墩設置對河川局部流場之影響，戴源宏(2003)應用有限體積法 模擬流場通過兩不同位置圓柱，流場所產生的變化。但非結構性格網

2 生成較為複雜，且計算時需耗費大量的記憶體及計算時間。 以有限差分法求解明渠流場具內部邊界問題時，內部邊界需以乾 床處理，謝德勇(2003)在判斷乾濕交界面時給予一小水深作為乾床的 分界，此法處理內部邊界時將該邊界底床拉高，使邊界內部計算點形 成乾點，此處理方式會造成數值不穩定的問題，導致程式不易收斂。 整體而言，淺水波模式在進行內部邊界處理時，以非結構性格網 處理會有格網生成不易的困難，以有限差分法配合結構性格網進行處 理時，需針對內部邊界額外進行複雜的數值處理，在乾濕交界處易造 成質量不守恆，常會面臨模式收斂不易的限制，使模式的實用性受到 較大的挑戰。 近年來發展更多的數值方法，如 Peskin(1972) 提出浸沒邊界法的 觀念，其為模擬血液(流體)與心臟瓣膜(彈性邊界)間的互制作用而提 出的數值方法。浸沒邊界法主要的精神為將固體邊界視為流體的一部 分作計算，以反饋力的概念來重新分配邊界周圍的物理量。 使用浸沒邊界法處理具內部邊界之流場問題時，是以座標點的方 式，將所有拉格蘭齊標記點(markers in Lagrangian coordinate)連結， 用 以 描 繪 內 部 邊 界 之 幾 何 形 狀 ， 並 以 尤 拉 格 網 點 (grid points in Eulerian coordinate)來描述整體流場分佈範圍。因此流場格網不需以 格網生成技術使其貼合內部邊界，在處理不規則內部邊界問題時變得 相對容易。且此法之計算格網點不需避開障礙物(即障礙物內可以存 在計算格網點)，故不需判斷乾濕分界點，不會有以乾床處理內部邊 界時，需判斷邊界點上流速為零，而造成數值不穩定的問題。 近期許多研究於 Peskin(1972)浸沒邊界法的架構下，將之擴展到 流體力學的計算。Lai(2000)延續 Peskin(1972)的觀念，利用面積權重 的觀念，以 4 點狄拉克脈衝函數法(Dirac Delta Function)(以下簡稱 4

3 點 法 ) 求 出 浸 沒 邊 界 上 每 一 標 記 點 的 拉 格 蘭 齊 虛 擬 力 (Lagrangian Force)給予流場的尤拉反饋力(Eulerian Force)，使其可以重新分配流 場中的流速分佈。並以流體通過圓柱體的案例進行流體力學計算。Su et al.(2007)以 2 點狄拉克脈衝函數法(以下簡稱 2 點法)對四種案例(包 括橋墩)進行流體力學計算，在低雷諾數(Re<150)的條件下進行模 擬。林慕塵(2008)在計算邊界力時應用 volume of solid 的概念，取得 該網格所佔的內部邊界權重值，作為計算內部邊界與流體間虛擬力 (virtual force)的權重值。 由上述回顧可知，在此之前，浸沒邊界法主要是應用於低雷諾數 之流體力學計算，本研究擬採用 2 點法與 4 點法進行探討，其最大差 異為在計算內部邊界與其鄰近格網點物理量傳遞情形時，計算格網點 範圍不同(詳述於 2.4 節)，因此，若標記點間距過小，在使用 2 點法 計算時，易因計算範圍較小，而產生奇異點。綜觀上述，若能適當配 置標記點，並將浸沒邊界法用於處理內部邊界問題，實具其可行性。 因此 Chung (2008)將浸沒邊界法應用於明渠流場的計算，用以求解流 場中具內部邊界的問題，克服以往在處理明渠流場內部邊界問題時遭 遇之困難，且在比對整體流場之分佈可獲致合理結果，顯示浸沒邊界 法應用於天然河川的模擬是可行的。

1.3 研究目的

Chung(2008)將浸沒邊界法植入謝德勇(2003)發展之 RESED2D 模 式的水理子模式中，本研究之主要目的為延續 Chung(2008)之研究， 考量自由液面及底床摩擦邊界，對明渠紊流流場具內部邊界的問題進 行水理參數適用性分析，藉此了解各項影響因子對流場分布的影響程 度。並以較佳之水理參數設定對其它橋墩、丁壩的定床案例進行驗 證。動床部分沿用謝德勇(2003)之研究成果，以橋墩、丁壩之動床案

4 例進行測試模擬，對模式在動床模擬時所造成之局部沖刷相對位置及 深度進行分析，並說明二維淺水波模式在動床模擬之限制性。

1.4 研究流程

本章節主要將本研究之研究目的以流程圖方式呈現，用以說明研 究之整體流程及架構，圖 1.1 為研究流程圖。

1.5 章節介紹

前四節已闡述本研究之研究動機與方向、文獻回顧、研究目的及 研究流程，以下將扼要說明本論文各章節之內容。 第一章為緒論，即闡述本研究之研究動機與方向、文獻回顧、研 究目的及研究流程，並介紹全文之架構。 第二章為理論基礎，本章將會說明以正交曲線座標系統所發展之 數值模式，其水理、沉滓運移控制方程式及邊界條件的處理。 第三章為數值架構，有關水理、沉滓運移控制方程及內部邊界處 理所採用的數值方法，均將於本章加以說明。 第四章為模式參數適用性分析，以橋墩及丁壩之定床實驗案例進 行參數適用性分析。 第五章為定床模擬驗證及動床案例測試，選用橋墩及丁壩其它定 床實驗案例進行水理參數驗證分析。最後以橋墩及丁壩之動床實驗案 例進行動床模式測試，針對演算結果進行解析。 第六章為結論與建議，針對本研究成果作綜合性之歸納說明，並 對不盡完善或往後尚須改進之處提出建議。

5

6

第二章 理論基礎

透過座標系統轉換將控制方程式轉換為正交曲線座標系統方程 式，再將此控制方程式作時間平均及水深平均後，即可推得水深平均 二維正交曲線座標模式所需之控制方程式。其中水理控制方程包含水 理連續方程式及動量方程式；沉滓運移部份則將輸砂通量分離為非均 勻質之底床載與懸浮載，所以其控制方程包含某一粒徑之懸浮載質量 守恆方程式、某一粒徑於作用層(active layer)之質量守恆方程式及整 體河床沉滓運移之質量守恆方程式。茲將水理及沉滓運移部分之理論 基礎敘述如下：

2.1 水理部分

2.1.1 水理控制方程式 為適度簡化複雜的控制方程式，需對數學模式作若干假設，分別 為(Ⅰ)不可壓縮牛頓流體(incompressible Newtonian fluid)；(Ⅱ)靜水壓 分布；(Ⅲ)忽略風剪力；(Ⅳ)忽略科氏力。則水深平均二維正交曲線 座標水理控制方程式可表示為： (1)水流連續方程式 1 2 2 1 d h h + (h ud)+ (h vd) t ξ η ∂ ∂ ∂ ∂ ∂ ∂ (2.1) (2)水流動量方程式

ξ

方向： 2 1 2 1 2 1 2 1 2 1 h 1 h u u u v u uv v t h h h h h h ∂ ∂ ∂ ∂ ∂ ∂ + ∂ξ + ∂η + ∂η − ∂ξ 2 11 1 12 1 1 2 1 2 1 1 ( _b ) ( ) ( ) g z d h T h T h h h d h h d ∂ ∂ ∂ ∂ξ ρ ∂ξ ρ ∂η = − + + + 1 1 2 12 22 1 2 1 2 1 h 1 h b T T h h d h h d d τ ∂ ∂ ρ ∂η ρ ∂ξ ρ + − −

7 2 11 2 11 1 12 1 12 1 2 1 ( ) s ( ) b ( ) s ( ) b s b s b z z z z h h h h h h d τ τ τ τ ρ ξ ξ η η ⎡ ∂ ∂ ∂ ∂ ⎤ + _⎢− + − + _⎥ ∂ ∂ ∂ ∂ ⎣ ⎦ (2.2)

η

方向： 2 2 1 1 2 1 2 1 2 1 h 1 h v u v v v uv u t h h h h h h ∂ ∂ ∂ ∂ ∂ ∂ + ∂ξ + ∂η + ∂ξ − ∂η 2 12 1 22 2 1 2 1 2 1 1 ( _b ) ( ) ( ) g z d h T h T h h h d h h d ∂ ∂ ∂ ∂η ρ ∂ξ ρ ∂η = − + + + 2 1 2 11 12 1 2 1 2 1 h 1 h b T T h h d h h d d τ ∂ ∂ ρ ∂η ρ ∂ξ ρ − + − 2 12 2 12 1 22 1 22 1 2 1 ( ) s ( ) b ( ) s ( ) b s b s b z z z z h h h h h h d τ τ τ τ ρ ξ ξ η η ⎡ ∂ ∂ ∂ ∂ ⎤ + _⎢− + − + _⎥ ∂ ∂ ∂ ∂ ⎣ ⎦ (2.3) 式中， 2 2 11 [ 11 ' ( ) ] s b z z T =

∫

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

∫

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

∫

τ −ρu v −ρ u−u v−v dz (2.6) 以上諸式中，ξ、η = 平面上兩正交曲線座標方向，其中ξ為縱方向， η為側方向；h₁ = ξ方向轉換係數；h₂ = η_{方向轉換係數；u =} ξ方 向速度； v = η方向速度；ρ = 流體密度；d _{= 水深；}g _{= 重力} 加速度；t _{= 時間；}z_b = 底床高程；z_s = 水面高程； i b τ = 底床剪 應力在ξ與η方向之分量；( ) = 時間平均；( ) = 水深平均；('_{) = 時} 間平均瞬時擾動量；下標s、b _{分別代表變數在水面與底床的值；} 11

T ,T₁₂,T₂₂ = 有效剪應力項(effective stress term)，包含層流剪應力、亂 流剪應力與延散剪應力(dispersion stresses)，但由於本研究所選用之案 例皆為直線渠道，故在此不考慮延散剪應力項。

2.1.2 輔助關係式

8

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

1 2 2 1/ 2 ( ) b Cf u u v τ = ρ + (2.7) 2 2 2 1/ 2 ( ) b Cf v u v τ = ρ + (2.8) 式中， 2 f C = g c = 摩擦係數；c = Chezy 係數 (2)層流與亂流剪應力 採用 Boussinesq 之渦流黏性理論，層流與亂流剪應力可合併表 示為： 2 11 1 1 1 2 1 ' 2 u v h u h h h τ _υ ρ ξ η ⎡ ∂ ∂ ⎤ − = _⎢ + _⎥ ∂ ∂ ⎣ ⎦ _(2.9) 2 22 2 2 1 2 1 ' 2 v u h v h h h τ _υ ρ η ξ ⎡ ∂ ∂ ⎤ − = _⎢ + _⎥ ∂ ∂ ⎣ ⎦ _(2.10) 12 2 1 1 2 2 1 ' ' 2 h v h u u v h h h h τ _υ ρ ξ η ⎡ ∂ ⎛ ⎞ ∂ ⎛ ⎞⎤ − = _{⎢ ⎜ ⎟}+ _{⎜ ⎟⎥} ∂ _{⎝ ⎠} ∂ _{⎝ ⎠} ⎣ ⎦ _(2.11) 式中，

υ υ υ

= +_{l t}；

υ

_{l = 層流黏滯係數；}

υ

_t = 亂流黏滯係數 = ku d_* / 6(Falcon 1979)； 1/ 2 * ( b/ ) u = τ ρ =剪力速度；k= von Karman’s 係數(約等於 4.0)。

2.2 沉滓運移部分

2.2.1 控制方程式 輸砂通量可分離為非均勻質之懸浮載與底床載，因此，沉滓運移 部分控制方程式包含某一粒徑懸浮載質量守恆方程式、某一粒徑於作 用層(active layer)之質量守恆方程式及整體河床沉滓運移之質量守恆 方程式 (1)某一粒徑懸浮載之質量守恆方程式 1 1 1 2 1 2 1 2 2 1 2 2 1 1 ( h ) h C u C v C C C S t h ξ h η h h d ξ ε h ξ h h d η ε h η ρd ⎡ ⎤ ∂ ∂ ∂ ∂ ∂ ∂ ∂ + + = + ⎢ ⎥+ ∂ ∂ ∂ ∂ ∂ ∂ _⎢⎣ ∂ _⎥⎦ (2.12)

9 (2)某一粒徑於作用層之質量守恆方程式 1 2 1 1 1 2 ( ) (1 ) m ( ) ( ) 0 s b b f E p h h h q h q S S t ∂ β ∂ ∂ ρ ∂ ∂ξ ∂η − + + + − = (2.13) (3)整體河床沉滓運移之質量守恆方程式 0 ) ( ) ( ) 1 ( 2 1 1 2 2 1 ⎥= ⎦ ⎤ ⎢ ⎣ ⎡ ₊ ∂ ∂ + ∂ ∂ + ∂ ∂ −

∑

h q hq S t z h h p b _{b b} s _{ξ η} ρ (2.14) 以上諸式中，C = 濃度；ε₁、ε₂ = ξ、η方向之亂流傳輸係數； ρ_s = 泥砂密度；β = 粒徑百分比；p = 孔隙率；Em= 作用層厚度；qbi = i 方向某一粒徑之底床載通量；S = 懸浮載源；Sf = 作用層源

(active-layer floor source)。

2.2.2 輔助關係式 (1)亂流傳輸 ξ 方向之亂流傳輸係數採用 Elder(1979)之經驗公式，可表示為 ε1=5.93u d* , ε2 =0.23u d* _(2.15) (2)底床載通量 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.16) 上式中， 1 3 * 2 ( 1) k k s g D D ν − ⎡ ⎤ = _{⎢ ⎥} ⎣ ⎦ = 無因次顆粒粒徑；Dk = 顆粒 k 之 粒徑； 2 2 * * 2 * ( ) ( ) c k k c k u u T u − = = 輸送參數；u*c = 1 u g c = 臨界剪力速度 ； 1 90 12 18 log( ) 3 d c D = = 顆粒蔡司係數； s s ρ ρ = = 砂比重。 假設底床載運移僅發生在作用層內，並考慮較小粒徑在水體中 會形成懸浮載及較細顆粒可能會被隱藏在較粗顆粒之間較不易 被水流帶動之機制，則某一粒徑之底床載通量可進一步修正為：

10 qbi =ζ βhk kqbi( )D k (2.17) 式中，ζh = 隱藏因子(hiding factor)(Karim 等 1987)。 (3)懸浮載源 S 懸浮載源係屬於河床質載之一部分，是由懸浮質向下之通量與 底床亂流剪力作用產生河床質向上之通量交互作用之結果。懸 浮質下移到河床表面，主要是受到重力的影響。對某一粒徑 k 之懸浮質而言，其向下通量可表為： Sdk = −ρw Cfk dk (2.18) 式中， * [3.25 0.55ln( )] k fk d k w C C u κ = + _{(Lin 1984)；}C_k：顆粒 k 之深度 平均濃度；w_fk = 顆粒 k 之沉降速度 其中wfk則在不同之粒徑下採用不同之經驗(VanRijn,1984b)，如 下三式所示 v gD s w k fk ) 1 ( 18 1 − = _；D_k <100μm _(2.19) ⎪⎭ ⎪ ⎬ ⎫ ⎪⎩ ⎪ ⎨ ⎧ − ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ₊ − =10 1 0.01( 1) 1 5 . 0 2 3 v gD s D v w k k fk ；100μm≤Dk ≤1000μm (2.20)

[

(

)

]

5 . 0 1 1 . 1 _k fk s gD w = − _；D_k ≥1000μm _(2.21) 另一方面，底床載成為懸浮質，主要受到底床之亂流作用所造 成。對某一粒徑 k 而言，底床載向上之通量可表為： Sek =ρ βwlk kCek (2.22) 式中， _0.3 * 5 . 1 015 . 0 k k D T a D C k k e = (Van Rijn,1984b)；a：沙丘高度之一半； w_lk = 顆粒 k 之躍起速度。 顆粒躍起速度定義為河床質發生跳躍(saltation)時，離開底床之 瞬間垂直速度。本研究採用 Hu and Hui(1996)提出之經驗公式

11 * 3.2 4.5log 12 3.1 1.2 lk w u − Θ Θ < ⎧ = ⎨ _{Θ >} ⎩ (2.23) 式中， ( ) b s gDk τ ρ ρ Θ = ₋ 。 故由(2.18)及(2.22)式知，對某一粒徑 k 之懸浮載源可表為： Sk =ρ(wlkβkCek −w Cfk dk) (2.24) (4)作用層源 Sf  作用層示意圖如圖 2.1 所示。作用層源由於母層(active stratum) 頂面之升降而產生，當其下降時。 Sf s(1 p) t[( s k) (zb Em)] ∂ ρ β ∂ = − − − (2.25) 式中，(β_{s k}) = 母層中某一顆粒 k 之粒徑百分比。 如母層之厚度增加，及其頂面上升時，(2.25)式中之(β_{s k}) 則改為 β_k。 (5)作用層厚度 Em

沖刷現象發生時，根據 Bennet and Nordin(1977)之研究，Em可

以下式表示： Em Cem(zbn 1 zbn) + = − − (2.26) 式中，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.27) 式中，Dm：不產生移動的最小顆粒粒徑。 另外，作用層在淤積期間可定義為： Emn 1 Emn (zbn 1 zbn) + _{= +} + ₋ (2.28)

12

2.3 外部邊界條件處理

2.3.1 水理邊界條件 本研究外部邊界分為開放式邊界與固體邊界。開放式邊界條件主 要為： 1. 渠道入流邊界條件設定為單位寬度入流量。 2. 渠道出流邊界條件則採用水位高程設定。

在固體邊界處應用側壁理論(Law of the wall)，設定靠近固體邊界 的邊界條件為： * 1 ln( ) w u Ey U k + = (2.29) 式中，u_w = 近固體邊界的水深平均速度；E = 糙度因子 = 9.0(Lien 等 1999b)；y+ =y U_{w *}/υ，y_w = 固體邊界與鄰近固體邊界格點的距離。 2.3.2 沉滓運移邊界條件 沉滓運移求解的未知變數為深度平均懸浮載濃度C、粒徑百分比 β 及河床高程z_b。 C邊界條件的處理，將渠道入流邊界條件設定為入流濃度分佈， 渠道出流處邊界條件為∂ ∂ =C ξ 0，在固體邊界處，邊界條件設定為 0 C η ∂ ∂ = 。β與zb部分，入流邊界條件可設定為∂ ∂ =β ξ 0與∂ ∂ =zb ξ 0， 在固體邊界處，其邊界條件則設定為∂ ∂ =β η 0與∂ ∂ =z_b η 0。

2.4 內部邊界條件處理

本研究為探討非浸沒式水工構造物(如橋墩、丁壩)對河川水流流 場的影響，在處理內部邊界時以反饋力的概念求解內部邊界周圍之流 速，故將浸沒邊界法植入水理子模式之動量方程式中，模擬水流與內 部邊界互制的現象，在求解沉滓運移部分，由於動床子模式並非採用 動量方程式，而是以質量守恆方程式求解，故動床子模式沿用謝德勇

13 (2003)之概念，不加以修改。浸沒邊界法的精神是藉由流場中的速 度，配合邊界不滑移(no-slip)的條件求出一虛擬力，進而解出浸沒邊 界給予流場的反饋力，並將兩方向的虛擬反饋力 f

( ) ( )

ξ、f η 置 入 (2.2)、(2.3)中右項，並改寫為式(2.30)、(2.31)，藉此推求新的流速分 佈情形。

ξ

方向： 2 1 2 1 2 1 2 1 2 1 h 1 h u u u v u uv v t h h h h h h ∂ ∂ ∂ ∂ ∂ ∂ + ∂ξ + ∂η + ∂η − ∂ξ 1 1 2 2 11 1 2 1 12 1 1 ( ) g ( _b ) ( ) ( ) f z d h T h T h h h d h h d ∂ ∂ ∂ ξ ∂ξ ρ ∂ξ ρ ∂η = − + + + 1 1 2 12 22 1 2 1 2 1 h 1 h b T T h h d h h d d τ ∂ ∂ ρ ∂η ρ ∂ξ ρ + − − 2 11 2 11 1 12 1 12 1 2 1 ( ) s ( ) b ( ) s ( ) b s b s b z z z z h h h h h h d τ τ τ τ ρ ξ ξ η η ⎡ ∂ ∂ ∂ ∂ ⎤ + _⎢− + − + _⎥ ∂ ∂ ∂ ∂ ⎣ ⎦ (2.30)

η

方向： 2 2 1 1 2 1 2 1 2 1 h 1 h v u v v v uv u t h h h h h h ∂ ∂ ∂ ∂ ∂ ∂ + ∂ξ + ∂η + ∂ξ − ∂η 2 1 2 2 12 1 2 1 22 1 1 ( ) g ( _b ) ( ) ( ) f z d h T h T h h h d h h d ∂ ∂ ∂ η ∂η ρ ∂ξ ρ ∂η = − + + + 2 1 2 11 12 1 2 1 2 1 h 1 h b T T h h d h h d d τ ∂ ∂ ρ ∂η ρ ∂ξ ρ − + − 2 11 2 11 1 12 1 12 1 2 1 ( ) s ( ) b ( ) s ( ) b s b s b z z z z h h h h h h d τ τ τ τ ρ ξ ξ η η ⎡ ∂ ∂ ∂ ∂ ⎤ + _⎢− + − + _⎥ ∂ ∂ ∂ ∂ ⎣ ⎦ (2.31) 如圖 2.2，圓點分佈所連結區域表示為內部邊界，稱之為拉格蘭 齊標記點，且標記點與分佈於流場範圍的尤拉格網點並不需要重合。 而虛擬力與反饋力的關係可表示為：

14 1 ( ) ( ) ( ) m m M i i h i i k f x F X δ x X s = =

∑

− Δ (2.32)

f x( )i 為 Eulerian Force，F X( im)為 Lagrangian Force，下標im為ξ方

向與η方向的變數，M 為標記點總數量，Δs為標記點的間距，δ_h =狄

拉克脈衝函數(Lai and Peskin 2000)，以(2.33)式求解標記點與格網點 之間物理量的權重關係。 δh(xi−Xi)=d xh( ξ −X d yξ) h( η −Yη) (2.33) , , (1 ) 0 i i i i h r h h r h d otherwise ⎧ − ≤ ⎪ = ⎨ ⎪⎩ (2.34) 2 , 2 , , 1 (3 2 1 4 4( ) ) 8 1 (3 2 7 12 4( ) ) 2 8 0 i i i i i i h i i i i i r r r r h h h h h r r r d r h h h h h otherwise ⎧ − + + − ≤ ⎪ ⎪ ⎪⎪ =⎨ − − − + − ≤ ⎪ ⎪ ⎪ ⎪⎩ (2.35) hi為尤拉正交格網點間距，(2.34)式為 2 點法，計算標記點與格 網之間的權重值時，在ξ 及η 方向往標記點左右及上下各抓取一個格 網點的距離做計算，如圖(2.3)，(2.35)式為 4 點法，計算標記點與格 網之間的權重值時，在ξ 及η 方向往標記點左右及上下各抓取二個格 網點的距離做計算，如圖(2.4)。 Eulerian 流速u x( )_i 與 Lagrangian 流速U X( _i)之關係如下式： ( im) ( )i h( i im) 1 2 i U X =

∑

u x δ x −X Δ Δh h (2.36) 以下列步驟求解 Lagrangian force(參考 Su et al. 2007)

* 1 2 ( ) ( ) ( ) ( ) m m m i i i h i i x U X U X f x x X h h t δ − = − Δ Δ Δ

∑

(2.37) 將(2.32)代入， * 1 2 1 ( ) ( ) ( ) ( ) ( ) m m m m m M i i i h i i h i i x k U X U X F X x X s x X h h t = δ δ − _{⎡ ⎤} = _⎢ − _⎥Δ − Δ Δ Δ

∑ ∑

⎣ ⎦ (2.38)

15 移項， * 1 2 1 ( ) ( ) ( ) ( ) ( ) m m m m m M i i i i h i i i k x U X U X x X x X s h h F X t = δ − _{⎡ ⎤} = _⎢ − − Δ Δ Δ _⎥ Δ

∑ ∑

⎣ ⎦ (2.39)

以高斯 – 喬登法(Gauss – Jordan method)求解(2.39)式中所形成的帶

狀矩陣，分別為： 2 4 ⇒ ⎧ ⎨ _⇒ ⎩ 點法 三階帶狀矩陣 點法 五階帶狀矩陣，可得 Lagrangian force，並將 其代回(2.32)式，求得 Eulerian force。

16 △L ß ß_s E_m Zb 作用層 底床第一層 底床第二層 圖 2.1 作用層示意圖

17 ( , ) m i X x y 圖 2.2 浸沒邊界示意圖

18

圖 2.3(a) 狄拉克 2 點脈衝函數法示意圖

19

第三章 數值架構

3.1 水理部份

3.1.1 分割操作趨近法 本研究基於分割操作之觀念，將動量方程式分割成兩個步驟(延 散步驟、傳播步驟)，並利用隱式數值方法求解。最後再以顯示數值 方法處理浸沒邊界求解步驟。延散步驟求解移流項和擴散項，傳播步 驟求解壓力項、底床剪應力項和連續方程式。據此，水理控制方程可 改寫為： 延散步驟 1 1 1 2 2 1 2 ( ) n n n n V V V T t ρ + + + ∂ ⎛ ⎞ _{= − ⋅∇ + ∇⋅} ⎜ _∂ ⎟ ⎝ ⎠ (3.1) 傳播步驟 * * 1 2 ( ) n n n b b V V g z d t t d τ ρ + ∂ ∂ ⎛ ⎞ ₋⎛ ⎞ _{= − ∇ + −} ⎜ _∂ ⎟ ⎜ _∂ ⎟ ⎝ ⎠ ⎝ ⎠ (3.2) * 0 n V ∇⋅ = (3.3) 浸沒邊界求解步驟 * 1 n n n V V f t t + ∂ ∂ ⎛ ⎞ ₋⎛ ⎞ ₌ ⎜ _∂ ⎟ ⎜ _∂ ⎟ ⎝ ⎠ ⎝ ⎠ (3.4) 式中，V表示速度向量；T 表示擴散及延散項；n+1表示(n+ Δ1) t時刻 之未知變數； n1 n t t + t Δ = − ；n表示n tΔ 時刻之已知變數；n*表示在(n+ Δ1) t 與n tΔ 間之未知變數。 (3.1)~(3.3)的一般式可表示成： 延散步驟 2 1 2 1 2 1 2 1 h h u u u v u uv v t h ξ h η h h η ξ ⎡∂ ∂ ⎤ ∂ _{= −} ∂ ₋ ∂ _{− −} ⎢ ⎥ ∂ ∂ ∂ _⎣∂ ∂ _⎦

20 2 11 1 12 1 2 12 22 1 2 1 2 1 2 1 2 ( ) ( ) 1 h T 1 h T 1 h 1 h T T h h d h h d h h d h h d ρ ξ ρ η ρ η ρ ξ ∂ ∂ ∂ ∂ + + + − ∂ ∂ ∂ ∂ 2 11 2 11 1 12 1 12 1 2 1 ( ) s ( ) b ( ) s ( ) b s b s b z z z z h h h h h h d τ τ τ τ ρ ξ ξ η η ⎡ ∂ ∂ ∂ ∂ ⎤ + _⎢− + − + _⎥ ∂ ∂ ∂ ∂ ⎣ ⎦ (3.5) 2 2 1 1 2 1 2 1 h h v u v v v uv u t h ξ h η h h ξ η ⎡∂ ∂ ⎤ ∂ _{= −} ∂ ₋ ∂ _{− +} ⎢ ⎥ ∂ ∂ ∂ _⎣∂ ∂ _⎦ 2 12 1 22 1 2 11 12 1 2 1 2 1 2 1 2 ( ) ( ) 1 h T 1 h T 1 h 1 h T T h h d h h d h h d h h d ρ ξ ρ η ρ η ρ ξ ∂ ∂ ∂ ∂ + + − + ∂ ∂ ∂ ∂ 2 12 2 12 1 22 1 22 1 2 1 ( ) s ( ) b ( ) s ( ) b s b s b z z z z h h h h h h d τ τ τ τ ρ ξ ξ η η ⎡ ∂ ∂ ∂ ∂ ⎤ + _⎢− + − + _⎥ ∂ ∂ ∂ ∂ ⎣ ⎦ (3.6) 傳播步驟 2 2 1 (zb d) C u uf v u g t h ξ d + ⎛∂ + ⎞ ∂ _{= − −} ⎜ ⎟ ∂ _⎝ ∂ _⎠ (3.7) 2 2 2 (z_b d) C v uf v v g t h η d + ⎛∂ + ⎞ ∂ = − _{⎜ ⎟}− ∂ _⎝ ∂ _⎠ (3.8) 和 2 1 1 2 ( ) ( ) 0 h ud h vd d h h t ξ η ∂ ∂ ∂ _{+ + =} ∂ ∂ ∂ (3.9) 針對n+1時的水深值(dn+1)做線性化處理，且僅保留一階項，(3.9)式可 改寫成 1 2 1 1 1 2 2 2 ( ) ( ) 0 d d d h h d d t ξ α ξ β γ η α η β γ ⎛ ⎞ ⎛ ⎞ ∂ ∂ ∂ Δ ∂ ∂ Δ + _⎜ + Δ + _⎟+ _⎜ + Δ + _⎟= ∂ ∂ _⎝ ∂ _⎠ ∂ _⎝ ∂ _⎠ (3.10) 式 中 ， 2 1 1 n h g t d C h_τ α = − Δ ； 1 2 1 2 2 1 1 [ ] n n n zb h h g t d u C_τ C h_τ ∂ ∂ β ∂ξ ∂ξ + + Δ = − + ; 1 1 n d γ = β ； 1 2 2 n h g t d C h_τ α = − Δ ; 12 1 1 1 2 2 [ ] n n n zb h h g t d v C_τ C h_τ ∂ ∂ β ∂η ∂η + + Δ = − + ; 2 2 n d γ =β ； 1 1 2 2 2 2 ( ) ( ) 1 n n f n C u v C t d τ + ₊ + = + Δ ； n1 n d d + d Δ = − 。 在求解延散及傳播步驟後，以(2.37)~(2.39)之求解流程求出(2.39)

21 式中之 ( i m ) F X ，並代回(2.32)式求解 f ，由於其求解流程相當複雜，故 在求解 f 時，以顯式數值方法求解 fn，再將之代入(3.4)式中，反覆疊 代求解n+1時刻的流速，直至收斂。 3.1.2 數值差分式 在數值差分方法選用的考量上，利用顯式數值方法求解時，計算 時距將會受到很大的限制，在模擬天然明渠水流問題時將耗費冗長的 演算時間與龐大的電腦計算量，在應用上有其困難存在。為解決這個 問題，本研究採用隱式數值方法求解。 本模式採用控制體積(control volume)法的觀念來離散控制方程 式，控制體積法的基本概念如圖 3.1 所示，其中(a)圖為實際區域，(b) 圖為計算區域，e、w、n、s 表控制面。在水理控制方程式中，除了 移流項採用一階精度混合型上風法(hybrid scheme)(Spalding 1972)差 分外，所有空間差分均採用二階精度的中央差分法。另外，時間項則 採用簡單的前項差分法。 中央差分法可表為 1 _{1 1} n _{n n} e w p ξ ξ + _{+ +} ⎛∂Ψ⎞ ₌Ψ − Ψ ⎜_∂ ⎟ _Δ ⎝ ⎠ (3.11) 1 _{1 1} n _{n n} n s p η η + _{+ +} ⎛∂Ψ⎞ ₌ Ψ − Ψ ⎜_∂ ⎟ _Δ ⎝ ⎠ (3.12) 式中， 1 1 1 1 1 1, , 0.5 ( ) 0.5 ( ) n n n n n e E P i j i j + + + + + + Ψ = ⋅ Ψ + Ψ = ⋅ Ψ + Ψ ； 1 1 1 1 1 , 1, 0.5 ( ) 0.5 ( ) n n n n n w P W i j i j + + + + + − Ψ = ⋅ Ψ + Ψ = ⋅ Ψ + Ψ ； 1 1 1 1 1 , 1 , 0.5 ( ) 0.5 ( ) n n n n n n N P i j i j + + + + + + Ψ = ⋅ Ψ + Ψ = ⋅ Ψ + Ψ ； 1 1 1 1 1 , , 1 0.5 ( ) 0.5 ( ) n n n n n s P S i j i j + + + + + − Ψ = ⋅ Ψ + Ψ = ⋅ Ψ + Ψ ； Ψ可表為 u ， v ，h1，h2，d，zs和zb。

22 混合型上風法為上風法(upwind scheme)與中央差分法組合而 成，當移流效應重要時，採用上風法；移流效應不重要時，則採用中 央差分法。至於移流效應重要性的判斷，則採用格網雷諾數(mesh Reynolds number) Rx、Ry做為判斷的因子，當 Rx 或 Ry 大於 2 時， 代表一流效應重要，差分方法採用能反應方向性的上風法；R_x 或 R_y 小於或等於 2 時，移流效應可視為不重要，差分方法採用中央差分 法。混合型上風法應用於本研究移流項的處理可表示成 , 1 1 1 1 1 , 1, , , 1, 1 1 0.5 (1 ) (1 ) i j n n n n n n n i j i j i j i j i j x x u u h ξ h α ξ α ξ + + + + + + − ⎡ ⎛Φ − Φ ⎞ ⎛Φ − Φ ⎞⎤ ⎛∂Φ ⎞_{= − + +} ⎢ ⎜ ⎟ ⎜ ⎟⎥ ⎜ _∂ ⎟ _{⎜ Δ ⎟ ⎜ Δ ⎟} ⎢ ⎥ ⎝ ⎠ ⎣ ⎝ ⎠ ⎝ ⎠⎦ (3.13) , 1 , , 1 , , , 1 2 2 0.5 (1 ) (1 ) i j n n n n n n n i j i j i j i j i j y y v v h η h α η α η + + − ⎡ ⎛Φ − Φ ⎞ ⎛Φ − Φ ⎞⎤ ⎛∂Φ ⎞_{= − + +} ⎢ ⎜ ⎟ ⎜ ⎟⎥ ⎜ _∂ ⎟ _{⎜ Δ ⎟ ⎜ Δ ⎟} ⎢ ⎥ ⎝ ⎠ ⎣ ⎝ ⎠ ⎝ ⎠⎦ (3.14) 其中 0 2 1 2 1 2 x x x x R R R α ⎧ _≤ ⎪ ⎪ =_⎨ > ⎪ ⎪− < − ⎩ ; 0 2 1 2 1 2 y y y y R R R α ⎧ _≤ ⎪ ⎪ =_⎨ > ⎪ ⎪− < − ⎩ (3.15) 上列諸式中， , 1 , n i j i j x u h R ξ μ ρ Δ = ； , 2 , n i j i j y v h R η μ ρ Δ = ；μ _{= 流體動力黏滯系數} (dynamic viscosity)；Φ可表成 u 或 v 。

3.2 沉滓運移部分

輸砂子模式控制方程中之(2.12)式為延散方程式，解此方程如同 解水流運動方程式，可分為下列二個步驟： (1)移流及反應(advection-reaction)步驟         1 2 ( C)a u C v C S t h h d ∂ ∂ ∂ ∂ + ∂ξ + ∂η = ρ (3.16) 式中( C)a t ∂ ∂ 為移流及懸浮載源 S  所引起的濃度變化。

23 (2)擴散(diffusion)步驟         1 1 1 2 2 1 2 2 1 1 ( C) ( C)a ( h C) h C t t h h d ξ h h h d η h ∂ ∂ _{ε ε} ∂ ∂ η η η η ⎡ ⎤ ∂ ∂ ∂ ∂ − = + ⎢ ⎥ ∂ ∂ ∂ _⎢⎣ ∂ _⎥⎦   (3.17)         經上述方法分解後之懸浮載方程式，在移流及反應步驟為一雙曲 線偏微分方程，而於擴散步驟則為一橢圓型偏微分方程。因此，先聯 立求解(2.13)、(2.14)、及(3.16)式；然後，所求得各變數之中間值再 與(3.17)式反覆疊代直至收斂為止。 本研究採用控制體積法差分離散控制方程，再利用 ADI 方法求 解。此種方法與水理延散步驟相似，在此不多贅述。

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

25

第四章 模式參數適用性分析

經由第二章與第三章在模式理論基礎與數值方法之說明介紹 後，本章節將選取定床案例(橋墩、丁壩)進行模式參數適用性分析， 用以檢定各項影響因子之敏感度及提供參數選定之依據。 本研究針對水流流經明渠流場具水工構造物的案例進行模擬，並 以不同的模擬參數配置進行分析。由於不同案例有其不同參數特性， 在此選取 2 點法及 4 點法，搭配不同之計算時距、不同標記點間距與 格網間距的比例、標記點間距非均勻性及質量守恆的影響，並對模擬 參數進行適用性分析。為展示各影響因子的敏感度，在此採用均方根 誤差(root mean square error)來表示對整體流場模擬值與實驗值的誤 差。 1 2 ( ) N i i i rms x y E N = ⎡ − ⎤ = ⎢ ⎥ ⎢ ⎥ ⎣ ⎦

∑

_(4.1) 式中xi為某一時刻的實驗值，yi為同一時刻的模擬值，N 表示比 較點的數量。E_rms為均方根誤差，E_rms越趨近於零，表示模擬結果與實 驗值的誤差越小。 由於本研究所選取之案例皆為直線渠道，因此在曲線座標系統中 之ξ、η方向可直接轉換為卡氏座標中的x 、y方向，因此以下案例皆 以x 、y 表示座標方向。

4.1 橋墩案例

4.1.1 案例介紹 本研究在橋墩案例部份，採用 Ahmed(1998)所做之定量流橋墩定 床實驗案例進行模擬。實驗渠道全長 18.4m，渠寬 1.22m，在距上游 13.1m 之渠道中心處設置一直徑 89mm 的非浸沒式圓柱橋墩，給定上

26 游均勻入流 0.065cms，水深為 0.182m，在橋墩部分設定為浸沒邊界， 圖 4.1 為橋墩實驗渠道配置圖。 在模擬 Ahmed(1998)的橋墩定床實驗案例時，固定格網數為 286 97× ，圖 4.2 為橋墩實驗格網佈置圖，在橋墩周圍鄰近實驗量測流 速處將格網局部加密為0.0075m×0.0075m之正方形格網，為節省計算時 間，在遠離橋墩處則將格網間距適度放大(最大處：0.5m×0.02m)。 實驗量測的範圍為距橋墩中心三倍半徑處縱斷面流速及距橋墩 中心上游五倍半徑處橫斷面流速，在與實驗值進行均方根誤差比對 時，範圍選取距橋墩中心三倍半徑處縱斷面流速除以距橋墩中心上游 十倍半徑處流速之無因次縱斷面流速，以及距橋墩中心上游五倍半徑 處橫斷面流速除以距橋墩中心上游十倍半徑處流速之無因次橫斷面 流速，將兩者之均方根誤差取平均值進行驗證。(斷面流速比較範圍 示意圖如圖 4.2)。 4.1.2 計算時距 在此僅改變不同計算時距Δt，探討Δt的不同對數值模擬產生的影 響為何？並設計 10 組不同Δt，利用流場中之最大可蘭數(Maximum Courant Number)與均方根誤差分別對 2 點法、4 點法進行測試，如表 4.1、4.2。 其中每一格網點之可蘭數(C_r)定義為： r t C u h Δ = ⋅ Δ (4.2) 式中， u =該格網點深度平均流速；Δt=計算時距；Δh=格網間距。 圖 4.3 為最大可蘭數(橫軸)，均方根誤差(縱軸)之關係圖，由圖 4.3 及表 4.1、4.2 可知，當 2 點法之最大可蘭數為 2.65(Δt為 0.06)時， 均方根誤差有最小值。4 點法之最大可蘭數為 6.87(Δt為 0.15)時，均 方根誤差有最小值。在均方根誤差最小的案例中，2 點法之平均可蘭

27 數為 0.98，4 點法之平均可蘭數為 2.46，故可得知，最大可蘭數反應 的是一局部流場行為。綜合上述，2 點法與 4 點法，在最大可蘭數為 2.65 及 6.87 時，為橋墩案例較佳之模擬參數。 4.1.3 標記點間距與格網間距比例 為瞭解格網與標記點配置對模擬結果產生的影響大小，在此定義 一無因次參數 s h Δ Φ = _Δ 為標記點間距與格網間距比例。並設計 8 組不 同Φ，設計的依據為固定格網點間距，並改變標記點間距，藉此比較 Φ與流速對 2 點法、4 點法之均方根誤差的影響關係，如表 4.3、4.4， 圖 4.4 為Φ(橫軸)，均方根誤差(縱軸)之關係圖，由表 4.3、4.4 及圖 4.4 可知，2 點法與 4 點法之均方根誤差相當接近，且以不同Φ進行模 擬時並無明顯差異，當 2 點法之Φ的比例為 2 時，有最小之均方根誤 差平均值為 0.02051，4 點法之Φ為 1.6，有最小之均方根誤差平均值 為 0.02100。綜合上述，2 點法與 4 點法，在Φ為 2 及 1.6 時，為橋墩 案例較佳之模擬參數。 4.1.4 標記點間距非均勻性的影響 由於 4.1.1.3 節為均勻分布標記點，故本章節設計四組案例，探 討非均勻分布之標記點對模擬結果產生之影響大小為何？ 在此選用兩種不同比例之Φ，分別為 1.6、2，搭配四組不同標記 點配置之案例進行模擬，如表 4.5，設計之依據分述如下： 案例一、二：以橋墩中心點為準，將橋墩縱切為正向水流沖擊的 沖擊面及非沖擊面。圖 4.5(a)、4.5(b)為案例一、二標記點配置示意圖。 案例三、四：以橋墩中心點為準，將橋墩橫切為對稱於渠道中心 線之右側面及左側面。圖 4.5(c)、4.5(d)為案例三、四標記點配置示意 圖。 非均勻分布標記點對流場之影響表如表 4.5，由表 4.5 可得知，

28 四個案例的均方根誤差皆與均勻分布標記點中Φ =2的均方根誤差較 為接近，由此可得一結論，如果標記點採兩種不同Φ的配置時，Φ較 大者對整體流場的均方根誤差影響將較為顯著。 4.1.5 質量守恆分析 本章節將對橋墩流場進行質量守恆的探討，圖 4.6 為流量比較斷 面之示意圖，比較範圍選取距橋墩上游 10 倍橋墩半徑處至距橋墩下 游 10 倍橋墩半徑處，觀察其流量變化，圖 4.7 為無因次縱斷面流量 比較圖，橫軸之無因次參數為 x 方向距橋墩中心點的距離與橋墩半徑 (r)之比值，以橋墩為中心點計，負值表上游區，正值表下游區。縱 軸之無因次參數為縱斷面流量Q除以上游未受橋墩擾動之流量Q₀，在 此 將 水 流 流 經 橋 墩 處 稱 為 束 縮 斷 面 ， 並 定 義 渠 道 的 束 縮 率 = 100% 構造物斷面積 渠道寬度 × ，此案例之最大束縮率為 7.3%。表 4.6 為上游入 流量與下游出流量及束縮斷面流量之比較關係表，由圖 4.7、表 4.6 可知，在流經束縮斷面時，流量有低估的現象，概因浸沒邊界法將內 部邊界視為流體的一部份作計算，因此在水流流經內部邊界時有些微 流量進入，造成流量在此處有低估的現象。

4.2 丁壩案例

4.2.1 案例介紹 在丁壩案例的模擬選用 Rajaratnam(1982)的丁壩定床實驗案例進 行模擬(A1 案例)，實驗渠道渠道全長 37m，渠寬 0.9m，丁壩長度設 定為 0.15m，給定上游均勻入流 0.045cms，水深為 0.189m，在丁壩 部分設定為浸沒邊界，圖 4.8 為丁壩實驗渠道配置圖。 在格網配置方面，根據鄭志斌(2000)研究，入、出流邊界需與丁 壩有足夠距離，以確保流況之完全發展，而不至受丁壩之擾動。因此 入流邊界設定在丁壩前方 25 倍丁壩長度處，出流邊界設定在丁壩後

29 方 36 倍丁壩長度處。模擬時固定格網數為246 76× ，在丁壩周圍鄰近 實驗量測流速處，將格網局部加密為0.01m×0.01m之正方形格網，並在 遠離丁壩處將格網適度放大(最大處：0.05m×0.02m)，以節省計算時間。 實驗量測流速處選取距丁壩上方 1、1.5、2、3、4 倍丁壩長度處 之無因次縱斷面流速，在與實驗值比較均方根誤差時，距丁壩上方 1、 1.5、2、3、4 倍丁壩長度處之無因次縱斷面流速除以上游未受擾動處 流速，取五者之均方根誤差平均值作為比較的基準。圖 4.9 為流速比 較範圍示意圖。 4.2.2 計算時距 在此僅改變不同計算時距Δt，設計 10 組不同Δt，如表 4.7、4.8， 利用最大可蘭數與均方根誤差分別對 2 點法、4 點法進行測試。 圖 4.10 為最大可蘭數(橫軸)，均方根誤差(縱軸)之關係圖，由表 4.7、4.8 及圖 4.10 可知，2 點法之Δt為 0.03 時，均方根誤差有最小值。 4 點法之Δt為 0.08 時，均方根誤差有最小值。在均方根誤差最小的案 例中，2 點法之平均可蘭數為 0.35，4 點法之平均可蘭數為 0.88，故 可得知，最大可蘭數反應的是一局部流場行為。綜合上述，2 點法與 4 點法，在最大可蘭數為 0.96 及 2.57 時為丁壩案例較佳之模擬參數。 4.2.3 標記點間距與格網間距比例 丁壩案例設計之依據為固定格網點間距，改變標記點間距，表 4.9、4.10 為 8 組不同Φ與流速對 2 點法及 4 點法之均方根誤差關係 表，圖 4.11 為橫軸Φ(橫軸)，均方根誤差(縱軸)之關係圖，由表 4.9、 4.10 及圖 4.11 可知，2 點法之Φ的比例為 1.2 時，有最小之均方根誤 差平均值為 0.07893。4 點法之Φ為 1.0 時，最小之均方根誤差平均值 為 0.07435。綜合上述，2 點法與 4 點法，在Φ為 1.2 以及 1.0，為丁 壩案例較佳之模擬參數。

30 4.2.4 標記點間距非均勻性的影響 在此選用兩種不同比例之Φ，分別為 1、1.5，搭配四組不同標記 點配置之案例進行模擬，如表 4.11，設計之依據分述如下： 案例一、二：以壩頭中心點為準，將丁壩縱切為正向水流沖擊的 沖擊面及非沖擊面。圖 4.12(a)、4.12(b)為案例一、二標記點配置示意 圖。 案例三、四：將丁壩橫切為上側及下側。圖 4.12(c)、4.12(d)為案 例三、四標記點配置示意圖。 非均勻分布標記點對流場之影響如表 4.11，由表 4.11 可得知， 四個案例的均方根誤差皆與均勻分布標記點中Φ =1.5的均方根誤差 較為接近，由此可得一結論，如果標記點採兩種不同Φ的配置時，Φ 較大者對整體流場的均方根誤差影響將較為顯著。 4.2.5 質量守恆分析 本章節將對丁壩流場進行質量守恆的探討，圖 4.13 為各比較斷 面流量之示意圖，比較範圍選取距丁壩上游 10 倍橋墩半徑處至距丁 壩下游 10 倍橋墩半徑處，觀察其流量變化，圖 4.14 為無因次縱斷面 流量比較圖，橫軸之無因次參數為 x 方向距丁壩中心點的距離與丁壩 長度(b)之比值，以丁壩為中心點計，負值表上游區，正值表下游區， 縱軸之無因次參數為縱斷面流量Q除以上游流量Q₀，在此將水流流經 丁壩處稱為束縮斷面，此案例之最大束縮率為 16.7%。表 4.12 為上游 入流量與下游出流量及最大束縮斷面流量之比較關係表，由圖 4.14、 表 4.12 可知，在流經丁壩周圍時，流量有低估的現象，概因浸沒邊 界法將內部邊界視為流體的一部份作計算，因此在水流流經內部邊界 時有些微流量進入，造成流量在此處有低估的現象。

31

4.3 參數分析彙整說明

經由上述分析，本節將彙整出影響模擬結果之參數，以較合適之 參數設定模擬，並於第五章展示模式對模擬流場具內部邊界案例之功 能。 在橋墩案例中，由表 4.1、4.2 及圖 4.3 可知，4 點法適用的最大 可蘭數範圍大於 2 點法。由表 4.3、4.4 及圖 4.4 可知，均勻分布之不 同Φ在以 2 點法及 4 點法模擬之均方根誤差相當接近，但 4 點法的模 擬結果則是更為一致，顯示在以不同Φ的配置時，4 點法穩定性比 2 點法好。對於標記點不均勻分布所造成的影響，由表 4.5 可知，2 點 法與 4 點法之均方根誤差都與Φ =2之均方根誤差較接近，由此可知， 如果標記點採兩種不同Φ的配置時，Φ較大者對整體流場的均方根誤 差影響將較為顯著。 在丁壩案例中，由表 4.7、4.8 及圖 4.10 可知，最大可蘭數適用 範圍在 4 點法中仍是優於 2 點法，由表 4.9、4.10 及圖 4.11 可知，2 點法Φ =2時已出現過大誤差，所以並未繼續進行測試，而 4 點法在 不同Φ時之模擬結果較為一致，顯示在以不同Φ的配置時，4 點法穩 定性較佳，且模擬結果之均方根誤差皆較 2 點法小。而標記點不均勻 分布結果與橋墩案例相同，在以兩種不同Φ的配置時，Φ較大者對整 體模擬結果影響將較為顯著，當模擬河道中具不規則形狀之內部邊 界，而無法以單一比例進行配置時，此結論能提供為標記點配置之參 考。 在質量守恆之分析，由表 4.6、4.12，圖 4.7、4.14 可知，橋墩及 丁壩案例在束縮斷面的流量出現誤差，主要原因是由於浸沒邊界法將 固體視為流體的一部分作計算，因此在內部邊界會有微小流量進入， 造成流量在束縮斷面有低估的現象。而橋墩和丁壩會造成不同的流量

32 誤差是由於兩案例在通過束縮斷面的通水面積不同，以兩案例之最大 束縮率與流量誤差比較，可知束縮率的大小為影響流量誤差之主要原 因。以浸沒邊界法進行模擬雖會造成質量不守恆的現象，但由 4.1.1 及 4.1.2 節的分析中得知，橋墩及丁壩案例與實驗值比對的均方根誤 差並未出現過大差距，顯示模式對模擬整體流場的流速分布仍具一定 程度之準確性，因此，流量在局部地區被低估的現象仍可接受。 綜觀上述，可知 2 點法與 4 點法在經參數適用性分析後，若以各 自較佳之水理模擬參數進行模擬，則兩者之誤差非常小，但經比較之 後，4 點法之均方根誤差仍略小於 2 點法，故本研究之模擬參數設定 擬採用 4 點法，橋墩與丁壩案例採用的各參數列於表 4.13。 最後以丁壩案例進行模擬，改變標記點數量的配置，探討模擬所

需之 CPU TIME。表 4.14 為無因次化的 CPU TIME，te為最長模擬時

間，t為每一案例之模擬時間，由表可知，不存在浸沒邊界的案例模

擬時間最短，而隨著標記點配置數量越多，所需的模擬時間隨之增 加，因此，標記點的數量和模擬時間的關係為正相關。

33 表 4.1 不同Δt對橋墩流場均方根誤差之關係表(2 點法) t Δ 最大可蘭數 縱斷面流速 均方根誤差 橫斷面流速 均分根誤差 均方根誤差 0.02 0.93 0.03361 0.01738 0.02549 0.04 1.85 0.02462 0.01860 0.02161 0.05 2.29 0.02229 0.01953 0.02091 0.06 2.65 0.02093 0.01995 0.02044 0.08 3.48 0.02085 0.02133 0.02109 0.1 4.32 0.02191 0.02245 0.02218 0.15 6.36 0.02572 0.02376 0.02474 0.2 8.55 0.02877 0.02624 0.02751 0.3 12.83 0.03179 0.02767 0.02973 0.4 17.16 0.03548 0.02694 0.03121 表 4.2 不同Δt對橋墩流場均方根誤差之關係表(4 點法) t Δ 最大可蘭數 縱斷面流速 均方根誤差 橫斷面流速 均分根誤差 均方根誤差 0.02 0.99 0.05115 0.02496 0.03806 0.04 1.98 0.04415 0.01775 0.03095 0.05 2.45 0.03964 0.01687 0.02826 0.06 2.92 0.03575 0.01664 0.02620 0.08 3.82 0.03151 0.01700 0.02426 0.1 4.68 0.02842 0.01758 0.02300 0.15 6.87 0.02198 0.01882 0.02040 0.2 8.66 0.02311 0.02006 0.02159 0.3 12.48 0.02378 0.02279 0.02328 0.4 17.62 0.02466 0.02324 0.02395

34 表 4.3 不同Φ對橋墩流場均方根誤差之關係表(2 點法) s h Δ Δ 縱斷面流速_{均方根誤差} 橫斷面流速_{均方根誤差} 均方根誤差 0.8 0.02892 0.01774 0.02333 1 0.02647 0.01818 0.02232 1.2 0.02632 0.01802 0.02217 1.4 0.02761 0.01768 0.02264 1.6 0.02301 0.01943 0.02122 1.8 0.02193 0.01936 0.02065 2 0.02055 0.02047 0.02051 2.2 0.02077 0.02161 0.02119 表 4.4 不同Φ對橋墩流場均方根誤差之關係表(4 點法) s h Δ Δ 縱斷面流速_{均方根誤差} 橫斷面流速_{均方根誤差} 均方根誤差 0.8 0.02400 0.01882 0.02141 1 0.02380 0.01897 0.02138 1.2 0.02396 0.01877 0.02136 1.4 0.02400 0.01875 0.02137 1.6 0.02297 0.01903 0.02100 1.8 0.02395 0.01878 0.02137 2 0.02400 0.01883 0.02141 2.2 0.02330 0.01902 0.02116

35 表 4.5 非均勻分布標記點對橋墩流場均方根誤差之關係表 案例名稱 標記點配置(Φ ) 均方根誤差 均勻分布標記點之均 方根誤差(1.6 與 2) 一(2 點法) 沖擊面 1.6 非沖擊面 2 0.02069 1.6 0.02122 2 0.02051 一(4 點法) 沖擊面 1.6 非沖擊面 2 0.02144 1.6 0.02100 2 0.02141 二(2 點法) 沖擊面 2 非沖擊面 1.6 0.02039 1.6 0.02122 2 0.02051 二(4 點法) 沖擊面 2 非沖擊面 1.6 0.02142 1.6 0.02100 2 0.02141 三(2 點法) 左側 1.6 右側 2 0.02028 1.6 0.02122 2 0.02051 三(4 點法) 左側 1.6 右側 2 0.02142 1.6 0.02100 2 0.02141 四(2 點法) 左側 2 右側 1.6 0.02055 1.6 0.02122 2 0.02051 四(4 點法) 左側 2 右側 1.6 0.02143 1.6 0.02100 2 0.02141 表 4.6 流量比較表(橋墩) 橋墩案例 上游-下游 上游-最大束縮斷面 2 點法 Q Δ _(cms) 0.0001 0.00195 2 點法 誤差百分比(%) 0.15 3 4 點法 Q Δ (cms) 0.00012 0.00176 4 點法 誤差百分比(%) 0.17 2.7

36 表 4.7 不同Δt對丁壩流場均方根誤差之關係表(2 點法) t Δ 最大可蘭 數 y/b=1 流 速均方根 誤差 y/b=1.5 流 速均方根 誤差 y/b=2 流 速均方根 誤差 y/b=3 流 速均方根 誤差 y/b=4 流 速均方根 誤差 均方根誤 差平均值 0.01 0.33 0.23385 0.04461 0.09179 0.07127 0.03619 0.09554 0.02 0.65 0.22810 0.05128 0.06395 0.04825 0.02893 0.08410 0.03 0.96 0.22217 0.04583 0.05741 0.04006 0.03861 0.08082 0.04 1.27 0.22006 0.04716 0.05262 0.04256 0.04532 0.08155 0.05 1.57 0.22337 0.03279 0.05293 0.04510 0.05363 0.08156 0.06 1.86 0.22719 0.03988 0.06395 0.04925 0.02993 0.08204 0.08 2.44 0.22146 0.03543 0.05671 0.04896 0.06217 0.08495 0.1 3.00 0.22498 0.02951 0.07754 0.07121 0.09072 0.09879 0.15 4.39 0.22062 0.03208 0.09326 0.08841 0.10824 0.10852 0.2 5.93 0.21144 0.03766 0.12264 0.07347 0.09047 0.10714 表 4.8 不同Δt對丁壩流場均方根誤差之關係表(4 點法) t Δ 最大可蘭 數 y/b=1 流 速均方根 誤差 y/b=1.5 流 速均方根 誤差 y/b=2 流 速均方根 誤差 y/b=3 流 速均方根 誤差 y/b=4 流 速均方根 誤差 均方根誤 差平均值 0.01 0.37 0.25443 0.06555 0.12172 0.10165 0.06114 0.12089 0.02 0.74 0.24321 0.07615 0.10017 0.08456 0.04493 0.09780 0.03 1.02 0.20557 0.06596 0.09525 0.07458 0.03693 0.09566 0.04 1.35 0.19852 0.06915 0.08126 0.07072 0.03568 0.09107 0.05 1.66 0.19410 0.06583 0.07035 0.06427 0.03017 0.08494 0.06 2.06 0.19091 0.06135 0.06580 0.05476 0.02653 0.07987 0.08 2.57 0.18798 0.04548 0.06226 0.04166 0.03436 0.07435 0.1 3.16 0.18409 0.04470 0.06290 0.04188 0.04728 0.07617 0.15 4.57 0.17832 0.03284 0.07389 0.05871 0.07413 0.08358 0.2 5.95 0.18021 0.03610 0.08488 0.07347 0.09047 0.09303

37 表 4.9 不同Φ對丁壩流場均方根誤差之關係表(2 點法) △s/△h y/b=1 流速 均方根誤差 y/b=1.5流速 均方根誤差 y/b=2 流速 均方根誤差 y/b=3 流速 均方根誤差 y/b=4 流速 均方根誤差 均方根誤差 平均值 1.0 0.22217 0.04583 0.05741 0.04006 0.03861 0.08082 1.1 0.22933 0.05104 0.06140 0.04732 0.02602 0.08302 1.2 0.21500 0.03976 0.05737 0.03932 0.04318 0.07893 1.4 0.22475 0.04665 0.05737 0.04053 0.03656 0.08117 1.5 0.22604 0.04904 0.05729 0.04135 0.03485 0.08171 1.6 0.22513 0.04285 0.05830 0.04207 0.04937 0.08355 1.8 0.22134 0.04069 0.05921 0.04428 0.05496 0.08410 2.0 0.22879 0.03203 0.10325 0.10651 0.12737 0.11959 表 4.10 不同Φ對丁壩流場均方根誤差之關係表(4 點法) △s/△h y/b=1 流速 均方根誤差 y/b=1.5流速 均方根誤差 y/b=2 流速 均方根誤差 y/b=3 流速 均方根誤差 y/b=4 流速 均方根誤差 均方根誤差 平均值 1.0 0.18798 0.04548 0.06226 0.04166 0.03436 0.07435 1.1 0.19884 0.05384 0.05717 0.04053 0.04004 0.07808 1.2 0.20071 0.04755 0.05624 0.04018 0.03948 0.07683 1.4 0.20119 0.04969 0.05560 0.04058 0.03809 0.07703 1.5 0.19828 0.0532 0.05659 0.04208 0.04687 0.07940 1.6 0.20074 0.04997 0.05527 0.04067 0.03882 0.07710 1.8 0.20399 0.04975 0.05555 0.04097 0.04078 0.07821 2.0 0.20056 0.04909 0.05590 0.04069 0.03968 0.07718