孔隙消能結構物之三維度流場數值研析

(1)

國立臺灣大學生物資源暨農學院生物環境系統工程學系碩士論文

Department of Bioenvironmental Systems Engineering College of Bioresources and Agriculture

National Taiwan University Master Thesis

孔隙消能結構物之三維度流場數值研析 Numerical Investigation of 3D Flows over

Energy-Dissipation Porous Media

張家銘

Chia-Ming Chang

指導教授：張倉榮教授 Advisor: Prof. Tsang-Jung Chang

中華民國 100 年 6 月 June, 2011

(2)

(3)

I

致謝

多根圓柱的研究總算告一段落，在這期間有許許多多的人成為我的幫助，就像柱子一樣支持著我。第一根要感謝的是我所信靠的神，

感謝祂帶我進入這間研究室，安排張倉榮老師當我的指導教授，多位照顧我的學長，以及一起奮鬥的同伴們，有趣的學弟妹們，在我最艱難的時候成為我的安息與力量。第二根要感謝的是我敬愛的老闆張倉榮老師，給我們自由的學習空間，適時的提攜、指點、督促與教導，

使我在研究上更能往前，幽默詼諧的個性也降低了學生對教授的距離感。第三根要感謝的是帶我研究的博班學長華哥，從該領域的入門基礎開始，研究的方法、作事的態度、團隊的合作，無一不讓我受益良多，雖然做研究時是個很嚴謹的學長，但是平常是相當平易近人的。

第四根要感謝的是夥伴加榮，從一進來就令人輕鬆的傢伙，考上高考又是老師所器重的愛將，相互的吐槽與熬夜寫論文、作研究，有你的陪伴使我走得更穩妥。還有其他的柱子：宏名學長、嘉和學長、龍哥、

利哥、均美學姐、家詮學長、菲比學長、浩維學長、丞瑋學長、小黑、

軍曹、渣渣維、史帥、小八、優哥、小風、讚讚、咩咩、195，我所親愛的家人、弟兄姊妹，我所喜愛的青少年們，若一一詳述，恐怕比論文本文還多…。

這篇論文的價值在於，它磨練了一個人，也把那個人連結到了其他的人，交織出一幅圖畫，嵌進了許多的回憶，這就像一個相框，又如一本回憶錄，記起曾經有過的年輕歲月…。最後將這篇論文完成的喜悅，歸給一切幫助過我的貴人！

(4)

II

摘要

本研究主要為探討三維自由液面流場受到筐網圓柱之影響。首先，以 k–ε 標準紊流模式搭配體積分率法進行單根筐網圓柱之流場模 擬，並利用孔隙介質流模式計算水體流經筐網圓柱時所受到的阻滯效應。經由模擬結果與物理試驗結果之比較，發現兩者趨勢頗為一致，

顯示數值模擬在單根筐網圓柱之流場變化上具有一定程度的準確性。

接著，在相同入流流況與相同圓柱直徑的情境下，本研究進行多組不同孔隙率的流場模擬，並定義臨界壓力去估算出不同孔隙率的筐網圓柱後方的保護範圍。我們將擁有最大保護範圍之孔隙率稱為最佳孔隙率，經由模擬結果發現最佳孔隙率會位於孔隙率0.43 左右。

最後，將單根筐網圓柱案例擴展到單排筐網圓柱之案例，我們進行單排實心圓柱與單排筐網圓柱之流場比較，以及不同根數的筐網圓柱對於流場之影響。由模擬結果得知，實心圓柱的阻水效應大，圓柱後方的流速分佈紊亂，相反的，筐網圓柱後方的流速分佈卻十分穩定，

顯示筐網圓柱比實心圓柱更能穩定流場。比較不同筐網圓柱根數的流場變化發現，根數越多之流場，後方的低速區範圍有明顯增大之趨勢。

關鍵字：筐網圓柱、最佳孔隙率、單排多根圓柱、k–ε 標準紊流模式、

體積分率法、孔隙介質。

(5)

III

Abstract

The main objective of this research is to investigate the effect of porous cylinders on flow patterns. The standard k–ε turbulence model and the volume of fluid method (VOF) are used to simulate 3-D free surface flows. The resistant force generated by porous cylinders is performed by porous media theory. The computed results are verified by the available experimental data and are in consistent with the experiment data. In addition, we are concerned with how a large area is protected behind the porous cylinder, so a critical pressure is defined to estimate the protection area. Through several kinds of porosity are used in simulations, it is found that there is the largest protection area when the porosity is around 0.43.

Thus, we call it as the optimal porosity in this study.

Finally, we apply our model to flows past two kinds of singled-row cylinders, porous cylinders and solid cylinders. From the numerical results, we can conclude that flows past porous cylinders are more stable than those past solid cylinders. The number of cylinders is in proportional to their protection area.

Keyword: Porous cylinder, Optimal porosity, Single-row cylinder, Standard k–ε turbulence model, Volume of Fluid, Porous media.

(6)

IV

目錄

致謝 ... I 摘要 ... II Abstract ... III 目錄 ... IV 表目錄 ... VI 圖目錄 ... VII 符號對照表 ... X

第一章緒論 ... 1

1.1 前言 ... 1

1.2 研究目的 ... 2

1.3 文獻回顧 ... 2

第二章理論模式 ... 5

2.1 流場控制方程式 ... 5

2.2 紊流模式 k–ε ... 6

2.3 邊牆函數 ... 7

2.4 孔隙介質流模式 ... 8

2.5 體積分率法 ... 9

第三章數值方法 ... 12

3.1 數值模擬架構 ... 12

3.2 有限體積法 ... 12

3.3 控制方程式的離散 ... 13

3.4 收斂條件 ... 14

(7)

V

3.5 邊界條件設定 ... 14

第四章模式驗證 ... 18

4.1 實心圓柱驗證 ... 18

4.2 孔隙圓柱驗證 ... 20

4.3 驗證結果與討論 ... 21

4.4 孔隙區域與孔隙跳躍之比較 ... 22

第五章模式應用 ... 42

5.1 最佳孔隙率的探討 ... 42

5.2 單排圓柱的探討 ... 44

5.2.1 不同根數之筐網圓柱的流場流速分佈比較 ... 45

5.2.2 筐網圓柱與實心圓柱的流場流速分佈比較 ... 46

第六章結論與建議 ... 64

6.1 結論 ... 64

6.2 建議 ... 65

參考文獻 ... 66

(8)

VI

表目錄

表3.1 鬆弛因子係數的設定 ... 17 表4.1 兩個驗證案例的水理參數 ... 24 表5.1 孔隙介質的參數設定 ... 49

(9)

VII

圖目錄

圖2.1 邊牆流場示意圖(節錄自 Fluent 6.1 manual ) ... 11

圖4.1Ahmed 等(1998)的試驗設計模型圖 ... 25

圖4.2 實心圓柱網格圖 ... 25

圖4.3 軸向速度沿 z 軸的變化 ... 26

圖4.4 垂向速度沿 z 軸的變化 ... 26

