行政院國家科學委員會專題研究計畫 成果報告

結合磁共振造影與計算流力技術之人體主動脈的數值模擬

計畫類別： 個別型計畫

計畫編號： NSC93-2212-E-216-001-

執行期間： 93 年 08 月 01 日至 94 年 09 月 30 日 執行單位： 中華大學機械工程學系

計畫主持人： 牛仰堯

報告類型： 精簡報告

報告附件： 出席國際會議研究心得報告及發表論文

處理方式： 本計畫涉及專利或其他智慧財產權，2 年後可公開查詢

中 華 民 國 94 年 11 月 17 日

行政院國家科學委員會補助專題研究計畫成果報告

※ ※※※※※※※※※※※※※※※※※※※※※※※

※ ※

※ 結合磁共振造影與計算流力技術之人體主動脈的數值模擬 ※

※ ※

※ ※

※ ※※※※※※※※※※※※※※※※※※※※※※※※

計畫類別：█個別型計畫 □整合型計畫 計畫編號：NSC

執行期間：93 年 8 月 1 日至 94 年 9 月 30 日

計畫主持人：牛仰堯

成果報告包括以下應繳交之附件：

□赴國外出差或研習心得報告一份

□赴大陸地區出差或研習心得報告一份

█出席國際學術會議心得報告及發表之論文各一份

□國際合作研究計畫國外研究報告書一份

執行單位：中華大學

中 華 民 國 94 年 11 月 17 日

結合磁共振造影與計算流力技術之人體主動脈的數值模擬 計畫編號：NSC

執行期限：93 年 8 月 1 日至 93 年 9 月 30 日 主持人：牛仰堯 中華大學機械工程學系 協同研究人員: 曾文毅、彭旭霞, 臺大醫院心臟外科

鄭守成,李隆正 國家高速電腦中心 參與人員: 伍邦銓

在本研究中藉由磁振造影（Magnetic Resonance Imaging，簡稱 MRI 由台大醫院提供）掃 描出ㄧ位 21 歲健康年輕人之主動脈外型、主動脈一週期 19 個時間點的血液流速、流量，並且 利用套裝軟體 CFX5 （AEA Technology）建構成計算用外型作一計算流體力學之分析與探討，

以及預測容易發生病變的區域。 在本研究中模擬結果如下：主動脈血液為牛頓流體時其壓力、

壁面剪力、內流場與二次流變化、速度變化。其中發現在主動脈 0.05m 至 0.1m 處因為容易產生 迴流區、剪力變化甚大且又有二次流的產生，故判定為較易發生病變的區域。且在其中發現主 動脈內流場變化為：血液進入上升主動脈入口為散射狀，爾後順時針旋轉然後在主動脈弓及下 降主動脈前端有著二次流的現象，在下降主動脈下游逆時針旋轉。

關鍵字: 主動脈，牛頓流體，剪應力，分歧管、二次流、MRI

前言

由於現今科技的進步，有越來越多的研究 是對於人體的各項器官去作數值模擬或是實 驗研究，但是絕大部分的研究報告或是研究論 文的外型（MODEL）或是邊界條件都是過度『理 想化』過後的產品，對於實際醫學或是臨床上 的幫助有限。所以此次研究就是要將『真實』

的外型以及邊界條件放入數值模擬中，以利將 最真實的現象以及情況計算出來。又由於主動 脈疾病致死率高且發病至死亡的時間短暫，所 以此次研究就鎖定主動脈部分來做深入的探 討。

在台灣，最常見的主動脈疾病就屬動脈剝 離以及動脈窄化了，其原因出自主動脈的異常 變化。為了對臨床研究有所貢獻本研究所採用 的外型以及邊界條件全為真實之人體並且利 用 MRI（Magnetic Resonance Imaging 核磁共 振造影）所測量出來，且為了預測病變的發生

所以採用健康的人來量測，但又因為真實活體 難尋所以就只有以本人來親自量測，MRI 所量 測出來的外型包括：主動脈（Aorta）、頭臂 動脈（Brachiocephalic artery）、左共通頸 部動脈（Left common carotid artery ）、

左鎖骨下動脈（Left subclavian artery），

以及其各方向（U、V、W）速度及流量。

在本研究中，所使用的血液物理參數是依 照一般人的平均值所設定的，其值如下表所 示：

Density

Dynamic viscosity 1080 kg/m^3 0.0067 Pa s

表一 主動脈血液參數

模擬外型如圖一所示：

第十一屆全國計算流體力學學術研討會 台東，中華民國九十三年八月 The 11^th National Computational Fluid Dynamics Conference Tai-Tung, August,2004

