利用沉浸邊界法於三維多球體沉降的數值模擬

(1)

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

Institute of Applied Mechanics College of Engineering

National Taiwan University Master Thesis

利用沉浸邊界法於三維多球體沉降的數值模擬

Development of an immersed boundary method to the numerical simulation of settling of multiple spheres

簡鈺霖 Yu-Lin Chien

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

中華民國 108 年 7 月

July, 2019

(2)

i

(3)

i

致謝

在碩士這兩年中，感謝周逸儒老師的細心指導，讓我在研究的這過程中不斷的成長，老師除了在研究上問題給予良好的建議，也在我研究犯錯的時候給予及時的改正和引導，並增進我面對問題時解決問題的能力，而最後讓我順利完成了碩士論文，並感謝在碩士論文口試時，口試委員牛仰堯教授和楊馥菱教授給予很好的意見且提出了不同的方法讓我參考並且改進。

在碩士研究的期間，在程式和研究上的問題，除了老師的指導，也感謝學長晨晏、佑任、哲榮和德耀的幫助及建議，且彼此相互討論而可以得到較好的解決辦法或結果，而在修課的期間，感謝同學國任、渙承、汶峰、良昇、忻融、裕愷、

威廷、宥羽、彥慶、德欣、毅倫和宗穎的陪伴，一起討論課業及互相幫忙彼此的不足，順利的完成修課階段。

感謝這兩年來一路陪伴我的師長和同學們，同時也感謝父母的支持，讓我能有這個機會來台大應力所完成碩士學位。

(4)

ii

中文摘要

本研究以拉格朗日-尤拉法的方式使用沉浸邊界法(immersed boundary method) 於直接數值模擬(direct numerical simulation)中，模擬多顆球體於流體中的運動行為。其目的主要是驗證並改善直接施力(direct forcing)沉浸邊界法[5]，並採用拉格朗日點向內收縮方法[8]，且同時考慮軟球及硬球碰撞模型，來模擬單球體和多顆球體於流場中的複雜運動行為。驗證和改善可區分為兩部分，第一部分是模擬穩態條件下靜止球體於均勻流場所受阻力並與前人模擬數據結果做比較。其中直接施力沉浸邊界法的力回授採用狄拉克正規化脈衝函數，此方法會增加阻力係數值，由本研究結果顯示球體直徑與網格大小的比值為8、雷諾數 Re = 50 時，模擬出來的阻力係數為1.805，比較文獻[24]模擬出來的阻力數值 1.574 來的大，故採用於拉格朗日點向內收縮的方法，此方法在不同流場網格解析度下，可以得到一Re、

內縮值和網格解析度的修正式，本研究結果指出在同樣解析度下，Re = 50 時可以得到的阻力係數為1.582，此值趨近於上述文獻模擬結果。而第二部分則是模擬多顆球體於流場中的運動行為，首先模擬和驗證單一球體自然沉降行為，而後模擬多顆球體於不同高度進行自由沉降，讓多個球體產生碰撞，最後參考前人文獻[8]

評估此不同條件下的運動行為。

關鍵字：直接數值模擬、沉浸邊界法、球體運動

(5)

iii

ABSTRACT

This study presents the direct numerical simulation with an immersed boundary method (IBM) to simulate the motion of multiple spheres in flow field. The purpose of this study is to verify and improve the direct forcing immersed boundary method [5]

and to adopt the inward contraction of Lagrangian points [8]; Simultaneously, we also consider the soft-sphere and hard-sphere collision model to simulate the collision behavior between sphere to sphere and sphere to wall respectively. The verification and improvement process consists of two parts. The first is to simulate the stationary sphere of steady state in the uniform flow and compare with the results of the simulation data of previous literature. The forced calculation of direct forcing immersed boundary method uses the regularized Dirac delta function which increases the drag coefficient value. When the results of this study show that the ratio of the sphere diameter to the grid size is 8 and Reynolds number (Re) is 50, the simulation is performed that the drag coefficient (Cd) is 1.805 which is larger than the literature [24]. Therefore, this study uses inward contraction of Lagrangian points. It can obtain a correction formula in form of Re, the contraction value and the grid resolution in different flow field. The simulating results of using this method indicate that the Cd can be obtained with Re = 50 is 1.582 at the same grid resolution. This value approaches the result of the literature [24]. The second part is to simulate the motion of multiple spheres in the flow field.

This study is to simulate and verify the natural settlement behavior of a single sphere firstly, and then it simulates the settling of multiple spheres at different heights. The process of settling generates collision of spheres by variation of the flow field and the spherical velocity.

Keywords: Direct numerical simulation, Immersed boundary method, the motion of spheres

(6)

iv

圖目錄

圖1.1 (a) 貼體法:結構性網格和非結構性網格結合應用 (b) 沉浸邊界法:流場網格

(藍色區塊)、結構體(橘色區塊)及其內插邊界點(黑色圓點) ... 2

圖1.2 1000 顆球體沉降瞬時位置圖[5] ... 5

圖1.3 (a) 圓柱與流場設置示意圖 (b) 在圓柱不同角度下的速度誤差[6] ... 5

圖1.4 Drafting-kissing-tumbling 現象由左到右不同時間位置圖[8] ... 5

圖2.1 三維卡式座標方向示意圖 ... 7

圖2.2 拉格朗日點各佔體積示意圖 ... 10

圖2.3 球體上拉格朗日點散佈圖(N =_L 200) (a)拉格朗日點初始位置 (b)拉格朗日點經由庫倫力平衡後的最終位置 ... 11

圖2.4 連續函數φ( )r ，橘色虛線即為拉格朗日點位置 ... 12

圖2.5 正規化狄拉克脈衝函數運作範圍示意圖，紅點、橘點和綠點為拉格朗日點，黑點表示流場網格 (a)表示單個拉格朗日點的內外插範圍，如橘點對應到橘框、綠點對應到綠框 (b)橘黃色區塊為固體於流場網格中內外插的整體運算範圍... 13

圖2.6 碰撞模型示意圖 ... 15

圖2.7 硬球碰撞模型示意圖 ... 18

圖2.8 同時多顆顆粒碰撞示意圖 ... 19

圖2.9 軟球碰撞模型示意圖[22, 23] ... 20

圖3.1 靜止球體在流場中的位置示意圖 ... 23

圖3.2 入流方向與球型顆粒x y平面角度示意圖 ... 23

圖3.3 在Re 50_D = 且不同迭代次數之拉格朗日點速度誤差點圖 ... 24

圖3.4 不同雷諾數經由不同迭代次數最大速度誤差圖 ... 24

圖3.5 不同權重對應其迭代次數 ... 26

