應用沉浸邊界法與移動網格模擬紊流流況下之結構物沖刷現象

國立臺灣大學工學院應用力學研究所 碩士論文

Institute of Applied Mechanics College of Engineering

National Taiwan University Master Thesis

應用沉浸邊界法與移動網格模擬紊流流況下之 結構物沖刷現象

Application of the immersed boundary method and arbitrary Lagrangian-Eulerian scheme to simulate local scour

in turbulent flow

楊善合 Shan-He Yang

指導教授：周逸儒 博士 Advisor: Yi-Ju Chou, Ph.D.

中華民國 103 年 7 月

July, 2014

致謝

時間過得很快，轉眼間人生最終的求學階段也要結束了，回想當初對自己的 未來感到懷疑，懵懵懂懂只想找件事情來學習，剛踏入校園內心的雀躍以及進入 實驗室擁有屬於自己的位置，就好像昨日一般。

學習的兩年間，首先感謝的是我的指導教授周逸儒老師，在這裡我學會到獨 立思考的重要性，能夠不受外在影響，擁有自我去判斷是非才能算是一個獨立的 人格，也謝謝我的口試委員，江允智教授與劉啟民教授，願意幫忙修正論文上的 疏失以及提供寶貴的建議。

也感謝實驗室的學長與助理們，宗緯、小光、宏志與秉憲，研究與修課中的 許多問題，謝謝你們不厭其煩的指導，還有一年級一起為必修奮鬥的好戰友，聖 堯、詩弘、科斯曼、阿罵、氣球、雨彤，沒有大家的幫忙，必修課一定會修得更 加辛苦吧！以及實驗室的學弟們，鈐鏞、郁誠，對我來說能認識大家真的是我來 應力所最大的收穫。

最後，謝謝我的父母，在我求學的過程中，從來沒有讓學習成為我的負擔，

一直都是讓我自由發展沒有壓力的過完每個階段，形式上的求學生涯已經結束，

但人生的學習才正要開始，希望我能牢記這段時間的訓練，在未來繼續追尋屬於 我的自由。

中文摘要

本研究使用大渦流模式結合沉浸邊界法與移動網格法，模擬流體經過一結構 物對底部泥沙所造成沖刷堆積的影響，大渦流模式(Large eddy simulation)為三維的 水動力模式，較能捕捉到在高雷諾數下的紊流現象，在高雷諾數的流場中特別重 要。結構物表面的複雜邊界主要使用沉浸邊界法(Immersed boundary method)來處 理，與傳統的貼體法(body-fitting)相比，沉浸邊界法假想力的形式寫進統御方程式 使其滿足邊界條件，在卡氏座標即可很有效率的處理複雜幾何或是移動網格的問 題，另外再底床網格的部分則使用移動網格法(Arbitrary Lagrangian-Eulerian scheme) 計算網格座標速度使其滿足泥沙傳輸，並且模擬紊流邊界層中的底床形貌變化。

本研究一開始使用沉浸邊界法模擬流體流經一圓柱，改變雷諾數由潛變流模 擬至紊流，觀察其尾流的變化，確認邊界與流場符合現實情況後，我們再依照 Roulund, et al.(2005)所做的實驗配置，設定邊界條件進行模擬，並且分別討論圓柱 前緣沖刷與圓柱下游沙漣(sand ripples)的形成過程，以及紊流動能在沖刷過程中所 扮演的腳色，結果發現沖刷最強的區域通常都受到一個持續穩定的上升流場給帶 動，如前緣的馬蹄形渦流(horseshoe vortex)。一開始前緣掏刷出來的泥沙會在圓柱 後方低壓區域產生堆積，堆積丘的下游處容易產生分離流動使得沙漣更加明顯。

而紊流動能大的區域通常為流場交界處，該區域由於其動能消散(disspation rate)十 分強，故會使部分懸浮泥沙失去動能而在此沉積。

關鍵字:大渦流模式、沉浸邊界法、移動網格法、局部沖刷、沙漣、馬蹄形渦流

ABSTRACT

In this study, we apply the large eddy simulation (LES) code that combines the immersed boundary method (IBM) and the arbitrary Lagrangian-Eulerian method (ALE) to simulate the evolution of the erodible bed around a structure. This code is a three dimensional computational fluid dynamics simulator, which is capable of resolving the detailed turbulent flow field. This is particularly important in the high Reynolds number flow. We employ the IBM to model the surface of the structure.

Compared to traditional body-fitting methods, IBM applies the body force to satisfy the desired boundary condition. It can efficiently handle the complex geometry and moving grids in Cartesian coordinate system. In addition, we apply the ALE method in our grid. It can calculate the grid velocity to guarantee conservation of sediment mass and simulate bed form dynamics in a turbulent boundary layer.

We simulate flow over a cylinder as a test case for IBM. Cases ranging from the regimes of creeping flow to turbulent flow are investigated. The present numerical model is then validated against the experimental results by Roulund et al.(2005), who investigated erosion around a cylinder in a laboratory setting. We discuss erosion in front of the cylinder, the development of sand ripples behind the cylinder, and turbulence kinetic energy (TKE). The numerical results show that the most dramatic erosion region is affected by the steady upflow associated with the horseshoe vortex.

Moreover, we observe that sediments in front of the edge pile up in the lower pressure area. Flow separation easily occurs behind the ripples, leading to the growth of the ripple amplitude. We found that high TKE occurs at regions of flow convergence, which is usually associated with strong dissipation, leading to the deposition of suspended sediments.

Keywords: large-eddy simulation, immersed boundary method, arbitrary Lagrangian -Eulerian method, local scour, sand ripples, horseshoe vortex

總目錄

致謝 ... i

中文摘要 ... ii

ABSTRACT ... iii

圖目錄 ... vi

表目錄 ... ix

Chapter 1 緒論 ... 1

1.1 前言 ... 1

1.2 研究動機 ... 2

1.3 文獻回顧 ... 3

1.3.1 沉浸邊界法 ... 3

1.3.2 底床沖刷 ... 4

Chapter 2 理論背景與方法 ... 7

2.1 大渦流模式 ... 8

2.1.1 統御方程式 ... 9

2.1.2 濾波統御方程式 ... 10

2.1.3 曲線座標格式 ... 12

2.1.4 流體離散化計算流程 ... 13

2.2 移動網格與泥沙傳輸方程式 ... 15

2.2.1 密度分層 ... 15

2.2.2 泥沙傳輸方程式 ... 16

2.2.3 移動網格法 ... 17

2.2.4 掏刷方程式 ... 18

2.2.5 底床高度方程式 ... 19

Chapter 3 沉浸邊界法介紹與驗證 ... 22

3.1 沉浸邊界法離散化 ... 22

3.2 流經圓柱之潛變流驗證 ... 27

3.3 流體流經圓柱(不同雷諾數) ... 31

Chapter 4 底床沖刷模擬與分析 ... 39

4.1 模擬配置 ... 39

4.2 流場分析 I(定床) ... 42

4.2.1 馬蹄形渦流 ... 42

4.2.2 尾跡渦流 ... 45

4.2.3 底床剪應力 ... 46

4.2.4 紊流動能 ... 46

4.3 流場分析 II(動床) ... 48

4.3.1 圓柱前緣分析 ... 48

4.3.2 圓柱後端分析 ... 53