【圖一 模擬主動脈外型及其各部位介紹】

其中，Ascending aorta 的直徑為 3.1 ㎝，

Descending aorta 的 直 徑 為 1.8 ㎝ ， Brachiocephalic artery 的直徑為 1.02 ㎝，

Left common carotid artery 的直徑為 0.67

㎝，Left subclavian artery 的直徑為 0.76

㎝。

數值模式

本 文 藉 由 商 用 計 算 流 力 軟 體 CFX5 來 求 解，其中採用三維的不可壓縮流，牛頓流體之 Navier-Stokes方程式，本文使用了CFX5之標

準k-

ε

Navier-Stokes方程式表示如下:

( )

•

∇

∂ +

∂ U

t ρ

ρ (1)

( )

( (

) )

T S

U U p U t U

U+∇• ⊗ =∇•− + ∇ +∇ +

∂

∂ρ ρ δ µ (2)

結果與討論

模擬外型如圖一所示，使用非結構性網 格，總格點數為 267876 個如圖五所示，主動 脈總長約 25 ㎝，壁面為剛體且假定血液為不 可壓縮流（incompressible flow）牛頓流體

（Newtonian flow）血液參數如表一，邊界條 件為利用 MRI 所掃出之速度 U、V、W 各分量如 圖六、圖七、圖八所示，雷諾數的巔峰值約為

4700 左右。從 MRI 掃描出的數據顯示，心臟跳 動一次週期為 0.855 秒，所以在程式模擬中的 一週期也同為 0.855 杪，為求精準所以總共模 擬五個週期，取第四至第五週期來進行討論。

又因為 MRI 掃描間隔為 0.045 秒，所以在結果 分析上是以 0.045 秒、0.090 秒、0.135 秒、

0.180 秒、0.225 秒、0.270 秒、0.315 秒、0.360 秒、0.405 秒、0.450 秒、0.495 秒、0.540 秒、

0.585 秒、0.630 秒、0.675 秒、0.720 秒、0.765 秒、0.810 秒、0.855 秒等 19 個時間點，以利 與 MRI 的數據進行比對，對於醫學上較有幫助。

在剪力部分，加速相時整體的剪力值升 高，而減速相所以剪力慢慢趨近在一定的範圍 內並無再明顯增加。值得一提的是在主動脈 0.05 m 至 0.1 m 處剪力曲線突然有巨大的起 伏，其原因是因為在此處幾何外型改變甚大，

且開始有歧管效應的作用產生使得速度變化 劇烈，在加上此處流體撞擊內壁或外壁，使得 外壁或內壁的流體量較內壁或外壁少而且流 動緩慢，而形成剪力巨大起伏，也造成迴流區 的產生。而在此區域因為剪力巨大起伏且加速 相與減速相的剪力顛峰值差距甚大，根據以往 文獻不管是高剪力區或是低剪力區或是剪力 變化較大的區域都是較易發生血管病變的地 方，所以在此區域都是在臨床上值得注意的方 向。而在流體通過此區域後，隨著流體的流 動，速度變化趨漸穩定，而剪力就不太有變化 了。

在內流場與二次流部份，其中 A、B、C 點 在上升主動脈（ascending aorta），D 點在主 動 脈 弓 頂 點 ， E 、 F 、 G 、 H 點 在 下 降 動 脈

（decending aorta）。在 1/19 T 時流體開始 加速，在上升主動脈流場變化並不明顯只有一 點像順時針方向（clockwise）旋轉的行為，

等到流體到了主動脈弓時因為流體尚在加速 流 場 並 不 穩 定 只 有 在 E 、 F 有 類 似 二 次 流

（secondary flow）出現，而流體到了下降動 脈由於幾何外型較無變化且速度穩定又由於 完全發展流的效應所以此處流場穩定呈逆時 針方向（counter-clockwise）旋轉。而到了 9/19 T 開始有了明顯的變化，流體在進入上升 主動脈時如 A 點為散射狀的進入，經 B 點及 C 點時為順時針旋轉，到了主動脈弓時開始產生 強烈的二次流（secondary flow），直至下降 動脈前端 E 點及 F 點都有強烈的二次流產生，

而流體到了下降動脈下游 G 點及 H 點為逆時針 方向旋轉。而在 10/19 T、11/19 T、12/19 T 中因為狀態與 9/19 T 差不多所以物理變化大 致相同。而在 13/19 T、14/19 T、15/19 T、