圖4.5 無因次化的速度強度沿 x 軸的變化 ... 27

圖4.6 無因次化的速度強度沿 y 軸的變化 ... 27

圖4.7 無因次化的速度強度沿 y 軸的變化 ... 27

圖4.8 石(2007)的試驗設計模型圖 ... 28

圖4.9 直徑 0.10 m 筐網圓柱網格圖 ... 28

圖4.10 直徑 0.16 m 筐網圓柱網格圖 ... 28

圖4.11 軸向速度沿 x 軸方向之變化（y = 0 m，z = 0.056 m） ... 29

(a)圓柱直徑 10 公分 (b)圓柱直徑 10 公分 ... 29

圖4.12 軸向速度沿 x 軸方向之變化（y = 0.15 m，z = 0.056 m） ... 30

圖4.13 軸向速度沿 x 軸方向之變化（y = 0.25 m，z = 0.056 m） ... 31

圖4.14 軸向速度沿 y 軸方向之變化（x = 0 m，z = 0.056 m） ... 32

圖4.15 軸向速度沿 y 軸方向之變化（x = 0.15 m，z = 0.056 m） ... 33

圖4.16 軸向速度沿 y 軸方向之變化（x = 0.25 m，z = 0.056 m） ... 34

(10)

VIII

圖4.17 軸向速度沿 z 軸方向之變化（x = 0.25 m，y = 0 m） ... 35

圖4.18 孔隙率 0.38 之 x–y 平面的流線圖 ... 36

(a)直徑 0.1 m (b)直徑 0.16 m ... 36

圖4.19 孔隙率 0.38 之 x–z 平面的流線圖 ... 37

(a)直徑 0.1 m (b)直徑 0.16 m ... 37

圖4.20 圓柱直徑 0.1 m 之 x–y 平面的流線圖... 38

(a)孔隙率 0.38 (b)孔隙率 0.62 ... 38

圖4.21 圓柱直徑 0.1 m 之 x–z 平面的流線圖 ... 39

(a)孔隙率 0.38 (b)孔隙率 0.62 ... 39

圖4.22 porous jump 與 porous zone ... 40

圖4.23 軸向速度沿 x 軸之變化(圓柱直徑 0.1 公尺，孔隙率 0.38) .... 40

圖4.24 試驗所使用的筐網圓柱圖 ... 41

圖5.1 最佳孔隙率案例使用的模型 ... 50

圖5.2 圓柱附近的網格示意圖 ... 50

圖5.3 不同孔隙率之相對低壓區的壓力分佈圖(z=0.032 m) ... 51

圖5.7 保護體積算法之示意圖 ... 55

圖5.8 不同孔隙率之自由液面圖 ... 56

圖5.9 不同孔隙率與影響區域範圍圖 ... 57

圖5.10 單排多根圓柱模擬用之模型 ... 58

(11)

IX

圖5.11 圓柱排列方式示意圖 ... 58

圖5.12 筐網圓柱的網格示意圖 ... 59

(a)單根 (b)兩根 (c)三根 (d)五根 ... 59

圖5.13 軸向速度沿 z 軸之變化(x=0.4m，y=0.2125m) ... 60

(a)筐網圓柱 (b)實心圓柱 ... 60

圖5.14 軸向速度沿 y 軸之變化(x=1m，z=0.064m) ... 61

(a)筐網圓柱 (b)實心圓柱 ... 61

圖5.15 筐網圓柱之 x-y 平面水平速度分佈圖(z=0.064m) ... 62

圖5.16 實心圓柱之 x-y 平面水平速度分佈圖(z=0.064m) ... 63

(12)

X

符號對照表

x 軸向座標 ^k_r 阻力係數

y 橫向座標 ^ 孔隙率

z 垂向座標 b 孔隙結構物厚度

U 軸向速度 F 體積分率

V 橫向速度 P^* 初始壓力值

W 垂向速度 U^* 初始軸向速度值

P 壓力 V^* 初始橫向速度值

 剪應力 W^* 初始垂向速度值

f 實體力 x 網格長度

U 平均速度 t 時間間距

u' 擾動速度 I 紊流強度

p' 擾動壓力 Re 流場特徵雷諾數

T 渦流係數 L 特徵長度

 渦流尺度 B 粗糙函數

k 紊流動能 ^K_s 實際的粗糙高度

 紊流動能消散率 Ks^ 無因次化後的粗糙高度

(13)

XI

c_ Prandtl 參數，為定值 0.09 Cs 粗糙常數

c1_ Prandtl 參數，為定值 1.44 ^U 流速強度 c2_ Prandtl 參數，為定值 1.92 ^r 圓柱半徑

ck Prandtl 參數，為定值 1 H 水深 c_ Prandtl 參數，為定值 1.3 D 圓柱直徑

 馮卡門常數，0.4~0.42

E 經驗常數，9.793

Up 在P 點之平均流速 Kp 在P 點之紊流動能 yp 邊牆到P 點的距離

 運動黏滯係數

U* 剪力速度

 滲透係數

C 慣性阻力因子

 密度

(14)

1

第一章緒論

1.1 前言

由於氣候的異常，豪大雨或颱風發生的機率增高，河川的水流量也亦趨不穩定，因此河川的整治與保護就顯得亦發重要。已往河川的保護工法有丁壩、蛇籠、箱籠、鼎塊、拋石及混凝土護岸等，不同的保護工法應用在不同的地形與環境，交錯搭配使用以達到保護河川的功效。目前河川的保護工法大約可以分為兩類：剛性保護工法與柔性保護工法。剛性工法指的是不透水結構物，其結構強度大保護能力也強，大水時不易被破壞，具有較高的安全性與耐用性。由於不透水的性質，對於導流和掛淤有相當好的效果，同時亦有造灘與護岸的功能，

並可束縮河槽，增加低水流槽水深，以利航運。在應用的地形方面適用於坡度陡，流勢較強、推移載較多的卵礫石河川。但是不透水結構物的投資成本相對較高，工程範圍較大，影響河道原有的通洪能力。

也因著水流無法通過結構物，故結構物所受之作用力較大，在其周圍較容易發生局部沖刷，影響結構體穩定性。對於原有河道形態與河川物理棲地也有著較大的衝擊，使河川的生物群聚有所改變，亦會阻止水中生物的縱向活動和水中物理及化學的傳輸，造成水中生態的負效應。

相對於剛性工法的透水性結構物之柔性保護工法，工程量相較於不透水結構物小，投資成本亦比較低。在工法應用的地形方面，透水性結構物的機動性高能適合不同地形變化，也比較不佔據河流的自然空間，對河道原有的通洪能力影響比不透水結構物為小。因著其可透水的性質，故結構物所受之作用力相對較小，周圍也較不容易產生局

(15)

2

部沖刷，結構物本身的安全性較高。透水性結構物也較易隨周邊外在條件調整，比不透水結構物容易與周邊環境融合，對生態環境的衝擊也比較小。透水性結構物也較不受原有河道排洪寬度、整治寬度等河道參數的影響，同時具有降低紊流強度、紊流動能、雷諾剪應力和渦度等作用，進而減少下游沖刷與破壞。