4.4 紊流動能分析 ... 57

4.5 侵蝕深度分析 ... 59

Chapter 5 結論與未來工作 ... 61

5.1 結論 ... 61

5.2 未來工作 ... 63

參考文獻 ... 64

圖目錄

圖 1-1 局部沖刷機制的示意圖(Roulund et al.,2005) ... 2

圖 1-2 沖刷實驗圓柱後緣高度變化(Dargahi, 1990) ... 6

圖 1-3 (a)圓柱沖刷實驗與(b)模擬結果(Roulund et al.,2005) ... 6

圖 1-4 底床形貌與濃度的演進(Chou & Fringer, 2010) ... 6

圖 2-1 卡氏座標方向 ... 9

圖 2-2 底床顆粒傳輸示意圖 (Chou & Fringer, 2010) ... 19

圖 2-3 靜止底床角度示意圖(Chou & Fringer, 2010) ... 21

圖 3-1 (a)網格分類表示圖與(b)內插速度示意圖 ... 23

圖 3-2 內插點示意圖 ... 24

圖 3-3 整體數值計算流程 ... 26

圖 3-4 圓柱表面無因次化壓力分佈圖 ... 27

圖 3-5 流場模擬配置圖 ... 28

圖 3-6 流場側視圖 x-y 剖面(z = 1.5 m) ... 28

圖 3-7 流場俯視圖 x-z 剖面(y = -0.5 m) ... 28

圖 3-8 x 方向速度(U)側視圖(z = 1.5 m) (Re = 0.1) ... 30

圖 3-9 x 方向速度(U)俯視圖(y = -0.5 m) (Re = 0.1) ... 30

圖 3-10 y 方向渦度(Vorticity) 俯視圖(y = -0.5 m) (Re = 0.1) ... 30

圖 3-11 無因次化壓力P/ (μU R/ )^在z 軸中心與圓柱表面上的分布 ... 30

圖 3-12 流場模擬配置圖 ... 31

圖 3-13 流場側視圖 x-y 剖面(z = 0.5 m) ... 31

圖 3-14 流場俯視圖 x-z 剖面(y = -0.5 m) ... 32

圖 3-15 x 方向速度(U)側視圖(z = 0.5 m) (Re = 10) ... 34

圖 3-16 y 方向渦度(Vorticity)俯視圖(y = -0.5 m) (Re = 10) ... 34

圖 3-17 y 方向渦度大小為 0.2 時，在三維座標下的等高面圖(Re = 10) ... 34

圖 3-18 x 方向速度(U)側視圖(z = 0.5 m) (Re = 30) ... 35

圖 3-19 y 方向渦度(Vorticity)俯視圖(y = -0.5 m) (Re = 30) ... 35

圖 3-20 y 方向渦度大小為 0.2 時，在三維座標下的等高面圖(Re = 30) ... 35

圖 3-21 x 方向速度(U)側視圖(z = 0.5 m) (Re = 100) ... 36

圖 3-22 y 方向渦度(Vorticity)俯視圖(y = -0.5 m) (Re = 100) ... 36

圖 3-23 y 方向渦度大小為 0.2 時，在三維座標下的等高面圖(Re = 100) ... 36

圖 3-24 x 方向速度(U)側視圖(z = 0.5 m)(Re = 1000) ... 37

圖 3-25 y 方向渦度(Vorticity)俯視圖(y = -0.5 m) (Re = 1000) ... 37

圖 3-26 y 方向渦度大小為 0.2 時，在三維座標下的等高面圖(Re = 1000) ... 37

圖 3-27 流體流經一光滑圓柱實驗結果(Blevins, 1990) ... 38

圖 4-1 流場模擬配置圖 ... 39

圖 4-2 展寬網格側視圖(x-y 方向) ... 39

圖 4-3 實際模擬網格分布側視圖 ... 40

圖 4-4 流場修正區示意圖 ... 40

圖 4-5 y 方向速度(V)側視圖(z = 0.3m)(時間平均流場) ... 42

圖 4-6 流線側視圖(z = 0.3m) (時間平均流場) ... 42

圖 4-7 x 方向速度(U)底床俯視圖(y = -1.999 m) (時間平均流場) ... 42

圖 4-8 z 方向渦度(Vorticity)側視圖(z = 0.3m) (時間平均流場) ... 43

圖 4-9 圓柱前緣渦度實驗結果 (Re_D =20000) (Dargahi,1990) ... 43

圖 4-10 q-criterion 三維示意圖(時間平均流場) ... 44

圖 4-11 壓力俯視圖(y = -0.1 m)(時間平均流場) ... 45

圖 4-12 y 方向渦度俯視圖(y = -0.1 m)(瞬時流場) ... 45

圖 4-13 底床剪應力俯視圖(y = -1.999 m) (時間平均流場) ... 46

圖 4-14 紊流動能(k)底床俯視圖(y = -1.999m) (時間平均流場) ... 47

圖 4-15 紊流動能(k)側視圖(z = 0.3m) (時間平均流場) ... 47

圖 4-16 紊流動能(k)前緣側視圖(z = 0.3m) (時間平均流場) ... 47

圖 4-17 圓柱前緣渦度與侵蝕變化側視圖(z = 0.3 m) ... 50

圖 4-18 (a)圓柱前緣實驗與(b)模擬結果(z = 0.3 m) (10 min) ... 50

圖 4-19 底床深度形貌俯視圖(0-300 s) ... 51

圖 4-20 底床深度形貌俯視圖(300-600 s) ... 52

圖 4-21 底床 x 方向速度(U)俯視圖(y = -1.999 m) (時間平均流場) ... 54

圖 4-22 底床 x 方向平均速度(U)側視圖(z = 0.3m) (時間平均流場) ... 54

圖 4-23 底床 y 方向平均速度(V)側視圖(z = 0.3m) (時間平均流場) ... 54

圖 4-24 泥沙濃度等高面三維座標示意圖 ... 55

圖 4-25 圓柱尾端沙漣演進側視圖(z = 0.3m) ... 56

圖 4-26 (a)圓柱後端實驗與(b)模擬結果(z = 0.3m) (10 min) ... 56

圖 4-27 底床紊流動能(k)俯視圖(y = -1.999 m) ... 58

圖 4-28 流場紊流動能(k)側視圖(z = 0.3 m) ... 58

圖 4-29 圓柱前緣與後端沖刷深度 ... 60

圖 4-30 圓柱前緣侵蝕深度變化 ... 60

圖 4-31 圓柱後端侵蝕深度變化 ... 60

表目錄

表 3-1 潛變流場(Re = 0.1)模擬配置(dt = 0.1 s) ... 29

表 3-2 流場(Re = 10)模擬配置(dt = 0.1 s) ... 34

表 3-3 流場(Re = 30)模擬配置(dt = 0.05 s) ... 35

表 3-4 流場(Re = 100)模擬配置(dt = 0.05 s) ... 36

表 3-5 流場(Re = 1000)模擬配置(dt = 0.1 s) ... 37

表 4-1 動床(Re = 46000)模擬配置(dt = 0.001 s) ... 41

Chapter 1 緒論

1.1 前言