16/19 T、17/19 T、18/19 T、19/19 T 等時間 點由於速度趨緩所以二次流的現象變弱，特別 是在 D 點以及 F 點之二次流最為明顯，但這些 時間的大多還是跟之前的時間點有著相同的 趨勢，流體進入上升主動脈入口為散射狀，爾 後順時針旋轉然後在主動脈弓及下降主動脈 前端有著二次流的現象，在下降主動脈下游逆 時針旋轉。二次流形成與速度、壓力與幾何外 型的轉變有著極大的關係，其主要的原因為高 低壓力差、速度差所構成的，也就是壓力梯度 與速度梯度在一特定的幾何外型下所形成出 來的特殊物理現象。二次流亦會對血管表面組 織造成損害，長期產生二次流的區域亦是容易 發生病變的所在，雖然不見得每一個人都會患 有心血管的病變，但只要心臟所給予的流量或 速度差異甚大時，或是血管壁上有破損時，或 是沉積在血管的雜質過多時，容易形成二次流 的區域都是會發生病變的所在。更值得一提的 是本研究發現二次流的時間點、發生位置與渦 流的運行方向，跟在 1993 年由 Kilner【13】

利用 MRI 所做出的實驗結果幾乎一樣，更確定 本研究的準確信與可信度。

在速度分布方面，其目的是觀察血液的流

動方向與速度的變化，在 1/19 T 時流體尚在 加速所以速度並未到達最高在 A、B、C、D 點 並無較特殊的行為發生，但在 E、F、G 點中因 為血液剛通過主動脈弓流體轉向而且有歧管 的效應存在，所以血液在外壁分佈的較多而內 壁分佈的較少。而在 2/19 T 時因為流量極速 度達到最高所以在主動脈弓 D 點的速度最大、

且在 E、F 點可以發現流體在內壁分部較多而 在外壁分部較少。在 8/19 T 時速度已趨近於 最低而且在 A、B、C、D 點有逆流的情形發生 而在 9/19 T、10/19 T、11/19 T 時除了在 A、

B、C 點有逆流的情形發生外在 E、F、G、H 點 由於幾何外型的急遽變化且為了使質通量守 恆而有迴流的情形發生。而 12/19 T、13/19 T、

14/19 T、15/19 T、16/19 T、17/19 T、18/19 T、19/19 T 等，由於這些時間點的流量都很小 甚至趨近於零（但大於零）所以速度都無太大 的變化。

結論

在剪力以及二次流的結果發現在主動脈 長度為0.05m至0.1m的區域，也就是主動脈弓 以及主動脈弓之後的部份區域，壁面剪力變化 落差很大，而且在此區域有二次流的形成，吾 人判斷此區域極有可能是容易發生主動脈病 變的所在。主動脈中內部血液的流動行為如 下；血液進入上升主動脈入口為散射狀，爾後 順時針旋轉然後在主動脈弓及下降主動脈前 端有著二次流的現象，在下降主動脈下游逆時 針旋轉。

參考文獻

13. P. J. Kilner, Z. Y. Yang, R. H. Mohiaddin, D.

N. Firmin, D. B. Longmore.,1993. Helical and retrograde secondary flow patterns in the aortic arch studied by three-dimensional magnetic

第十一屆全國計算流體力學學術研討會 台東，中華民國九十三年八月 The 11^th National Computational Fluid Dynamics Conference Tai-Tung, August,2004

resonance velocity mapping. Circulation 88, 2235-2247.

【圖二 MRI掃出之主動脈外型】

【圖三 簡化過後之MRI主動脈外型】

【圖四 於CFX5重建之主動脈外型】

【圖五 模擬網格】

【圖六 入口流量圖 】

【圖七 入口速度 U 分量波形圖 】

【圖八 入口速度 V 分量波形圖 】

【圖九 入口速度 W 分量波形圖 】

【圖十 加速相剪力分布於 1/19 T】

【圖十一 加速相剪力分布於 2/19 T】

【圖十二 減速相剪力分布於 8/19 T】

【圖十三 減速相剪力分布於 9/19 T】

【圖十四 減速相剪力分布於 10/19 T】

第十一屆全國計算流體力學學術研討會 台東，中華民國九十三年八月 The 11^th National Computational Fluid Dynamics Conference Tai-Tung, August,2004

【圖十五 內流場截面位置示意圖】

【圖十六 入口內流場於 1/19 T】

【圖十七 加速相內流場於 2/19 T】

【圖十八 減速相內流場於 8/19 T】

【圖十九 減速相內流場於 9/19 T】

【圖二十 減速相內流場於 10/19 T】

【圖二十一 入口速度立體圖於 1/19 T】

【圖二十二 加速相速度立體圖於 2/19 T】

【圖二十三 減速相速度立體圖於 8/19 T】

【圖二十四 減速相速度立體圖於 9/19 T】