圖3.6 Re 50_D = 不同解析度下的阻力係數值C_d，而應用Clift et al(1978)[24]經由實驗所得公式計算Re 50_D = 阻力係數實驗值為1.539，顯示模擬結果皆大於實驗值。 ... 27 圖3.7 球體厚度及拉格朗日點向內收縮值r_L示意圖，黑線為顆粒真實半徑，深藍

(9)

vii

色虛線為因厚度向外延伸的球殼，紅色虛線為經由收縮後的球體半徑，黑點為初始設定之拉格朗日點位置，綠點為向球體內收縮後的拉格朗日點位置 ... 27 圖3.8 Re 50_D = 且r h_L = 0之流線圖

(a) h = D 8 (b) h = D 12 (c) h = D 16 ... 29 圖3.9

Re 100

_D

=

且r h_L = 0^之流線圖

(a) h = D 8 (b) h = D 12 (c) h = D 16 ... 30 圖3.10

Re 50

_D

=

不同向內收縮值之流線圖(h = D 16) (a) r h_L = 0

(b) r h_L = 0.25 (c) r h_L = 0.35 (d) r h_L = 0.4 (e) r h_L = 0.5

(f) r h_L = 0.75 ... 32 圖3.11

Re 100

_D

=

不同向內收縮值之不同向內收縮值流線圖(h = D 16)

(a) r h_L = 0 (b) r h_L = 0.25 (c) r h_L = 0.35 (d) r h_L = 0.4 (e) r h_L = 0.5 (f) r h_L = 0.75 ... 34 圖 3.12 應用線性擬合在不同解析度大小的阻力係數值(藍點)，得到h ∼ D / ∞的阻力係數(紅點) (a)

Re = 50

_D (b)

Re 100

_D

=

... 36 圖3.13 雷諾數和阻力係數值的模擬結果與實驗結果[24]比對圖 ... 37 圖3.14 (a)

Re = 50

_D (b)

Re 100

_D

=

在不同解析度下及不同向內收縮值與相應的阻力係數值圖(實線)，而向內收縮值r h_L = 0^且h ∼ D / ∞所得到的阻力係數值 (紅色點)且紅色虛線為此值的延伸於不同向內收縮值以利於觀察。 ... 38 圖3.15 (a)

Re = 50

_D (b)

Re 100

_D

=

不同解析度及固定直徑不同向內收縮值下與相應阻力係數圖，不同解析度對應_{h ∼}_{D /}∞時的阻力係數值即可得到該解析度下該使用的向內收縮值。 ... 39 圖3.16 雷諾數

Re

_D 、向内收縮值r_L 和網格解析度h之擬合曲面圖 ... 40 圖4.1 球體沉降速度與時間關係圖，沒有使用拉格朗日向內收縮方法為虛線，而有使用該方法的為實線 ... 42 圖4.2 靜止球體於均勻入流流場前後擺放示意圖... 44 圖4.3 阻力係數值對應球體間擺放之距離圖[8,25,27]，

C

_do為單顆球體於均勻流場中且網格解析度極高_{h ∼}_{D /}∞的阻力係數值。 ... 45

(10)

viii

圖4.4 靜止球體於均勻入流流場並排擺放示意圖... 46

圖4.5 阻力係數值對應球體間擺放之距離圖[8,26] ... 47

圖4.6 不同時間點的兩顆球體沉降位置圖 ... 50

圖4.7 兩顆球體隨時間變化的沉降速度圖 ... 51

圖4.8 隨時間變化的兩顆球體間之距離圖 ... 51

(11)

ix

表目錄

表3.1 驗證邊界條件之模擬參數 ... 22

表3.2 模擬參數 ... 28

表3.3 模擬參數—雷諾數、入流出流速度和時間步長 ... 35

表4.1 模擬球體沉降速度之模擬參數 ... 41

表4.2 球體終端速度與文獻比較 ... 43

表4.3 模擬兩個靜止球體前後擺放於均勻入流流場之模擬參數 ... 44

表4.4 模擬兩個靜止球體並排擺放於均勻入流流場之模擬參數 ... 46

表4.5 模擬兩個球體上下擺放於流場中自然沉降之模擬參數 ... 48

(12)

1

Chapter 1 緒論

1.1 研究背景

在自然界中，固體大多存在於流體中，並兩者相互作用，如下雨天空氣中落下的水滴、沙塵暴中沙粒隨著風飄移和經由河流沖刷岸堤所產生的沙粒被流體帶往出海口而沉降等自然現象;在工程和生醫應用上，如風力發電機的發電效率、離岸水下洋流發電機、化學沉降法分離水與廢水和人體血液中血球的運動方式等問題，皆是人們想要研究和了解的問題及現象，因此對固體和液體交互作用的行為進行有效且準確的數值模擬，可對於基礎研究、了解複雜流動問題是有很大的幫助。

目前模擬固體於流場中的方式主要可以分為兩種，分別為貼體法(Body-fitted) 和沉浸邊界法(Immersed boundary method)。貼體法的流場網格主要由結構體邊界來進行繪製，流場網格可分為結構性網格和非結構性網格，並且為了模擬更好的邊界層流動現象，會在結構體周圍網格進行加密的動作(如圖 1.1(a))，並求解流場;

沉浸邊界法則不需要依照結構體來繪製網格，可直接使用求解流場的結構性網格，並將結構體放置於流場網格中(如圖 1.1(b))，而後計算其中受力，藉以模擬出結構體於流場中的運動行為及結構體對流場作用所產生的變化。若結構體為複雜幾何形狀，則應用貼體法結構性網格會難以繪製，而使用非結構性網格，倘若結構體是會隨流場進行移動，如同顆粒沉降或懸浮於水中，網格必須不斷地重新更新，會造成計算量的增加，而沉浸邊界法與貼體法相比會較容易模擬結構體的運動行為。

(13)

2

(a)

(b)

圖 1.1 (a) 貼體法:結構性網格和非結構性網格結合應用 (b) 沉浸邊界法:流場網格 (藍色區塊)、結構體(橘色區塊)及其內插邊界點(黑色圓點)

(14)

3

1.2 文獻回顧

沉浸邊界法由 Peskin (1977)為了研究人類心臟中彈性瓣膜周圍血液流動模式而提出的方法，基本的方法設計是在確定邊界(瓣膜)的位置上標記拉格朗日點，依據瓣膜隨血液流動變化，給予瓣膜邊界(拉格朗日點)一具有變形和彈性參數的外力函數，並使用正規化狄拉克脈衝函數(regularized Dirac delta function) 沉浸邊界周圍 2 個網格進行速度和受力內插，將結構體邊界的分佈外力應用於納維爾史托克方程式(Navier-Stokes equations)上進行計算，用以模擬出瓣膜邊界對於血液流動的影響 [1]，且 Peskin (2002)對沉浸邊界法的進行整理及撰寫完整的推導[2]，其上述過程將結構固體和流體分別應用尤拉(Eulerian)和拉格朗日(Lagrangian)描述法進行推導。