近年來，國人的環保意識逐漸抬頭，自日本發生 311 大地震所引起的核災之 後，國內反核的聲浪日漸嚴重，台灣地狹人稠，廣設核電廠對安全來說的確是一 個很大的隱憂，但若要減少核電，就必須額外增加其他生產能源的方式，現今能 源主要由化石燃料與核能發電提供，化石燃料有限且產生能量的火力發電過程會 排放大量二氧化碳，加劇溫室效應，故乾淨的可再生性能源近年來逐漸受到人們 的重視。

台灣位於亞熱帶區域，夏季受颱風與西南氣流吹拂，冬天則受東北季風的影 響，風能資源十分豐富，非常適合發展風力發電，但風力發電也是有其缺點，像 是附近會產生大量噪音，而且靠近內陸的部分由於地表摩擦力的關係，所能擷取 的風能有限，故近幾年，設置在較外海的離岸(off shore)風力發電有著顯著的進展，

我國政府也在民國99 年開始了離岸風電廠的初步規劃，預計會在民國 104 年於苗 栗竹南、彰化芳苑完成示範組機台，在海上建立離岸風力發電機組有很多，像是 不佔用土地，且由於擾流小所以風能更加穩定豐富，對生態的衝擊也相對較小，

離岸風力發電已是現在已開發國家的趨勢，若能善用台灣這得天獨厚的天然資 源，在化石能源日漸縮減的未來，相信會扮演很重要的角色。

發展離岸風力發電，並不是一個簡單的工程，在建立之前，必須先考慮不同 地點的風能密度以及穩定性，以尋找最適合的風場，台灣位於板塊交界地帶，故 其底座結構強度要能負荷未來不可知的強震，海底電纜的整體配置也很重要，生 態方面，必須考慮保育類動物的棲息地(中華白海豚)，以及須考慮風力發電所產生 的各種不同頻率的噪音與震動，對當地的浮游生物所或鳥類造成的影響 。台灣所 設立的離岸風力發電機組位於西部沙岸，容易受到沿岸海流及波浪的影響，故若 能在建立前先做一個數值模擬，了解各流況下的底床沖刷現象，便能延長風力發 電機組的壽命，減少額外的成本支出。

1.2 研究動機

當一穩態流流經一垂直圓柱，流況會產生大幅改變(如圖 1-1)，流體上游處，

流體由於受到圓柱減速的影響，會在圓柱表面產生一停滯壓力，且因流速的不同 導致壓力由上到下逐漸減小，故會在在圓柱上游方向會產生一個向下射流(down flow)，此向下射流沿著結構物往下移動碰至底床而改變方向，在底部形成一個三 維渦流，由於此三維渦流結構與馬蹄十分相像，我們就稱之為馬蹄形漩渦(horseshoe vortex)，若底床為砂石等可被侵蝕的物質，則會加速漂沙的傳輸作用，導致底床的 前緣部分遭到沖刷，馬蹄形渦流是造成局部沖刷的主因。另外由兩側流經圓柱後 的流場會產生渦流分離(vortex shedding)，分離出來後的不規則細部渦流可能會與 上游的馬蹄形窩流交互作用，在尾部產生一低壓中心，將部分的沖刷的泥沙帶至 下游處堆積，稱為尾跡渦流(lee wake vorticities)。

近年來，台灣開始了離岸風力發電計畫，海流流經離岸風力發電基座也會造 成上述所說的局部侵蝕現象，導致海底地形變動，故在建立離岸風力發電機組之 前，最好先利用數值模式做個事前評估，但在過去，許多模擬採用商業用的套裝 數值模式，其中內部設定的參數、假設都無法掌握且難以修改，無法模擬出較真 實的流況，而且較沒有著重在底部侵蝕堆積的模擬，僅以初步的渦度帶過，只能 用經驗公式做一個概略性的預估，故我們希望可以建立一套數值模式，能完整的 模擬底部砂石與流場之間交互作用的現象，並且探討局部侵蝕的細部機制，了解 沖刷的整體物理流程後，在未來才能依照不同情況，去改變配置並且做一些應用。

圖 1-1 局部沖刷機制的示意圖(Roulund et al.,2005)

1.3 文獻回顧

1.3.1 沉浸邊界法

沉浸邊界法(immersed boundary method)最早是由 Peskin(1977)為了解決模擬血 液流進人體心臟時所發表的。在模擬心臟跳動時，彈性的血管壁可能會產生形變 而造成上邊界設定變得十分複雜，使用傳統的方式會十分損耗資源，故 Peskin 則 在統御方程式中加上了一個假想的人工外力項去使邊界速度近似於該邊界之表面 速度，達到不用修改邊界條件也能使結果近似於加入邊界的效果，之後依循著 Peskin 的理論，各種不同類型的沉浸邊界法相繼被發表出來。

Goldstein et al.(1993)發展了虛擬邊界法，模擬二維環繞圓柱流場、三維平板與 紊流肋管。

Saiki & Biringen(1996)使用了力回饋(feedback-forcing)沉浸邊界法去計算低雷 諾數流(Re≤ 400)經靜止、轉動的水平震盪圓柱，該方法改良了 Goldstein 的虛擬邊 界法，透過內插流體在網格上的速度，並適當的分佈力在邊界點上，以消去使用 虛擬邊界法中回饋力項所造成的錯誤震盪，他們發現力回饋沉浸邊界法還能應用 在會移動的固體邊界上，但力回饋沉浸邊界法應用在複雜三維的情況下常常會發 生一些穩定性的問題，另一種新的沉浸邊界法因此而生。

Fadlun et al.(2000)提出了 Direct-forcing 沉浸邊界法，成功地運用在複雜三維 的流體模擬上，該方法是將沉浸邊界法裡的假想力項，以物理邊界速度與附近的 點以速度的形式，直接作內插求得，這個方法的好處是比較不會受到穩定性限制 的影響，而且不需要像力回饋沉浸邊界法使用一些經驗公式求得的常數作最佳化。

Tseng & Ferziger(2003)延伸 Fadlun(2000)的構想創造出了虛擬網格沉浸邊界法 (ghost cell immersed boundary method)，定義出虛擬點(ghost point)，希望可以藉由 較高階的線性內差得到更高的精度，並與使用 body-fitted 建立出來的波浪狀底部 網格做驗證，結果速度分布與使用body-fitted 的方式十分吻合，另外也將 GCIBM 應用在用在較大尺度的三維海洋模式MIT Global Circulation Model(MITGCM)中，

與傳統的階梯法(stair-step)做比較，GCIBM 較不容易受到網格解析度變化的影響。

1.3.2 底床沖刷

橋墩局部沖刷機制的研究實驗在過去已有相當多的人在研究，不過早期大多 都是做實驗，如 Dargahi(1989)提出了影響局部沖刷的主因是由於圓柱前緣的馬蹄 形渦流。並且將沖刷分為三個階段，第一個階段為初始期，並無明顯的沖刷型態，

底床由定床逐漸轉變成動床，第二個階段為主要期，沖刷作用十分劇烈，圓柱前 緣侵蝕深度隨時間有個週期震盪向上的趨勢，第三個階段為衰退期，整體沖刷逐 漸趨緩，底床的侵蝕型態達到一個穩態平衡。另外還探討了影響圓柱前緣局部沖 刷現象的各種不同因子，並且提出了一些減少局部沖刷現象的方式。