【圖二十五 減速相速度立體圖 10/19 T】

ISTP-16, 2005, PRAGUE 16^TH INTERNATIONAL SYMPOSIUM ON TRANSPORT PHENOMENA

The Sixteenth International Symposium on Transport Phenomena (ISTP-16)學術研討會心得報告

中華大學機械系 報告人: 牛仰堯

前言：

參加國際性會議不僅可以了解國外在相關領域的最新研究發展狀況， 同時也是與國外 學者專家技術交流最直接最具經濟效益的方法之一。 因為一個成功的研討會，可以在 短短的幾天中，聚集了來自各地達數百，甚至數千位學者專家於同一地點，做面對面直 接的討論與溝通。 因此為求得在最短時間內獲取最新的資訊，參加謹慎選擇的國際性 研討會，是最值得的。 本次會議 ISTP-16(The Sixteenth International Symposium on Transport Phenomen) 於 8 月 29 日至 9 月 1 日舉行為 2 年一度的且有十餘年歷史之 熱流工程研討會，有來自各地的相關專家與學者，包括大學教授、研究人員、各國家實 驗室的專職人員、工業界相關研究人員出席或發表論文。在會場上並有出版商提供各種 熱流工程相關專業書刊與資料，以及最新 Flow 電腦模擬所需的軟體，實驗測量,儀器所 需的各項器材之展覽。

過程

本次會議的地點在 Prague, Czech Republic 會議於 8 月 28 日開始辦理報到註冊，主 辦單位同時展示一些會議相關論文集及專業書刊，並有部份工作人員舉行會議前之準備 集會。接下來的四天為正式的會議，包括專題演講、8 個邀請演講、專題討論及論文發表 等議程。 由大會依專業細分為 10 大項，論文發表以 3 個場地同時進行，共計有 24 場，

每場約有 5-6 篇論文，會場外的大廳有軟體廠商、出版商與實驗器材廠商的展示及研究生 之論文成果展。大會的子題包括了：

1. Turbulence and Flow Instabilities 2. Boundary Layer and Free Shear Flows

3. Experimental and Computational Fluid Dynamics - Gases 4. Experimental and Computational Fluid Dynamics - Liquids 5. Heat and Mass Transfer

6. Rheologically Complex Systems and Biofluid Dynamics 7. Two Phase Flow

8. Chemical Process Systems 9. Combustion and Reacting Flows

10. Industrial Aerodynamics and Wind Engineering

11. Heat Exchangers 12. Environmental Systems 13. Electronic Equipment Cooling 14. Fluid Dynamics in Micro-systems 15. Special Topics

由於場次眾多，個人僅能選擇性的參與 cavitation& Liquid-Gas flows 論文發表。大本次會 議台灣除筆者外，另外有臺灣、清華、交通、成功、逢甲大學、虎尾科大學教授學生近 五十多人參與論文發表可謂盛況空前。

心得

. 由本次研討會所發表的 200 餘篇論文及張貼海報可以看出，Multi-Phase flow 領域本 世紀發展之重要方向，有相當數量牽涉與奈微流體相關。例如 DNA 或蛋白質分子於奈微 流體之動態特性，往往影響此類生醫分子之輸送及流體物理，為困難度極高但卻急需輸 送解決之問題。奈微尺度生醫流體的混合（mixing）與分離（separation）再微生醫機械 系統內，有時，生醫流體需要配合，例如藥劑的混合。有時，生醫檢體需要分離，以便 進 行 進一步的分析檢驗。但應 奈 微 尺度 的 流 動都 將變為層 流 ，且是低雷諾數 （low Reynolds number）、Multi-Phase flow 的流動，此時混合便得非常困難。分離也是困難且 耗時的，極需新技術的研發。從 1950 年到 1960 年代，隨著半導體工業的蓬勃發展，製 程技術的提昇，在為電子的世界裡產生了微型化及自動化的革命。電腦的尺寸由房間般 的大小，一直縮到膝上型電腦。在最近一些切割技術已被需求，應用到醫藥診斷實驗 上，不過還是許多實驗室的工作依然是效率低，辛苦耗時的。 然而，在過去五年間，

研究及發展臨床診斷的系統於 lab-on-a-chip 的技術已相當驚人。大約有 20 篇的文章主題 是研究 lab-on-a-chip， Multi-phase flow 是一個發展快速的研究領域，無論在日常生活 上，或是在工業上應用範圍都相當廣泛。本次會議除了論文發表外，更有許多來自世界 各地從事此一領域的學者專家齊聚一堂交換研究心得，探討各種 multi-phase flow 問題，