Roma et al. (1999)延續 Peskin(1978)的沉浸邊界法應用，且採用適應網格加密的方式，在固體沉浸邊界周圍局部加密流場網格，進而增加固體邊界條件的計算精度，並應用二維基本模型模擬問題成功驗證了此方法的可行性，另外提出一沉浸邊界周圍1.5 個網格正規化狄拉克脈衝函數進行內插速度和受力[3]。

Fadlun et al. (2000)提出直接施力(direct-forcing)沉浸邊界法，將其用以模擬三維具有複雜結構的流體上，其使用方法為在每個時間步長於邊界點上得到期望速度(desired velocity)且計算其受力，而固體邊界的速度為利用流場網格速度進行內插計算藉以描繪出固體邊界[4]。

Uhlmann (2005)應用 Fadlun(2000)直接施力沉浸邊界法模擬二維的圓柱運動和加入碰撞模型用以三維多顆球體於流場中的沉降行為(如圖 1.2)，其中應用尤拉和拉格朗日法將固體與流體分開描述，並使用Roma et al. (1999)的正規化狄拉克脈衝函數進行內插結構體拉格朗日點與流場網格的速度和受力，並且將受力回授於流場的更新流場網格速度，且應用牛頓-尤拉方程式(Newton–Euler equations)將其受力用以計算結構體的速度及角速度，藉以模擬結構體於流場中的運動行為[5]。

由於 Uhlmann (2005)的直接施力沉浸邊界方法在固體邊界上無法滿足無滑移 (no-slip)、無穿透(no-penetration)條件，Kempe & Fröhlich (2012) 發現時間步長會影響其邊界計算受力的結果，因此其應用一簡易步驟改善其方法，其為沉浸邊界受力的進行迭代計算，使結構體邊界會被描述的更準確，並對此也運用二維圓柱問題進行了驗證(如圖 1.3)，並改善 Uhlmann (2005)牛頓-尤拉方程式的數值計算方

(15)

4

法，應用體積分率的方式將結構體中的流體內能計算出來，讓顆粒質量密度比所產生數值穩定限制

ρ ρ

_p

/

_f

≥ 1.2

_{[5]的降低為}

ρ ρ

_p

/

_f

≥ 0.3

，使模擬問題的固體密度

ρ

p可小於流體密度

ρ

_f[6]，其後加入了碰撞模型，模擬兩顆球體的碰撞行為。

Breugem (2012)同樣提出 Uhlmann (2005)的數值方法所產生的邊界問題，並認為應用正規化狄拉克脈衝函數會使得球體邊界的擴大，所以用向內收縮和受力迭代的方式去計算，且應用數值模擬測試，測試出收縮值為0.3 個流場網格時，此數值方法會達到二階精度 [7]，同時也採用 Kempe & Fröhlich (2009)[13]的牛頓-尤拉方程式的計算方法，可使質量密度比降低並對其進行了驗證。

Luo et al. (2019)延續直接施力沉浸邊界法進行解析度測試，並應用 Breugem (2012)的固體邊界向內縮方法，但不同的是其應用不同解析度進行模擬測試，經由外插可得到該雷諾數下的阻力係數值，並與向內縮值構成一函數用以模擬在不同雷諾數和不同解析度下滿足阻力係數值的向內縮值，並用於模擬球體沉降問題，

成功模擬驗證單顆顆粒沉降達終端速度，且模擬兩顆顆粒上下不同位置自然沉降，模擬出Drafting-kissing-tumbling 現象(如圖 1.4)並與他人進行比對[8]。

(16)

5

圖1.2 1000 顆球體沉降瞬時位置圖[5]

圖1.3 (a) 圓柱與流場設置示意圖 (b) 在圓柱不同角度下的速度誤差[6]

圖1.4 Drafting-kissing-tumbling 現象由左到右不同時間位置圖[8]

(17)

6

1.3 研究動機

沉浸邊界法的應用越來越廣泛，除了基本球體的驗證，近年來對於複雜幾何的應用也非常多，可幫助人們了解在流體中泥沙沉降或動物的行為、也可幫助人們發展科技，如機翼或螺槳的翼型[9,10]、動物(如蜻蜓和魚等)運動於水中或空氣的造成的流場[11,12]等。在科技的應用上，像是近年來綠能產業的蓬勃發展，如同沿岸的風力發電機機和離岸發電機，其依靠大自然現象產生的電力供人們使用，

而在發電機中重要的是它的發電效率會隨著位置方向和結構外型而有所不同，為了提升其效率，則必須測試不同結構外形才能達成，若以實驗來測試則會耗費大量資金和資源，而只要得到實地流場流向及速度數據，應用模擬可不必於實地中做太多的實驗，就可得到發電效率，在其他方面的應用也是如此，模擬能節省成本並且加快產業發展的速度。

沉浸邊界法不斷被改進和發展，目前大多使用迭代方法或是在結構體邊界網格加密，藉以達成近似無滑移、無穿透邊界條件，但是應用迭代方法雖然可以讓邊界的誤差降低，但不能界定其所造成的誤差量;而網格加密方法對於移動的結構體，則需要在每個時間重新進行網格加密的動作，此動作會增加程式的運算量。

為此本研究希望發展一在固定流場網格中有效描述結構邊界，並且能精確模擬出結構體在流體中隨流場移動或是讓結構體在流體中運動的數值方法，且藉由此方法來應用顆粒沉降問題來進行模擬驗證。

1.4 論文內容概述

本研究包含本章節分為五章，在第二章節論述基本理論、數值方法及相關模型應用，主要提及求解流場及沉浸邊界法的結合應用，並加入碰撞模型模擬多顆球體及球體與牆壁的碰撞行為，而第三章節為改善沉浸邊界法並與其他文獻進行驗證，利用迭代法降低沉浸邊界法求解時所產生的邊界誤差，且應用拉格朗日點向內收縮方法使球體阻力計算更貼近前人文獻之實驗值[24]，第四章節為應用前述理論、改善方法及模型進行模擬單顆球體沉降和兩顆球體於流場中的相互作用，

而第五章為本研究數值模擬驗證及模擬結果的總結，並針對本研究在某些情況下的不足之處，提出未來工作及發展方向。

(18)

7

Chapter 2 理論和數值方法

.Equation Chapter 2 Section 2

本研究為了模擬顆粒沉降問題，必須了解描述流體與固體的耦合方式，而在模擬上要知道其在數值方法上運作模式，並假設模擬之顆粒形狀為球體，以利於與其他相關文獻模擬或實驗作驗證和比對，本章節會先提到流體統御方程式及沉浸邊界法，並以數值離散化和向量的形式來呈現，而後為了模擬多顆於流場中運動的顆粒，顆粒間會彼此碰撞，所以會加入碰撞模型以模擬此現象。