Olsen & Melaaen(1993)是第一個做出三維圓柱的沖刷模式模擬，底床的材質是 使用黏性泥沙，流體計算的部分是使用雷諾平均納維爾史托克模型(RANS)，漂沙 的傳輸方程式合併在沖刷計算模型當中，在圓柱前緣發現馬蹄形漩渦的現象，且 模擬出來的沖刷結果與實驗也大致相同，Olsen & Kjellesvig(1998)延續了之前的模 擬，用了相同的方程式，做出來的模擬不僅限於初始沖刷階段，還包含了完整的 沖刷過程，得到的沖刷深度計算與經驗公式得到的沖刷深度吻合。

Melville & Chiew(1999)針對橋梁底部沖刷做了許多實驗去探討平衡沖刷深度 (定義為達至穩態後沖刷的最大深度)與給其定的物理量之關係，其中歸納出幾個結 論，像是(a)僅需 10%的平衡時間，沖刷深度即可達到平衡深度的 50%至 80%，這 主要取決於流場的強度。(b)清水沖刷橋墩至掏刷平衡，總共所需的平衡時間主要 受流場強度( / )V V 、流場深度_c ( / )y D 與泥沙粗糙度( /D d₅₀)給影響。(c)歸納出平 衡時間的的經驗公式如下，在水深超過一定範圍之後，即不受水深的影響。

( ) 48.26 ( 0.4) 6

e

c

D V y

t days

V V D

= − > (1.1)

( ) 30.89 ( 0.4)( )0.25 6

e

c

D V y y

t days

V V D D

= − ≤ (1.2)

Richardson & Panchang(1998)也做了流體流經一垂直圓柱的三維模擬，他們做 了三種不同情況的測試，分別是平板底床、沖刷過程與平衡的沖刷過程，沖刷過 程的邊界是採用Melville & Raudkivi’s (1977)實驗所測量出邊界條件，計算紊流的 部分使用了普朗托混合長度理論(Prandtl’s mixing-length theory)、黏滯渦流模型

(eddy-viscosity model)兩種不同方程式的k−ε 模型，結果成功的模擬出穩態流與馬 蹄形渦流的現象。

Tseng & Song(2000)是第三個做出三維結構物沖刷模式的團隊，他們使用大渦 流 模 式(LES)成功的模擬出圓柱前的馬蹄流與圓柱之後的渦流分離，並且與 Dargaji(1989)的實驗做驗證，雖然有討論流體在侵蝕過程所受到的影響，但是卻沒 有模擬底層高度受沖刷所產生的變化。

Roulund, Sumer, Fredsoe & Michelsen(2005)建立了一套以 RANS 模型為基礎的 沖刷模式，去探討圓柱前緣的馬蹄形渦流與其他的物理量之關係，結果發現馬蹄 形渦流的規模會隨著邊界層厚度與直徑的比值( / )δ D 大小產生關係，且當雷諾數 (Re )_D 愈大時，底床剪應力以及馬蹄形渦流的大小也會隨著變大，並與實際實驗做 比較如下圖1-3(a)，成功的模擬出圓柱前緣與後端的侵蝕堆積型態圖 1-3(b)。

Paik, Escauriaza & Sotiropoulos(2007) 使 用 分 離 渦 流 模 式 (detached-eddy simulation)模擬高雷諾數翼體結構前緣的馬蹄形渦流分離現象，該模式混合了 RANS 與 LES 的特點，模擬出來的結果則使用 POD 法去分析，並且以 q-criterion 表示圓柱前緣渦流結構隨時間的變化，觀察出渦流結構主要有兩種基本狀態，第 一種是圓柱前緣上方主要的馬蹄形渦流，另一個就是包覆著馬蹄形渦流的外部混 亂髮夾渦流(hairpin vortices)，最後認為髮夾渦流的成因，是由於馬蹄形漩渦所造成 不穩定性的結果。

Chou & Fringer(2010)利用質量守恆與傳輸方程在曲面座標系統上建立出移動 網格模式(Arbitrary Eulerian-Lagrangian scheme)，該模式的特點在於移動網格上的 物理量仍能保持高階精度，且質量守恆與傳輸方程式的離散化精度必須相同，否 則會不滿足濃度守恆導致結果發散。

Chou & Fringer(2010)應用了先前的移動網格，以大渦流模式為基礎建立了一 套模擬紊流邊界層底床演進的三維模式，與過去模擬底床形變的數值模式不同之 處是對其臨界剪應力作一個流速與邊界的角度修正，以及包含所受到重力的影 響，模擬出來一個平坦的底床，在經過像波浪一樣的週期震盪速度邊界，演進至 沙漣的過程，與實驗比較結果十分吻合。

圖 1-2 沖刷實驗圓柱後緣高度變化(Dargahi, 1990)

圖 1-3 (a)圓柱沖刷實驗與(b)模擬結果(Roulund et al.,2005)

圖 1-4 底床形貌與濃度的演進(Chou & Fringer, 2010) (a) (b)

Chapter 2 理論背景與方法

我們想了解沖刷過程中的細部現象，但風機底座的實驗規模實在太大，所以 想 先 建 立 數 值 模 式 來 進 行 模 擬 ， 流 體 計 算 的 部 分 我 們 使 用 的 是 大 渦 流 模 式 (Large-Eddy Simulation)，該 LES 模式原是由 Zang(1993)在美國史丹福大學環境流 體力學實驗室所發展出來，此模型為一個三維的水動力模式，支援平行運算，且 可以模擬漂沙懸浮在水中的濃度與泥沙傳輸，我們會再2-1 節介紹大渦流模式的特 點以及這套模式所使用的離散方法與流程。

在底床邊界的部分使用移動網格法(Arbitrary Lagrangian-Eulerian)，該模式是由 Chou & Fringer(2010)將濃度守恆與傳輸方程式結合大渦流模式，能計算出網格的 相對速度，模擬底床受到流體沖刷或是堆積所造成形貌變化的過程，我們會再2-2 節說明泥沙傳輸方程式和掏刷機制，以及所使用的經驗公式。

在 處 理 較 複 雜 橋 墩 邊 界 的 部 份 是 使 用 了 沉 浸 邊 界 法(immersed boundary method)，好處是能夠處理複雜邊界且不需要增加太多的計算量，我們在第 3 章會 詳細介紹所使用的沉浸邊界法，其中離散化和整體數值模式的計算流程，並且在 後面做一個穩態流流經圓柱的範例，藉由改變邊界速度造成雷諾數的變化，去觀 察圓柱後方流場的現象，與實際的實驗比較作驗證。

2.1 大渦流模式

一般來說，計算流體力學隨著離散方法主要分為三種不同的模型。

第一種為直接數值模擬(Direct Numerical Simulation)，此方法是直接捕捉所有 相關尺度的亂流運動，故不須對最小的尺度建立其他模型，結果也是最準確的，

但是由於流場中所有微小擾動都會做計算，故極耗費運算資源，大多只能做到實 驗室尺度的分析。

第二種為雷諾平均納維爾史托克方程(Reynolds Averaged Navier-Stokes)，該方 式是假設流場中的物理量(B)可以表示為一個時間平均值與變動量的和如下式 (2.1)，將此物理量帶入納維爾史托克方程式，可得到統計平均物理量的統御方程式 (2.2)、(2.3)、(2.4)，即為 RANS 方程式，其中動量方程式裡的非線性項會產生一個 雷諾應力(Reynolds Stress)，為了要模擬這個二階張量，通常是使用k−ε 模型求 解，但由於時間平均假設的關係，RANS 模型求得的為時間平均物理量，當雷諾數 變大，流場中的擾動性變強時，許多紊流現象RANS 模型並不能捕捉到。