不僅使參加者有機會發表自己的研究心得，更可了解他人的研究現況。藉著論文資料的 蒐集，未來尚可提供國內相關研究人員會議狀況，對於促進國內 Multi-phase 科技研究，

提升學術水準有莫大的助益。

ISTP-16, 2005, PRAGUE 16^TH INTERNATIONAL SYMPOSIUM ON TRANSPORT PHENOMENA

NUMERICAL SIMULATION OF TRANSIENT SLURRY-CAVITATED MULTIPHASE FLOWS

Yang-Yao Niu and Yee-Ming Lin

Institute of Mechanical and Aerospace Engineering, Chung-Hua University Hsin-Chu, Taiwan, ROC

Corresponding author: Yang-Yao Niu (yniu@chu.edu.tw)

Keywords: slurry-cavitation, water-vapor flow, equilibrium saturation model, upwind methods

Abstract

In this work, a multi-scale modeling of cavitation-erosion phenomena performed in an Eulerian–Lagrangian two-way coupling approach is developed. Cavitation-erosion phenomena are widely seen in multi-phase flow with different sizes of time and length scale. Based on the multi-scale concepts, the information obtained at the level of small length scales can be used to provide closure information at the level of larger length scales.

First of all in our implementation, a two-phase mixture flow model and a temperature dependent homogeneous equilibrium saturation models are used simulate phase transitions and steam-water interfaces.

Subsequently, solid particle trajectories caused by the translation and rotation and particle-wall collisions are evaluated by a Lagrangian approach based on the continuum flow information. Then, an Eulerian–

Lagrangian two way coupling are used to evaluate the interactions between the gas-liquid mixture and dispersive solid particles such as the drag, virtual mass and lift forces. An erosion model is also implemented in this Eulerian-Lagrangian frame work.

1 Introduction n

I the power industry, surface material-loss process caused by the erosion and cavitation substantially shortens the lives of pipelines, heat exchanger systems, and turbomachinery surfaces. In this work, the erosion we are concerned is the so-called slurry erosion which occurs at a surface impinged upon by solid particles suspended in a liquid stream. The progressive loss of original material from a solid surface due to mechanical interaction between that solid surfaces and impinging liquid or solid particles in a fluid or, a multi-component fluid. The corresponding particle-laden flow motions with erosive wear in industrial applications [1-4] can be categorized as (I) impingement flows with application to metal cutting, boundary layer flows occurring in turbomachinery (III) ;obstructed flows passing over heat exchanger confined flows applied to pipeline systems. In our previous work [4, 5], numerical prediction of erosion by particle impact has been studied the confined flows in pipeline systems. Computational results reveal that turbulence intensity, particle size and inlet liquid or gas flow velocity have significant influences on erosion. However, our work only considered the solid erosive wear influenced by the pure gas flow. The influence of cavitated mixture flows on erosion has not been done yet. To our knowledge, the combined effects of cavitation-erosion are not widely seen in the literatures. As noted in [7,8], cavitation is a physical phenomenon associated with three aspects: formation, growth and collapse of bubbles within the body of a liquid due to the process of nucleation in a liquid flow where the pressure falls below the vapor pressure. The formation and collapse of vapor bubbles or cavities in a fluid can produce extremely high pressures transmitted from the collapsing bubbles to the surface either in the form of a shock wave or by microjets, depending on the distance from the surface. They collapse with enough force to erode microscopic pieces of metal if they are close to the tube wall. The cycle of formation and collapse of the bubbles occurs at a high frequency and dynamic stress can cause the material loss by fatigue. If considering the erosion-cavitated bubbly flow, solid particle with bubbly type turbulent liquid-gas regimes, where all types of interfaces exist, not only transport equations and closure models need to be developed for each phase, but also each phase should accommodate the existence of that phase indifferent forms and shapes. For example, cavitated bubble flows can be spherical, distorted, confined and/or elongated. In this regard, several flow or phase phenomena co-exist at the different sizes. This type of representation would encompass the various flow regimes which are required to deal with the multi-scale approach based on their transport equations and closure conditions.

The Eulerian-Lagrangian approach could be one of the state-of the-art and major step towards providing the possibility of the multi-scale numerical simulation.

In the Eulerian-Lagrangian framework, the phenomena at a certain length scale are allowed to be analyzed by different approaches. The idea of this study is to use different levels of modeling, each developed to study phenomena at a certain length scale. Information obtained at the level of small length scales can be used to provide closure information at the level of larger length scales. In this work, the mixture model with an equilibrium phase change model is used to simulate the interfaces of liquid-gas and gas bubble morphology in a cavitated flow field. The liquid-gas model can be used in Eulerian approach such as a two-fluid model or mixture model to simulate the larger scale two-phase flow phenomena; also can accurately capture smaller scale sizes as the interfaces and the formation of gas bubble simultaneously.