2.1 統御方程式

描述流體動力的方程式為質量守恆方程式和納維爾斯托克方程式 (Navier-Stokes equation)，在三維卡式座標系統(Cartesian coordinate system)和不可壓縮流體(incompressible flow)的假設下，並用索引表示法(index notation)的形式可表示為:

0

i i

u x

∂ =

∂ ^(2.1)

( ) 2 + f

i i

j i i

j i j j

u u u p u

t x x ν x x

∂ + ∂ = − ∂ + ∂

∂ ∂ ∂ ∂ ∂ (2.2)

其中i = 1, 2, 3分別代表 x、y、z 三個不同方向，u 為流體速度，_i p^{為壓力除以} 流體密度，

t

^為時間，ν為運動黏滯係數和f_i為單位流體密度受力。方程式(2.2)左式依序分別為加速度項和對流項，右式依序分別為壓力項、黏滯項和外力項，而外力項也是加入沉浸邊界法模擬結構物體重要的一項。

圖2.1 三維卡式座標方向示意圖 x

y

z Domain

(19)

8

為了求解流場，先將前述方程式其進行數值離散化，在動量方程式上採用 Chorin (1967)[14]提出的求解方法—投影法(Projection method)，也稱為分步法 (fractional step method)，並由 Zang et al (1994)[15]將其方法改為使用有限體積法並應用於曲線座標的計算，最後由Cui (1999)[16]使其程式能以平行化運算，由於本研究不會使用到曲線座標，以下方程式都以卡式座標呈現。

分步法將動量方程分成兩部分求解，第一部分先解對流項和黏滯項，為預測式(predictor) ，第二部分應用壓力項並計算出下一步流場流速，為修正式 (corrector)，其過程如下:

第一步: 對時間上採取半隱式求解的方式使用 Adams-Bashforth method 於對流項，

應用Crank-Nicolson method 於黏滯項，而在空間上使用 QUICK 求解對流項，應用二階中央差分(central difference)計算黏滯項，而後將速度更新為預測速度u_i^*^，其方程式如下列所示:

* 2 * 2

3 1 1 1 (u) (u)

( (u u ) ) ( (u u ) ) ( )

2 2 2

n n

i i

j i j i

j j j j j j

u u

t x x ⁻ ν x x ν x x

− = ∂ − ∂ + ∂ + ∂

∆ ∂ ∂ ∂ ∂ ∂ ∂ (2.3)

第二步: 利用壓力項和u_i^*來更新下步速度u_iⁿ⁺¹，如下式(2.4)

1 *

n

i

u u p

t x

+ − ∂

∆ = −∂ (2.4)

但由於壓力項尚未被計算，所以應用質量守恆方程式(2.1)求解壓力，所以對(2.4) 式取散度，可以得到

1 * 2

( ⁱⁿ ⁱ ) =

i i i

u u p

x t x x

+ −

∂ − ∂

∂ ∆ ∂ ∂ (2.5)

1 * 2

1 ( ⁱⁿ ⁱ ) =

i i i i

u u p

t x x x x

∂ + −∂ − ∂

∆ ∂ ∂ ∂ ∂ (2.6)

由連續方程式(2.1) ⁱⁿ ¹ = 0

i

u x

∂ +

∂ 可將(2.6)式化簡為

* 2

1 ⁱ =

i i i

u p

t x x x

∂ ∂

∆ ∂ ∂ ∂ (2.7)

方程式(2.7)為波松方程式(Poisson equation)，並應用 multigrid method 對(2.7)式的壓力進行求解，得到壓力項後代入(2.4)式即可計算出下一步的流場速度。

(20)

9

2.2 直接施力沉浸邊界法

本研究採用Uhlmann (2005)[5]的直接施力沉浸邊界法(direct-forcing immersed boundary method)處理流場中固體邊界，可應用於分步法中，並使用流場和固體間的速度差異來進行計算，從而改變固體周圍流場用以描述在流場中固體邊界，但這方法會無法完全滿足固體邊界條件(無滑移、無穿透)，其改進辦法與相關驗證會在第三章提到。

2.2.1 流固耦合方法

直接施力沉浸邊界法主要是利用拉格朗日(固體)和尤拉(流體)描述法來運作，

以下符號小寫代表尤拉(流體)描述法、大寫代表拉格朗日(固體)描述法，運作方式如下列方程式(2.8) ~ (2.11)，首先將流場預測步方程式(2.3)求解出的預測速度_u^*應用正規化狄拉克脈衝函數內插於顆粒上各個拉格朗日點上得到_U^Γ，此過程即是先判斷在拉格朗日點周圍的流場網格，並使用流場網格上的速度內插於拉格朗日點上，且與固體邊界上的速度_U^d進行受力F運算，同樣地在將計算完的受力F^應用正規化狄拉克脈衝函數外插於拉格朗日點周圍流場網格上得到f^{，隨後在用}f^以更新速度_u^**，在方程式(2.7)以_u^**取代_u^*計算求解壓力，並最後運算方程式(2.4)得出

un+1，其流固耦合運作方程式如下所示

Nx Ny Nz

Γ * 3

l i,j,k h i,j,k l

i=1 j=1 z=1

( ) =

∑∑∑

( ) (δ ) h−

U X u x x X (2.8)

d Γ

l l

l

( ) ( ) ( ) =

t

−

∆

U X U X

F X (2.9)

Np

i,j,k l i,j,k l

p=1 1

( ) = ^N^L ( ) (_h - )

l

V δ

=

∑∑

∆

f x F X x X (2.10)

** = * + ∆t

u u f (2.11)

其中拉格朗日點上的速度U 上標

Γ

和d 分別表示為流場和球體本身，x為流場網格的位置，X^{為拉格朗日點的位置，}h 為網格寬度， V∆ 為拉格朗日點所代表體積，

N

x、

N

_y、

N

_z分別為流體三個方向的網格數量，N_p為顆粒的數量，N_L^為拉格朗日點的數量，而

δ

_h為正規化狄拉克脈衝函數，其運作方式在後續小節會詳細提到。

(21)

10

2.2.2 球體均勻分佈點

應用直接施力沉浸邊界法於球體上，此方法假設球體表面有一層外殼，給予條件為一個拉格朗日點代表的體積∆ 約略等於一個流場網格的體積大小V _h³，利用這個假設和條件，可以使用球體半徑_r計算拉格朗日點數量的推導如下方程式 (2.12)~(2.16) [5]:

應用球體半徑r和假設球殼厚度為固定寬度h ，也等同於網格寬度，以定義出外圈 半徑

r

₂、內圈半徑

r