近年來有學者利用透水性結構物的概念，發展透水性的筐網圓柱，用作橋墩的保護工以及河川的導流工；2006 年洪勝榮等更利用筐網圓柱群，於八掌溪作河道導流工程，用以保護河道的凹岸，發現有泥砂淤積且穩定河川底床的功效(許，2010)。

1.2 研究目的

目前對於筐網圓柱的研究多以現地實驗或是室內的水工試驗為主，需要繁複的步驟，通常耗時、耗費也耗人力；而隨著數值模擬計算方法的成熟，電腦計算效率的提升，用模擬搭配實驗是一種輔相成的方式，且模擬比試驗更具有易調整的特性，可以重複多組相似的模擬，亦可針對實際情況作邊界以及參數的修改；同時實驗不易量測的地方或物理量，亦可藉由數值模擬來補足(黃，2007)。

因此，本研究將以數值模擬真實流場，配合室內水工實驗以調整成適合的參數，進一步探討筐網圓柱對週遭流場的影響，並與不透水的實心圓柱作比較，以分析可透水筐網圓柱的特性，以期將來能成為現地試驗或規畫的一個參考依據。

1.3 文獻回顧

Ahmed 和 Rajaratnam(1998)在室內水槽進行單根實心圓柱體的流場試驗，底床型態分為動床與定床。文獻中針對各位置的流速以及柱

(16)

3

前壓力，進行分析並探討平滑底床與粗糙底床等情形對流況的影響。

Salaheldin 等(2004)根據 Ahmed 和 Rajaratnam (1998)的水槽試驗，利用數值模擬軟體 Fluent 成功模擬出水槽內的流況，其以不同型態的 k–ε 紊流模組搭配VOF 模擬水槽流況，其中以標準 k–ε 紊流模式與 RNG k–ε 紊流模式對受實心圓柱影響的流場有較為準確的呈現。

Wilson(1985)發表，因孔隙結構物的幾何難以描述，故在流體流經一孔隙結構物區域時，在那維爾–史托克方程式(Navier–Stokes equation) 中加入一動量損失項，以呈現孔隙結構物所產生的阻水效應。Jue(2003) 以有限元素法模擬受孔隙方塊影響之流場的變化，並以不同的達西數、雷諾數以及孔隙率作分析，發現在低雷諾數時(R_e 250)，相較於孔隙率，達西數為流場的主要控制因子。Sobera(2008)以 RNG k–ε模式求解三維杜雷諾平均的那維爾–史托克方程式，探討一實心圓柱在外圍包上一層孔隙介質影響流場的情形，並以直接數值解進行模式的驗證，發現該模式非常適合此種情況的流場 (  R_e 80000) 。 Bhattacharyya(2006)則是以數值方法模擬一孔隙圓柱周圍的流場變化，探討雷諾數、達西數以及孔隙率對流場流況之影響，以阻力係數、

分離角度、回復長度、流線圖和速度圖來作分析。阻力係數隨著雷諾數的增加和達西數的減少而降低。分離角度和尾流長度也隨著達西數的增加而減少。

Santiago 等在 2007 年以實驗以及數值模擬的方法，來探討風場流過風柵後流況的改變，文獻中使用三種形式的k–ε紊流模式，模擬結果與試驗值相當穩和，並且該文獻試驗了不同孔隙率的風柵，發現最佳孔隙率為 0.35。Dong 等(2007)發表的文獻說明孔隙率為影響風柵效

(17)

4

能的主要因子，因此藉由風洞試驗，以不同的孔隙率搭配不同的風速，

來探討最佳孔隙率，從流速、流線分佈、回復距離的多項數據來分析，

發現最佳孔隙率落在0.2 至 0.3。

石(2007)在室內水槽試驗單根透水性筐網圓柱之流場變化情形，並利用流速分佈定義出受筐網圓柱影響而導致速度折減的低流速區域，

在穩定低流速區範圍內，由於流速小且流場的擾動程度低，故若設置於動床河道上將有利於柱後渠床的穩定。梁(2009)在室內水槽針對單排透水性筐網圓柱的流場進行探討，以固定間距及角度排列圓柱的情境下執行試驗，並發現水流會隨著筐網擺設角度之方向行成一股導流作用且當沿岸壁架設筐網圓柱群時，其內側區受到遮蔽，使得框群內側下游段行程流速緩慢之低流速區域。

(18)

5

第二章理論模式

2.1 流場控制方程式

本研究之流場必須滿足不可壓縮流連續方程式與三條動量守恆方程式，以及一條紊流動能方程式和一條紊流消散率方程式：

U V W 0

x y z

  

  

   (2.1)

1 1 xx yx zx

x

DU P

Dt x x y z f

  

 

   

           (2.2)

1 1 xy yy zy

y

DV P

Dt y x y z f

  

 

  

 

           (2.3)

1 1 xz yz zz

z

DW P

Dt z x y z f

  

 

   

           (2.4)

其中，三個維度 x、y、z 的速度分別由 U、V、W 代表，P 為壓力，τ 為剪應力，f 為各種的實體力(如重力、科氏力...等)。

本研究的流體為水，視為牛頓流體，即流體的應變率和所受的剪應力成線性關係，亦即流體之黏滯係數不受剪應力的影響，故上述的三個動量方程式可簡化為：

1 2

i i i

i i

j i j j

U U P U

U f

t x x  x x



        

     (2.5)

此即為那維爾–史托克方程式(Navier–Stokes equation)，式 2-8 為張量形式的寫法。

英國學者雷諾(Osborne Reynolds)於 1985 年提出建議，紊流流場中任何物理量的瞬時變化皆可以分為時間平均量與擾動量，此即為雷諾分解(Reynolds decomposition)：

U U u (2.6)

P P p  (2.7)

(19)

6

利用雷諾分解帶入連續方程式以及那維爾–史托克方程式，並對整個方程式取時間平均可以得到，時間平均的連續方程式以及時間平均的那維爾–史托克方程式：

0

i i

U x

 

 (2.8)

1 ⁱ ^j 2

i i i

i j j i

DU P u u U

Dt x x  x x



  

  

   

   

    (2.9)

式中的u ′u ′項是為雷諾剪應力，為紊流所造成的動量通量。

從上述方程式可以發現未知數共有十個(U、V、W、P、U′V′、

U′W′、V′W′、U′U′、V′V′、W′W′)多過於方程式的數目，故無法直接求解，須先建構能以解出雷諾剪應力的模式，才能進一步解出流場的時間平均速度與壓力，這又稱為紊流的閉合問題(朱，2003)。

2.2 紊流模式 k–ε

為了解決上述的問題，Boussinesq 在 1877 年提出紊流通量與速度梯度成一個正比的關係：

i j

i j ij

j i

dU dU

U U^{ }  ^dx dx ^

    

 

  (2.10)

