光滑粒子流體動力學在微流道中紅血球變形之模擬與應用
87
0
0
全文
(2)
(3) 中文摘要 本文發展了光滑粒子流體動力學法(SPH)應用在二維紅血球於微流道 中的變形，紅血球的細胞膜以一群粒子離散並與相鄰的細胞膜粒子以線 性彈簧連接，而血漿及細胞質由一群不可壓縮牛頓流體粒子所表示。為 驗證 SPH 算法，模擬紅血球於壓力驅動流中的變形，紅血球往下游移動 並由原先的雙凹透鏡形變成特徵形狀降落傘形(parachute shape)，與文獻 結果一致，為模擬紅血球變形提供了另一種可行的途徑。較早的研究， 大多假設紅血球的細胞內液和血漿的性質是相同的，因此本文模擬不同 密度比、黏滯係數比對紅血球在微流道中的變形的影響，結果顯示改變 彈簧常數、密度比、黏滯係數比皆影響其變形和速度，其中以改變細胞 膜的彈簧常數影響最大。最後模擬了紅血球進入微流道的變形過程，兩 翼向中心線靠攏並向後延伸，中間增厚。通過流道的時間隨彈簧常數的 增大而增加，但此表現在常數變小較不明顯。. 關鍵字 : 光滑粒子流體動力學法、紅血球. I.
(4) Abstract This thesis presents a smooth particle hydrodynamics (SPH) model for simulation of two dimensional deformation of a red blood cell (RBC) in micro-channel. In the model, the RBC membrane is discretized by membrane particles that are connected to neighboring membrane particles by linear spring. The plasma and cytoplasm are modeled as incompressible homogeneous Newtonian fluids, discretized by SPH particles. In order to verify the model, the deformation of a RBC in capillary, with the RBC moving downstream and processing a characteristic parachute shape in the steady state, results is well-agreement of previously reported. Then the effects of parameter were also investigated, which include the spring modulus, density and viscosity ratios of the fluid to RBC. The simulation results reveals that all of the parameters have effects on the dynamic behavior of the RBC. Finally, the dynamic behavior of an RBC entering a contraction micro-channel was presented. As the RBC moves, two wings fold toward each other, and bending of the membrane is increased. The result agrees well with the experiment results of reported previously in the literature, and the traveling time increases as the spring modulus is increased. Key words: smooth particle hydrodynamics、red blood cell. II.
(5) 誌謝 首先誠摯的感謝指導教授陳朝光教授及賴新一教授，老師悉心的教導， 不時的指點我正確方向，使我在這兩年中獲益匪淺。老師對學問的嚴謹 更是我學習的典範。感謝口委周煥銘教授給予我許多寶貴意見。 兩年的日子裡，實驗室共同的生活點滴，學術上的討論、趕作業的革 命情感……，感謝宜謙、修逸、韋豪學長們給予的照顧，同學詩璇的幫 忙，恭喜我們順利走過這兩年。以及實驗室可愛的學弟妹泉泰、睿群、 祥宇、季培、姿菱你們的幫忙我銘記於心。 最後，謹以此文獻給我摯愛的家人，在求學過程中默默支持、體諒、 包容更是我前進的動力。. III.
(6) 目錄. 中文摘要 ............................................................................................................I Abstract............................................................................................................. II 誌謝 ................................................................................................................. III 目錄 ................................................................................................................. IV 表目錄 .......................................................................................................... VIII 圖目錄 .......................................................................................................... VIII 符號目錄 ......................................................................................................... XI 第1章. 緒論 .................................................................................................. 1. 1.1. 研究動機 .......................................................................................... 1. 1.2. 研究目的 .......................................................................................... 2. 1.3. 章節瀏覽 .......................................................................................... 2. 第2章. 文獻回顧與基本假設 ...................................................................... 4. 2.1. 文獻回顧 .......................................................................................... 4 2.1.1 以粒子法模擬血液之文獻回顧 ........................................... 4 IV.
(7) 2.1.2 光滑粒子流體動力學法之文獻回顧 ................................... 5 2.2. 血液模型的基本假設與控制方程.................................................. 7 2.2.1 基本假設 ............................................................................... 7 2.2.2 流體控制方程 ....................................................................... 7 2.2.3 紅血球的形狀方程 ............................................................... 8 2.2.4 紅血球的粒子模型 ............................................................... 8. 第3章. 模型的建構與模擬流程 ................................................................ 10. 3.1. 光滑粒子流體動力學法的基本原理............................................ 10. 3.2. 光滑粒子流體動力學法的基本方程............................................ 11 3.2.1 函數的積分表示法 ............................................................. 11 3.2.2 粒子近似法 ......................................................................... 13. 3.3. 支持域與影響域 ............................................................................ 14. 3.4. 光滑函數 ........................................................................................ 16 3.4.1 光滑函數的求解步驟 ......................................................... 17 3.4.2 常用的光滑函數 ................................................................. 21. 3.5. 控制方程式之離散 ........................................................................ 24 V.
(8) 3.5.1 密度計算 ............................................................................. 24 3.5.2 動量方程 ............................................................................. 26 3.5.3 狀態方程 ............................................................................. 27 3.6. 粒子的位移修正法 ........................................................................ 27. 3.7. 邊界修正法 .................................................................................... 28 3.7.1 邊界力法 ............................................................................. 28 3.7.2 耦合邊界法 ......................................................................... 29 3.7.3 鏡像邊界法 ......................................................................... 29. 3.8. 相關程式原理與步驟 .................................................................... 29 3.8.1 粒子間相互作用的對稱化 ................................................. 29 3.8.2 相鄰粒子搜索法 ................................................................. 30 3.8.3 週期性邊界條件 ................................................................. 33 3.8.4 最小映像法則 ..................................................................... 34 3.8.5 紅血球變形指數(DI 值) ..................................................... 35. 3.9. 程式流程 ........................................................................................ 35. 第4章. 數值模擬與結果討論 .................................................................... 48 VI.
(9) 4.1. 基本流體算例驗證 ........................................................................ 48 4.1.1 Couette 流動 ........................................................................ 48 4.1.2 Poiseuille 流動..................................................................... 49. 4.2. 紅血球於微流道中變形模擬與結果討論.................................... 50 4.2.1 在壓力驅動流中二維紅血球的變形模擬 ......................... 50 4.2.2 紅血球進入微流道中二維的變形模擬 ............................. 53. 第5章. 總結與展望 .................................................................................... 68. 5.1. 總結 ................................................................................................ 68. 5.2. 展望 ................................................................................................ 69. 參考文獻 ......................................................................................................... 70. VII.
(10) 表目錄 表 4-1 COUETTE 流資料.................................................................................. 55 表 4-2 POISEUILLE 流資料............................................................................... 57. 圖目錄 圖 2-1 紅血球原始比例 ................................................................................... 9 圖 2-2 紅血球粒子模型示意圖 ....................................................................... 9 圖 3-1 支持域位於問題域內 ......................................................................... 37 圖 3-2 支持域與問題域交錯 ......................................................................... 37 圖 3-3 粒子法求解示意圖 ............................................................................. 38 圖 3-4 支持域示意圖 ..................................................................................... 38 圖 3-5 影響域示意圖 ..................................................................................... 39 圖 3-6 影響域不同示意圖 ............................................................................. 39 圖 3-7 高斯函數及其一階導數 ..................................................................... 40 圖 3-8 三次樣條函數及其一階導數 ............................................................. 40 圖 3-9 邊界力法示意圖 ................................................................................. 41 圖 3-10 耦合邊界法示意圖 ........................................................................... 41 圖 3-11 鏡像邊界法示意圖 ........................................................................... 42 VIII.
(11) 圖 3-12 在二維空間中應用全配對搜索法 ................................................... 42 圖 3-13 在二維空間中應用鏈表搜索法 ....................................................... 43 圖 3-14 週期性邊界示意圖 利用週期重複以組成真實系統 ..................... 43 圖 3-15 一維週期性邊界示意圖( xi > L ) ....................................................... 44 圖 3-16 一維週期性邊界示意圖( xi < L ) ....................................................... 44 圖 3-17 一維最小映像法則示意圖( xij <. L ) .................................................. 44 2 L 2. 圖 3-18 一維最小映像法則示意圖( xij < − ) ............................................... 45 L 2. 圖 3-19 一維最小映像法則示意圖( xij > ) .................................................. 45 圖 3-20 程式流程圖 ....................................................................................... 46 圖 3-21 程式結構圖 ....................................................................................... 47 圖 4-1COUETTE 流粒子配置圖 ...................................................................... 55 圖 4-2 COUETTE 流在 1 秒時粒子的位置速度分布圖 ................................. 56 圖 4-3 COUETTE 流動的速度分布圖 ............................................................. 56 圖 4-4 POISEUILLE 流粒子配置圖................................................................... 57 圖 4-5 POISEUILLE 流動在 1 秒時粒子的位置速度分布圖 .......................... 58 圖 4-6 POISEUILLE 流動的速度分布圖 .......................................................... 58 圖 4-7 紅血球在壓力驅動流的 SPH 粒子配置圖 ....................................... 59 圖 4-8 紅血球在壓力驅動流中隨時間變化的形變圖 ................................. 59. IX.
(12) 圖 4-9 紅血球達穩定狀態的比較圖 ............................................................. 60 圖 4-10 不同彈簧常數紅血球變形比較圖 ................................................... 60 圖 4-11 不同彈簧常數速度比較圖 ............................................................... 61 圖 4-12 不同密度比紅血球變形比較圖 ....................................................... 61 圖 4-13 不同密度比速度比較圖 ................................................................... 62 圖 4-14 不同黏滯係數比紅血球變形比較圖 ............................................... 62 圖 4-15 不同黏滯係數比速度比較圖 ........................................................... 63 圖 4-16 紅血球進入微流道的粒子配置圖 ................................................... 63 圖 4-17 紅血球進入微流道( 5µ m )的形變圖................................................. 65 圖 4-18 不同彈簧常數通過微流道示意圖 ................................................... 66 圖 4-19 不同彈簧常數通過的時間圖 ........................................................... 66 圖 4-20 通過的時間圖 ................................................................................... 67. X.
(13) 符號目錄 羅馬字母 c. 人工音速. D. 邊界力模型中邊界力公式之係數. dim. 問題維度. F. 單位質量所受的力. g. 任意函數. h. 光滑長度. m. 質點質量. n. 外法線. N. 整個計算區域內的總粒子數. p. 壓力. r. 質點座標向量 r= xiˆˆ+ yj. S. 面域. 希臘字母. δ. Dirac 函數. ∆t. 真實的物理時間步長. ∆u. XSPH 修正項. ∆x. 質點的平均距離 XI.
(14) ε. 剪切應變. ε. XSPH 修正項中的可調參數. κ. 光滑函數的緊支域比例係數. µ. 動黏滯係數. ρ. 質點的密度. σ. 應力張量. τ. 黏滯應力. Ω. 體積域. 下標 0. 初始值. i, j. 質點序號. XII.
(15) 第1章. 緒論. 1.1 研究動機 眾所皆知，在血液循環中，是由心臟加壓，經動脈到微血管，最 後再從靜脈回到心臟，藉由此過程達到交換養分、帶走廢棄物的功能。 然，血液在血管內有一定的阻力，而紅血球的血容比( hematocrit )深 深的影響血液黏度的高低，當血容比越大，血液就越黏越不易流動。 正常人的血容比大約在 40%~45%，在如此高濃度的情況下，欲維持 暢通，主因為紅血球具強大的變形能力。 健康的紅細胞藉由變形擠過毛細血管輸送氧氣至身體的各個部 分，而研究發現許多疾病中，紅血球的變形能力異常、細胞的機械性 質改變，往往加速了病理症狀的惡化。舉例來說，常見的糖尿病、心 肌梗塞、末梢血管疾病，都被證實紅血球的變形能力異常，使得微循 環時間增長，導致了疾病的惡化。當紅血球被瘧原蟲感染時，會使得 楊氏模數增為原來的十倍左右[24]，使它不能通過毛細血管、進而阻 塞血管，導致患者昏迷甚至死亡。又如，癌細胞的楊氏模數為健康細 胞的十分之一左右[12]，因此癌細胞的變形能力比正常細胞優越，使 其更容易轉移至其他組織、器官。由以上例子，研究血細胞的變形，. 1.
(16) 可了解病情惡化的背後機制。更進一步，將有可能設計出新的診斷設 備和技術來檢測疾病的早期階段。. 1.2 研究目的 本研究的主要目的為利用光滑粒子流體動力學法(SPH)來模擬二 維紅血球的變形，並驗證此法的可靠性和可行性。 本文主要工作為以下幾點: 1. 對研究區域採用粒子方法來描述、並推導 SPH 方法粒子表達的公 式，熟悉相關的離散過程。 2. 瞭解 Navier-Stokes 方程的推導過程,採用 SPH 方法對 NavierStokes 方程的每一項進行離散, 得到 SPH 控制方程組，按照 SPH 的計算步驟，編寫程式。 3. 通過計算紅血球於壓力驅動流中的變形與文獻比對，驗證該模型， 並探討細胞膜彈性、密度比、黏滯係數比對其產生的影響。 4. 最後模擬血細胞進入微流道的變形。. 1.3 章節瀏覽 本論文共分為五章，各章安排如下: 第一章為緒論，闡述研究動機與目的。 2.
(17) 第二章為文獻回顧與基本假設，首先介紹粒子法用於模擬血流的 研究概況，接著是光滑粒子流體動力學法的文獻回顧，最後說明本文 的基本假設。 第三章為光滑粒子流體動力學法的基本理論，首先介紹光滑粒子 流體動力學法的基本方程、光滑函數和邊界修正法，接著對控制方程 進行離散，最後說明求解步驟及程式流程，並對所編寫的程式進行驗 證。 第四章為數值模擬與討論，本文模擬紅血球於壓力驅動流中的變 形與文獻比對驗證模型正確性並探討各參數的影響。最後模擬紅血球 進入微流道的變形，探討細胞膜的彈性與通過流道時間的影響。 第五章為總結，整合前述之討論與分析結果，說明本文之研究成 果並對未來研究方向提出改進與建議。. 3.
(18) 第2章. 文獻回顧與基本假設. 2.1 文獻回顧 2.1.1 以粒子法模擬血液之文獻回顧 血液在宏觀上可被當作均質流體，因此常使用 FEM、FDM、FVM 作為分析方法。然而，微循環系統是指血管直徑小於 180μm，其中 微動脈和微靜脈 10μm~180μm、毛細血管 3μm~10μm，故不能被 視為均質流體需視為懸浮溶液，且此時的脈動性已消失。所以，使用 粒子法來模擬是一個自然的選擇。同時，粒子法用於研究結構與疾病 的耦合與傳統網格法相比，有著獨特的潛力。粒子法在內部結構的建 模提供了更多的自由，如不同的粒子可代表不同的胞器、具不同性質， 並且在處理大變形及複雜流場具有較高的靈活性，可使模型更為精確， 故本文使用粒子法來模擬血液。 模擬紅血球在毛細血管的力學行為有多種數值方法，近期最常用 的是半隱式法(MPS)，此法亦為拉格朗日法。Tsubota 等人於 2005 年 使用此法研究血栓的形成和破壞導致血小板的聚集及紅血球的軸向 遷移運動[26]，隔年 Tsubota、Wada、Yamaguchi 使用半隱式法分別分 析了紅血球於毛細血管的運動[28]和血球容積比(Htc)對血壓、流速的. 4.
(19) 影響[29]。Pivkin 和 Karniadakis 更於 2008 年使用耗散分子動力學(DPD) 模擬三維紅血球的運動情形[1]。但離散粒子模型用於微觀血液的模 擬尚處於起步階段，例如:細胞膜用簡單的彈簧系統模擬，未考慮到 胞器、Tsubota 使用的模型未考慮到細胞質及約束條件會限制膨脹收 縮的總面積、顆粒級模型的參數未有生理測量的數據，這些問題都會 造成模擬結果與現實有落差。 近幾年來，許多研究團體使用半隱式法(MPS)來模擬紅血球在毛 細血管中的力學行為。而 SPH 與 MPS 相比，SPH 方法優點在易於給 出更為清晰光滑的粒子分佈形狀，且具有較高的計算效率，適合求解 大型複雜流動問題[32]，因此本文使用 SPH 方法方法來模擬二維紅血 球的變形。 2.1.2 光滑粒子流體動力學法之文獻回顧 本論文的中心是 SPH 方法，SPH 方法是一種無網格的拉格朗日質 點方法，具體概念簡單、容易編寫程式且易於添加新的物理模型。它 是由 Lucy、Gingold 和 Monaghan 在 1977 年分別獨立提出的。SPH 法最初用於解決三維開放空間中的天體物理問題，隨後在天體物理的 諸多領域中成功應用，如 Benz、Monaghan 等應用 SPH 方法模擬雙 子星及恆星的碰撞問題，Monaghan 甚至用 SPH 法研究了宇宙形成的 5.
(20) 過程。而 Monaghan 於 1992 年在 Annual Reviews of Astronomical and Astrophysics 發表的 SPH 綜述被認為是此法走向成熟的指標[14]。 現今，SPH 方法已被應用於許多不同的領域，Swegle 用 SPH 方 法 研 究彈 性流 動問 題 ， Morris 用 此 法 研究 磁流 體動 力 學問 題， Monaghan 和 Morris 等人將 SPH 方法應用於不可壓縮流問題的研究， ； Benz 和 Asphaug 則將 SPH 應用於易碎固體的斷裂模擬，Libersky、 Randles 及 Johnson 等人研究固體高速碰撞問題上做出了突破性的貢 獻，而 Lui 等人應用 SPH 方法模擬了一系列的爆炸現象[11]。 SPH 方法也有一些不足之處。早期的 SPH 方法存在穩定性和精度 方面的缺陷以及沒有通用的邊界處理方法等問題。因此諸多學者 Monaghan、Morris、Balsara、Ben、Swegle 等人都對 SPH 方法進行理 論探討和數值分析。Monaghan 等人於 1983 提出了適用於 SPH 方法 的人工黏滯項、1989 給出了粒子運動的平均速度修正項。Balsara、 Swegle 和 Hicks 等分別對 SPH 進行了穩定性分析。通過研究，逐漸 了解了 SPH 方法的穩定性、精度及收斂性質等。然而，大多數研究 都假設質點均勻分布，且大多數還是在一為情況下分析，因此，所得 的結果大多都限制於理想狀態。但在大多數情況下，特別是有大變形 或是有沖擊負載作用的情況，粒子分布都是不規則的，由於目前還很 6.
(21) 難釐清粒子的不規則分布是如何影響模擬結果的精度，因此，這些分 析所得到的結果並不十分可靠[11]。. 2.2 血液模型的基本假設與控制方程 2.2.1 基本假設 1. 本文模擬空間為二維空間。 2. 設血液只由血漿及紅血球組成。 3. 血液流動假設為穩態、不可壓縮牛頓流體，遵守 Navier-Stoke 方 程式。 4. 紅血球可簡化為細胞膜和細胞內液的組合。 5. 紅血球視為彈性薄膜包圍細胞內液所製成的膠囊。 2.2.2 流體控制方程 1. 連續性方程: dρ = − ρ∇ ⋅ v dt. (2.1). dv 1 = ∇ ⋅σ dt ρ. (2.2). 2. 動量方程:. 忽略重力作用， σ 表總應力張量，由壓力 p 和剪應力 τ 組成 σ αβ = − pδ αβ + τ αβ 7. (2.3).
(22) 式 ρ 表密度， v 表速度向量， p 表壓力， t 表時間。 2.2.3 紅血球的形狀方程 正常紅血球懸浮於等滲透壓溶液時，是一個直徑 8.1µ m 、體積約 90 µ m 2 、表面積約 140 µ m 2 ，呈扁平雙凹透鏡形，而與其相同體積的球. 形細胞表面積卻只有 97 µ m 2 ，可以見得，如此大的剩餘表面積提供了 紅血球優越的變形能力。 本文使用 1972 年 Evans 和 Fung 教授所提出的狀態方程[4]來描述 二維紅血球靜止時的形狀，如下所示: x = ac cos φ ac 2 4 y= 2 sin φ (c1 + c2 cos φ + c3 cos φ ). (2.4). 其中 c = 1.3858189 ; c1 = 0.207 ; c2 = 2.003 ; c3 = −1.123 ; a ≈ 2.82µ m 為紅血球的 半徑。圖 2-1 2.2.4 紅血球的粒子模型 見圖 2-2，使用彈簧系統來模擬彈性細胞膜。彈性薄膜為固體粒 子，相鄰粒子間以線性彈簧串聯[27]。 = F k (l − l0 ). 8. (2.5).
(23) 圖 2-1 紅血球原始比例. 圖 2-2 紅血球粒子模型示意圖. 9.
(24) 第3章. 模型的建構與模擬流程. 3.1 光滑粒子流體動力學法的基本原理 流體動力學問題主要是基於解速度、密度、能量等變量場的偏微 分方程組。除了一些簡單問題可得偏微分方程組的解析解外，大部分 的流體動力學問題只能求其數值解。 光滑粒子動力學的原理如下: 一、. 使用 SPH 方法求解流體動力學問題時，問題域由一群任意分. 佈的粒子來表示，粒子間無任何連接，因此 SPH 方法是無網格性 質的。 二、. 場函數用積分法來近似。(SPH 中稱核近似法). 三、. 接著使用支持域內的相鄰粒子對應值疊加求和來取代場函數. 的積分表示式。(SPH 中稱粒子近似法) 四、. 將所有偏微分方程組以粒子近似法進行離散。. 五、. 最後應用顯示積分法得到所有粒子的場變量。. 由以上分析可知，SPH 方法是一種純拉格朗日且具無網格、自適 應性質的流體動力學解法。 10.
(25) 3.2 光滑粒子流體動力學法的基本方程 SPH 的基本方程常按兩個關鍵步驟進行。第一步為積分表示法， 又稱場函數核近似法，第二步為粒子近似法。[11] 3.2.1 函數的積分表示法 SPH 方法是採用核近似方法將微分形式的方程式轉為積分形式。 SPH 最初想法來自於 Dirac 函數，函數 g ( x) 的積分表達式為. ∫ g ( x ')δ ( x − x ')dx '. = g ( x). (3.1). Ω. 其中 x ∈ Ω ， δ ( x − x ') 為 Dirac 函數。使用光滑函數(亦稱核函數)代替上 述的 Dirac 函數，可得到下式做為 g ( x) 的近似: g ( x) ≈ ∫ g ( x ')W ( x − x ', h)dx '. (3.2). Ω. 其中 W ( x − x ', h) 稱為光滑函數， h 為光滑長度，它是光滑核函數影響區 域的半徑，直接影響了方法的精確度。在 δ ( x − x ') 被 W ( x − x ', h) 取代後， 式(3.2)中的積分表達式只是近似相等，以< >表達近似式。 = < g ( x) >. ∫ g ( x ')W ( x − x ', h)dx '. (3.3). Ω. 由 Taylor 級數分析可得，當粒子均勻分布在點 x 周圍時，式(3.3)具有 二階精度。說明如下，給定光滑函數的支持域為 || x − x ' ||≤ kh ， g ( x) 是 可微的，將式(3.3)在 x 處泰勒級數展開可得. 11.
(26) <= g ( x) >. ∫ [ g ( x) + g '( x)( x − x ') + r (( x '− x) )]W ( x − x ', h)dx ' 2. Ω. = g ( x) ∫ W ( x − x ', h)dx ' + g '( x) ∫ ( x − x ')W ( x − x ', h)dx ' + r (h 2 ) Ω. (3.4). Ω. 其中 r 為殘餘值。因為 W 為與 x 相關的偶函數， ( x '− x)W ( x − x ', h) 應為 奇函數，故可得. ∫ ( x − x ')W ( x − x ', h)dx '=0. (3.5). Ω. 利用式(3.3)和式(3.5)，式(3.4)可改寫為 < g (= x) > g ( x) + r (h 2 ). (3.6). 由上述推導可知，在 SPH 方法中，函數的數值積分法和光滑函數 近似法為二階精度。若函數非偶函數，或不滿足歸一性，則光滑函數 近似法不一定有二階精度。函數的導數積分表示法( ∇ ⋅ g ( x) )可應用式 (3.3)，將 ∇ ⋅ g ( x) 取代式中的 g ( x) 得 < ∇ ⋅ g ( x) >=. ∫ [∇ ⋅ g ( x ')]W ( x − x ', h)dx '. (3.7). Ω. 對上述積分內可改寫為 [∇ ⋅ g ( x ')]W ( x − x ', h) = ∇ ⋅ [ g ( x ')W ( x − x ', h)] − g ( x ') ⋅ ∇W ( x − x ', h) (3.8). 將式(3.8)代入式(3.7)可得 < ∇ ⋅ g ( x) >=. ∫ ∇ ⋅ [ g ( x ')W ( x − x ', h)]dx ' − ∫ g ( x ') ⋅∇W ( x − x ', h)dx ' (3.9). Ω. Ω. 方程(3.9)右邊第一項應用散度定理將在體域 Ω 上的積分變換為在面 12.
(27) 域 S 上的積分 <= ∇ ⋅ g ( x) >. . ∫ g ( x ')W ( x − x ', h) ⋅ n dS − ∫ g ( x ') ⋅∇W ( x − x ', h)dx ' (3.10) Ω. S. 這裡 n 為表面 S 上的單位法向量，方向朝外。 由於光滑函數 W 是定義在緊支域內，故當支持域位於問題域內 部時(圖 3-1)，在式(3.10)等號右邊面積分項為零。若支持域與問題 域交錯則光滑函數 W 被邊界截斷(圖 3-2)，這種情況下，面積分項 還用零來處理，則邊界效應就須用特殊方法來彌補。 因此對支持域在問題域內的問題式(3.10)可化簡為 < ∇ ⋅ g ( x) >= − ∫ g ( x ') ⋅∇W ( x − x ', h)dx '. (3.11). Ω. 由式(3.11)可知，在函數 g ( x) 導數的核近似方法中，可化轉成光滑函 數 W ( x − x ', h) 的導數積分形式，即函數的導數可以用場函數值和光滑 函數導數來計算。 3.2.2 粒子近似法 SPH 計算式建立在有限個離散粒子上，每個粒子攜帶一定訊息， 如:質量、密度、壓力、速度… …等等。粒子近似法就是把連續積分 方程轉換為支持域內所有粒子疊加求和的離散化形式。如(圖 3-3)所 示 Ω 表支持域的體積， S 表支持域外表面面積。 13.
(28) 粒子的質量 m j 可表示為: m j = ∆V j ρ j. (3.12). 其中 ρ j 表粒子j的密度，N表粒子j支持域內的粒子數。 g ( x) 的連續SPH積分表示式可改寫成以下的離散化粒子近似式:. = < g ( xi ) >. N. x j )Wij ∆V j ∫ g ( x ')W ( x − x ', h)dx ' ≈ ∑ g (=. Ω. N. mj. ∑ρ. =j 1 =j 1. g ( x j )Wij (3.13). j. 其中 = Wij W ( xi − x j , h). (3.14). 式(3.13)說明了粒子i處任意函數值可通過光滑函數對其緊支域內所 有粒子相對應的函數值進行加權平均近似。同理可得函數空間導數的 粒子近似式為: = < ∇ ⋅ g ( x) >. N. mj. ∑ρ j =1. g ( x j ) ⋅∇Wij. (3.15). j. 其中 xi − x j ∂Wij ∇iWij = ∂rij rij. (3.16). 3.3 支持域與影響域 在前面的敘述中，使用了支持域的概念。另外影響域的概念也被 14.
(29) 廣泛的應用，上述兩種概念的使用，程式編碼不同，然而，這兩個觀 念十分容易混淆，故在此做進一步的區分。 i.. 支持域:給定一點 x ， x 的支持域定義為此範圍內所有點的訊 息都被用於決定 x 點的訊息。見圖 3-4， x 的支持域為虛線圓 所圍起的範圍，點 1、2、3 在 x 的支持域內，而點 4、5 在 x 的 支持域外，故計算時使用點 1、2、3。. ii.. 影響域: x 的影響域定義為此範圍內所有點的訊息都將被 x 點 的所影響，在計算 x 點訊息時需先判斷其落於那些點的影響 域內。見圖 3-5， x 點的訊息將由點 1、3、4 所決定。. 而 SPH 方法中，粒子同時具有支持域與影響域，且由光滑長度 h 乘以因子 κ 所決定。 在 SPH 的歷史中，這兩種不同的概念產生了兩種不同的模型，分 別為分散模型與集中模型。 i.. 分散模型:通過應用那些影響域覆蓋到粒子 x 的粒子來對 x 進 行近似計算。. ii.. 集中模型:通過應用粒子 x 支持域內的粒子來對 x 進行近似運 算。. 15.
(30) 由於兩個粒子的光滑長度不一定相同，故有可能發生以下情況，粒子 i 位於粒子 j 的影響域內，但粒子 j 卻不位於粒子 i 的影響域內，這種. 非對稱性的現象會造成模擬失真[14]。見圖 3-6 粒子 j 對粒子 i 施加了 一作用力，但明顯的粒子 j 並未落於粒子 i 的影響域內，所以粒子 i 未 對粒子 j 施加作用力，明顯違反了牛頓第三定律。為了解決此種問題， 對相互作用的粒子對取其光滑函數的平均值來定義粒子的支持域與 影響域。而在這種情況下，分散模型和集中模型是一樣的。本文採用 此法，故粒子的支持域與影響域是一樣的。. 3.4 光滑函數 無網格方法中的一個關鍵問題是如何在一系列不需要預先定義 網格來提供節點連接的任意分佈的離散點上有效的進行函數近似。函 數的近似方法主要有三類:積分表示法、級數表示法、微分表示法。 SPH 是以光滑函數引進積分式，因此光滑函數非常重要，它不但決定 了函數近似的形式、定義粒子支持域的尺寸，還決定了核近似和粒子 近似的一致性跟精度。光滑函數 W ( x − x ', h) 需具有以下特性: 1. 光滑函數需滿足正則化(歸一性)。 1 ∫ W ( x − x ', h)dx ' =. Ω. 2. 當光滑半徑趨近於零時具有 Dirac 函數的性質。 16. (3.17).
(31) lim W ( x − x ', h) =δ ( x − x ') h →0. (3.18). 3. 為節省計算量、提高程式效率，使 SPH 近似式由全區計算變為局 部計算，光滑函數需滿足緊支性，即在支持域內才有作用。 W ( x −= x ', h) 0. 當 || x − x ' ||> κ h. (3.19). κ 為與光滑半徑有關的係數。 4. 光滑函數隨粒子間距單調遞減。 5. 在點 x 上的粒子的支持域內任意點 x ' 處存在 W ( x − x ', h) ≥ 0 (非負 性)。 6. 光滑函數為偶函數(對稱性質)。 7. 為了獲得更好的近似結果，光滑函數需充分光滑。因光滑的光滑 函數對質點的不規則分佈不敏感。 3.4.1 光滑函數的求解步驟 場函數的近似 對 g 場函數進行近似，得積分描述形式如下所示: = < g ( x) >. ∫ g ( x ')W ( x − x ', h)dx '. Ω. 若 g ( x) 充分光滑，在 x 點附近對進行泰勒級數展開可得:. 17. (3.20).
(32) 1 g ''( x)( x '− x 2 ) + ... 2 n (−1) k h k g k ( x) x − x ' k x − x' ( ) + rn ( ) = ∑ k! h h k =0 g ( x ')= g ( x) + g '( x)( x '− x) +. 其中 rn (. (3.21). x − x' ) 是餘數項。 h. 將式(3.21)代入式(3.20)可得 = g ( x). x − x' (−1) k h k g k ( x) x − x ' k ( ) W ( x − x ', h)dx ' + rn ( ) ∫Ω ∑ k! h h k =0 n. x − x' (−1) k h k g k ( x) x − x ' k ( ) W ( x − x ', h)dx ' + rn ( ) (3.22) = ∑ ∫ k! h h k =0 Ω n. =. n. ∑A g k =0. k. k. ( x) + rn (. x − x' ) h. 其中 = Ak. (−1) k h k x − x ' k ( ) W ( x − x ', h)dx ' k ! Ω∫ h. (3.23). 比較式(3.22)的左右兩邊，知要達 n 接近似須滿足下述條件: A0 =W ( x − x ', h)dx ' = 1 ∫ Ω x − x' A = )W ( x − x ', h)dx ' = 0 −h ∫ ( 1 h Ω h2 x − x ' 2 A ', h)dx ' 0 = 2 ∫Ω ( h ) W ( x − x= 2! (−1) n h n x − x ' n A ( ) W= ( x − x ', h)dx ' 0 = n n ! Ω∫ h . 式(3.24)又可進一步地化簡為 18. (3.24).
(33) M 0 =W ( x − x ', h)dx ' = 1 ∫ Ω 0 M1 = ∫ ( x − x ')W ( x − x ', h)dx ' = Ω 2 0 M 2 =− ∫Ω ( x x ') W ( x − x ', h)dx ' = M n =− ( x x ') n W ( x − x ', h)dx ' = 0 ∫ Ω . (3.25). 其中式(3.25)的第一項即為歸一性條件，第二項為對稱性條件。滿足 此二項即可確保函數的 SPH 核近似的一階連續性。 場函數導數的近似 計算流體力學中，偏微分控制方程中的場函數最高階導數是二階 的，而彈性力學中，最高階導數項是四階，在此推導一階導數近似， 高階依此類推。通過(3.20)微分得 = < g '( x) >. ∫ g '( x ')W ( x − x ', h)dx '. (3.26). Ω. 接著使用部分積分 = < g '( x) >. . ∫ g ( x ')W ( x − x ', h) ⋅ ndS − ∫ g ( x ')W ( x − x ', h)dx ' '. Ω. S. 將式(3.21)代入(3.27)等號右邊第二項即得. 19. (3.27).
(34) . ∫ g ( x ')W ( x − x ', h) ⋅ ndS. '( x) g=. S. (−1) k h k g ( k ) ( x) x − x ' k x − x' ( ) + rn ( ))W ' ( x − x ', h)]dx ' − ∫ [(∑ k! h h Ω k =0 (3.28) = ∫ g ( x ')W ( x − x ', h) ⋅ ndS n. S. (−1) k h k g ( k ) ( x) x − x ' k ' x − x' −[∑ ( ) W ( x − x ', h)dx ' + rn ( )] ∫ k! h h k =0 Ω n. =. . ∫ g ( x ')W ( x − x ', h) ⋅ ndS + A. ' k. g ( k ) ( x) + rn (. S. x − x' ) h. 其中 = A'k. (−1) k h k g ( k ) ( x) x − x ' k ' ∫ ( h ) W ( x − x ', h)dx ' k! Ω. (3.29). 顯然若以下方程能滿足，則 g '( x) 達 n 接近似: M 0 = W ' ( x − x ', h)dx ' = 1 ∫ Ω ( x − x ')W ' ( x − x ', h)dx ' = 0 M1 = ∫ Ω 2 ' 0 M 2 =− ∫Ω ( x x ') W ( x − x ', h)dx ' = M n =− ( x x ') n W ' ( x − x ', h)dx ' = 0 ∫ Ω . (3.30). 和 W ( x − x ', h) |s = 0. (3.31). 按照同樣的推導流程，可知一個 k 階導數欲達到 n 階近似精度， 其緊支條件為: 20.
(35) W ( x − x ', h) |s = 0 W ' ( x − x ', h) |s = 0 . (3.32). W k −1 ( x − x ', h) |s = 0. 若光滑函數僅與光滑半徑有關，則光滑函數可寫成下式: W ( x − x ', h) = W ( R) = a0 + a1 R + a2 R 2 + ... + an Rn. (3.33). 對式(3.33)，若 R = 0 處光滑函數的第 i 階導數存在可得 i !ai ; x − x ' → 0+ W (0) = i − (−1) i !ai ; x − x ' → 0 i = 1, 2,..., n (i ). (3.34). 由上式可知只有下述成立時， W (i +1) ( x − x ', h) 才在 R = 0 處存在: ai 0;= i 奇數 = ai 0;= i 偶數 =. (3.35). 若光滑函數的二階導數存在，則 a1 = 0 ，故由(3.35)式(3.33)可改寫成 W ( x − x ', h) = W ( R) = a0 + a2 R 2 + ... + an R n. (3.36). 再由式(3.30)和式(3.32)即可求得係數。 3.4.2 常用的光滑函數 常用的有三種，分別為高斯型光滑函數、分段三次樣條光滑函數、 高次樣條光滑條函數。而本文採用的為分段三次樣條光滑函數。 高斯型光滑函數 21.
(36) Gingold 和 Monaghan(1977[14])最早使用以下高斯型光滑函數來 仿真非球型星體，如下所示(圖 3-7): W ( R, h) = α d e − R. r h. 其中 R= =. 2. (3.37). 1 | x − x'| 1 ，在一維空間中 α d = 1 ，在二維空間中 α d = 2 ， h πh π 2h. 在三維空間中 α d =. 1 3 2. π h. 。 3. 高斯型光滑函數是充分光滑，穩定且精度很高，特別是對於不規 則粒子分佈的情況，所以它是一種較好的選擇。但是，高斯型光滑函 數計算量要求很大，因其趨於零時要經歷的距離較大，特別是對於光 滑函數的高階導數。這樣就會產生一個很大的支持區域，區域中域於 粒子近似的粒子也會很多，導致計算時間較長。因此 Morris 對高斯 型光滑函數進行了修正，修正後的高斯光滑函數比原始高斯型光滑函 數更穩定、具有更高的效率，以下是二維修正高斯光滑函數。 e − R − e − (δ / h ) 2. W ( R , h, δ ) =. δ. 2π ∫ r (e. − R2. −e. 2. − (δ / h ) 2. (3.38) )dr. 0. r h. 其中 R= =. | x − x'| ， δ 是截斷半徑，在二維空間中一般取 δ = 3h 。 h. 分段三次樣條光滑函數(Cubic spline). 22.
(37) Monaghan 和 Lattanzio(1985)提出分段三次樣條光滑函數，如下所 示: (圖 3-8) 1 3 2 2 3 − R + 2 R 0 ≤ R<1 1 3 1 ≤ R<2 ) α d × (2-R) W ( R, h= 6 0 R≥2 . r h. 其中 R= =. (3.39). | x − x'| 1 15 ，在一維空間中 α d = ，在二維空間中 α d = ， h 7π h 2 h. 在三維空間中 α d =. 3 。 2π h3. 到目前為止，三次樣條光滑函數是應用最廣泛的光滑函數，因為 在狹窄的緊支域中，三次樣條光滑函數與高斯型光滑函數類似。但由 於其光滑函數是分段的，使用起來會比只有一段的光滑函數困難一 些。 分段五次樣條光滑函數 Morris(1997[17])提出了更加近似於高斯型光滑函數並更加穩定 的高次樣條函數(四次和五次)，其中分段五次樣條光滑函數如下所示: (3 − R)5 − 6(2 − R)5 +15(1 − R)5 5 5 (3 − R) − 6(2 − R) W ( R, h= ) αd × 5 (3 − R) 0 . 23. 0 ≤ R ≤1 1≤ R ≤ 2 2 ≤ R<3 R ≥3. (3.40).
(38) r h. 其中 R= =. 7 | x − x'| 120 ，在一維空間中 α d = ，在二維空間中 α d = ， 478π h 2 h h. 在三維空間中 α d =. 3 。 359π h3. 3.5 控制方程式之離散 在第二章中已經建立了本文的控制方程式(2.1)-(2.3)，在本節將使 用 SPH 法得到離散式，接著就可以編譯程式進行模擬。 3.5.1 密度計算 由於 SPH 方法中粒子運動是由相鄰粒子間的密度差與壓力差產 生，故在 SPH 法中密度近似法非常重要。本文使用在 SPH 法中對密 度進行展開的常用方法，正則密度求和法。對任意粒子 i ，由式(3.13) 可寫為: N. mj. j =1. ρj. < ρ ( xi ) >= ∑. ρ ( x j )Wij. (3.41). 整理過，其密度可被表達成下式 N. ρi = ∑ m jWij. (3.42). j =1. 其 中 N 為 粒 子 i 支 持 域 內 粒 子 的 總 數 ， mj 為 粒 子 j 的 質 量 ， = Wij W ( xi − x j , h) 為粒子 j 對粒子 i 產生影響的光滑函數。. 接著將式(3.42)等號右端正則化得: 24.
(39) N. ρi =. ∑m W j =1 N. j. mj. ∑ρ j =1. ij. (3.43) Wij. j. 簡而言之，式(3.42)稱密度求和法，意義為粒子密度可通過其支持 域內所有粒子密度的加權平均近似得到。在 SPH 密度近似法中，應 用較為廣泛，因為它體現了 SPH 近似法的本質。優點是在整個積分 區域中，積分嚴格的遵循值量守恆定律，但當粒子位於流體區域邊界 上時，會產生邊界粒子缺陷，導致結果出錯，且計算量大，因為在計 算其他參數前需先計算密度和自身的光滑函數。現今有不少修正的方 法，其一為正則密度求合法式(3.43)可提高自由邊界處和相同材料粒 子密度不連續交界面處的精度。此式適合模擬非連續廣義流體的流動 問題。 除了上述算法，亦可將連續方程(3.44)以 SPH 方法進行離散得式 (3.45)，稱為連續密度近似法。 dρ = − ρ∇ ⋅ v dt. d ρi = dt. (3.44). ∂Wij. N. ∑ m v β ⋅ ∂x β j =1. j ij. (3.45). i. 連續性密度近似法的優點在於計算效率較高，並適用於平行運算法編 程，主要原因為其不需先計算密度，因此可能會不遵守質量守恆定 25.
(40) 律。 3.5.2 動量方程 拉格朗日型的動量守恆方程式可備表達為: dvα 1 ∂σ αβ = ρ ∂x β dt. (3.46). 代入 SPH 公式後得 αβ σ iαβ σ j ∂Wij m j ( 2 + 2 ) β + Fext ∑ ρi ρ j ∂xi j =1. dviα = dt. N. (3.47). − pδ αβ + τ αβ 代入式(3.47)得如下: 將 σ αβ = αβ N µiε iαβ µ jε j ∂Wij dviα pi p j ∂Wij N = −∑ m j ( 2 + 2 ) α + ∑ m j ( 2 + + Fext (3.48) ) ρi ρ j ∂xi ρi ρ 2j ∂xiβ dt j 1= j 1 =. 上式中，等號右邊第一項為壓力近似式而第二項為黏滯力包含物理黏 滯係數，其中應力張量又可表達為 ∂v β ∂vα 2 + − (∇ ⋅ v)δ αβ ∂xα ∂x β 3. αβ ε=. (3.49). 將式(3.49)代入 SPH 公式，可得 = ε iαβ. N. mj. ∑ρ. v βji. ∂Wij. N. +∑. ∂xα. =j 1 = j 1 j i. mj. ∂Wij. 2 N mj ∑ v ji ⋅∇iWij )δ αβ (3.50) ∂xiβ 3 j 1 ρj ρ j=. α. 使用(3.48)和(3.50)即可求得 dvi 。 dt. 26. vαji. −(.
(41) 3.5.3 狀態方程 在求解 SPH 方法中，粒子的運動是由於壓力梯度的作用產生，而 粒子的壓力式通過狀態方程由粒子自身的密度和內能計算而來。在此 使用狀態方程: p = c2 ρ. (3.51). 其中 c = 0.1 。Morris[17]使用此人工狀態方程通過 SPH 方法數值模擬 了不可壓縮流動。Zhu 等使用此人工狀態方程通過 SPH 方法用於微孔 尺度的數值模型，用於數值模擬穿越多孔介質的流動。. 3.6 粒子的位移修正法 為了使粒子保持整齊，不絮亂，在此使用粒子位移修正法，由 Monaghan[14]提出，應用於不可壓縮流體時非常有效。該方法中，流 體粒子的運動方式如下所式: N m dxi j = vi − ε ∑ xijWij dt j =1 ρ j. (3.52). 式中， ε 為常數，取值範圍在 0 ≤ ε ≤ 1.0 。 由式(3.52)可知，該公式考慮了相鄰粒子的影響，可使粒子 i 的運 動速度與相鄰粒子的平均速度相近，可使粒子保持較為整齊。. 27.
(42) 3.7 邊界修正法 前面提到，當支持域完全落在研究區域內時，由面積分項等於零 可得式(3.11)。然而，欲探討的情況往往支持域不完全落在研究區域 內，必須使用補償方法，稱邊界處理。 固壁邊界的處理辦法，最常用的三種方法分別為 Monaghan 提出 的邊界力法，Libersky 等人提出的鏡像邊界法，還有 Dalrymple 等人 提出的耦合邊界法。至今，仍未有一種較好且較通用的邊界處理法。 [13][15] 3.7.1 邊界力法 此法最早由 Monaghan 提出，在邊界附近放入邊界虛粒子，當流 體粒子接近邊界時，邊界虛粒子將會提供斥力用以防止流體粒子的非 物理性穿透(圖 3-9)。最早使用的邊界力函數是參考分子動力學中 Lennard-Jones 是能函數， r r0 n1 r0 n2 D ( ) − ( ) as r < r0 g (r ) = r 2 r r 0 as r ≥ r0 . (3.53). 其中 D 是常數， r0 是支持域半徑， n1 、 n2 常使用 12 和 4。上式可清楚 得到當粒子間距趨近於零，斥力將趨近於無窮大。此法 D 、 r0 、 n1 、 n2 的選擇對計算結果影響非常大，如 r0 選擇過大，會使得初始粒子受到 28.
(43) 不必要的干擾，甚至爆破。而選擇過小，會使流體粒子直接穿透邊界。 3.7.2 耦合邊界法 由 Dalrymple 等人提出，此法是把固壁邊界離散成一系列的邊界 點(圖 3-10)，令這些點具有流體質點的部分屬性(質量、密度、壓力 等)，在求解方程組時，一併計算，而邊界虛粒子的厚度由光滑長度 決定。耦合邊界法的優點在於無須對邊界進行特殊處理，邊界對流體 質點的力已經耦合在求解方程的過程中，缺點為得到的邊界壓力並非 真實壓力。 3.7.3 鏡像邊界法 此法在在邊界粒子的另一側產生鏡像粒子，其速度大小相同方向相反， 鏡像粒子可在第一部生成後固定不動，也可一段時間更新一次，視需 求選擇，但此方法較適合用於線形邊界(圖 3-11)。. 3.8 相關程式原理與步驟 3.8.1 粒子間相互作用的對稱化 在前面已經提到了支持域與影響域的觀念，若光滑長度是隨著時 間、空間變化，則每個粒子都具有獨立的光滑長度。假設 hi ≠ h j ，那 麼粒子 i 的影響域可能包含粒子 j ，但反之則不一定，這與牛頓第三 29.
(44) 運動定律矛盾。欲解決此問題，必須設法保持粒子間相互作用的對稱 性。 對光滑長度進行修正為一可行的方法。因此，本文對光滑長度取 幾何平均數進行修正 hij =. 2hi h j hi + h j. (3.54). 3.8.2 相鄰粒子搜索法 在粒子 i 的計算過程中，只需使用被緊支域所包含的粒子(有限個)。 因此如何找尋這些粒子是一個重要課題。一般將這些被包含於緊支域 的粒子稱為粒子 i 的最近相鄰粒子(NNP)。而尋找最近相鄰粒子的過程 稱為相鄰粒子搜索法(NNPS)。下面介紹兩種常用的方法: 3.8.2.1 全配對搜索法(all-pair search) 全配對搜索法是一個直接又簡單的相鄰粒子搜索法方法。如圖 3-12 所示，給定粒子 i，計算粒子 i 到粒子 j 的距離 rij，其中 j = 1, 2..., N 。 若 rij < κ h ，則粒子 j 為粒子 i 的最近相鄰粒子。假使粒子 i 的緊支域半 徑同於粒子 j ，粒子 i 亦為粒子 j 的最近相鄰粒子。顯然，此法複雜度 的階數精度為 O( N 2 ) ，且每個時間步的計算都要進行 NNPS，故所耗 費的計算時間長，僅適用於粒子數量少的問題。 30.
(45) 3.8.2.2 鏈表搜索法(linked-list search algorithm) 鏈表搜索法相較於全配對搜索法的效率較高。首先，在問題域上 鋪設臨時網格，網格大小與支持域的空間大小一致(支持域尺度為 κ h ， 則網格尺度亦為 κ h )，見圖 3-13。對給定的粒子 i 而言，其相鄰粒子 必在同一網格或相鄰的網格內，所以在一維、二維、三維空間的搜尋 範圍分別是 3、9、27 個單元內，接著在計算粒子 i 到粒子 j 的距離 rij ， 其中 j ∈ 搜尋單元。若 rij < κ h ，則粒子 j 為粒子 i 的最近相鄰粒子。因 粒子 i 的緊支域半徑同於粒子 j ，粒子 i 亦為粒子 j 的最近相鄰粒子。 若每個單元內的平均粒子數足夠少，則此法複雜度的階數為 O( N ) 。 但當每個粒子的光滑長度不一致，抑或光滑長度可變時，網格空間就 不適用於所有粒子，故此法適用於每個粒子的光滑半徑為相同常數的 條件下。 3.8.2.3 粒子對的相互作用 前述了找尋相互作用粒子對的方法，在 SPH 模擬中，定義並記錄 所有相互作用粒子對是非常重要且耗時的一部分。下面敘述程序: //初始化 ninteraction=0. //初始化總交互作用對數. 31.
(46) do i=1,ntotal count(i)=0 //初始化每個粒子的相鄰粒子數 end do //開始記錄 do i=1,ntotal-1 do j=i+1,ntotal if(abs(x(i)-x(j))<=mhsml)then. //nnps 檢驗條件. ninteraction=ninteraction+1 //記錄總交互作用對數 pair_i(i)=i. //儲存相互作用粒子對 i 裡的第一個粒子標號. pair_j(i)=j. //儲存相互作用粒子對 i 裡的第二個粒子標號. count(i)=count(i)+1 //紀錄粒子 i 的相鄰粒子數 count(j)=count(j)+1 //紀錄粒子 j 的相鄰粒子數 endif enddo enddo 上述循環中，ntotal 為粒子總數;ninteraction 為交互作用粒子對數; 數列 pair_i 和 pair_j 為儲存交互作用粒子對的第一及第二粒子標 號;count 用來記錄粒子 i 的相鄰粒子數。如此通過一個循環即可記錄 32.
(47) 所有的粒子對，在使用時將數列讀入即可。 為舉例如何使用，下面給出密度求和法的程序: do i=1,ntotal density(i)=mass(i)*W(0) enddo do k=1,ninteraction //粒子對 i~ninteraction i=pair_i(k) //讀取交互作用粒子對 k 的第一個粒子 j=pair_j(k) //讀取交互作用粒子對 k 的第二個粒子 density(i)= density(i)+mass(j)*W(k) density(i)= density(i)+mass(j)*W(k) enddo 粒子近似法的其他變量的求和法皆可應用此程序實現。使用成對 粒子法可有效的結省儲存量和提高計算效率。 3.8.3 週期性邊界條件 本文使用的數值方法為粒子法，離散的粒子數目是影響運算效率 的主要因素，如何以有限顆粒子模擬無限長流道，在此介紹週期性邊 界條件(Periodic Boundary Conditions, PBC)，其為邊界條件的一種， 即在合理的範圍內假設模擬系統的重複性，利用邊界條件替代所選擇 33.
(48) 的系統受周邊環境的影響。在粒子運動過程中，當粒子即將離開指定 區域穿越邊界時，粒子就會受到週期性邊界作用回到指定區域，因此 在進行積分和粒子間相互作用時都需考慮週期性邊界的影響。以管流 為例，若使用週期性邊界條件，即表示所模擬的系統是真實管道的一 小段，忽略了出入口等效應在模擬所造成的影響。 見圖 3-14，灰色區塊為主模擬空間(Primary Cell)，而四周區塊為 主模擬空間所複製出來的映象空間(Image Space)。對一維模型有 2 個 映像空間，二維模型有 8 個，三維模型則有 26 個映像空間。映像空 間內的粒子和主模擬空見內粒子只有座標位置不同，速度的大小、方 向皆相同，即 v = v ' 。在固定時間步階檢查粒子是否跑出模型外，若跑 出模型外則處理法如下: i.. 若 xi > L ，則實際位置 x p= xi − L 。圖 3-15. ii.. 若 xi < L ，則實際位置 x p= xi + L 。圖 3-16. 此法可大幅減少模擬所需的粒子數及增加計算效率。 3.8.4 最小映像法則 對一個採用週期性邊界的系統而言，除了計算模型內的真實粒子 與真實粒子作用力外，也須考慮模型內的真實粒子與映像空間內之假. 34.
(49) 想粒子的作用力。計算作用力時，模型內的粒子與映像空間內的粒子 之間距必須找到最小作用距離，此法稱最小映像法則(minimum image convention)。欲使用此法，主模擬空間的尺寸必大於兩倍的支持域半 徑 2κ h 。舉一維空間為例，以下為處理程序: i.. L 2. 若 xij < ，則 xij < xij ' 且 xij < xij '' 。粒子 i 和粒子 j 的距離取 xij 。見 圖 3-17. ii.. L 2. 若 xij < − ，則 xij '' > xij > xij ' 。粒子 i 和粒子 j 的距離取 xij ' 。見圖 3-18. iii.. 若 xij >. L ，則 xij '' < xij < xij ' 。粒子 i 和粒子 j 的距離取 xij '' 。見圖 2. 3-19 二維空間和三維空間處理法與上述相同。 3.8.5 紅血球變形指數(DI 值) 2001 年 Tsukada 等人定義 Deformation Index，為現今測量紅血球 變形能力的指標之一。其定義如下:[25] deformation index( DE ) =. 其中 l 為長; d 為寬。. 3.9 程式流程 35. l d. (3.55).
(50) 前文已給出數學模型及其數值離散過程。在程式模擬基本上可分 為以下幾點(具體程式流程圖見圖 3-20、程式結構見圖 3-21): 1. 先離散整個問題域，並且初始化流場，給定質點初始量，如:位置、 速度、質量、密度等。 2. 細胞膜以彈簧粒子模型表達，並計算所受的彈簧力 F(當作外力項 計算)。 3. 求解式(3.42)-(3.50)。 4. 時間積分(跳蛙演算法)，計算下一個時間步階個質點的物理量。 5. 重複 2-4 步驟，直至結束。. 36.
(51) 圖 3-1 支持域位於問題域內. 圖 3-2 支持域與問題域交錯 37.
(52) 圖 3-3 粒子法求解示意圖 求解粒子 i 可通過使用光滑函數 W 支持域內的粒子來進行，其中支持 域為半徑 κ h 的圓. 圖 3-4 支持域示意圖 虛線為 x 的支持域，計算 x 訊息時所使用的點為 1 、 2 、 3 點 38.
(53) 圖 3-5 影響域示意圖 虛線為各點的支持域，計算 x 訊息時所使用的點為 1 、 3 、 4 點. 圖 3-6 影響域不同示意圖 實線為為粒子 j 的影響域範圍虛線微粒子 i 的影響域範圍. 39.
(54) 圖 3-7 高斯函數及其一階導數. 圖 3-8 三次樣條函數及其一階導數 40.
(55) 圖 3-9 邊界力法示意圖. 圖 3-10 耦合邊界法示意圖. 41.
(56) 圖 3-11 鏡像邊界法示意圖. 圖 3-12 在二維空間中應用全配對搜索法 此圖中灰色粒子具 4 個相鄰粒子。. 42.
(57) 圖 3-13 在二維空間中應用鏈表搜索法 此圖中對黑色粒子而言，搜尋範圍為灰色的九個網格。. 圖 3-14 週期性邊界示意圖 利用週期重複以組成真實系統 43.
(58) 圖 3-15 一維週期性邊界示意圖( xi > L ). 圖 3-16 一維週期性邊界示意圖( xi < L ). 圖 3-17 一維最小映像法則示意圖( xij <. 44. L ) 2.
(59) L 2. 圖 3-18 一維最小映像法則示意圖( xij < − ). L 2. 圖 3-19 一維最小映像法則示意圖( xij > ). 45.
(60) 開始. 以粒子離散問題域. 計算薄膜所受的彈簧力. No. 以SPH法計算控制方程. 模擬下一個步階. Yes. 結束. 圖 3-20 程式流程圖. 46.
(61) 初始化模組 產生邊界虛 粒子 搜尋相鄰粒子 A. 全配對搜索法 B.鏈表搜索法 計算光滑函數 A.高斯函數 B.三次樣條函數 C.五次樣條函數 下一個時間步階. 密度求和法 計算彈簧力 計算內力. 計算人工黏滯力. 計算動量改變率 更新物理量 (X、V…) 輸出. 圖 3-21 程式結構圖. 47. 計算外力.
(62) 第4章. 數值模擬與結果討論. 4.1 基本流體算例驗證 以上推導了 SPH 形式的流體動力學模型，並說明其數值求解過程。 本文以模擬二維 Couette 流動和 Poiseuille 流動為基礎，對所編寫的程 式進行驗證。本算列中問題域的大小及粒子的分布和 G.R.Lui[11]所提 供的數據一致。 4.1.1 Couette 流動 本算例模擬二維 Couette 流動，流體在兩塊固定並水平放置在 y = 0 和 y = l 處的無限大平板中運動，y = l 處之平板突然以速度 v0 向右運動。 根據 Morris[17]，可以得到 Couette 流動隨時間的變化之解析解 Vx1 ( x2 , t ) 。. V0 x2 l ∞ 2V nπ n 2π 2ν z ) exp(− 2 t ) + ∑ 0 (−1) n sin( l l n = 0 nπ. Vx1 ( x2 , t ) =. (3.56). 本節中， l = 10−3 m ， ρ = 103 kg / m3 ， υ = 10−6 m 2 / s ， v= 2.5 ×10−5 m / s ， 0 相對應的雷諾數 Re = 2.5 × 10−2 。圖 4-1 為 SPH 粒子配置圖，問題域為 0.0005m × 0.001m 的長方形，如圖所示，長方形的上下邊分別放置邊界. 粒子，在計算過程中，上邊界粒子向右運動，下邊界粒子保持不動， 48.
(63) 左右兩邊界施加週期性邊界。圓形代表流體粒子，長 20 顆，寬 40 顆， 共 800 顆流體粒子，方形代表邊界粒子，上下邊界各 60 顆，共 120 顆邊界粒子，總粒子數為 920 顆。 當達穩態時，粒子的分佈圖如圖 4-2 所示。由圖可見，因週期性 邊界作用，有些本應在下游的粒子會運動到上游。而圖 4-3 為不同時 刻流體速度分布圖，SPH 算法所得的結果與級數求解法相當符合。 4.1.2 Poiseuille 流動 本例模擬二維 Poiseuille 流動。流體在兩塊固定並水平放置在 y = 0 和 y = l 處的無限大平板中運動，兩板固定不動。 根據 Morris[17]，可以得到 Poiseuille 流動隨時間的變化之解析解 Vx1 ( x2 , t ) 。. Vx1= ( x2 , t ). F x2 (l − x2 ) 2ν (3.57) ∞ 4 Fl 2 (2n + 1) 2 π 2ν πz t) sin( (2n + 1)) exp(− +∑ 3 3 l l2 n = 0 νπ (2n + 1). 在 本 例 中 ， l = 10−3 m ， ρ = 103 kg / m3 ， υ = 10−6 m 2 / s ， 傳 動 力 Fext = 2 ×10−4 m / s 2 ，由穩態時的速度解析解式(3.58) 1 dp l 2 vx ( y ) = ( − y2 ) − 2ν dx 4. (3.58). 可得流體前沿速度為 2.5 ×10−5 m / s ，相對應的雷諾數 Re = 2.5 × 10−2 。 49.
(64) 圖 4-4 為 SPH 粒子配置圖，問題域為 0.0005m × 0.001m 的長方形，如圖 所示，長方形的上下邊分別放置邊界粒子，在計算過程中，邊界粒子 保持不動，左右兩邊界施加週期性邊界。當達穩態時，粒子的分佈圖 為圖 4-5 所示。由圖可見，因週期性邊界作用，有些本應在下游的粒 子會運動到上游。而圖 4-6 為不同時刻流體速度分布圖，SPH 算法 所得的結果與級數求解法符合。. 4.2 紅血球於微流道中變形模擬與結果討論 紅血球是血液的主要成分，占總血量的 45%，它在血流中的變形 可以做為判斷正常細胞和病變細胞的一個標準，故其變形的特性是研 究臨床醫學和生物力學實驗的重要問題。 在本章中，先模擬細胞於壓力驅動流中的變形問題來驗證以粒子 為基礎之模型的基礎能力，並分析不同密度比及黏滯係數比對其變形 的影響。在驗證了以粒子為基礎的模型後，本文進一步模擬細胞進入 微流道的運動變形並探討不同剛度和黏滯力對其影響。 4.2.1 在壓力驅動流中二維紅血球的變形模擬 紅血球在壓力驅動流中的變形運動與紅血球在毛細血管中的變 形相似。其中最重要的發現便是紅血球會變形為特有的降落傘形. 50.
(65) (parachute shape)[22][25]以便通過毛細血管。此問題不單只是驗證數 值模型正確性的基本問題，同時對微循環的研究也很重要。在較早的 理論模型和數值模擬，大多假設紅血球的細胞內液和血漿的密度、黏 滯係數是相同的，因此本節在驗證模型正確性後，將模擬不同密度比 和黏度比的對變形的影響。 圖 4-7 為 SPH 粒子配置圖，模擬二維紅血球於壓力驅動流中的變 形，流道道寬度為 10µ m ，略大於紅血球長度，紅血球置於具入口 5µ m 處，長方形的上下邊分別放置邊界粒子，在計算過程中，邊界粒子保 持靜止，左右兩邊界施加週期性邊界。粒子間距為 0.5µ m ，粒子總數 共 2154 顆。 圖 4-8 記錄了紅血球於不同時間的形狀，其由上游往下游運動， 紅血球的形狀一開始為雙凹透鏡型，最終變為前凹後凸透鏡型，此形 狀上的差異是由血液的剪切應力和兩側管壁摩擦力產生的拋物線型 速度分布所造成，由白努力定律得知流道中央速度快、壓力小，因此 細胞膜中央較厚兩翼較薄，此形狀對應實驗觀察，即所謂的降落傘形 狀(parachute shape)，此種特性使紅血球能通過比它更窄的管道。更進 一步比較，圖 4-9 為本文和 Secomb 的比較，兩者特徵一致，皆發展 成降落傘型，但本文模型結果較不對稱，其原因在於 Secomb 使用軸 51.
(66) 對稱模型，故其上下對稱，而本文則是使用平面模型，而形狀和曲度 差異則由於兩者所使用的細胞膜模型不同。 接著模擬紅血球在不同彈簧常數下的變形，分別為 k1 = 8 ×102 N / m、 k2 = 8 ×103 N / m 、 k3 = 8 ×104 N / m ，圖 4-10 為不同彈簧常數紅血球的變. 形比較圖結果與預期相同，彈簧常數越大、抗拉伸的能力就強，因此， 兩翼較短，紅血球變形較不易。改變彈簧常數，可視為改變外部流場 的速度，直接影響細胞膜的變形，圖 4-11 為不同彈簧常數速度比較 圖，當彈簧常數增為原來的十倍時，速度降為原來的 0.43 倍，當彈 簧常數降為原來的十分之一倍時，速度增為原來的 2.3 倍，與預期相 同，彈簧常數越大，速度越小。較早的理論模型和數值模擬，大多假 設紅血球的細胞內液和血漿的密度是相同的，因此本文模擬不同密度 比對紅血球在微流道中的變形的影響。本節檢驗了密度比為 1:0.7、 1:1、和 1:1.5 的情況，如圖 4-12 所示，由於紅血球的變形阻力隨著 密度比的增加而減少，因此紅血球也變得越來越細長，拉伸的幅度也 較大。圖 4-13 為速度比較圖，如預期的紅血球的密度比越大，速度 越快。最後，模擬不同黏滯係數比對紅血球變形的影響，圖 4-14， 紅血球黏滯係數增加，意味著變形阻力增加，故由圖可見其變形較小。 圖 4-15 為速度比較圖，如預期的紅血球的黏滯係數比越大，速度越 52.
(67) 快。 4.2.2 紅血球進入微流道中二維的變形模擬 毛細血管管徑平均為 3µ m ~ 9µ m 。其數量很大，遍部全身，主要功 能為利於血液和組織之間進行物質交換。當血流受阻時，會造成高血 壓、肺動脈高壓等疾病，造成此情況的其中一個原因為紅血球異常， 所以一個廉價且快速的檢驗方法是必要的。透過模擬紅血球通過流道 的時間，來探討各個因子的影響。 圖 4-16 為 SPH 粒子配置圖，進出口流道道寬度為 15µ m ，在距出 入口處 5µ m 時開始變窄直至流道寬度變為 5µ m ，流道的上下邊分別放 置 3 排的邊界粒子，在計算過程中，邊界粒子保持靜止，左右兩邊界 施加週期性邊界。紅血球中心置於 x 軸上距入口 2.5µ m 處，粒子間距 為 0.5µ m ，光滑長度設為 2 倍的粒子間距，粒子總數共 2406 顆。 圖 4-17 記錄了紅血球進入微流道的變形狀況，當紅血球進入微 流道時，兩翼向中心線靠攏並向後延伸，由於流道中央速度快、壓力 小因此中間增厚，最終成一穩定狀態通過毛細血管。而此模擬的變形 過程與文獻結果一致[1]。圖 4-18 分別為 k = 8 ×103 N / m 和 k = 5 ×104 N / m 在 0.001s 、 0.002s 時的位置圖，與預期相同彈簧常數越大，越不易變 形，進入微流道前的變形過程較慢，同時移動速度較慢，可知通過微 53.
(68) 流道的時間較長。圖 4-19 為不同彈簧常數通過的時間圖，明顯可見， 彈簧常數越大，所花的時間越長，但在彈簧常數小於正常值時的變化 差異不大，而當彈簧常數增至約 10 倍時，時間變化大。圖 4-20 為黏 滯係數比 0.67 和密度比 0.67 的比較圖，時間變化符合預期，紅血球 的密度越大、黏滯係數越大則穿越流道的時間越長，但可明顯發現其 變化遠小於彈簧常數的改變。因此，此裝置較適合用來做紅血球剛度 變硬的快篩。. 54.
(69) 圖 4-1Couette 流粒子配置圖 表 4-1 Couette 流資料 空間維度. 2. 長寬. 0.0005m × 0.001m. 動黏滯係數. 1× 10−3 Pa ⋅ s. 密度. 103 kg / m3. v0. 2.5 ×10−5 m / s. 55.
(70) 圖 4-2 Couette 流在 1 秒時粒子的位置速度分布圖. 圖 4-3 Couette 流動的速度分布圖. 56.
(71) 圖 4-4 poiseuille 流粒子配置圖 表 4-2 poiseuille 流資料 空間維度. 2. 長寬. 0.0005m × 0.001m. 動黏滯係數. 1× 10−3 Pa ⋅ s. 密度. 103 kg / m3. 傳動力. Fext = 2 ×10−4 m / s 2. 57.
(72) 圖 4-5 poiseuille 流動在 1 秒時粒子的位置速度分布圖. 圖 4-6 poiseuille 流動的速度分布圖. 58.
(73) 圖 4-7 紅血球在壓力驅動流的 SPH 粒子配置圖. 圖 4-8 紅血球在壓力驅動流中隨時間變化的形變圖 (a) 0.001s (b) 0.005s (c) 0. 01s (d) 0. 015s (e) 0.02s(f) 0.025s 59.
(74) 圖 4-9 紅血球達穩定狀態的比較圖[1]. 圖 4-10 不同彈簧常數紅血球變形比較圖. 60.
(75) 2.5. 速度比. 2 1.5 1 0.5 0 0. 2. 4. 6. 8. 10. 彈簧常數 圖 4-11 不同彈簧常數速度比較圖. 圖 4-12 不同密度比紅血球變形比較圖. 61. 12.
(76) 2.5 2. 速度比. 1.5 1 0.5 0 0. 0.5. 1. 1.5. 密度比 圖 4-13 不同密度比速度比較圖. 圖 4-14 不同黏滯係數比紅血球變形比較圖. 62. 2.
(77) 2.50. 速度比. 2.00 1.50 1.00 0.50 0.00 0. 0.5. 1. 1.5. 2. 黏滯係數比 圖 4-15 不同黏滯係數比速度比較圖. 圖 4-16 紅血球進入微流道的粒子配置圖. 63. 2.5.
(78) 64.
(79) 圖 4-17 紅血球進入微流道( 5µ m )的形變圖. 65.
(80) 圖 4-18 不同彈簧常數通過微流道示意圖. 2.50 2.00. 時間比. 1.50 1.00 0.50 0.00 0.00. 2.00. 4.00. 6.00. 8.00. 10.00. 彈簧常數比 圖 4-19 不同彈簧常數通過的時間圖. 66. 12.00.
(81) 2.50 2.00. 時間比. 1.50 黏滯係數比(0.67) 1.00. 不改變 密度比(0.67). 0.50 0.00 0. 2. 4. 6. 8. 10. 彈簧常數比. 圖 4-20 通過的時間圖. 67. 12.
(82) 第5章. 總結與展望. 5.1 總結 光滑粒子流體動力學法是最早出現的無網格算法，此法應用範圍 極廣，由小尺度的分子模擬至大尺度的天體力學都可使用，且其在於 大變形的問題優於傳統網格。 本文利用此法來模擬二維紅血球在血流中的變形，證明了此法應 用的可靠性及可行性，綜合前面的討論及分析，內容總結為以下幾點: 1. 本文應用光滑粒子流體動力學法的基礎理論(核近似法、粒子近似 法)推導出函數及其空間導數的粒子近似表達式，轉變成 Navier-Stokes 方程組的離散形式。建立紅血球在血流中的二維 SPH 數值模型，並完成電腦程式的建構與基礎流力的印證。 2. 本文進行了二維紅血於壓力驅動流中的變形模擬紅血球一開始 成雙凹透鏡型，受血流影響最終呈降落傘形狀(Parachute-shaped)。 改變細胞膜的彈性速度變化約在 0.4~2.3 倍、密度比和黏滯係數 比速度變化約在 0.8~1.2 倍，皆會影響紅血球的變形，影響最大 的為細胞的膜彈性。 3. 本文也估算出紅血球進入微流道時，兩翼向中心線靠攏並向後延 68.
(83) 伸，中間增厚。通過流道的時間隨彈簧常數的增大而增加之趨勢， 但此表現在常數變小較不明顯。 4. 以紅血球通過微流道之時間作為篩檢方法，較適用於檢測紅血球 變硬的情況。. 5.2 展望 縱使光滑粒子流體動力學法在很多方面都取得不錯的成果，然而 模型和邊界處理的發展仍存在未臻成熟之處，有些根本性的問題尚未 得到理想的解決，尚有很多需要修正及發展的空間。以下提出一些建 議及可探討之方向： 1. 模擬紅血球於彈性管之變形。 2. 將模型發展至三維空間，可得到更真實的模擬。 3. 欲模擬真實的物理現象，需要更多的粒子進行計算，會使得計算 量十分龐大，因此在今後的研究中可以考慮使用平行算法，或與 有限元素等方法耦合計算來提高效率、節省時間。 4. 紅血球的粒子模型可更加精細，加入細胞的微觀結構，如離子通 道、胞器等，進而可以分析細胞病變所造成的病理變化。. 69.
(84) 參考文獻 [1]. Constantine Pozrikidis. (2010). Computational hydrodynamics of capsules and biological cells.. [2]. Fischer, T. M. and Sch¨onbein, H. S. (1977)., Tank tread motion of red cell membranes in viscometric flow: Behavior of intracellular and extracellular markers. Blood Cells, 3 : 351–365.. [3]. Fisher, T. M. (1980)., On the energy dissipation in a tank-tread human red blood cell. Biophysical J., 32 : 863–868.. [4]. Fung, Y. C. (1984) Biodynamics: Circulation. Springer–Verlag, New York.. [5]. Goldsmith, H. L. (1971). Red cell motions and wall interactions in tube flow. Fed. Proc.,30 : 1578–1588.. [6]. Hosseini, S.M., Amanifard, N. (2007). Presenting a modified SPH algorithm for numerical studies of fluid-structure interaction problems. Int. J. Eng. Trans. B 20, 167–178.. [7]. Hosseini, S.M., Manzari, M.T., Hannani, S.K. (2007). A fully explicit three-step SPH algorithm for simulation of non-Newtonian fluid flow. Int. J. Num. Methods Heat Fluid Flow 7, 715–735.. [8]. Hashemi, M. R., Fatehi, R., Manzari, M.T., (2011) SPH simulation of interacting solid bodies suspended in a shear flow of an Oldroy-B fluid. J. Non-Newtonian Fluid Mechanics, 166, 1239-1252.. [9]. Hashemi, M.R., Fatehi, R., Manzari, M.T. (2012) A modified SPH 70.
(85) method for simulating motion of rigid bodies in Newtonian fluid flows. J. Non-Linear Mechanics, 47, 626–638. [10]. Karnis, A., Goldsmith, H. L. and Mason, S. G. (1966). The flow of suspensions through tubes. Canadian J. Chem. Eng., 44 : 181–193.. [11]. Liu, G.R., Liu, M.B., (2005).Smoothed Particle Hydrodynamics: A Meshfree Particle Method.. [12]. Lee, G.Y.H., Lim, C.T., (2007). Biomechanics approaches to studying human diseases. Trends Biotech. 25, 111-118.. [13]. Liu, MouBin., Shao, JiaRu., Chang, JianZong., (2012). On the treatment of solid boundary in particle hydrodynamics. J. Science China Technological Sciences, 55, 244-254.. [14]. Monaghan, J.J., (1992). Smoothed particle hydrodynamics. Annual Review of Astronomical and Astrophysics, 30, 543-574.. [15]. Monaghan, J. J. (1994). Simulating free surfaceflows with SPH. J. Comput. Phys., 110 :399–406.. [16]. Morris, J.P. (1996). Analysis of Smoothed Particle Hydrodynamics with Applications, Department of Mathematics, Monash University.. [17]. Morris, J.P., Fox, P.J., Zhu, Y.J., (1997). Modelling low Reynolds number incompressible flows using SPH. J. Computational Physics, 136, 214–226.. [18]. Monaghan, J.J., Kajtar, J. (2008). SPH simulations of swimming linked bodies. J. Comput. Phys.,, 227, 8568–8587. 71.
(86) [19]. Monaghan, J.J., Kajtar, J.B., (2009). SPH particle boundary forces for arbitrary boundaries. Computer Physics Communications 180 (2009) 1811–1820.. [20]. Pozrikidis, C. (2003). Numerical simulation of the flow. -induced. deformation of red blood cells. Ann. Biomed. 31, 1194–1205. [21]. Pivkin, I.V., Karniadakis, G.E., (2008). Accurate coarse-grained modeling of red blood cells. J. Phys. Rev. Lett. 101, 118-105.. [22]. Skalak, R., Branemark, P.-I. (1969)., Deformation of red blood cells in capillaries.Science 164, 717–719.. [23]. Skalak, R., ¨Ozkaya, N., Skalak, T.C. (1989). Biofluid mechanics. Ann. Rev. Fluid Mech. 21, 167–204.. [24]. Suresh, S., (2006). Mechanical response of human red blood cells in health and disease: some structure–property–function relationships. J. Mater. Res. 21, 1871-1877.. [25]. Tsubota, K. et al , (2001), Direct measurement of erythrocyte deformability in diabetesmellitus with a transparent microchannel capillary model and high-speed video camera system. J. Microvasc. Res. 61, 231–239.. [26]. Tsubota, K. et al , (2005), A particle method for blood flow simulation application to flowing red blood cell and platelets. J. Earth Sim. 5, 2-7.. [27]. Tanaka, N., Takano, T. (2005).Microscopic-scale simulation of blood flow using SPH method. International Journal of Computational 72.
(87) Methods, 2(04), 555-568. [28]. Tsubota, K., Wada, S., Yamaguchi, T., (2006), Particle method for computer simulation of red blood cell motion in blood flow, Computer Methods and Programs in Biomedicine, 83, 139-146.. [29]. Tsubota, K., Wada, S., Yamaguchi, T., (2006), Simulation study on effects of hematocrit on blood flow properties using particle method, J. Biomechanical Science and Engineering, 6, 159-170.. [30]. Ye T., Li H., Lam K.Y.,(2010), Modeling and simulation of microfluid effects on deformation behavior of a red blood cell in a capillary. Microvasc Res. 80(3):453-63.. [31]. Yap, B., Kamm, R.D., (2005b) Mechanical deformation of neutrophils into narrow channels induces pseudopod projection and changes in biomechanical properties. J. Appl. Physiol. 98, 1930–1939.. [32]. Zhang, C., Zhang, Y. X., Wan, D. C., (2011). Comparative study of SPH and MPS methods for numerical simulations of dam breaking problems. J. Hydrodynamics Ser. 26, 736-746.. [33]. 韓旭、楊綱、強洪夫.(2005).光滑粒子流體動力學-一種無網格粒 子法.. 73.
(88)
數據
+7
相關文件
3: Calculated ratio of dynamic structure factor S(k, ω) to static structure factor S(k) for "-Ge at T = 1250K for several values of k, plotted as a function of ω, calculated
We compare the results of analytical and numerical studies of lattice 2D quantum gravity, where the internal quantum metric is described by random (dynamical)
All steps, except Step 3 below for computing the residual vector r (k) , of Iterative Refinement are performed in the t-digit arithmetic... of precision t.. OUTPUT approx. exceeded’
唇音 b巴 p趴 m媽 f花 舌尖音 d打 t它 n拿 l啦.. 舌葉音 z渣 c茶 s沙 j也 舌根音 g家
2 System modeling and problem formulation 8 3 Adaptive Minimum Variance Control of T-S Fuzzy Model 12 3.1 Stability of Stochastic T-S Fuzzy
void SetZTypeFunc(string FuncName, double zero, double one); //修正 ZType 歸屬函數 void SetSTypeFunc(string FuncName, double zero, double one); //修正 SType 歸屬函數
It is found that pressure increased gradually from fan inlet to fan outlet through the maximum flow rate, operation point, down to the cut-off point.. This implies the
(5)The Direction-Giving Language and the Empathetic Language of the principal have reach to the outstanding level of anticipa t i on f or t he t e a c he r ’ s j ob