₁、(如圖 2.2)，外圈體積減去內圈體積即為外殼體積，利用顆粒所代表的體積∆ 就可以計算出顆粒的數量V

N

_L

2 2

r = r + h (2.12)

1 2

r = r − h (2.13)

3 3

2 1

4 4

3 3

N VL∆ = πr − πr (2.14) 將(2.12)和(2.13)式代入(2.14)式

2 2

(12 )

L 3 h

N r h

V

= π +

∆ (2.15)

並利用∆_V≈ _h³可得到

2 2

( 12 1)

L

3

r

N h

= π +

_(2.16)

圖2.2 拉格朗日點各佔體積示意圖

(22)

11

為了代表整個球體，必須先將拉格朗日點均勻分配於球體表面上，才能在流場中模擬出結構體邊界，在本研究採用庫倫斥力均勻分散方法[5]。先在球體表面上隨意散佈給定數量的拉格朗日點，並給予每個拉格朗日點一同號庫倫電荷q^，讓在球體表面上各個相近拉格朗日點進行同號電荷產生相互排斥力如下式(2.17)，而後受到排斥力的作用，讓拉格朗日點產生加速度，並更新拉格朗日點位置，最終達到力平衡的狀態，此結果(如圖 2.3)即為球體表面均勻分布的拉格朗日點。

3

'( ')

= '

E ke qq −

− F d d

d d (2.17)

其中

F

_E^{為電荷互斥力，}q^、q'為拉格朗日點帶有的電荷，d^、d 為拉格朗日點位' 置，

k

_e^{為庫倫常數。}

圖2.3 球體上拉格朗日點散佈圖(

N =

_L

200

) (a)拉格朗日點初始位置 (b)拉格朗日點經由庫倫力平衡後的最終位置

2.2.3 正規化狄拉克脈衝函數

在應用沉浸邊界法模擬固體於流體中，若不使用正規化狄拉克脈衝函數，而將拉格朗日點受力回授於該點所在的流體網格中，會因為受力使得該網格所更新的流體速度和周圍網格的速度差異太大，在數值計算上會有非物理現象的震盪產生，因此如2.2.1 小節所示，將受力應用正規化狄拉克脈衝函數傳遞於周圍網格中，

用以消除震盪現象[2]。在本文應用 Roma (1999)[3]發展的脈衝函數方程式，運作方式為利用拉格朗日點周圍1.5 個流場網格來內外插速度和受力，其運作寬度為3h，

也被稱作使用3 點脈衝函數，Uhlmann (2005)[5]比較使用 3 點脈衝函數[3]和使用 4

(23)

12

點脈衝函數[2]，在相同條件下利用阻力係數來證明使用 4 點脈衝函數會使固體於

流場中的平均阻力上升，故本研究使用 3 點脈衝函數，其脈衝函數方程式如

(2.18)~(2.20)式所示，此方法亦可應用於一維至三維問題上。

1 1 1

h h h h

δ ( ) = δ ( ) δ ( ) δ ( )

x X

−

^D x X

−

^D y

−

Y ^D z

−

Z (2.18)

h1

δ ( ) ^D x X 1 ( ) hφ α

− = (2.19)

正規化狄拉克脈衝函數由φ 這一連續函數所組成如下式:

( )

1 2 6 1 2 3

5 3 3( 1 ) 1 , 0.5 1.5 ( ) 1 3 1 , 0.5

0, otherwise,

α α α

φ α α α

 − − − − + ≤ ≤



=  + − + ≤





(2.20)

其中α ( - ) /= x X h。而連續函數φ α 分布圖如下圖_{( )} 2.4 所示

圖2.4 連續函數φ α ，橘色虛線即為拉格朗日點位置 ( )

(24)

13

將其運用於球體模擬上可得到運作的範圍如下圖2.5 所示(以二維的方式呈現)

圖2.5 正規化狄拉克脈衝函數運作範圍示意圖，紅點、橘點和綠點為拉格朗日點，

黑點表示流場網格 (a)表示單個拉格朗日點的內外插範圍，如橘點對應到橘框、綠點對應到綠框 (b)橘黃色區塊為固體於流場網格中內外插的整體運算範圍

2.3 牛頓—尤拉方程式(Newton-Euler equations)

描述顆粒運動行為的統御方程式為牛頓-尤拉方程式[5,31,32]，並透過相應於固體周圍流體域的動量守恆來計算顆粒於流體裡的運動，方程式分別為動量及角動量方程式如下式(2.21)和(2.22)。

d ds + ( )

d

^p

p

p f p p f

m V

t

ρ ρ ρ

= ∫

Γ

−

u



τ n



g ^(2.21)

I d ( ) ds

d

^p

p

p t

ρ

f

= ∫

Γ

×

ω



r τ n



^(2.22)

其中

Γ

_p為顆粒的表面，m_p

、、 ρ

_p V_p分別為顆粒的質量、密度和體積，而其關係為

= V

p p p

m

ρ

，I_p

、、

u ω_p _p為顆粒的慣性舉、質心速度和質心角速度，r X x

=

_l

−

_p為表面拉格朗日點位置

Xl^{和顆粒質心位置}xp^{的位置向量，而}τ為流體剪應力張量，g 為重力加速度。方程式(2.21)等號右式第一項和方程式(2.22)右式可經由柯西定理 (Cauchy’s principle)[30]進而應用動量守恆和角動量守恆描述顆粒內部受力(internal force)和顆粒內部扭矩(internal torque)，如(2.23)和(2.24)式將此兩項拆解為外力項 (external force)、外扭矩項(external torque)、動量變化率項和角動量變化率項來計算。

(25)

14

ds d d d d

f p f V f V

ρ ρ t ρ

Γ = − Ω + Ω

∫

^{τ n}^

∫

^f

∫

^u



^(2.23)

( ) ds = d d d d

f p f V f V

ρ ρ t ρ

Γ × − Ω × + Ω ×

∫

^r ^{τ n}^

∫

^r ^f

∫

^r ^u



^(2.24)

其中Ω為顆粒的體積，f^為實體力(body force)。在此 Uhlmann (2005)[5]將方程式 (2.23)和(2.24)右式的第一項—外力項、外部扭矩，可用(2.9)式算出的總受力來計算。

1 l

d ^N^L ( )

l

V V

Ω =

=

∑

∆

∫

^f ^{F X} ^(2.25)

1 l

d ( ) ^N^L

l

V V

Ω =

× =

∑

× ∆

∫

^{r f} ^{r F X} ^(2.26)

而(2.25)及(2.26)右式的第二項—動量及角動量變化率項，基於不可壓縮流體並應用沉浸邊界法，可視顆粒內部為一封閉體積，動量和角動量變化率等同於剛體運動的方程式，如(2.27)和(2.28)式。

d d = d

