• 沒有找到結果。

脈衝流於非等截徑彈性管(動脈血管)中之流場數值模擬

N/A
N/A
Protected

Academic year: 2021

Share "脈衝流於非等截徑彈性管(動脈血管)中之流場數值模擬"

Copied!
7
0
0

加載中.... (立即查看全文)

全文

(1)

脈衝流於非等截徑彈性管(動脈血管)中之流場數值模擬

Numerical Simulation of Pulsatile Flow in Elastic Tubes (Arteries) with Nonuniform

Cross-Sections

計劃編號:NSC 89-2212-E-002-144 執行期限:89 年 8 月1 日至 90 年 7 月 31日 主持人:伍次寅 國立台灣大學機械工程學系 共同主持人:王淑音 中國文化大學畜產學系 計劃參與人員:馮建忠、蘇致遠 國立台灣大學機械工程學系 一、中文摘要 本計劃運用數值法來計算脈衝流(pulsatile flow)於彈性管(elastic tube)內之流場變化。將 血液視為一不可壓縮之 Newtonian 流體,描述 流場部份之統御方程式為不可壓縮形式之非 定常(unsteady) Navier-Stokes 方程組。在薄管 壁之假設下,運用彈性力學中應力及應變之 力平衡式可推導出管壁之運動方程式。管內 流 體 與 管 壁 運 動 之 交 互 作 用 (fluid-structure interaction)乃經由流體施加於管壁上之壓力以 及 管 壁 位 移 速 度 作 動 態 之 耦 合 (dynamical coupling)。在數值方法上,我們採用雙時步 (dual time-step) 虛 擬 壓 縮 策 略 (pseudo-compressibility)配合有限體 積法(finite-volume formulation)來求得 Navier-Stokes 方程式之數 值解,而後將所求得之流場於管壁面之壓力 值代入管壁運動方程式以計算出管壁之位移 和速度,作為下一時刻計算流場速度與壓力 分布用。本計劃中我們將探討管截面變化及 管徑收縮對管內流場之影響,藉以瞭解血管 因狹窄(stenosis)等因素所造成之潛在危機。本 計劃之初期目標為培養及累積計算與結構動 力耦合之內流場之能力,並結合彈性力學和 血液動力學(Hemodynamics)方面的知識,作為 日後發展建立一套計算複雜血管內血液流動 情形之數值程式的基礎。 (關鍵詞:脈衝流、彈性管、血管流、數值計 算) 英文摘要

In this project, we apply a numerical scheme to calculate the motion of a pulsatile flow in an

elastic tube (blood vessel) with variable cross-section. The blood is assumed to be an incompressible Newtonian fluid, and the flow is governed by a set of unsteady incompressible Navier-Stokes equations. Under the assumption of a thin-wall tube, the transverse motion of the tube is obtained from a force-balance equation using an elastic stress-strain relation. Dynamic coupling of the fluid flow and the tube-wall motion through the pressure at the tube wall is considered in the present formulation. In the fluid part, the time-dependent incompressible Navier-Stokes equations are first written into a conservative form and then solved by using a dual time-step plus pseudo-compressibility strategy. A finite-volume formulation is adopted, and an upwind difference combined with the TVD scheme is applied to evaluate the flux terms. Calculated flow pressure at tube wall is then fed into the equation governing the tube wall motion. The obtained tube-wall displacement and transverse velocity are then used in advancing the internal flow field to the next time level. In this project, the effect of nonuniform cross section and constriction in tube diameter on the general features of the pulsatile flow will be studied. The short-term aim of this study is to build up the knowledge and capability of calculating fluid-structure interaction problems, and to set firm foundation for treating more complicated blood flow dynamics in the future.

(keywords: pulsatile flow, elastic tube, blood vessel, numerical calculation)

(2)