ν 一般寫為ν_T稱為渦流黏滯係數(eddy viscosity)；而紊流傳輸的速率遠大於層流的傳輸速率，所以那維爾–史托克方程式可以表示成下式：

i 1 i

T i

i j j

DU P U

Dt x x  x f



  

 

       (2.11)

在 1974 年，Launder 和 Spalding 提出了 k–ε 紊流模式，建議在紊 流流場中，大渦流的尺度l 應該與紊流動能 k 以及動能消散率 ε 有關：

(20)

7

3/2/

l k  (2.12)

渦流係數ν_T則為紊流動能k 以及動能消散率 ε 的函數：

2 T

c_ k

 ^  (2.13)

此外，k–ε 紊流模式加入了兩條控制方程式，分別為紊流動能方程 式以及紊流動能消散率方程式：

i j i

T

i i j i j

U U U

Dk k

Dt x c_ x x x x

  ^ ^ ^ 

   

 

          (2.14)

2

1 2

i j i

T

i i j i j

U U U

D c c

Dt x c_ x ^ k x x x ^ k



 ^ ^ ^ ^  ^^ ^ ^^   (2.15)

其中，模式常數c_μ 0.09，c _ε 1.44，c _ε 1.92，c 1.0，c_ε 1.3皆是藉由大量實驗結果來作驗證而得到的，因此k–ε 紊流模式被稱為半經驗 的紊流模式。藉由連續方程式、三個動量方程式，再加上紊流動能跟紊流消散率以及紊流黏滯係數方程式，共七個方程式就可以求解七個未知數 U、V、W、P、k、ε 以及ν_T (朱，2003)。

2.3 邊牆函數

在靠近底床或者是邊牆的地方，於有會有粗糙度以及不滑動邊界的影響，使得垂直於邊牆的紊流速度受到阻滯變小，水平於邊牆的紊流速度則較不受影響，進而使得紊流呈現非等向性，不適用以等向性紊流為假設的k–ε 紊流模式來進行求解；因此本研究搭配使用邊牆函數 (wall function)，以進行接近邊壁流場的求解。

在接近邊牆的流場中，大致可分為三個部分，如圖2.1 示，黏性層 (viscous sublayer)，此層中幾乎都是層流狀態，故黏滯力為主要影響力；

最外層為完全發展紊流層(fully-turbulent layer)，慣性力為主要影響力；

(21)

8

而介於兩層中間的則為過度層(buffer layer)，兩力對此層則有相同的影響力。使用邊牆函數時，不直接對層流層以及過度層作計算，藉由半經驗的公式來連結邊牆到完全發展紊流層；由於不需在邊牆加密網格，此法相當節省計算時間，精確度也相當不錯，故常跟紊流模式搭配使用使用。

本研究採用標準邊牆函數(standard wall functions)是根據 Launder 和 Spalding 所提出的，對接近邊牆的區域進行運算。對平均流速則修正為下式：

1/4 1/2 1/4 1/2

2

*

1ln

p p p p

U c k c k y

U E

  

 

 

   (2.16)

其中，κ 為馮卡門常數，通常為 0.4~0.42 之間；E 為經驗常數，其值為 9.793；U 為流體在 P 點之平均流速；k 為在 P 點之紊流動能；y 為邊牆到P 點的距離；μ 為流體的運動黏滯係數；U 則為剪力速度 U_* 

  

 

 

(Fluent6.1 manual )。

2.4 孔隙介質流模式

為了探討透水性圓柱體對於流場之影響，本研究將透水性圓柱體視為一孔隙介質結構物(porous media)，因此在原本的動量方程式中增加一動量損失項，由流體經過孔隙結構物所受之阻滯效應–壓降 (pressure drop)所產生的；此動量損失項包含了兩個部分–黏滯項損失以及慣性項損失：

1

i i 2 i i

S U C U U



 

   

  (2.17)