d d

p

V Vp

t

∫

^Ω^u ^ut ^(2.27)

d d = d

d ^p_p d ^p

V I

t

∫

^Ω^r ^× ^u ρ ^ωt ^(2.28) 最後整理(2.21)~(2.28)式可得到計算顆粒速度及角速度的方程式

1 l

d ( ) +

d ( )

NL

p f

p p f l

t V V

ρ

ρ ρ ₌

= − ∆

−

∑

u F X g (2.29)

1 l

d = ( )

d ( )

NL

p p f

p p f l

t I V

ρ ρ

ρ ρ =

− × ∆

−

∑

ω r F X (2.30)

將(2.29)和(2.30)式速度與角速度的微分式應用前向尤拉法(forward Euler method)離散，如下所示

1

1 l

( ) +

( )

n n NL

p p f

p p f l

t V V

ρ

ρ ρ

+

=

− = − ∆

∆ −

∑

u u

F X g (2.31)

1

1 l

= ( )

( )

n n NL

p p p f

p p f l

t I V

ρ ρ

+

=

− − × ∆

∆ −

∑

ω ω

r F X (2.32) 得到顆粒的質心速度u 可以進行顆粒的位置計算 _p

1

(1 1 + ) 2

n n

p p n n

p p

t

+

− +

∆ =

x x

u u (2.33)

最後運用剛體運動方程式將顆粒表面上拉格朗日點的速度進行更新

(26)

15

d( ) ul = _p + ( )ω_p × _l − _p

U X X x (2.34)

2.4 顆粒間、顆粒與牆壁的碰撞模型

在模擬多顆顆粒沉降的問題時，由於會受到邊界和流場的影響，會產生顆粒間彼此相當靠近甚至產生碰撞行為，且球體向邊壁移動與邊壁進行碰撞，而顆粒沉降於底部時會產生的堆積現象，此時除了顆粒間的碰撞同時也要考慮顆粒與牆壁間也會產生的碰撞行為，本研究應用三種模型去模擬碰撞問題，分別為潤滑模型(lubrication model)、軟球模型(soft-sphere model)和硬球模型(hard-sphere model)，

運作範圍如圖2.6 所示。

圖2.6 碰撞模型示意圖

2.4.1 潤滑模型(lubrication model)

當顆粒彼此間相當靠近時，其兩顆粒間的間距

ξ

非常的細小，因顆粒間相對速度靠近會擠壓間隙中的流體，使其向外被擠出，而在顆粒因碰撞而分開時，流體又會被推回間隙中，在此時主要的受力行為由黏滯力所主導，即為潤滑效應。但由於顆粒間的間隙相當狹小，流場網格大小不夠細小並不足以涵蓋此範圍中的黏滯力運算，所以必須應用潤滑模型來計算，而本文參考Kempe & Fröhlich (2012)[19]

的潤滑模型，其模型可以分成兩部分:

(27)

16

2.4.1.1 顆粒間的潤滑模型

顆粒間的潤滑受力(lubrication force)為

2

lub min max

6 ,

ij i j

ij ij ij

ij i j

c r r r r

πνρ ξ ξ ξ

ξ

 

= −  +  ≤ ≤

F n (2.35)

其中

ξ

_min和

ξ

_max分別為進入潤滑模型的最小和最大距離，

ξ

_min假設為自然表面粗糙度ξ_min ≈ 10 D⁻⁴ ，而 ξ_max ≈ 10 ³ ξ_min，n 為兩個顆粒間的單位向量，_ij c 為兩個顆_ij 粒的法向相對速度(純量)，ξ_ij為兩個顆粒之間的距離，各計算方程式如下

j i

ij

j i

= −

−

x x

n x x (2.36)

( )

ij i j ij

c = u − u n (2.37)

( r )

ij j i ri j

ξ = x − x − + (2.38)

2.4.1.2 顆粒與牆壁間的潤滑模型 顆粒與牆壁間的潤滑受力為

lub 2

min max

6 ^w ,

w w w

w

πνρ

c r

ξ ξ ξ

= −

ξ

≤ ≤

F n (2.39)

其中

ξ

_min和

ξ

_max為顆粒與牆並進入潤滑模型的最小和最大距離，其計算範圍同顆粒間的潤滑模型為ξ_min ≈ 10 D⁻⁴ 和 ξ_max ≈ 10 ³ ξ_min，n_w為顆粒與牆壁間的單位向量，c_w為顆粒與牆壁的相對速度，ξ_w為顆粒與牆壁之間的距離，各計算方程式如下

w i

w

w i

= −

−

x x

n x x (2.40)

w i w

c = u  n ^(2.41)

w w i ri

ξ = x − x − (2.42)

其中(2.40)和(2.42)式的x 為牆壁的位置。_w

(28)

17

2.4.2 硬球模型(hard-sphere model):顆粒與牆壁碰撞

本研究將硬球模型應用於顆粒與牆壁的碰撞行為，其方法優點為計算的快速，可應用於較少的碰撞行為，如較稀濃度的顆粒問題，而在模擬的設定上，由於邊界並不會移動，顆粒與牆壁的碰撞基本上是一對一碰撞問題，符合硬球碰撞模型的計算方式，計算方法如下所示

應用顆粒原始速度即可判定顆粒是否會撞到牆壁，而顆粒若是被判定為會進入硬球碰撞模式，則會將其過程分成兩個階段，第一階段為顆粒經由原始速度移動到牆面的過程，利用這個過程可以得到第一階段的歷程時間

= ^w ^p

p

t r

δ ^x ⁻ ^x ⁻

u (2.43)

第二階段則利用第一階段的時間可以得到反彈的時間，並且利用恢復係數 (coefficient of restitution)e^{進行反彈速度的計算}

'_p = − e _p

u u (2.44)

並利用反彈速度可以計算出顆粒反彈後的位置

( )

' p = p ∆ −t δt + w − r

x u' x (2.45)

而e^{為恢復係數}(coefficient of restitution)，其值介於 0 ~ 1之間， = 0e 為完全非彈性碰撞，顆粒碰觸到牆面即黏於牆面上，而e = 1為完全彈性碰撞，顆粒與牆面進行碰撞後，反彈速度會與原始速度相等，同樣地應用在顆粒間，e = 0即顆粒間碰撞會黏在一起，而e = 1顆粒間碰撞後，則動量會完全傳遞，若同樣質量和大小的顆粒彼此相撞，其兩者速度會進行互換。

(29)

18

圖2.7 硬球碰撞模型示意圖

2.4.3 軟球模型(soft-sphere model):顆粒間碰撞