Regarding the small scale forces as the drag, virtual mass and lift forces of dispersive phase can be simulated based on the Lagrangian model; the Eulerian-Lagrangian method can consider the coupling between

ISTP-16, 2005, PRAGUE 16^TH INTERNATIONAL SYMPOSIUM ON TRANSPORT PHENOMENA

dispersive phase and continuum phase concurrently. Therefore, in this study, the multi-scale model for slurry solid particle loaded in liquid-vapor cavitated flow problems are developed based on our existing AUSM solvers with the mixture model [5, 6] and the Euerian-Lagragian approach [2-4]. In our previous works [5,6], The AUSMDV with a two-phase mixture model was presented its simplicity to maintain pressure equilibrium over contact discontinuities on the resolution of interfaces, contact surfaces and shock waves under a high-density ratio up to 1:1000 like air to water. To continue our efforts, an Eulerian type mixture model with a phase transition model is proposed here to simulate steam-water mixture flow problems. In order to enhance convergence for low-speed flow regime, the mixture model incorporating a preconditioning algorithm [9, 10] is used to enhance the convergence and accuracy for solving the low-speed fluid flow phase. For the particulate flow phase, the fourth-order Runge-Kutta method is adopted to solve the Lagrangian trajectory equations. Also, a renormalization group theory (RNG) based k-ε turbulence model is used to evaluate the turbulence effects for liquid-vapor mixture flows. Subsequently, solid particle trajectories caused by the translation and rotation and particle-wall collisions are evaluated by a Lagrangian approach based on the continuum flow information. Test cases in this work include 1-D cavitation and 2-D transient slurry-cavitated U tube flows.

2 GOVERNING EQUATIONS

To continue the previous efforts on solving the mixture model [6] a two-dimensional two phase Navier- Stokes equations with its mixture model with a saturation model based on is used.

2.1 Saturation model

In the simulation of the phase change process among the steam-water mixture, three equations of state for water, vapor and saturation states modified from [11-12] based on local phase equilibrium are selected together to yield a smooth link at the so-called spinodal points, the extreme points on the both sides of the saturation zone in the temperature-pressure diagram, also accurately cover the range (273.15K≤T≤623.15K

& ). They are briefly described as the following section and the details can be found in [11-12].

( )

ps ≤ ≤100

Liquid Phase EOS

An established EOS for modeling liquid water as a compressible fluid is given as

Specific volume

( )

πγ

v , (1) Specific internal energy ( )

π

τ πγ

τγ −

RT = p T

u , (2)

Specific enthalpy

( )

τγ

RT = p T h ,

(3) Speed of sound

( )

( )

2 2

2 2

, a T p

RT

π

π πτ

ππ ττ

γ γ τγ τ γ γ

= −

−

(4)

Inverse reduced temperature

T 1386K

τ

Gas phase EOS

The corresponding water vapor EOS is

Specific volume

( ) (

)

RT p p T

v , =πγ_π⁰−γ_π (5) Specific internal energy

( ) (

) (

)

RT p T u

π π π

π γ π γ γ

γ

τ + − +

= ^{0 0}

, (6) Specific enthalpy

( ) (

R

)

T p T h

τ

τ γ

γ

τ −

= ⁰

, (7) Speed of sound

( )

( )

( ) ( ) ( )

2 2 2

2 2

2 0

, 1 2 1 1

r r

r r

r

r

a T p RT

π π

π πτ

ππ

ττ ττ

πγ π γ πγ τπγ

π γ τ γ γ

+ +

= + −

− +

+

(8)

Derivatives of the residual part are also shown in [12].

Saturation EOS

The saturation region are assumed to be in equilibrium: _{Pl = Pg = P}and , where the pressure equals the saturation pressure

T T T_l = _g = )

(T P

P = _sat , and the temperature is equal to the saturation temperature:T =T_sat . The saturation pressure and density can be obtained respectively as

4 5 . 0

2 4 )

( 2 1

)