其中 α 是滲透係數，C 是慣性阻力因子；等號右邊第一項即從達西定律(Darcy's law)得來的。

(22)

9

由於本研究的流場為高雷諾數，省略黏滯損失項，故只定義慣性損失項的慣性因子C：



²



2 0

1.04 1 kr Cdx C b 



 

    ^(2.18)

其中，k 為阻力係數，ψ為孔隙率，∆ 為孔隙結構物的厚度；慣性阻力係數 C，利用阻力係數k 與孔隙率ψ的關係式，此關係式由 Pinker 和Herbert (1967)及 Reynolds(1969)所提出的(林，2005)。

2.5 體積分率法

為了使模擬更為精確，自由液面是需要考量的，必須同時計算空氣與水的流場，故加入體積分率法(volume of fluid)，以求解流場。

在體積分率法中，引入一個變數F(volume fraction)，來定義一個控制體積內的物理量，而該控制體積內的體積分率和為1，亦即當水的體積分率為F 時，空氣的體積分率則為 1–F，故 F=1 時，則整個控制體積為水，若F=0 時，則該控制體積整個為空氣，而 0 < F < 1 時，控制體積內則有水與空氣，亦即有自由液面存在。

控制體積內的體積分率，是藉由連續方程式計算而得到的，同時兩相流體是不能相互混合的，體積分率計算如下：

i F

S

F F

t U x 

 

 

  (2.19)

一般等號右邊為0，經過上式計算出體積分率後，重新定義控制體積內的其它物理量：



¹



water air

F F

     (2.20)



¹



water air

F F

     (2.21)

(23)

10

接著用上式算出的密度與黏滯係數等，去計算動量方程式、能量方程式，把其它物理變量計算出來(Fluent6.1 manual)。

(24)

11

圖 2.1 邊牆流場示意圖(節錄自 Fluent 6.1 manual )

(25)

12

第三章數值方法

3.1 數值模擬架構

本研究的數值模擬主要分為兩個部分：(1)Gambit 用來建構模擬的範圍，再對建構好的範圍進行網格的劃分，最後再設置適當的邊界條件。(2)Fluent 為計算流場的模擬軟體，讀入前一個步驟所輸出的網格檔以進行運算，並可將所求得的數據保存並輸出(陳，2006)。

主要的執行流程如下，先在Gambit 畫出所需的物理模型，接著在預計流場變化大的地方，或者是想要有細部探討的地方加密網格，選擇三角網格或者是四邊網格，並在相對應的面或體積進行適當的邊界條件設定，儲存後並輸出網格檔。第二部分則為將輸出的網格檔讀入 Fluent 裡，首先對網格進行確認，是否有體積出現負值的情形；接著選取時間項：是否為穩定流或非穩定流(steady flow 或 unsteady flow)；並選取所要使用的紊流模式，本研究使用標準k-ε 紊流模式以及標準邊牆 函數；之後選取所需模擬的介質；設定模擬環境，主要是重力項的設定；再者是邊界條件，入流、出流以及邊牆的粗糙度等；選擇適合的離散方法，再給定初始條件後，就可以開始進行流場的模擬計算。

3.2 有限體積法

目前常用的離散方式有：有限差分法(finite difference method)、有限元素法(finite element method) 以及有限體積法 (finite volume method)。有限元素法一般使用在固體力學為多，也有應用在流體力學，

但大多是在求解層流流況；有限差分法則為最具歷史的解法，現在多已不再使用；而有限體積法則廣泛的在商用軟體以及各個研究領域中被使用。

(26)

13

有限體積法是基於物理量守恆而提出的方法，主要的概念是將整個計算區域，畫分成有限個且互不重疊的控制體積，並將待解的微分方程式對控制體積進行積分，而將得出離散的方程式進行求解；其離散方程式在物理上的意義為控制體積內的通量(flux)平衡，即使是非結構網格、存在不連續解或網格較疏的地方，仍能保持守恆的狀態，對解的穩定及準確度大為提升(Versteeg 和 Malalasekera，2007)。

3.3 控制方程式的離散

本研究採用二階上風法(second order upwind)對動量、體積分率、

紊流動能以及紊流消散率進行離散，壓力則使用body force weighted 離散；壓力–速度的耦合，則用 Patankar 和 Spalding 在 1972 年所提出的 SIMPLE 法(semi-implicit method for linked equation)，為不斷地對壓力進行猜測與修正(guess-and-correct)直到收斂為止(謝，2004)。

SIMPLE 法一開始的第一次迭代需要先給定一個假設的初始壓力值P ，算出一個近似的速度場(U 、V 、W )，將這三個近似的速度場帶入連續方程式求解，得到P^，壓力修正值，將壓力修正為P P^，，速度亦因壓力修正而得到修正的速度場U U^，，將修正的壓力值與速度值，當作第二次迭代的初始值，直到連續方程式滿足收斂條件後停止迭代；這就是所謂的分段計算(segregated)，先求解速度場，再代入連續方程式後，修正壓力與速度，而非同時計算壓力與速度場，通常此法較耦合計算(couple)省時(Versteeg 和 Malalasekera，2007)。

(27)

14

3.4 收斂條件

為了使流場能快速的收斂以減少每一時步的迭代次數，並達到合理的數值解，以下幾點需要注意：

(1) 需要適當網格的劃分，特別是在流場變化較大的部分、阻礙物的附近以及自由液面附近，較為細密的網格是必須的；相對的，流場變化較為小的或是較不需要觀察的部位，則可使用較為稀疏的網格，以減少計算的時間，此即為疏密網格。

(2) 各項計算量殘差值的設定，是影響數值解正確與否的重要因素，

因此本研究對於三個維度的速度、紊流動能、紊流消散率以及連續方程式設定10 算為收斂，才進到下一個時步的計算。

(3) 為了加快解的穩定與速度，並控制每次迭代各個變數值的修正，

需給予適當的鬆弛因子(under–relaxation factor)，各個鬆弛因子的係數給定如表3.1。

(4) 由於流場會隨時間變化，時間間距(time step size)的選取是相當重要的，根據庫倫數(courant number)的計算而得出，如下式：

Courant nu t ber x

U

m 

  (3.1)

其中，Δx 為最小網格長，U 為流場流速，本研究是給定為入流速度(顧，

2007)。

3.5 邊界條件設定

為了使得控制方程式得以被求解，需要設置邊界條件，故本研究在入流處使用速度入口作為邊界條件，頂端使用對稱性邊界條件，出流處則以自由出流或是速度入口作為邊界條件，其它地方則為不滑動

(28)

15

邊界條件(Dalpe，2009)。

(1) 速度入口(Velocity inlet)：

給定入流的速度以及方向，由於是紊流流場，必須給定合適的紊流量值，一般是給紊流動能k 以及紊流動能消散率 ε，但由於此二值不容易量測，故通常以紊流強度 I (turbulent intensity)以及長度尺度 l (length scale)作為替代，如下所示。

 

¹⁸

0.16 Re

I  ^ (3.2)

0.07L

 (3.3)

其中，Re 為特徵雷諾數，並以圓柱直徑為特徵長度；L 跟渦流尺度有關，以水深或圓柱直徑代入。

(2) 自由出流(Outflow)：

當出口的流速或壓力不清楚時，可以使用自由出流邊界，因此該邊界上不需要定義任何的物理量，而是由內部的模擬值外插到邊界上，也就是任何物理量的通量在邊界上皆為零，但在計算可壓縮流或密度會變動的流體時，不可使用自由出流邊界，且須確定到達自由出流邊界時為完全發展流，才不會影響計算。

(3) 對稱性邊界(Symmetry)：

對稱性邊界不會有任何流體通過，垂直該邊界的速度設為零，且任何物理量的梯度都為零，其變數值皆相同於最靠近該邊界的結點值。

(4) 邊牆(Wall)：

為不滑動(no slip)邊界，亦即在該邊界上速度皆為零；在此邊界上可以設定粗糙高度(roughness high)，因此邊牆函數 2-19 式須加入一修正項ΔB(為一粗糙函數)：

(29)

16

1/4 1/2 1/4 1/2

2

*

1ln

p p p p

U C k C k y

E B

U

  

 

 

   

  (3.4)

而ΔB 決定與無因次化後的粗糙高度K_s ρK_sC_μ^1/4κ_p^1/2/μ有關，一般分為三個區域討論：

當K_s^ 2.25時：

0

 B (3.5)

當2.25K_s^90時：

 

 

2.25

1ln sin 0.4258 ln 0.811

87.75

s

s s s

B K C K K



  

  

     

  (3.6)

當90K_s^時：

 

1ln 1 _s _s

B C K



    (3.7)

其中，K 為實際的粗糙高度(roughness height)，C 是粗糙常數(roughness constant)。

(30)

17

表3.1 鬆弛因子係數的設定

鬆弛因子係數

Pressure 0.3 Density 1

Body forces 1

Momentum 0.7 Turbulent Kinetic Energy 0.8

Turbulent Dissipation Rate 0.8 Turbulent Viscosity 1

(31)

18

第四章模式驗證

4.1 實心圓柱驗證

為了建置適當的邊界條件以及相關紊流參數，本研究藉由 Ahmed 和Rajaratnam 在 1998 年進行之單根實心圓柱流場試驗，此試驗是在長 20 公尺、寬 1.22 公尺，以及水深維持 0.18 公尺的水槽中所進行的，並將圓柱擺放在距離入流處13 公尺遠的地方，使水流在達到圓柱前能完全發展，圓柱直徑為0.089 公尺並突出水面，底床鋪上一層細砂，細砂的平均粒徑大小D 為 0.00184 公尺。並針對動床、定床，平滑或是粗糙底床等行情形進行流速分佈、圓柱前壓力分佈之探討。從實驗結果發現，粗糙底床會增大壓力梯度，使得圓柱前的向下流(downflow)變強；動床而出現淘刷的坑洞(scour hole)亦會增加向下流的強度。不論底床情形為何，不透水的實心圓柱都會產生極為複雜的渦流流況，像圓柱後的卡門渦街和馬蹄型渦流(Rahman 等，2007)。

圖4.1 為驗證 Ahmed 和 Rajaratnam(1998)的實心圓柱試驗所建構的模型，圖4.2 為模擬時所使用的網格圖，為了能精確的捕捉到實心圓柱週遭流況的變化，在圓柱四周予以適當的加密，網格最小體積為 1.08 10 立方公尺，最大體積則為2.81 10 立方公尺，總網格數為76076 個，時間間距為 0.005 s。同時搭配 Salaheldin 等(2004)的模擬情況描述，使用k–ε 紊流模組搭配 VOF 二相流模式來計算整個流域，

以期能更真實的模擬出水流狀況外，底床鋪細砂所造成的影響，亦在底床加設粗糙高度，反映出接近底床的模擬情形。入流邊界為速度入口，頂端為對稱性邊界，出流邊界分為兩部分，空氣層為自由出流，

水體則為速度入口，其它水理條件如表4.1。

(32)

19

在圓柱前方2.25 倍半徑(x=-0.1，y=0)及 10 倍半徑(x=-0.445，y=0) 處，觀察軸向速度(U)和垂向速度(W)沿 Z 軸之變化。由圖 4.3 可知，軸向速度受底床不滑動邊界的影響，隨著距離底床之距離增加，軸向速度亦隨之上升，在接近自由液面時速度達到最大值；而距圓柱之距離越近速度越慢，是由於水體受到圓柱的阻滯效應。圖 4.4 為 W 在 z 軸之分佈，底床為不滑動邊界，故速度為0，隨著距離底床之距離增加，

垂向速度先減小再上升，這是由於水流通過圓柱等鈍形體(blunt body) 時，柱前停滯點壓力甚大，而產生二股分別向上與向下之水流，其中向下水流(downflow)會與橫向來流結合而在柱前形成「馬蹄型渦流」，

故W 呈現出負值的情形(石，2007)。

在圓柱側邊三個位置作進一步的驗證，第一個位置為距離圓柱側邊3 倍半徑的地方(y = 0.178 m)，水深 0.1456 公尺處，我們將觀測沿 x 軸的速度變化(圖 4.5)。由圖可以發現，在未到達圓柱時，速度緩慢的增加，經過圓柱時速度有著較為顯著的增加，隨後速度逐漸降低，這是因為圓柱減少了通水面積，使得在圓柱側邊的速度增快。第二個位置為距離圓柱上5 倍半徑(即 x=-0.267 m)，水深 0.1456 公尺的地方，我們將觀測沿 y 軸方向的速度變化(即圖 4.6)。中心部分速度最慢，是由於受到圓柱阻滯效應的關係，使得速度沿 y 軸漸漸增加。第三個位置為水深0.1456 公尺，圓柱中心線，沿 y 軸的速度變化圖(即圖 4.7)。由於通水斷面積受到圓柱的阻礙而減小，故速度在經過圓柱時明顯的增加。

(33)

20

4.2 孔隙圓柱驗證

石(2007)使用長 18 公尺、寬 1 公尺和高 0.4 公尺的水槽，進行室內水工試驗，在渠道尾端有蓄水槽以及馬達，將水送至渠道前方，經由整流網流入試驗渠道，以完成一循環水路系統。在距離入流 8 公尺處擺放圓柱體，使水流在到達筐網圓柱時能形成完全發展流，圓柱直徑分別為0.1 公尺與 0.16 公尺，孔隙率則為 0.38 和 0.62，筐網圓柱厚 0.0006 公尺，底床以及邊壁均為光滑表面，其中孔隙率的定義為表面孔洞的面積總和佔筐網圓柱表面積總和之比例，其它的水理條件如表 4.1。

完成k–ε 紊流模組以及 VOF 模式的建置後，引入孔隙介質流模式，

以呈現筐網圓柱的阻水效應，並藉由石(2007)的筐網圓柱室內水工試驗，作為孔隙介質流的模式驗證。圖4.8 為驗證石(2007)的筐網圓柱試驗所建構的模型，圖 4.9 與圖 4.10 為直徑 0.1 公尺和 0.16 公尺筐網圓柱模擬時所使用的網格圖，在筐網圓柱四周加密，是為了捕捉圓柱附近的流場變化。直徑 0.1 公尺筐網圓柱的案例，網格最小體積為 4.78 10 立方公尺，最大體積為6.54 10 立方公尺，總網格數為 61030；直徑 0.16 公尺筐網圓柱的案例，網格最小體積為8.96 10 立方公尺，最大體積為1.14 10 立方公尺，總網格數為 51700 個，時間間距皆為0.01 s。邊界條件部分，入流邊界為速度入口，頂端為對稱性邊界，出流邊界則分兩部分，空氣層為自由出流，水體為速度入口。

先觀察軸向速度(U)沿 x 軸方向的變化，取圓柱中心線，水深 0.064 公尺深的地方(即圖 4.11)。當水體接近筐網圓柱時，因受到圓柱阻滯，

速度開始下降，但由於筐網圓柱會透水，且從筐網圓柱旁流過的側流

(34)

21

會匯入，故速度降低到一定之後，就會開始回升，但不會回升到原本的入流速度；相反地，因著實心圓柱的不透水性，圓柱後有著強烈的負壓區，進而產生向下之渦流，使得水體回流向圓柱(Laws and Livesey，1978)。另外觀察圓柱側邊 1.5 倍圓柱(y=0.15，z=0.056)與 2.5 倍圓柱(y=0.25，z=0.056)處之流場變化(即圖 4.12 與圖 4.13)。在水流流經筐網圓柱時，由於通水斷面積減小，流速會增快，孔隙率小(ψ=0.38) 又比孔隙率大(ψ=0.62)的流速增加的多，主要是流體穿越孔隙率大的 (ψ=0.62)較多，從圓柱兩旁流過的流體相對的減少，故流速增加量也就相對的低。

接下來探討軸向速度(U)沿 y 軸方向的變化，分別在圓柱中心 (x=0，z=0.056)、下游距圓柱 1.5 倍圓柱(x=0.15，z=0.056)和下游距圓柱 2.5 倍圓柱(x=0.25，z=0.056)三個位置，如圖 4.14 至圖 4.16 所示。

由於通水斷面積減小而產生的束縮效應，使得流速加快；而孔隙率小的流速增加幅度比孔隙率大的來的多，亦符合物理現象。

至於軸向速度(U)沿 z 軸方向的變化情形，在圓柱中心線、下游方向距圓柱2.5 倍圓柱處(x=0.25，y=0)(即圖 4.17)。底床為不滑動邊界，

故速度仍然為 0，但是因為筐網圓柱可以透水，所以不論孔隙率大或小，流速皆隨著距離底床的距離增加而變快。

4.3 驗證結果與討論

在實心圓柱的案例中，執行了k–ε 紊流模組和 VOF 自由液面模式 的參數建置；在可透水筐網圓柱的案例中，進行測試孔隙介質流模式是否可以呈現受筐網圓柱影響而產生之阻水效應。模擬結果與實驗結果經對照後相當吻合，顯示模組的正確性與適用性。

(35)

22

在筐網圓柱的案例中，孔隙率因子比圓柱直徑因子對於流場有著更大的影響。圓柱直徑大小對於流場之影響，在軸向速度沿 y 軸方向的變化上較為明顯，主要原因為室內試驗水槽的渠寬為固定的，當筐網圓柱直徑變大時，通水斷面積減小更多，流體受擠壓的程度更大加速趨勢也較明顯，特別在孔隙率小(ψ=0.38)時尤其顯著。相反地，在沿 z 方向之水平速度的變化上，幾乎沒有任何的差異，顯示筐網圓柱直徑大小的變化，對沿 z 方向之水平速度的變化影響不大(Franssona 等，2004)。

由流場流線圖觀察流場變化的情形(即圖 4.18 與圖 4.20)，實心圓柱的下游方會有週期性的卡門渦街；相反地，筐網圓柱的下游方無卡門渦街產生，這是由於筐網圓柱可透水，讓圓柱前後壓力差值減小，

使得流體分離點延後。另外，由圖4.19 及圖 4.21 可發現，孔隙率 0.38 的案例中，當流體穿透筐網圓柱後，流線會呈現上揚的現象，其原因應為流體在流經圓柱兩側後，流向水槽中央並抬升穿透筐網之流體，

造成底床附近流線密集而且圓柱後方流線上揚之結果。

4.4 孔隙區域與孔隙跳躍之比較

針對孔隙介質流模式的模擬方式有兩種，一為孔隙區域(porous zone)，另一為孔隙跳躍(porous jump)邊界，如圖 4.22。孔隙區域，需要設定孔隙率以及三個維度的黏滯係數與阻力係數；孔隙跳躍邊界需給定黏滯係數與阻力係數以及孔隙結構物的厚度，因為三個維度的黏滯係數與阻力係數相同，所以動量扣減為等向性的，且當孔隙結構物的厚度非常薄時，適用孔隙跳躍邊界。本研究針對進行模擬，並與石 (2007)試驗結果比對(圖 4.23)。比較結果發現，使用孔隙跳耀邊界時，

(36)

23

水體經過筐網圓柱後速度的扣減與試驗結果吻合，但使用孔隙區域時，

經過筐網圓柱後水體速度有低估的趨勢。兩者的結果差異在於邊界條件的設置是否與物理試驗裝置一致，物理試驗的筐網圓柱為一不鏽鋼薄片所圍成之透水性圓柱(圖 4.24)，厚度只有6 10 ^⁴公尺且筐網圓柱中間皆為空心，與一般充滿物質的孔隙結構物不相同。因此，本研究之後的應用案例，皆選用孔隙跳躍邊界作為呈現阻水效應之孔隙介質的模擬方式。

(37)

24

表4.1 兩個驗證案例的水理參數

Ahmed 等(1998) 石武融(2007)

流量Q (l/s) 65.0 38.1

平均水深 H(m) 0.18 0.14

平均流速Uo(m/s) 0.293 0.272

福祿數Fr 0.22 0.26

底床坡度So 0 0.02

(38)

25

圖 4.1Ahmed 等(1998)的試驗設計模型圖

圖4.2 實心圓柱網格圖

(39)

26

圖4.3 軸向速度沿 z 軸的變化

圖4.4 垂向速度沿 z 軸的變化

(40)

27

圖4.5 無因次化的速度強度沿 x 軸的變化

圖4.6 無因次化的速度強度沿 y 軸的變化

圖4.7 無因次化的速度強度沿 y 軸的變化

(41)

28

圖4.8 石(2007)的試驗設計模型圖

圖4.9 直徑 0.10 m 筐網圓柱網格圖

圖 4.10 直徑 0.16 m 筐網圓柱網格圖

(42)

29

圖 4.11 軸向速度沿 x 軸方向之變化（y = 0 m，z = 0.056 m）

(a)圓柱直徑 10 公分 (b)圓柱直徑 10 公分

(43)

30

圖4.12 軸向速度沿 x 軸方向之變化（y = 0.15 m，z = 0.056 m）

(44)

31

圖4.13 軸向速度沿 x 軸方向之變化（y = 0.25 m，z = 0.056 m）

(45)

32

圖4.14 軸向速度沿 y 軸方向之變化（x = 0 m，z = 0.056 m）

(46)

33

圖4.15 軸向速度沿 y 軸方向之變化（x = 0.15 m，z = 0.056 m）

(47)

34

圖4.16 軸向速度沿 y 軸方向之變化（x = 0.25 m，z = 0.056 m）

(48)

35

圖4.17 軸向速度沿 z 軸方向之變化（x = 0.25 m，y = 0 m）

(49)

36

圖4.18 孔隙率 0.38 之 x–y 平面的流線圖 (a)直徑 0.1 m (b)直徑 0.16 m

(50)

37

圖 4.19 孔隙率 0.38 之 x–z 平面的流線圖 (a)直徑 0.1 m (b)直徑 0.16 m

(51)

38

圖4.20 圓柱直徑 0.1 m 之 x–y 平面的流線圖 (a)孔隙率 0.38 (b)孔隙率 0.62

(52)

39

圖 4.21 圓柱直徑 0.1 m 之 x–z 平面的流線圖 (a)孔隙率 0.38 (b)孔隙率 0.62

(53)

40

圖4.22 porous jump 與 porous zone

圖4.23 軸向速度沿 x 軸之變化(圓柱直徑 0.1 公尺，孔隙率 0.38)

(54)

41

圖 4.24 試驗所使用的筐網圓柱圖

(55)

42

第五章模式應用

5.1 最佳孔隙率的探討

在筐網圓柱對流場影響的因子中，孔隙率因子比直徑大小因子更為顯著，故本節將針對孔隙率作進一步的分析與探討。因孔隙率小，

阻水面積大，穿越筐網圓柱的流量少，流場較為不穩定，與實心圓柱的流況也較為相似。反之，當孔隙率越大，阻水面積小，穿越筐網圓柱的流量就多，流場就相對穩定，比較像穿越無阻礙的流場，但流速折減的少，後方的穩定低速區也小。為此，找到一個能適當降低流速，

又盡可能使後方穩定低速流場的區域大的孔隙率，是本節所要探討的重點(Chan 和 Lien， 2005)。

模擬使用長11 公尺、寬 1.1 公尺、高 0.24 公尺的模型，筐網圓柱直徑 0.1 公尺，且放置在距離入流口處 5 公尺的地方並突出水面(即圖 5.1)。流況方面，流量為 43.52 公升/秒，水深維持 0.16 公尺，福祿數為 0.22。並將筐網圓柱視為一孔隙結構物，使水流經過筐網圓柱而產生之阻水效應，當作一動量損失項加在動量方程式中，並利用孔隙率與孔隙結構物的厚度來設定阻力係數。除筐網圓柱的厚度固定不變外，將分為九組孔隙率進行模擬分析(即表 5.1)，分別為低孔隙率 0.1、0.2、

0.3，中孔隙率 0.35、0.4、0.45 和高孔隙率 0.5、0.55、0.6。為了兼顧計算精確性以及計算效率，網格的劃分採用不均勻網格，並在筐網圓柱附近加密(即圖 5.2)，並在遠離筐網圓柱後網格逐步變疏。其中，網格最小體積為2.31 10 ⁶立方公尺，最大體積則為2.64 10 ⁴立方公尺，總網格數為28164 個，經由 courant number 計算所得的時間間距為 0.01 秒。模式則使用 k–ε 標準紊流模式以及 VOF 模擬流場流況，並以

(56)

43

孔隙介質模式呈現筐網圓柱所產生之阻水效應。選用速度入口作為入流邊界，頂端部分使用對稱性邊界，出流口分為兩個部分，空氣的出流選用自由出流邊界條件，水的出流則使用速度入口。

藉由總壓力的分佈，來圈選出筐網圓柱後方保護區域的大小，將總壓力小於一參考壓力值(0.8 倍的總壓)視為受到筐網圓柱保護的區域 (即圖 5.3 至圖 5.6)，並將水深分成四個等份，深度分別為由底床向上水深0.2、0.4、0.6 以及 0.8 倍，以此來探討最佳孔隙率。

由相對低壓區的壓力分佈圖(即圖 5.3 至圖 5.6)得知，當孔隙率小於0.3 時，受筐網圓柱保護之區域較小，這是因為筐網圓柱的孔隙率過小，其流況類似於實心不透水圓柱的流況，筐網圓柱後方的速度大幅降低，而經由筐網圓柱兩側流過的流體流速則大幅增加，並在筐網圓柱後方快速的擠壓進來，因此受到筐網圓柱保護之區域較小。孔隙率大於0.3 時，受到筐網圓柱保護之區域顯著地增加，這是由於穿過筐網圓柱的流體(一般稱穿越流，bleed flow)增多，並且開始降低從筐網圓柱兩側流過之流體的速度，此時流經筐網圓柱兩側之流體需較遠的距離，才會再次擠入從筐網圓柱穿越出來的流體，此現象即為分離點延後，因此受到筐網圓柱保護的區域增加。但當孔隙率超過 0.45 時，受筐網圓柱保護的區域開始縮小，這現象是因為孔隙率過大，穿越筐網圓柱的流體流速並無減慢多少，與從筐網圓柱兩側流過的流體流速相差無多，因此兩者很快就匯合在一起，受筐網圓柱保護的區域也隨之減小。

圖5.3 至圖 5.6 為不同深度的筐網圓柱保護區域。當深度越深時，

筐網圓柱保護的區域越小；反之，當深度越淺時，筐網圓柱保護的區域越大。不過，筐網圓柱的保護區域大小與孔隙率之間的關係並無改

(57)

44

變，仍然介於孔隙率0.35 到 0.45 時有最大的保護區域。

本研究利用四個不同深度所求得的相對低壓區面積，並假設不同深度間的面積為線性關係(圖 5.7)，所以就將保護體積近似為梯形體，

始得孔隙率與保護體積之關係圖(圖 5.9)，規定擁有最大保護體積的孔隙率為最佳孔隙率。從孔隙率與保護體積的關係圖可看出，隨著孔隙率增大保護體積隨之增加，尤其是在孔隙率0.3 以後，保護體積大幅的上升，這是因為孔隙率0.3 之後，穿越框網圓的穿越流變成流場的主控因子，其流速強度開始增大並能影響從筐網圓柱兩側經過的側流。保護體積之峰值出現在孔隙率 0.4~0.45 左右，之後呈現遞減的趨勢，原因為穿越筐網圓柱的流體與經過筐圓柱兩旁的流體間的速差已無太大差異，受到筐網圓柱保護之區域因而開始縮減。綜觀以上結果，最佳孔隙率約為0.43 的位置 Dong(2007)。此外，從不同孔隙率之自由液面圖(圖 5.8)也可看出，當孔隙率為 0.4 之後，水體受圓柱阻擋而堆高的情形已經趨於平緩，亦印證了前述所定義的最佳孔隙率。

5.2 單排圓柱的探討

本研究將進一步探討單排筐網圓柱對流場的影響。由於實際應用於現地時，不會只使用單根的筐網圓柱，而是使用多排的型式，在此先以探討單排型式對流場的影響情形為主。以固定的角度及固定的圓柱間距來排列，案例有一根、兩根、三根以及五根筐網圓柱，由渠道中向渠道邊壁依序增加圓柱擺放的數量，以期能觀察出不同圓柱根數對流場的影響。

模擬將使用長11 公尺、寬 1 公尺、高 0.24 公尺的模型，維持水深 0.16 公尺，筐網圓柱直徑 0.025 公尺，並放置在距入流口 5 公尺的處且

(58)

45

突出水面(圖 5.10)，流量為 37.12 公升/秒，福祿數為 0.19，圓柱間相距0.1 公尺，並以 30 度角排列(圖 5.11)(梁，2009)。另外，以相同根數與相同排列的實心圓柱作為比較。為了兼顧計算的精確性以及計算效率，採用不均勻網格，在圓柱附近加密，如圖 5.12 所示，在遠離圓柱後網格逐漸變疏，網格最小體積為4.40 10 立方公尺，最大體積則為5.96 10 立方公尺，總網格數為 82884 個格點，時間間距為 0.005 秒。模式則使用k–ε 標準紊流模式以及 VOF 模擬流場流況，並以孔隙 介質模式呈現筐網圓柱所產生的阻水效應。入流邊界使用速度入口，

頂端使用對稱邊界，出流邊界分兩部分，上層空氣使用自由出流，下層水體則使用速度入口。

以下將分別討論不同根數之筐網圓柱對流場的影響，以及實心圓柱與透水性筐網圓柱，在相同根數與相同排列下對流場的影響。

5.2.1 不同根數之筐網圓柱的流場流速分佈比較

在第五號圓柱後方約兩倍圓柱直徑(x=0.4 m，y=0.2125 m)處，來觀察軸向速度沿z 軸的變化(即圖 5.13(a))。穿越五號圓柱而流出的水體速度，在近底床處受到底床不滑動邊界之影響，隨著高度增加流速也隨之增加。在只有單根筐網圓柱時流速最快。當有兩根筐網圓柱時，因著四號筐網圓柱的導流效應，部分流體從四號筐網圓柱兩側流去，在加上五號筐網圓柱本身的阻水效應以及導流效應，水體進入五號筐網圓柱的速度已經有減緩的趨勢，所以流體流出五號筐網圓柱的速度比只擺放單根筐網圓柱時緩慢。但是當有三根筐網圓柱時，水體流出五號筐網圓柱的速度反而有稍微回升的趨勢，比擺放兩根筐網圓柱時的速度快，這是由於三號筐網圓柱導流後，使得流體進入五號圓柱的速