模擬多顆顆粒沉降問題時，於在底部堆積的情況會使顆粒濃度上升，會在同一時間內多顆顆粒進行碰撞行為(如圖 2.8)，若使用硬球碰撞模型，由於此模型是顆粒一對一進行碰撞的計算，若要達到此時間內的平衡，必須使這過程不斷迭代計算，才能計算下一步，會需要大量的計算時間，而軟球模型則是可以同時計算多顆顆粒給予的總受力計算，進而更新速度和顆粒位置。本研究採用由Cundall &

Strack (1979)[20]發展的軟球碰撞模型，其使用質量—阻尼—彈簧系統(MCK system)，假設顆粒碰撞行為如受壓及回彈狀態的彈簧，並考慮在碰撞的過程中所產生的能量損失進而加入了阻尼項進行運算(如圖 2.9)，且由此系統可以計算出碰撞時的受力，其計算受力方式如下:

col_ij = − kξ_{ij ij} − η _ij, d_ij < r_i + r_j

F n c (2.46)

其中k 為彈性係數，

η

為阻尼係數，d 為兩顆顆粒的質心距離，_ij c 為兩顆粒間的法_ij 向相對速度(向量)，ξ_ij為兩顆球體重疊(overlap)距離，即壓縮量。

ij j i

d = x − x (2.47)

ij =cij ij

c n (2.48)

(30)

19

= ( ) ξ_ij r_i + r_j − x_j − x_i (2.49) 而後經由求解MCK system 的二階常微分方程

0

mij

ξ η ξ

+ +k

ξ

= _(2.50) 其中m 等於約化質量(reduced mass) _ij

=

i j

ij

i j

m m m

m + m (2.51)

可以得到彈性係數k 和阻尼係數

η

2 2

2ln

ln m kij

e e

η ^{= −} + π (2.52)

(

² ²

)

2 ln (N )

mij

k e

t π

= +

∆ (2.53)

其中N 為一固定常數，其數值可為 8[21]、10[19]和 15[23]，其存在目的於碰撞時間遠小於求解流體時間步長，為避免彈係係數過大而影響碰撞的準確性，在Brändle de Motta, et al. (2013)[21]經由數值測試，在 N < 8 時在高史托克斯數(Stokes number)無法很好地呈現碰撞行為，而且較高的 N 並不會有更好的碰撞結果產生只會增加運算時間，而本文研究應用van der Hoef, et al (2006)[22]測試出來的範圍，

設定N = 15去進行計算。

圖2.8 同時多顆顆粒碰撞示意圖

(31)

20

圖2.9 軟球碰撞模型示意圖[22, 23]

(32)

21

Chapter 3 沉浸邊界法驗證和改善

Equation Chapter (Next) Section 1

為了模擬顆粒沉降問題，本章節使用第二章節所提到的直接施力沉浸邊界法[5]

應用於三維流場中單顆靜止的球型顆粒，參考 Uhlmann (2005)[5]和 Kempe &

Fröhlich (2012)[6]的數值方法，並應用固體邊界上的流場速度以驗證其方法準確性，而後驗證顆粒所受阻力，且使用 Breugem (2012)[7]的拉格朗日點向內收縮方法(inward retraction of Lagrangian points)，並且利用 Luo et al. (2019)[8]的建立修正函數的方式加以改進沉浸邊界法。

3.1 固體邊界速度驗證及改善:迭代法

在流場中滿足固體的邊界條件即為(2.10)式中U X^Γ( ) = _l U X ，流場網格速度^d( )_l 內插在固體邊界上的速度要等同於固體拉格朗日點上速度，對於靜止固體

Γ( ) = 0l

U X ，但是直接施力沉浸邊界法並無法一次準確計算出此條件，而後Kempe

& Fröhlich (2012)[6]進行迭代的方法來使固體邊界的速度更趨近於無滑移及無穿透條件，應用第二章的方程式(2.8)至(2.11)可得到基本的直接施力沉浸邊界法運算方式[5]，而在這運算方式加入迭代的過程，本文使用靜止球體在均勻入流流場中進行模擬，並運用迭代方法對球體邊界誤差進行驗證。

其運作方法為重複計算在方程式(2.8)至(2.11)，重複計算次數為n_f，利用重複地計算流場網格內插於拉格朗日點上的速度U X 與拉格朗日點本身的速度^Γ( ) _l

d( )l

U X 差值並除以時間步長，從而再次計算球體受力，且讓受力外插於流場網格 中進行速度更新，此舉能讓流場網格內插於拉格朗日點的速度更趨近於拉格朗日點本身的速度，即為U X^Γ( ) = _l U X ，而達穩態後利用(2.9)式計算出來的^d( )_l x^方向受力進行加總即為球體所受阻力。

(33)

22

3.1.1 模擬設置

模擬實驗的流場大小參考Luo et al. (2019)[8]的模擬設置，如表 3.1 和圖 3.1 所示，在流場邊界條件上，_x方向的邊界條件為給定入流速度和出流速度，並使其相等以確保在計算上的進出流量相等，其餘的邊界條件皆設為無應力(stress free)邊界條件。(Re_D ^{U D}^b

= ν )

表3.1 驗證邊界條件之模擬參數

流場大小 ₍_{X Y Z}× × ₎ 1.92 10 m 0. 96 10 m 0. 96 10 m× ⁻² × × ⁻² × × ⁻² 網格解析度 ₍_N_x_×_N_y_×_N_z₎ 192 96 96× ×

入流速度U_b( inflow )

1.25 10 m s (× ⁻² ⁻¹ Re_D =10 ) 6.25 10 m s (× ⁻² ⁻¹ Re_D =50 ) 1.25 10 m s (R× ⁻¹ ⁻¹ e_D =1 0 )0

出流速度 ( outflow )

1.25 10 m s (Re× ⁻² ⁻¹ _D =1 0 ) 6.25 10 m s (Re× ⁻² ⁻¹ _D=5 0 ) 1.25 10 m s (R× ⁻¹ ⁻¹ e_D =1 0 )0

流體密度_ρ_f 1000 kg m^-3

流體運動黏滯係數

ν

_{10 m s}⁻⁶ ² ^-1

顆粒直徑D 0.8 10 m× ⁻³

顆粒質心位置 _{(x, y,z)}

(

0.64 10 m, 0.48 10 m, 0.× ⁻² × ⁻² 48 10 m× ⁻²

)

時間步長∆t

4 D

5 10 s (Re× ⁻ =10 )

4 D

1 10 s (Re× ⁻ =50 )

5 D

5 10 s (Re× ⁻ =100 )

拉格朗日點數目N_l 202

迭代次數 n_f 1 , 3 , 5 , 7 , 10 , 12 , 15 , 20

(34)

23

圖3.1 靜止球體在流場中的位置示意圖

3.1.2 模擬結果