由於電腦硬軟體在過去十幾年來快速的 發展大幅度的提昇了其計算的能力,計算流 力(computational fluid dynamics)之方法和技術 已被廣泛應用於許多領域中;運用數值計算 來模擬血管內血液之流動即為一例。近年來 之文獻中指出一般常見之心血疾病譬如狹心 症 (cardiac stenosis) 、 動 脈 硬 化 (atherosclerosis)、血栓(thrombosis)、主動脈塌 陷(artery collapse)等與血管彈性和血管內流場 速度及壓力之分布有著密切的關連,因而更 突顯了以計算流力來模擬研究血管內流場之 重要性與迫切性。 鑒 於 完 整 心 血 循 環 系 統 (cardiovascular circulation system)的複雜性,一般在計算模擬 血管流時均作了一些簡化問題的假設而專注 於某特定現象之探討。例如最近的一些文獻 裡有將血管視為一剛體而針對分歧血管之流 場所作之計算([1-4])。另有文獻探討單一彈性 管管壁變形對管內流場及流體作用於管壁剪 應力之影響([5-8]),然而此類計算所使用之流 體與管壁交互作用模式為一靜態之耦合(static coupling),且多半假設流場為一不隨時間變化 之 穩 態 流 (steady flow) 。 亦 有 文 獻 探 討 Non-Newtonian 流 體 於 剛 性 擴 張 管 中 之 流 動 情 形 [9]。文獻[10-11]則是模擬具滲透性管壁之管 流,然其亦假設管壁為一剛體。文獻[12-13] 藉由流場計算來探討血管因狹窄(stenosis)而形 成局部低壓所引發管壁塌陷(collapse)症狀之問 題,同樣的在這些研究裡亦僅考慮穩態流場 及靜態之流體與管壁應變耦合模式。 在初步的規劃裡,本研究亦針對所欲模 擬之血管流作如下之假設: (1) 假設血液為一單相不可壓縮 Newtonian 流 體; (2) 僅考慮單一管路,管壁具彈性但不具滲透 性; (3) 假設流場為軸對稱(axisymmetric)形式; (4) 假 設 薄 管 壁 , 忽 略 管 壁 沿 徑 向 之 應 力 變 化。 為 模 擬 因 心 臟 壓 縮 舒 張 所 產 生 之 脈 衝 血 流 (pulsatile flow),流場採用非定場(unsteady)之 Navier-Stokes 運動方程式。管壁運動與流場間 之交互作用亦採取較精確之動態耦合模式。 結合後之流場及管壁運動方程式則利用數值 計算法來求得其解。 本計劃中將針對漸縮、頸縮等不同管截 面變化對脈衝流振盪所造成之影響來探討。 同時,管壁應力與應變物理性質之改變效應 亦可經由流場計算而得知。各種不同狀況下 所計算出之流場分布將有助於醫學上對一些 常見心血疾病從血液流動學的觀點提出合理 之現象解釋和研究出對策。 三、研究方法 1.數學模式 流場部份 流場之統御方程為一組軸對稱座標(z,r)上 之非定場不可壓縮 Navier-Stokes 方程式(參見 圖一) 0 ) ( ) ( = ∂ ∂ + ∂ ∂ J rV J rU η ξ (1) s f f e e t qˆ (ˆ ˆv) (ˆ ˆv)= ˆ ∂ − ∂ + ∂ − ∂ + ∂ ∂ η ξ (2a,b) 其中     = rvru J qˆ 1 ,         − = r v Re p J s 2 0 1 ˆ ,       + + + + = t r t z rv rp rvU ru rp ruU J eˆ 1 ξξ ξξ ,      + + + + = t r t z rv rp rvV ru rp ruV J fˆ 1 ηη ηη     + + + + = r r z z r r z r z z v r u v rv v u r ru JRe e ξ ξ ξξ 2 ) ( ) ( 2 1 ˆ ,     + + + + = r r z z r r z r z z v r u v rv v u r ru JRe f η η ηη 2 ) ( ) ( 2 1 ˆ , Re:雷諾數, r z v u U = ξ + ξ , V=z+r, J =zξrηzηrξξ =ξ(z,r,t),η=η(z,r,t)為軸對稱座標(z,r,t)與 計算座標(ξ,η,t)間之轉換式。由於流場之壓力 項並未出現於不可壓縮 Navier-Stokes 方程組 中 之 連 續 方 程 式 內 , 在 數 值 處 理 上 較 為 困 難 ; 因 此 引 介 一 虛 擬 壓 縮 因 子 (pseudo-compressibility factor)β將流場之壓力與連續方 程式結合而形成一虛擬可壓縮系統,流場之 運動方程式則可寫成傳統可壓縮流場之保守 形式(conservative form): S F F E E t Qˆ (ˆ ˆv) (ˆ ˆv) = ˆ ∂ − ∂ + ∂ − ∂ + ∂ ∂ η ξ (3a,b,c) 其中

(3)

        = rv ru rp J Qˆ 1 ,             − = r v Re p J S 2 0 0 1 ˆ ,           + + + + = t r t z rv rp rvU ru rp ruU ru J E ξ ξ ξ ξ β 1 ˆ ,           + + + + = t r t z rv rp rvV ru rp ruV rv J F η η η η β 1 ˆ           + + + + = r r z z r r z r z z v rv v u r v u r ru JRe E ξ ξ ξ ξ 2 ) ( ) ( 2 0 1 ˆ ,           + + + + = r r z z r r z r z z v rv v u r v u r ru JRe F η η η η 2 ) ( ) ( 2 0 1 ˆ 適當的選取β 值及採取下一節中所述之雙時步 法可求得同時滿足不可壓縮連續方程式(1)與 動量方程式(2a,b)之流場解(p,u,v)。 管壁運動部份 將管壁面上所受壓力p和管壁本身因變形 所產生之應力代入力平衡式,在小變形之假 設下可推導出管壁徑向位移Rw(z,t)之運動方 程式: ) ( 2 o o w o o w w w R R R R Yh p p t R h = − − − ∂ ∂ ρ (4) 其中 Y為管壁之楊氏係數,ρwh為管壁之密 度與厚度, Ropo為管半徑和壓力之參考 值。詳細之推導可參見文獻[14]。 邊界條件 藉由下述之座標轉換: ) , ( , t z R r z = = η ξ (5) 我們可將在物理座標(z,r)下隨時間移動之管 壁位置r =R(z,t)轉換成計算座標(ξ,η)中η=1之 固定邊界位置形式。在壁面上(η=1),結合管 壁與流場運動之動態邊界條件為: 0 = u (6a) t R v=∂ w ∂ (6b) 在管中心軸上(r =0),軸對稱之邊界條件為: 0 = v (7a) 0 = ∂ ∂u r (7b) 由於我們僅計算一段開放式之管路,因此需 給定管流入口和出口處之流動條件。在上游 處我們採用 Ling 和 Atabek[15]以局部脈衝流 近似理論所推導出之流速與壓力分布作為入 口處之條件。在下游出口處則假設流體不再 受拘束而滿足所謂之 traction free 條件。 2.數值方法 由於在流場方程組中之連續方程式(3a)裡 我們引入了一虛擬壓縮之時間微分項,而真 正不可壓縮流場之連續方程式中並無此項, 因此在求解過程中我們採用所謂之雙時步法 則(dual time-step)[16]。亦即,我們首先僅對 動量方程式(3b,c)作真實物理時步∆t之差分, 然而此差分式所求得之流場(p,u,v)並不滿足不 可壓縮流場之連續方程式。因此在此真實物 理時間步驟內另引介一虛擬時間步驟∆τ ,迭 代流場之(p,u,v)值使其在虛擬時間τ→∞時滿 足虛擬連續方程式(3a)之穩態解;亦即迭代 0 ) ( ) ( ) ( 1 = ∂ ∂ + ∂ ∂ + ∂ ∂ J rV J rU J rp η ξ τ β (8) 式直至 ( ) J rp τ ∂ ∂ 項趨近於 0,而後再步進至下一 個真實物理時刻。 在數值差分法方面,我們先以有限體積 法(finite-volume method)來離散流場方程式。 在空間差分上,我們採用 upwind difference 配 合 TVD 法則[17-18]來求取通量(flux)之差分, 在時間差分上則採兼具隱性穩定和節省計算 量之 LUSGS 法[19-20]。將每一物理時刻所計 算出之流場壓力代入管壁運動方程式(4)則可 算出管壁之位移量Rw和位移速度∂Rwt 。整 個計算流程如下: (1) 給定流場之初始和邊界條件; (2) 求取流場方程式(3)之數值解; (3) 將步驟 2 所解得之流場數據代入管壁運動 方程式(4)中計算出管壁位移和速度; (4) 前進至下一時刻,以步驟 3 之結果作為新 的流場邊界條件,重複步驟 2-4 直至流場 循環數個週期為止。 四、研究結果 本計算範例中所使用有關血液及管壁物 理性質之參數值如下: 流場最大參考流速 umax =100cmsec 流場雷諾數 Re=3000 流場參考壓力 po =100mmHg 血管參考半徑 Ro =0.6cm

(4)

血管管長 Lo =10Do =12cm 管壁厚度 h=0.06cm 管壁密度 1.3g cm3 w= ρ 管壁楊氏係數 8 106 2 cm dyne Y= × 收縮-舒張一週期時間 T=0.8sec 我們選用127×17之網格(圖二)分別計算了等管 徑 、 漸 縮 管 徑 (taper =1o ) 、 頸 縮 管 徑 (Dmin =0.7Do)三種情況下之流場,各個不同管 徑變化例子中管入口處之瞬時流速分布(如圖 三所示)皆是由積分 Ling 和 Atabek[15]所提出 之局部近似模式而得。圖四與圖五比較了位 於等管徑和頸縮管徑上、中、下游處管中央 軸上之壓力隨時間一週期之變化,由圖中可 知各處壓力差異不大。圖六則比較了等管徑 及漸縮管徑之瞬時軸向壓力分布,圖中顯示 了漸縮管徑需較大的壓力差來推動血流,因 此心臟所需承擔之負荷亦較大。頸縮管徑最 顯著之效應在於頸縮處管壁之彈性位移較等 管徑情形下小了許多(圖七),此意味著血管頸 縮處管壁儲存和釋放流體動能之能力明顯下 降。圖八亦顯示了在管徑頸縮處管壁所受之 剪應力大增,造成血液流經此處時紅血球較 易產生破裂現象。同時,流經此處之血小板 亦較易活化而於頸縮處後方之迴流區形成血 栓(thrombosis)情形。 數值模擬最大的優點在於其不受現實技 術之限制而能任意變更流場之物理參數值, 使我們得以便利的了解與掌握各個參數對流 場之影響效應。在繁雜之心血疾病發生過程 中往往存在著許多與流體及流場息息相關之 現象。日後冀望藉由更精確與實際之數值模 擬來幫助我們探討心血疾病形成的原因和提 供醫學臨床診斷上的一個有利的工具。 五、參考文獻

[1] Perktold, K., Peter, R. “ Numerical 3D-simulation of pulsatile wall shear stress in an arterial T-bifurcation model, ” J. Biomed. Eng., Vol. 12, pp. 2-12, 1990.

[2] Perktold, K., Peter, R. O., Resch, M., Langs, G. “Pulsatile non-Newtonian blood flow in three-dimensional carotid bifurcation models: a numerical study of flow phenomena under

different bifurcation angles, ” J. Biomed. Eng., Vol. 13, pp. 507-515, 1991.

[3] Paz, D., Einav, S., Elad, D., Avidor, J. M. “Numerical model of blood flow through a bifurcation: interaction between an artery and a small prosthesis, ” Medical & Biological Engineering & Computing, Vol.

30, pp. 543-550, 1992.

[4] Perktold, K., Thurner, E., Kenner, Th. “ Flow and stress characteristics in rigid walled and compliant carotid artery bifurcation models,” Medical & Biological Engineering & Computing, Vol. 32, pp.

19-26, 1994.

[5] Tandon, P. N., Rana, U. V. “A new model for blood flow through an artery with axisymmetric stenosis,” Int. J. Bio-Medical Computing, Vol. 38, pp. 257-267, 1995.

[6] Tang, D., Yang, C., Huang, Y., Ku, D. N. “ Wall stress and strain analysis using a three-dimensional thick-wall model with fluid-structure interactions for blood flow in carotid arteries with stenoses,” Computers & Structures, Vol. 72, pp. 341-356, 1999.

[7] Tang, D., Yang, C., Ku, D. N. “A 3-D thin-wall model with fluid-structure interactions foe blood flow in carotid arteries with symmetric and asymmetric stenoses, ”

Computers & Structures, Vol. 72, pp.

357-377, 1999.

[8] Tang, D., Yang, J., Yang, C., Ku, D. N. “A nonlinear axisymmetric model with fluid-wall interactions for steady viscous flow in stenotic elastic tubes,” J. Biomech. Eng.,

Vol. 121, pp. 494-501, 1999.

[9] Vradis, G. C., Otugen, M. V. “ The axisymmetric sudden expansion flow of a non-Newtonian viscoplastic fluid,”J. Fluids Eng., Vol. 118, pp. 193-200, 1997.

[10] Rappitsch, G., Perktold, K. “ Computer simulation of convective diffusion processes in large arteries,”J. Biomechanics, Vol. 29,

pp. 207-215, 1996.

[11] Rappitsch, G., Perktold, K. “ Pulsatile albumin transport in large arteries: a numerical simulation study,” J. Biomech. Eng., Vol. 118, pp. 511-519, 1996.

(5)

[12] Eler, D. F., Blackketter, D. M., Budwig, R. S., Johansen, K. H. “The influence of shape on the stresses in model abdominal aortic aneurysms,”J. Biomech. Eng., Vol. 118, pp.

326-332, 1996.

[13] Downing, J. M., Ku, D. N. “ Effects of frictional losses and pulsatile flow on the collapse of stenotic arteries,” J. Biomech. Eng., Vol. 119, pp. 317-324, 1997.

[14] Perktold, K., Rappitsch, G. “Mathematical modeling of local arterial flow and vessel mechanics,”In Computational Methods for Fluid-Structure Interactions, eds. J. M.

Crolet & R. Ohayon, Longman, New York, pp. 230-245, 1994.

[15] Lin, S. C., Atabek, H. B. “ A nonlinear analysis of pulsatile flow in arteries,” J. Fluid Mech., Vol. 55, pp. 493-511, 1972.

[16] Roers, S. E., Kwak, D. “ Upwind differencin scheme for the time-accurate incompressible Navier-Stokes equations, ”

AIAA J., Vol. 28, pp. 253-262, 1990.

[17] Van Leer, B. “ Towards the ultimate conservative difference scheme V, A second-order sequel to Godunov’s method,” J. Comput. Phys., Vol. 32, pp.

234-245, 1979.

[18] Van Leer, B. “Upwind-difference methods for aerodynamic problems governed by the Euler equations,”Lett. Appl. Math., Vol. 22,

pp. 327-336, 1985.

[19] Yoon, S., Jameson, A. “ An LU-SSOR scheme for the Euler and Navier-Stokes equations,”AIAA Paper no. 87-0600, 1987.

[20] Yoon, S. Kwak, D., Chang, L. “LU-SGS implicit algorithm for three-dimensional incompressible Navier-Stokes equations with source term,”AIAA Paper no. 89-1964,

1989. 六、圖表

r

z

圖一、座標示意圖。 4 5 6 0 0.25 0.5 0 2 4 6 8 10 0 0.5 r z 圖二、127×17 網格配置圖。 0 0.2 0.4 0.6 0.8 u / umax 0 0.2 0.4 0.6 0.8 1 r /R τ= 0.05 τ= 0.10 τ= 0.20 τ= 0.40 τ= 0.50 τ= 0.60 圖三、管入口瞬時流速分布圖。

(6)

0 0.2 0.4 0.6 0.8 1 t / T 40 60 80 100 120 140 in m out ( m m H g ) p p p p 圖四、等管徑上、中、下游處於軸對稱位置 上壓力隨時間之變化情形。 0 0.2 0.4 0.6 0.8 1 t / T 40 60 80 100 120 140 in m out ( m m H g ) p p p p 圖五、頸縮管徑上、中、下游處於軸對稱位 置上壓力隨時間之變化情形。 constant 0 2 4 6 8 10 z 70 72 74 76 78 80 82 84 p taper = 1o ( m m H g ) 圖六、 t=0.05時等管徑與漸縮管徑軸向瞬時 壓力分佈圖。 0 0.2 0.4 0.6 0.8 1 t / T -6 -4 -2 0 2 4 6 in m outR x 1 0 ( m m ) ∆RRR 圖七、頸縮管上、中、下游處管徑位移量隨 時間之變化。 0 0.2 0.4 0.6 0.8 1 t / T -2.0 -1.0 0.0 1.0 2.0 in m out τw τw τw τw x 1 0 -3 ( d y n e /c m 2 ) 圖八、 頸縮管上、中、下游處管壁壁面上剪 應力隨時間之變化。

(7)

參考文獻

相關文件

Stress and energy distribution in quark-anti-quark systems using gradient flow.. Ryosuke Yanagihara

 Replace the wall in observation room with the projected image of the remote room...

In flow field boarding up for the smooth wall surface, gets down the board is the Sinusoidal wavy wall, but under board wavy wall oscillation amplitude A=2.54(mm), the wave length λ

“Model Tests on Excavation Problems with Different Wall Friction and Wall Stiffness,” Proceedings, 32th Conference of Japanese Society of Soil Mechanics and Foundation

Followed by the use of an important degree of satisfaction with the service quality attributes, by Kano two-dimensional quality model, IPA analysis and

The result indicated that the constitutive relationship among the flow stress, deformation temperature, and strain rate for IN 600 during hot deformation satisfied the

The analysis of biomechanics after the total hip replacement can be reached by comparing with the stress and strain distribution of the intact femur and “Tripot Stemless Total

From those examples, it indicates that the efflorescence will reoccur after three to six months using the waterproofing membrane painted on the inside of an exterior wall and it