( ⎥

⎦

⎢ ⎤

⎣

⎡

− +

= −

AC B

B

C MPa

T

P_sat (9)

where

2 1

2 n n

A

= ϑ + ϑ +

ϑ

ϑ

8 7 2

6 n n

n

C=

ϑ

ϑ

and the liquid and gas densities along the saturation curve are also chosen from [11] like

( )

3 6 43 3 5

16 3 2 1

1 1θ θ θ θ

ρ

ρ T b b b b

c

tsat = + + + + (10)

( ) ³⁷₆

6 5 18 6 4 8 6 3 4 6 2 2

1θ θ θ θ θ

ρ

ρ T c c c c c

n

c

gsat ⎥= + + + +

⎦

⎢ ⎤

⎣

A ⎡ ⁶

71 6θ

+c (11)

where

10 9

1

1 n

K T

n K T

− +

ϑ= and the remained coefficients are seen in [12].

ISTP-16, 2005, PRAGUE 16^TH INTERNATIONAL SYMPOSIUM ON TRANSPORT PHENOMENA

Based on the above equation of state, the mixture flow sound speed and the volume fraction of gas phase diagram for temperature=300k and 400k is depicted in Figure 1.

Figure 1 Sound speed--gas volume fraction diagram

3 Results and Discussion

1-D cavitation tube

The first case [11] considers a tube filled with a constant specific-heat perfect gas, initially at uniform pressure and density. Initially, the fluid occupying the right-half of the tube is set into motion to the right, and the fluid on the left-half of the tube is set into motion in the opposite direction. The initial conditions chosen for this example are temperature 300 K,density 1 kg/m3, velocity in left-hand tube is 3000 m/s, whereas the fluid in right-hand tube is -3000 m/s. The solution to this test case corresponds to two rarefaction waves running in opposite directions, resulting in vaporization of the liquid in the central region.

(A)

(B)

(C)

Figure 2 One-dimensional cavitation tube plots ((A)-gas volume fraction, (B)-density, (C)- pressure )

Results of various properties are shown at 10,30,50, and 70 µs in Fig 2. Fig 2(A) shows the velocity profiles. The fluid velocity decreases across the rarefaction waves reaching zero velocity at the center of the domain. Fig 2(B) obtains strong rarefaction waves (right and left running) induce strong decrease in density from the initial value of 1000 kg/m3 to pure vapor values~3 orders of magnitude less dense. Vapor is being produced at a constant temperature because of the assumption of local equilibrium during phase change.

When the fluid becomes a pure gas, the rarefaction effects continue. The pressure decreases abruptly in the pure liquid, then reaches the saturation pressure and remains constant. It the decreases again in the pure gas phase near the center of the domain. For the initial conditions studied, Fig. 2 (C) rarefaction waves are sufficiently strong to evaporate all of the liquid. Across these waves, the gas volume fraction starts at α = 0 and reaches α = 1 progressively.

.

2-D Slurry-Cavitation

In this study, numerical evaluation of slurry-cavitation in a two-way coupling system is performed for the turbulent water-vapor flow through a U type tube. First of all, computational domain and grid mesh distribution used for a U tube bend is 150x59 grids. The flow conditions used in [2, 4] are assumed as a dilute particle-laden gas flow. They are the bulk velocity is temperature =375k, water density=1000 ㎏

/m³; U=52.19 m/s corresponding to the Reynolds number 3.47×10⁵; the inlet turbulence intensity of the water phase is 1﹪; the particle diameter size 25, 50 and 75 µm; the particle material density ρ_s=2,990 ㎏

ISTP-16, 2005, PRAGUE 16^TH INTERNATIONAL SYMPOSIUM ON TRANSPORT PHENOMENA

/m³; the inlet particulate bulk density ρ_p,in=1.8×10^-4 ㎏/m³ ; The corresponding particle loading and volumetric ratios are 1.5×10^-4 and 6×10^-8.

t=260 µS

t=320 _µS

t=350 _µS

t=600 µS

Figure 7 The gas volume fraction evolution of cavitated vapor-water flow in an U tube( T=370 K、t=810

µS )

Figure 8 the particle traces in the cavitated fluid in an U tube (T=370K, particle size=50

µ

ISTP-16, 2005, PRAGUE 16^TH INTERNATIONAL SYMPOSIUM ON TRANSPORT PHENOMENA

Figure 9 the particle traces in the cavitated fluid in an U tube (T=370K, particle size=50

µ

Figure 10 the particle traces in the cavitated fluid in an U tube (T=370K, particle size=75

µ

Inlet flow condition is assumed subsonic flow condition (fixed velocity and temperature and extrapolated pressure) and subsonic outflow (fixed pressure and extrapolated velocity and temperature) are used in the calculations. The RNG k-εturbulence is used to evaluate the turbulence effects. First, figure 7 illustrates that cavitated bubbles first emerging from the turning corner, then getting growth and moving into the flow separation flow region near the upper wall after the turning section due to an adverse pressure gradient, then convects into the downstream. After the computation of cavitated water-vapor flow, the transient solutions from 0-810 µs are used for the initial conditions for the calculation of particle- water-vapor mixture in a two-way coupling way. First, to decide the inlet particle number used in solving the particle trajectory equation, the particle-wall frequency, defined in [4] as the ratio of the particle number colliding with the wall to the particle number initially released at the inlet, is used to study the effect of particle size and number on the particle impact distribution. 20 particle trajectories are chosen to compute per each of 50 starting locations uniformly distributed over the entrance of the bend. I It is found that a secondary flow region near the inner wall as previous work [4] was shown to be a particle-free and erosion-free region. Here the particle trajectories of different particle sizes loaded in cavitated mixture flows at different time are demonstrated in Figure 8-10. Different particle size in dm＝,25, 50 and 75 µm are selected for discussing the influence of particle size. It is shown that particles larger than 50 µm collide with the outer wall. It is also noted that the effect of particle collision is not significant near the inner wall. In Figure 9, it is seen that the trajectories of the particles with 50 µm deviate considerably from the cavitated bubble regions when the time is 810 µs. However, the traces of the particles with 50 µm influenced by cavitation effects are not observed during the period early than 810 µs such as depicted in Figure 8. As shown in Figure 10, the particle size larger than 50 µm due to larger inertial momentum, the particle traces are shown not be changed by cavitation effect in the simulation. Regarding the study on the influence of cavitation on the particle motion by turbulence and flow conditions will be required further studies in the near future.

Concluded Remarks

In this work, an Eulerian-Lagrangian approach based on the AUSMDV scheme with preconditioning technique and a modified steam-water phase transition model is developed for the multi-scale modeling of the slurry-cavitating flow simulation. The cavitation test cases are verified on 1D cavitation tubes.. In addition, the cavitation effects from water-vapor mixture flow on the solid particle motion have been simulated for a 2D U tube flow. A further work is required to study physical phenomena contain many sizes of length and time scales such as the turbulence effects and closures among different phases. Also a more sophisticated multi-scale modeling relied on the current CFD techniques will be highly required.

Acknowledgments

We acknowledge financial support from the National Science Council under NSC 93-2623-7-216-001-NU.

Also we thanks the computer facility supported by the National Center for High-Performance Computing, Hsin-Chu, Taiwan.

2. References

1. C. Humphrey, “Fundamentals of Fluid Motion in Erosion by Solid Particle Impact”, Intl. J. Heat and Fluid Flow, 11, No. 3,pp. 170-194 (1990)

2. J. Y. Tu and C. A. J. Fletcher, “Numerical Computation of Turbulent Gas-Solid Particle Flow in 90°

Bend”, AIChE J., 41, 2187~2196, (1994).

3. S. Naik and I. G. Bryden, “Prediction of Turbulent Gas-Solid Flow in Curved Ducts Using the Eulerian-Lagrangian Method”, Int. J. for Numer. Met. Fluids, 31,579-600 (1999).

4. Yang-Yao Niu and Jia-Chih Tsao, “Numerical Evalution of Erosion in Curved Ducts”, Numerical Heat Transfer, Vol. 41:341-356, March, 2002.

5. Yang-Yao Niu, “On Some Simple and Robust Riemann Solvers for Compressible Two-Phase Flows”, AIAA CFD Conference, AIAA paper, 2001–2647, June 11-June 13, 2001

6. Yang-Yao Niu, “Simple Conservative Flux Splitting for Multi-component Flow Calculations”, Numerical Heat Transfer, Part B, Vol.38 No. 2, pp.203-222, September, 2000

7. Jiangong Qin,”Numerical Simulation of cavitating Flows by the Space-Time Conservation Element and Solution Element Method”, PhD Thesis, Wayne State University, Michigan, 2000.

8. David P. Schmidt, Christopher J. Rutland, and M. L. Corradini, “A Fully Compressible Two- Dimensional Model of Small High-speed Cavitating Nozzle”, Atomization and Sprays, Vol.9, pp255- 275,1999

9. Jack R. Edwards, Randall K Franklin and Meng-Sing Liou, “Low-Diffusion Flux Splitting Methods for Real Fluid Flows at All Speeds, AIAA Journal, pp1624–1633, 38, 2000.

10. Merkle, C. L., Sullivan, J. Y., Buelow, P. E. O., and Venkateswaran, S.,“Computation of Flows with Arbitrary Equations of State,” AIAA Journal, Vol. 36, No. 4, 1998, pp. 515–521,1998

11. Richard Saurel and Jean Pierre Cocchi,”Numerical Study of Cavitation in the Wake of a Hypervelocity Underwater projectile, Journal of Propulsion and Power, Vol. 15., No.4, July-August,1999