經由模擬測試，應用將流場內插於拉格朗日點的速度誤差投影在經過球體中心x y平面，並使用角度θ 為x軸，即為顆粒指向出流方向的角度為0 度，速度誤差為 y^{軸進行繪製}(如圖 3.3)，且將同個迭代值得的所有誤差點運用多項式曲線擬合，比較紅、黑和藍線，迭代次數的增加讓在顆粒邊界上的速度誤差降低。而後利用最大速度誤差與迭代次數來繪製圖3.4，同樣可以得到迭代次數高而顆粒邊界速度誤差降低的結論，但也隨著迭代次數的增加，誤差與迭代次數的斜率逐漸趨緩，不過在顆粒邊界上的速度誤差仍然會降低。

圖3.2 入流方向與球型顆粒x y平面角度示意圖

(35)

24

圖3.3 在Re 50_D = 且不同迭代次數之拉格朗日點速度誤差點圖

圖3.4 不同雷諾數經由不同迭代次數最大速度誤差圖

3.1.3 改善方法

經由模擬測試，可得到迭代次數多誤差越低的結果，而 Kempe & Fröhlich (2012)[6]基本使用的迭代次數為 3，但由圖 3.4 顯示在不同雷諾數且相同迭代數的情況所造成的誤差不同，並且隨著迭代數的上升能得到更趨近於顆粒表面邊界條件的結果，因此本研究利用設定誤差的方式，讓此計算受力的過程不斷的重複運算，直到低於設定誤差值才計算下一步，並在更新速度的方程式裡，以增加一權重(weight)的方式加速迭代過程，計算過程的虛擬碼(pseudo code)如下所示:

(36)

25

1 while U X^d( ) ( ) Error_l − U X^Γ _l >

判定拉格朗日點本身的速度與流場網格內插在拉格朗日點上的速度差值，若大於給定誤差時則需重複執行下列運算。

2 ^if U X^d( ) ( ) Error^l − U X^Γ ^l >

^u^*⁼^u^**

若上述條件滿足則更新的速度會在計算一次。

3 ^Γ l ^Nx ^Ny ^Nz ^* i,j,k h i,j,k l ³ i=1 j=1 z=1

( ) =

∑∑∑

( ) (δ ) h−

U X u x x X

將流場網格速度內插在拉格朗日點上

Γ( ) l

U X 。 4

d Γ

l l

l

( ) ( ) ( ) =

t

−

∆

U X U X

F X ^{計算拉格朗日點上的}_受力。

5 _i,j,k ^N^p _l _i,j,k _l

p=1 1

( ) = ^N^L ( ) (_h - )

l

δ V

=

∑∑

∆

f x F X x X 將受力外插於流場網

格上。

6 u_temp = u^* + f ∆t

u^** = weight (u^* + × u_temp - )u^*

利用權重(weight) 加速收斂方法，進行流場網格上的速度更新。

其中權重(weight)經由模擬實驗測試，給定一均勻流場，並放置一靜止球體於流場中，測量其迭代次數。模擬結果如圖 3.5 所示，x 軸為權重，y 軸為迭代次數，圖中藍色點為不同權重值時運算達穩態時的迭代次數，顯示當權重值為2.5 時會使模擬運算達到穩態時有最低迭代次數8，所以使用此權重能使每步運算有最低迭代次數，能讓計算速度提升使整體使用時間下降。

(37)

26

21

14 12 11 10

9 8 10

0 5 10 15 20 25

1 1.5 2 2.5 3

iteration numbernf

weight

weight_iteration test

圖3.5 不同權重對應其迭代次數

3.2 邊界範圍

應用3.1.3 節改善的方式，可將顆粒邊界速度誤差降至設定的範圍內，然而在進行阻力計算時，如圖 3.6 所示，x 軸為解析度，y 軸為阻力係數值，圖上各藍色的點表示不同解析度下所模擬出球體的阻力係數值，可以發現測量出的阻力係數值比實驗值1.539[24]高，此問題在 Breugem (2012)[7]和 Luo et al. (2019)[8]的期刊論文裡被討論及以其他的方式進行改善，原因由於使用第二章所論述的正規化狄拉克脈衝函數，會使得在球型顆粒表面上增加了一層多孔性的厚度進而增加了顆粒實際大小，造成顆粒於流場中所受阻力上升，其改善方法為使用拉格朗日點向顆粒內部收縮(inward retraction of Lagranian points, IRL)[7]的方式(如圖 3.7)，

Breugem (2012)[7]應用此方法達成數值方法上的二階精度，而 Luo et al. (2019)[8]

參考 Breugem (2012)的方法利用不同解析度計算的阻力係數值外插得到的解析度極高時的阻力係數值，且利用拉格朗日向內收縮法去對應該解析度收縮值應為多少可達到正確的阻力係數值，而最後可以得到一收縮值、解析度和雷諾數的函數來進行模擬應用，本文對此方法進行了驗證且將其實際應用於模擬上。

(38)

27

圖3.6 Re 50_D = 不同解析度下的阻力係數值C_d^，而應用Clift et al. (1978)[24]經由實驗所得公式計算Re 50_D = 阻力係數實驗值為 1.539，顯示模擬結果皆大於實驗值。

圖 3.7 球體厚度及拉格朗日點向內收縮值r_L示意圖，黑線為顆粒真實半徑，深藍色虛線為因厚度向外延伸的球殼，紅色虛線為經由收縮後的球體半徑，黑點為初始設定之拉格朗日點位置，綠點為向球體內收縮後的拉格朗日點位置

利用沉浸邊界法於三維多球體沉降的數值模擬

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

Institute of Applied Mechanics College of Engineering

National Taiwan University Master Thesis

利用沉浸邊界法於三維多球體沉降的數值模擬

Development of an immersed boundary method to the numerical simulation of settling of multiple spheres

簡鈺霖 Yu-Lin Chien

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

中華民國 108 年 7 月

July, 2019

致謝

中文摘要

ABSTRACT

目錄

圖目錄

Re 100

=

Re 50

=

Re 100

=

Re = 50

Re 100

=

Re = 50

Re 100

=

Re = 50

Re 100

=

Re

C

表目錄

Chapter 1 緒論

1.1 研究背景

1.2 文獻回顧

ρ ρ

/

≥ 1.2

ρ ρ

/

≥ 0.3

ρ

ρ

1.3 研究動機

1.4 論文內容概述

Chapter 2 理論和數值方法

2.1 統御方程式

t

2.2 直接施力沉浸邊界法

∑∑∑

( ) ( ) ( ) =

−

∆

∑∑

Γ

N

N

N

δ

r

r

N

( 12 1)

3

= π +

F

k

N =

200

δ ( ) = δ ( ) δ ( ) δ ( )

−

−

−

−

( )

( )

2.3 牛頓—尤拉方程式(Newton-Euler equations)

d ds + ( )

d

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

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

、、 ρ

、、