第三種為大渦流模式(Large-eddy Simulation)，一開始是以所設的網格尺度對渦 流去做過濾，分成主要網格與次網格(subgrid scale)，LES 模式計算主要網格，次 網格的部分是使用其他模型去做模擬，在帶入其方程式中求解，由於 LES 並沒有 作時間平均，故在高雷諾數中，其較細微的紊流現象較RANS 模型來的準確。

前面有提到，局部侵蝕的主要原因是由於結構物前面的馬蹄形漩渦，故紊流 現象在我們的模擬中十分重要，所以我們使用大渦流模式去計算流場。

( , , , ) ( , , ) '( , , , )

B x y z t =B x y z +B x y z t (2.1)

2 2

1 _x [ ( ') ( ' ') ( ' ')]

Du P

g u u u v u w

Dt x ν x y z

ρ

∂ ∂ ∂ ∂

= − + + ∇ − + +

∂ ∂ ∂ ∂ (2.2)

2 2

1 _y [ ( ' ') ( ') ( ' ')]

Dv P

g v v u v v w

Dt y ν x y z

ρ

∂ ∂ ∂ ∂

= − + + ∇ − + +

∂ ∂ ∂ ∂ (2.3)

2 2

1 _z [ ( ' ') ( ' ') ( ') ]

Dw P

g w w u w v w

Dt z ν x y z

ρ

∂ ∂ ∂ ∂

= − + + ∇ − + +

∂ ∂ ∂ ∂ (2.4)

2.1.1 統御方程式

為了要滿足三維、非穩態的不可壓縮流，我們列出了不可壓縮流的質量守恆 方程式、納維爾-史托克方程式，以及濃度傳輸方程式，其中式(2.6)中的非線性項 已將連續方程式帶入化簡。

j 0

j

u x

∂ =

∂ (2.5)

2 2

( ) 1

i i

j i i

j i j j

u p u

u u g

t x x δ ν x x

ρ

∂ + ∂ = − ∂ − + ∂

∂ ∂ ∂ ∂ ∂ (2.6)

[( _{j s i}2) ] 0

j

C u w C

t x δ

∂ + ∂ − =

∂ ∂ (2.7) 其中 i,j=1,2,3

u 為卡氏座標速度 i

p 為流場壓力

ν = μ

ρ為運動黏滯係數 C為泥沙濃度

w 為顆粒的下沉終端速度 s

重力加速度的方向為x (-y)的方向，如下圖所示 ₂

圖 2-1 卡氏座標方向 x y

g

z

為了簡化統御方程式，包氏近似(Boussinesq approximation)常用於簡化流體動 量守恆的密度變化量，流體密度ρ可分成一個假想的靜止密度ρ₀^{以及相對擾動密} 度ρ'^。

0 '

ρ ρ ρ^{= + (2.8)}

壓力p 可分成一個假想的參考靜壓p 和相對擾動壓力 '₀ ρ

0 '

p= p +p (2.9)

此外靜壓梯度與參考密度的關係為

0 0

p ρ g

∇ = 

(2.10)

將式(2.8)、式(2.9)與式(2.10)帶入連續方程式(2.5)與動量守恆式(2.6)中，可得到包 氏近似假設的統御方程式如下

j 0

j

u x

∂ =

∂ (2.11)

2

0 2

0 0

( ) 1 ( )

i i

j i i

j i j j

u p g u

t x u u x ρ ρ δ ν x x

ρ ρ

∂ + ∂ = − ∂ − − + ∂

∂ ∂ ∂ ∂ ∂ (2.12)

2.1.2 濾波統御方程式

在大渦流模式中，我們將物理量對網格尺度過濾分成主要網格(large-scale) 與次要網格(subgrid-scale)兩個不同項去做處理

φ φ φ= + (2.13) '

其中經過濾波函數G 過濾後的主要網格物理量定義為

3

1 2 3 1 1 2 3 1 2 3

( , , ) _{D i}( , ') ( ',_{i i} ', ') ' ' ' x x x i G x x x x x dx dx dx

φ φ

=



j 0

j

u x

∂ =

∂ (2.15)

2

0 2

0 0

1 ( ) ^ij

i i i

j i

j i j j j

u u p g u

t u x x x x x

ρ ρ δ ν τ

ρ ρ

∂ + ∂ = − ∂ − − + ∂ +∂

∂ ∂ ∂ ∂ ∂ ∂ (2.16)

( _{j s i}2 ) ^j

j j

C u C w C

t x x

δ ^∂χ

∂ + ∂ − = −

∂ ∂ ∂ (2.17)

其中 u_j 表示在 j 方向的過濾後速度分量 p 為修正後的壓力項

從濾波後的統御方程式，會發現多了兩項額外的物理量，分別為次網格應力(SGS) 與次濾波通量(SFS)

i j ij u ui j u u

τ = − 為次網格尺度的雷諾應力(subgrid-scale stress)

j u C u Cj j

χ = − 為次濾波通量(subfilter-scale flux)

在我們的研究裡，次網格尺度的雷諾應力項中，包含了被過濾掉的局部小尺 度渦流，會根據局部速度場所調整其所受到的影響，我們是使用動態混合模型 (dynamic mixed model)去模擬次網格應力(SGS)，與傳統的 Smagorinsky model 相 比，動態混合模型利用了紊流尺度中自我相似(self-similarity)的特性，較不會產生 過度消散(dissusive)的效應。

2.1.3 曲線座標格式

在處理複雜幾何邊界的條件下，曲線座標( , , )ξ ξ ξ_{1 2 3} ^{相較卡氏座標}( , , )x x x 具1 2 3

有較好的適應性，我們可以對卡氏曲標與曲線座標做一個簡單的連鎖率座標轉換 如下

m

j j m

x x

ξ ξ

∂ = ∂ ∂

∂ ∂ ∂ (2.18) 如果流體為不可壓縮流，且運動黏滯係數為常數，我們將座標轉換方程式帶 入質量守恆(2.15)、動量守恆(2.16)與漂沙傳輸方程式(2.17)可得

0

m m

U ξ

∂ =

∂ (2.19)

1 1 1 0

2

0 0

,

( ) ( ) 1

( )

( )

i m i

ij i

m m m

mn i

i m

m n

J u U u p

J J g

t

G u T

δ ρ ρ δ

ξ ξ ρ ξ ρ

ξ ν ξ

− − −

∂ +∂ = − ∂ ∂ + −

∂ ∂ ∂ ∂

∂ ∂

+ −

∂ ∂

(2.20)

[ ( _{m m i}2)]

m m

C U W

C F

t

δ

ξ ξ

∂ −

∂ + = − ∂

∂ ∂ ∂ (2.21) 其中 根據愛因斯坦求和約定 i, j, k, m, n 1, 2,3=

1 m

m j

j

U J u

x ξ

− ∂

= ∂ 為轉換體積通量

1 det( ⁱ )

j

J x

ξ

− = ∂

∂ Jacobian 轉換，可視為網格的體積

1

mn m n

j j

G J

x x

− ∂ξ ξ∂

= ∂ ∂ 為目偏張量

1

, ^m( )

i m i j i j

j

T J u u u u

x ξ δ

− ∂

= − 為為座標轉換的次網格應力

1

2

m m s

W J w

x ξ

− ∂

= ∂ 為座標轉換後的終端沉降速度

1 ^m( )

m j j

j

F J Cu Cu

x ξ δ

− ∂

= − 為座標轉換的次網格應力漂沙通量

ρ0 ^{為純水密度}

2.1.4 流體離散化計算流程

我們解流體所使用的數值方法為分步法(fraction step method)，首先將前面推得的包 氏近似統御方程式，分成兩個部分，簡單離散如下兩式

** 2

( )

n

i i i

j i

j j j

u u u

t x u u ν x x

− = − ∂ + ∂

Δ ∂ ∂ ∂ (2.22)

1 **

0 n 1

i i

i

u u p

t ρ x

+ − = − ∂

Δ ∂ (2.23) 其中式(2.22)，由於還沒有滿足其質量守恆，稱之為預估式(predictor)，先求出預估 的速度，再將其帶入校正式(corrector) (2.23)，求出下一個時階的速度。

我們將式(2.22)也分解成兩步來做，先解完對流項(convective term)如式(2.24)，再去 解黏滯項(viscous term)如下式(2.25)

* n

i i i

j j

u u u

t u x

− = − ∂

Δ ∂ (2.24)

** * 2

i i i

j j

u u u

t ν x x

− = ∂

Δ ∂ ∂ (2.25) 式(2.24)我們對其使用 QUICK 作空間離散，為一個三階精度的離散方法，時間離 散則是使用Adams-Bashforth，為一個二階精度的離散方法

*

3 1 1

( ) ( )

2 2

n

n n

i i

j i j i

j j

u u

u u u u

t x x

− = − ∂ + ∂ −

Δ ∂ ∂ (2.26) 式(2.25)使用中央差分作空間離散，時間離散則是使用 Crank-Nicolson

** *

1 1

2

i j 2 i j i j

i i

j

u u u

u u

t ν ^{+ −} x ^{+ −}

− =

Δ Δ (2.27)

* * * ** ** **

** *

1 1 1 1 1 1

2 2

2 2

1 ( )

2

i j i j i j i j i j i j

i i

j j

u u u u u u

u u

t ν ^{+ −} x ^{+ − + + −} x ^{+ − −}

− = +

Δ Δ Δ (2.28)

其中u_i 代表的是速度的方向

1

ui j₊ 代表的是u 在座標_i j+1的速度

經上面離散計算後可以求出二階精度的u ，但是要帶入式(2.23)時，卻發現還沒有_i^**

壓力的值，故我們對上式(2.23)做散度，可得

1 **

2

0

( ) 1

n

i i

u u

t ρ P

+ −

∇ • = − ∇

Δ (2.29)

且為了要滿足質量守恆定理 即∇ •u_iⁿ⁺¹= ，帶入式(2.29)可得 0

** 2

0

1 1

ui P

t ρ

− ∇ • = − ∇

Δ (2.30)

可改寫成

2P 0 ( u_i**) t

∇ = ρ ∇ •

Δ (2.31)

其中右邊的∇ •u_i^**可由前面求得之預估式算出

式(2.31)為一波松方程式(Poisson Equation)，計算其壓力分佈後代入式(2.23)，即可 求得下一步的速度。

2.2 移動網格與泥沙傳輸方程式

我們所使用的泥沙傳輸模式是由 Chou(2010)所建立，要模擬出整體泥沙傳輸 的現象，分好幾個步驟去執行。2.2.1 節定義出泥沙濃度，流場的密度可以泥沙濃 度的形式表示。2.2.2 節是在介紹泥沙的濃度傳輸方程式，要滿足濃度質量守恆以 及泥沙顆粒重力的沉降作用。2.2.3 節是在介紹處理底床邊界的移動網格法，在掏 刷模擬中，底床的邊界會隨著泥沙濃度的變化而移動，在這小節我們會說明處理 底床的統御方程式，給予底床網格速度，使其網格滿足泥沙濃度的移動，模擬出 底床形貌的變化。2.2.4 節介紹的是掏刷方程式，底床掏刷量是一個非常複雜的問 題，現今仍沒有完整的解析解，故我們是使用經驗公式去描述底床掏刷量，並且 依照沖刷的角度去對原先的經驗公式進行修正。2.25 節利用上一小節所求得的底 床掏刷量，推導出高度方程式，在這裡為了要滿足顆粒所受到的重力影響，再對 底床高度方程式中增加了一個擴散項做修正，使整體模型更符合實際現象。

2.2.1 密度分層

在模擬懸浮泥沙的部分，為了簡化計算，我們將懸浮泥沙視為流體，以尤拉 法(Eulerian formulation)的觀點去處理顆粒濃度，整體的密度可以流體與顆粒的體 積百分濃度表示為

(1 )

liquid C sC

ρ ρ= − +ρ (2.32)

其中 ρ_liquid 為液體之密度

C 為懸浮泥沙的體積百分濃度

ρs 為泥沙密度(通常為 2650 kg-m ) ^-3

2.2.2 泥沙傳輸方程式

泥沙傳輸方程式(Sediment transport equation)通常是使用上述的尤拉法，以濃 度的形式去表示，但尤拉法在實際應用上，僅能適用於體積百分濃度尺度為O(0.01) 左右的濃度模擬，在底層區域高泥沙濃度的部分卻不適用，實際上高濃度泥沙顆 粒與水的二相耦合僅能使用拉格朗日描述法(lagrangian description)去做模擬，但在 過去並沒有人做過較大尺度的泥沙懸浮模式，而且泥沙與水之間的交互作用機制 實在太過複雜，一直都是個未解決的難題，所以我們要先解出底床的高度方程式 套用在移動網格上(Arbitrary Eulerian-Lagrangian)，盡量使泥沙傳輸方程求解的流場 大部分的區域泥沙濃度皆小於O(0.01)，以便達成我們一開始的假設。

在經過空間濾網過濾後的懸浮濃度場的泥沙傳輸方程式為

[ ( _{m m i}2)]

m m

C U W

C F

t

δ

ξ ξ

∂ −

∂ + = − ∂

∂ ∂ ∂ (2.33)

其中 ¹

2

m m s

W J w

x

− ∂ξ

= ∂ 為座標轉換後的終端沉降速度

1 ^m( )

m j j

j

F J Cu Cu

x ξ δ

− ∂

= − 為座標轉換的次網格應力漂沙通量

2.2.3 移動網格法

在掏刷問題中，底部的邊界並不是固定而會隨著泥沙的濃度而改變，上一章 有提到我們所使用的尤拉法泥沙濃度假設，故我們希望所模擬的流場裡的網格能 自由移動且明確的區分出底部邊界，故我們使用由 Chou(2010)建立的移動網格法 (Arbitrary Lagrangian-Eulerian scheme)去給予座標中的網格速度，使其會隨著濃度 變化而做出位移或形變，其中推導如下

首先我們先列出曲線座標的濃度守恆方程式

( _j) _c

j

C Cu S

t ξ

∂ + ∂ =

∂ ∂ (2.34)

若以尤拉法的觀點去考慮網格的移動則可將上式化為

1

( ) ( _m) ( _{g m}, ) _c

m m

J C CU CU S

t ⁻ ξ ξ

∂ + ∂ − ∂ =

∂ ∂ ∂ (2.35)

其中 ¹ det( ⁱ )

j

J x

ξ

− = ∂

∂ jacobian 的倒數，物理意義為體積轉換率

1 m

m j

j

U J u

x ξ

− ∂

= ∂ 代表座標轉換後的絕對體積通量

1 1 ,

, ,

m m g j

g m g j

j j

U J u J dx

x x dt

ξ ξ

− ∂ − ∂

= =

∂ ∂ 代表座標轉換後的網格體積通量

S 為濃度質量所受到壓力項、黏滯項、重力項等其他因素的影響 c

若流體為不可壓縮流( ^m 0)

m

U ξ

∂ =

∂ ，則式(2.32)可化簡為

1 ,

( ) ^{g m} 0

m

J U

t ⁻ ξ

∂ −∂ =

∂ ∂ (2.36)

由式(2.35)以及式(2.36)，即可計算出其中網格的體積通量所對應的濃度傳輸方 程式。

2.2.4 掏刷方程式

影響掏刷的因素有很多，但最主要的因素還是由於流場帶給底床顆粒的剪應 力，Shields(1936)提出希爾斯參數來量化底部砂石所受到水捲起來的難易程度

希爾斯參數

0

( )

( 1) shields parameter b

s gd θ τ

= ρ

− (2.37) 其中τ_b =ρC U_{D h}² 為底床剪應力

2.65

s= 為泥沙的比重

d 為泥沙顆粒直徑 0

事實上希爾斯參數即為無因次化的底床剪應力，表示為底床剪應力與泥沙顆 粒所受到的重力比值，底床剪應力愈大時，顆粒愈容易被水給捲起來，故希爾斯 參數大於臨界希爾斯參數時，底層的砂石會被帶動捲起至流場中，我們使用 Van Rijn (1993)所提出的掏刷方程式 (pickup function)如下來描述上面的現象。

* *

0 0

( 1)

k D T

P s gd

β γ

α 

=  

−   ^c otherwise

θ θ>

(2.38)

θc 為臨界希爾斯參數

2 2

1 2

Uh = u +u 為水平速度之大小(u u₁, ₂為x x 方向之速度) ₁, ₂

1 0 2

0

[ ln(1 )]

D

z z

C κ z ⁻

= + 為底部紊流邊界層的阻力係數

1

* 2 3

0[( 1) / ]

D =d s− g ν 為無因次粒徑

* ( _c) / _c

T = −θ θ θ

κ =0.41為von Karman 常數

z 為底層粗糙度 0

z 為底床與最底層網格中心的距離 1

0.00033

α = 0.3β = 1.5γ =

2.2.5 底床高度方程式

我們將注意力放在底床的形變上，底床高度的變化主要是因為流場中帶來的 泥沙沉積以及剪應力侵蝕所導致，以上的現象可以寫出最簡易的底床高度方程式 如下

1 1

(1 ') _B h _B ( _{s b k}) p J J w C P

t

− ∂ −

− = −

∂ (2.39) 其中w 為泥沙的沉降速度 _s

P 為掏刷方程式(pickup function)所求出來的泥沙濃度通量 k

'

p 為泥沙的孔隙率

但底層形貌不單只受剪應力的影響，泥沙顆粒被侵蝕或堆積後，所呈現的形 貌變化，必須考慮重力對其所產生的影響(gravity effect)，原始所求得的臨界希爾 斯參數並沒有考慮重力影響，故我們先列出 Soulsby(1997)做實驗所求得的臨界希 爾斯參數(θ_c,0)^{經驗公式如下}

*

,0 *

0.3 0.055[1 exp( 0.02 )]

1 1.2

c D

θ = D + − −

+ (2.40)

再對原始的臨界希爾斯參數做一個角度的修正，求得修正後的臨界希爾斯參數( )θ_c

,0

sin( ) sin( )

c rp

c rp

φ φ θ

θ φ

= + (2.41)

圖 2-2 底床顆粒傳輸示意圖 (Chou & Fringer, 2010)

其中φ_rp 為靜止底床角度 φ ^{為底床角度}

在二維座標下的底床角度可表示為 tan ¹dh φ = ⁻ dx

但在三維座標下，考慮底部剪應力的方向影響的底床角度可表示為

2

2,

1, 1 2, 2

2 1,

sin sin sin( )

bed

bed bed

bed

u u

u u

φ φ

φ = ⁺

+

其中φ₁^和φ2分別為水平方向速度的斜角

另外還要考慮底床角度( )φ ^{大於靜止底床角度(}φrp⁾時，一部分的泥沙會受重力 作用導致堆積在附近，我們可以用以下的方程式去表示

,

mn

g m B

n

Q kG h ξ

= ∂

∂ (2.42) 其中Q_{g m,} 為重力所引起的體積流率

1

, , 1,2

mn m n |

B B m n j

j j

G J

x x ξ ξ

− ∂ ∂ =

= ∂ ∂

將原本推導出來的底床高度方程式(2.34)，加入上面對重力項所做的修正，我 們可以得到最終的底床高度方程式

1 1

, 1,2

(1 ') _b ( _b^mn ) |_{m n b} ( _{s b k})

m n

h h

p J kG J w C P

t ξ ξ

− −

=

∂ ∂ ∂

− = + −

∂ ∂ ∂ (2.43)

其中當底床角度( )φ ^{大於靜止底床角度(}φrp⁾時，這是一個擴散方程式(diffusion equation，主要的效果是會平滑整體侵蝕與堆積的高度。

0

rp k

k rp

α φ φ φ

−

 

 

=  

 

 

^rp otherwise

φ φ>

(2.44)

我們考慮重力項去修正臨界希爾斯參數( )θ_c ，使得原本的經驗公式更貼近了實 際流場流動時，顆粒堆積角度所帶給掏刷方程式的影響。另外在底床高度方式中 加入了重力項的修正，定義出靜止底床角度(φ_rp)，描述當侵蝕角度過大時，會由 於重力的作用導致部份的泥沙傳輸至附近，平滑較極端的侵蝕曲線，較能模擬出 真實底床的侵蝕現象。

圖 2-3 靜止底床角度示意圖(Chou & Fringer, 2010)

Chapter 3 沉浸邊界法介紹與驗證

在計算流體力學中，複雜幾何邊界的處理一直是個很大的問題，傳統方式是 使 用 貼 體 法(body-fitting) ， 產 生 符 合 複 雜 幾 何 之 曲 線 網 格 或 非 結 構 性 網 格 (unstructured grid)，或是卡氏座標網格藉由切割成較細的解析度去貼近複雜幾何邊 界，再修改其邊界條件使其滿足，但缺點主要是當邊界過於複雜，貼體法往往需 要在網格的精細度上耗費大量資源，且若邊界會隨著時間改變位置時，每次的計 算量更是驚人，故我們使用了沉浸邊界法(immersed boundary method)去處理較複雜 的邊界問題，沉浸邊界法主要是在統御方程式中假想出一個人工外力項，去對邊 界的動量方程式做修正而使其滿足邊界條件，故無需增加網格或是對網格作座標 轉換，即可在卡氏座標中模擬複雜邊界或移動邊界之流場問題。

我們所使用的沉浸邊界法，是來自Tseng & Ferziger(2003)提出的虛擬網格沉浸 邊界法(ghost cell immersed boundary method)(GCIBM)，主要的優點就是較簡單使 用，且已經驗證過可以應用在各種不同尺度的流體計算模式之中，處理邊界速度 的部分是使用多項式作反矩陣找出虛擬點(ghost point)上的速度內插方程，在帶入 該點所定義的座標求出速度，若要增加所求的速度精度階數，也較容易執行，在 本研究中所使用的是一階線性內插模型。

3.1 沉浸邊界法離散化

首先我們先列不可壓縮流的連續方程式(3.1)，以及經過包氏近似化簡後的納維 爾-史托克方程式(3.2)，為了減少計算量又要符合結構物的邊界條件，我們在固液 邊界網格中的統御方程式加入了假想力項去對原先的速度做平衡如式(3.3)。

j 0

j

u x

∂ =

∂ (3.1)

2

0 2

0 0

( ) 1 ( ) in other point

i i

j i ij i

j i j j

u p g u

t x u u x δ ρ ρ δ ν x x

ρ ρ

∂ + ∂ = − ∂ − − + ∂

∂ ∂ ∂ ∂ ∂ (3.2)

2

0 2

0 0

( ) 1 ( ) in ghost point

i i

j i ij i i

j i j j

u p g u

u u f

t x x δ ρ ρ δ ν x x

ρ ρ

∂ + ∂ = − ∂ − − + ∂ +

∂ ∂ ∂ ∂ ∂ (3.3)

假想力為一個位置與速度的函數，主要是用來平衡邊界速度，為了定義出所 須施加的網格座標點，我們將網格分成三類如圖 3-1(a)，第一類為固體網格(inner point)，第二類為固液網格(ghost point)，第三類為液體網格(outer point)，在固液網 格中我們加入假想力去對其速度做修正，使其邊界滿足無滑移條件如圖3-1(b)，離 散化後的動量方程式如下

1

n n

i i

i i

u u

RHS f t

+ − = +

Δ in ghost point (3.4) 其中RHS 包含對流項、黏滯項與壓力梯度項 _i

故我們可將假想力視為外插速度的函數，反推出假想力，表示如下

1

n n

i i

V u f RHS

t

+ −

= − +

Δ in ghost point (3.5)

其中Vⁿ⁺¹為固液網格中利用附近流場速度外插出的值

(a) (b)

圖 3-1 (a)網格分類表示圖與(b)內插速度示意圖

(a)點 1 (咖啡色區域)代表為固體網格(inner point) ，點 2 也就是黃色(固體)與淺藍色(液體)區域代表 為固液網格(ghost point)，點 3(藍色區域)代表為液體網格(outer point)

(b)表示為流體流經一平板，我們藉由修正 ghost point 的速度(下方黃點內的速度)使其邊界(紅點)滿 足無滑移條件(no-slip condition)

定義出所需外插速度的網格後，我們可將座標上的物理量視為一階線性分布，

表示為

0 1 2

a a x a y

φ = + + (3.6) 其中φ可代表座標點上的任意物理量

Ba= (3.7) φ a B^{= −}1φ (3.8) 在這裡我們是使用最簡易的一階線性內插方式，並且以二維的形式表示，先 找出ghost point 最鄰近的兩個液態網格座標點(不同方向)，以及該網格固體液體邊 界切點，如圖3-2 的x 、₁ x 與₂ x ，即可列出矩陣 B 之值。 ₀

0 0

1 1

2 2

1 1 1

x y

B x y

x y

 

 

=  

 

 

(3.9)

圖 3-2 內插點示意圖

x1為離ghost point 最近的 y 方向座標液體網格點(網格中心)，x₂為離ghost point 最近的 x 方向座 標液體網格點(網格中心)，x₀為該ghost point 上的邊界座標切點(紅點)，有了這三個點後即可內插 出ghost point 上的速度(網格中心)

對 B 取反矩陣後帶入式(3.8)，求得該流場中的線性內插參數 a，帶回式(3.7) 即可求得該網格中心所代表的速度V_iⁿ⁺¹(由於我們的網格是屬於非交錯性網格 (non-staggered)，故速度是分布在網格中心，即圖 3-2 的中心 G 點)。

我們所使用的計算流程，跟傳統的 GCIBM 的流程稍微有點不同，傳統的方 式，壓力項是分兩步再做，求完預估式(predictor)之後，直接內插出速度(V ) ，ⁿ⁺¹ 改變ghost point 的速度後，再帶入修正式(2.31)解波松方程式取得壓力差( )φ ^以滿足 質量守恆定理，但是我們的模型中壓力項是一次求解，為了要滿足沉浸邊界法的 核心，也就是用下一秒的速度資料內插進ghost point 上的速度，將 fraction step 分 兩次求解，計算第一次的速度資料內插出 ghost point 速度再將其儲存，重新返回 該時階的初始速度，但 ghost point 上的速度則用先前求出的內插速度取代，改變 了 ghost point 上的速度後，再繼續對其對流項、黏滯項與壓力項做平衡，第二次 求得的才是真正的速度，兩次解完後才進入下一個時階。

沉沒邊界法計算流程

1. 將網格分成三類，找出需施加假想力，也就是固體與液體交會的的網格區域 (ghost point)。

2. 計算出第一次求出的流場速度 (corrector)。

3. 取得 ghost point 外部的網格座標、座標上的速度以及邊界切點座標資料，再將 每個ghost point 上的內插速度給求出。

4. 回到一開始的速度，將第 3 步求出的內插速度帶進 ghost point 中。

5. 再次計算流場速度(corrector)

6. 更新速度場，進行下一個時階(u_iⁿ⁺¹)

(1)

(2)

(3) (14)

(4)

(5)

(6)

(7)

(8)

(9) (10) (11)

(12)

(13)

圖 3-3 整體數值計算流程 初始化網格、流場

將網格分成三類

(inner point,outer point,ghost point)

計算流程開始

更新速度

1

n n

u =u ⁺

求解掏刷方程式、底床高度方程式 並且更新網格

* *

0 0

( 1)

k D T

P s gd

β γ

α 

=  

−  

c

otherwise θ θ>

1

, 1,2 1

(1 ') ( ) |

( )

mn

b b m n

m n

b s b k

h h

p J kG

t

J w C P

ξ ξ

−

=

−

∂ ∂ ∂

− =

∂ ∂ ∂

+ −

0

rp k

k rp

α φ φ φ

−

 

 

=  

 

 

^rp otherwise

φ φ>

求解出校正項之速度

1 **

0 n 1

i i

ij i

u u p

t x δ

ρ

+ − = − ∂

Δ ∂

求解出ghost point 之內差速度

0 0

1 1

2 2

1 1 1

x y

B x y

x y

 

 

=  

 

 

a B= ⁻1φ

n+1

g g

X a=V X 為 ghost point _g

將初始速度與ghost point 速度重新代入

n save

u =u in inner point & outer point

1

n n

u =V ⁺ in ghost point 求解預估項之速度

** 2

( )

n

i i i

j i

j j j

u u u

t x u u ν x x

− = − ∂ + ∂

Δ ∂ ∂ ∂

將初始速度儲存

save n

u = u

帶入預估項之速度求解壓力場

2P 0 ( u_i**) t

∇ = ρ ∇ • Δ