報告內容
(1) 前言: 血管是人體重要得組織之ㄧ，負責運送血液至身體各處，進行養分、能量的供應，廢 物排出及氣體交換，維持身體機能的運作，而主動脈則是人體最大口徑的動脈，直接承接 左心室所博出之充氧血液，以循環至全身。但隨著時代的進步，生活飲食習慣的改變，心 血管疾病已成為危害人類生命健康最嚴重的疾病之ㄧ，民國九十七年的十大死因中，心臟 血管疾病排名第二，相較於前幾年的第三名提升了不少[4]。 造成心血管疾病的誘導因素有很多，冠心病、主動脈剝離、高血壓症和高血脂症的發 病，除和遺傳基因相關以外，也和個人的生活習慣與生活環境有密切關係。工作壓力繁重， 嗜食高鹽分和多脂肪、高膽固醇食物，吸煙和不運動，都會讓病情惡化，而導致心血管疾 病和中風甚至死亡。且近年來在臨床上發現心血管疾病年齡層有越來越年輕化的趨向，醫 師認為這和現代社會型態及生活改變有關。所以心血管疾病發生原因之探討及疾病發生後 對身體機能所造成的影響已成為現代人不容忽視的問題。 近來許多研究以血液動力學(Hemodynamic)來研究血液在血管內的流體運動，與心血管 疾病發病機制的相對關係。例如異常的血流動力因數(Hemodynamic Factors)：血液流速、 壁面剪應力(Wall Shear Stress)及壁面壓力的變異，是形成主動脈剝離(Aortic Dissection)與動 脈粥樣硬化(Atherosclerotic)等心血管疾病的危險因子[5]。其中主動脈剝離係因異常的血管 壁面剪應力與壁面壓力所導致，過高的壁面剪應力會對血管內皮細胞造成傷害，形成主動 脈剝離的現象；而過低的壁面剪應力亦會導致動脈粥樣硬化或是血栓。另外，血液作用於 血管壁面的壓力與血液流速突然升高，血管壁面容易於較薄處開始破裂或是產生剝離的現 象。 現代醫學影像學進步快速，已發展出多套非侵入方式的影像技術，例如：放射線閃爍 造影術(Radionuclide Scintigraphy)、超音波(Echocardiography)、電腦斷層掃描(Computed Tomography，CT)以及核磁共振影像(Magnetic Resonance Imaging，MRI)。其中核磁共振影 像不同於一般的X光攝影或電腦斷層掃描，它並不是利用X光來形成影像。其成像原理簡單 來說是利用原子「核」自轉及在「磁」場中「共振」的現象及特性，經由電磁波來形成影 像的。除了不具游離輻射之外，核磁共振影像還有許多優點，包括：優良的軟組織對比、 高的解析度、可同時取得不同切面的影像，且若需要注射顯影劑時，不僅注射量少，而且 人體幾乎不會產生過敏反應。核磁共振影像除了顯示詳細的解剖構造之外，更可以提供生 理及生化的訊息，正因核磁共振影像具有這些優點，而且對人體幾乎無任何傷害性，因此 成為臨床診斷以及醫學研究的最佳利器之一。 (2) 研究目的: 主動脈剝離分類中最廣泛引用的為Debakey et al.[1]採用血管內膜撕裂的所在位置與主 動脈剝離的範圍作為分類要點，如圖1 所示。第一型的剝離發生在主動脈瓣膜以上，並擴 展到降主動脈；第二型的剝離發生只限於昇主動脈內；而第三型的剝離發生在左鎖股下動 脈近端的降主動脈中並延伸到腹部的主動脈。Erbel et al.[6]分類則認為Debakey一型和二型 是為近側的或上升段的主動脈剝離，又Debakey三型則為遠端的或下降段的主動脈剝離。Barakat et al.[7]利用兔子的主動脈為一迴路實驗，說明了主動脈內複雜幾何形狀所產生的流 場特性，和確立異常血液動力學行為與動脈粥樣硬化的病理學現象間之相互關係。然而， 未治療的主動脈剝離病患在24小內時有21%的死亡率，48小時內有37%死亡率，2週內有74% 的死亡率，3個月內有90%死亡率[8]。病理學指出動脈呈現僵硬、血管內壁的異常增厚，將 導致心臟收縮壓升高，造成周邊的血管壓力增加，左心室肥大，血壓異常增高，甚至嚴重 狀況下，因異常的血流動力因數(Hemodynamic Factors)如血液流速、壁面剪應力及壁面壓 力的變異，更是形成主動脈剝離(Aortic Dissection)與動脈粥樣硬化(Atherosclerotic)等心血管 疾病的危險因子[5]。其中主動脈剝離係因異常的血管壁面剪應力與壁面壓力所導致，過高 的壁面剪應力會對血管內皮細胞造成傷害，形成主動脈剝離的現象；而過低的壁面剪應力 亦會導致動脈粥樣硬化或是血栓。另外，血液作用於血管壁面的壓力與血液流速突然升高， 血管壁面容易於較薄處開始破裂或是產生剝離的現象。又以胸主動脈(Thoracic Aorta)和它 的三個主要分岔處內的血液流場結構較為複雜，因此值得深入研究。
本研究藉由 PC-MRI (Phase-Contrast Magnetic Resonance Imaging, PC-MRI)量測真實人 體胸主動脈弓幾何形狀和內部血流速度分佈，並利用快速成型(Rapid Prototyping, RP)製作 剛體主動脈模型進行體外實驗，並且將實驗的結果與剛體的模擬結果互相比較，做為驗證 計算程式正確性的依據，過程中使用計算流體力學軟體 ACE+®，模擬暫態、三維的胸主動 脈流場，而後由於血管的變形以及升胸主動脈入口為旋轉流場這兩個條件對整個胸主動脈 內流場的影響不可輕忽，所以加入了模擬變形的模組以及利用程式設計軟體 fortran®撰寫升 胸主動脈入口的旋轉流場的副程式，並掛載在計算流體力學軟體 ACE+®中，以期能更精確 的獲得人體胸主動脈內部的資訊。所呈現出來的結果例如三維的流線、截面的軸向速度分 布、壁面剪應力及壓力、壁面剪應力震盪值等，另一方面，也嘗詴將臨床人體以及數值模 擬的脈波傳導速度(Pulse Wave Velocity, PWV)做比較探討，希望能彌補實驗量測以及臨床 觀察上的限制，更進一步的了解人體胸主動脈。 (3) 文獻探討: 3.1 心血管模擬與實驗 過往在血液流場之研究中，已有許多文獻藉由數值模擬的方式，探討因複雜的血液流 動特性，如血液流速、壁面剪應力及壁面壓力的變異。1997 年，Barakat et al.[7]利用兔子 的主動脈完成一迴路實驗，其結果說明了主動脈內複雜幾何形狀所產生的流場特性，和確 立異常血液動力學行為與動脈粥樣硬化的病理學現象間之相互關係。同年，Perktold, et al.[9] 使用雷射都普勒測速儀的實驗方式，量測體外冠狀動脈分枝的入口速度與內部的軸向速度 分佈，並將其實驗量測得到的數據作為數值模擬的邊界條件與驗證依據。體外冠狀動脈分 枝的實驗與數值計算之幾何模型則依真實冠狀動脈的幾何形狀，以透明且不具彈性的圓管 來近似。研究過程中先說明了冠狀動脈的軸向速度隨心臟脈波的週期性變化分布，爾後更 從截面二次流的變化探討壁面剪應力的方向隨週期性的變化關係。1998 年，Gijsen, et al.[10] 以實驗與數值模擬 90 度彎曲的血管中，於暫態的條件下比較血液依據牛頓及非牛頓流體的 兩種假設條件，輸入一給定的流量波形為入口條件，以觀察其內部血液的流動與二次流的 變化。影響血管流場的主要因素為雷諾數，數值計算顯示，模擬中採用牛頓或非牛頓流體 條件對所得結果無明顯差異。1999 年，Hugo, et al.[11]利用 MRI 所得到的三維流場數據加 上時間解析，觀察到血液的流動方式，在健康者的升胸主動脈到主動脈弓(Aortic Arch)之
間，是呈現順時針螺旋方向前進，在主動脈弓之降胸主動脈之間，則是以逆時針螺旋方向 前進，且此流動特徵在年輕者要比年長者更為明顯。2000 年，Zabielski, et al.[12]藉由一近 似主動脈弓的螺旋漸縮圓管，以數值計算方式，探討在不同種活體體重下，其螺旋漸縮圓 管會有不同的管徑與血液流量，進而改變 Womersley numbers 與 Reynold numbers；導致血 管截面二次流的速度分佈與壁面剪應力的變化會有不同的結果，對於動脈粥樣硬化提供了 另一種研究觀點。同年，Zhao, et al.[13]利用 MRI 掃描出頸動脈的外型並且加入彈性流體模 式觀察其壁面剪應力、速度與二次流的變化，並且於流場型態中發現，動脈粥樣硬化好發 於流場分離與停滯流的區域。2002 年，Shahcheraghi, et al.[14]利用 CAT SCAN Image 掃描 出年輕人的主動脈外型，並且利用一假定波型為入口邊界條件配合程式模擬，以探討主動 脈內速度、壓力、剪應力的變化，其中指出了無論高剪應力區或是低剪應力區或是剪應力 變化較大的區域，皆易導致局部的動脈粥樣硬化。2003 年，Buchanan, et al.[15]採用數值計 算方式，探討兔子腹動脈內部血流動力因素對動脈粥樣硬化的關係，研究過程中先預測出 血流動力學的各項壁面參數；在數值計算過程中使用非牛頓流體並導入單核白血球，觀察 其沉積位置與血流動力學的各項壁面參數的相對關係，提供在醫學研究上，藉此了解血管 內壁病變的好發位置與原因。同年，Jin, et al.[16]利用 MRI 掃描出人體胸主動脈的模型以及 流場的情況，並找出胸升主動脈入口端晃動的軌跡，並以此當做模擬的邊界條件，加入胸 升主動脈部分壁面晃動以及主動脹縮等條件，結果顯示出了考慮血液與壁面的流固耦合 (Fluid-Structure Interaction)後的壁面剪應力的變化。2006 年，Nakamura, et al.[17]利用加入 左心室的主動脈數值模型，取代均勻流的胸升主動脈入口條件，結果顯現出心室對主動脈 壁面剪應力的影響需要被列入考慮。 國內相關研究：國防大學丁大為教授[18-19]，清華大學劉通敏教授[20]等人皆使用數 值計算的方式，各別針對冠狀動脈、肺動脈、顱內之動脈瘤與血管內部植入枝架後的血液 流場加以分析。而中華大學牛仰堯教授等人[21]使用 MRI 掃描分析出一位 21 歲健康年輕 人之主動脈相關幾何形狀與流場數據，並利用計算流體力學加以分析探討，但其使用 MRI 分析過程時，並未完整的定義出胸主動脈壁面邊界位置，亦未對胸主動脈的病變進行分析。 因此，對於胸主動脈的壁面剪應力分析尚有探討空間。國內研究另有國家高速網路與計算 中心的李明龍博士等人[22]，係使用數值模擬冠狀動脈繞道血管的血液流場分析，其中在管 壁變型的模擬過程，因計算流體力學和固體力學分析過程做分開運算，此可能造成在數值 分析流場的時間解析度方面會有誤差產生，另外成功大學溫志湧教授、台北科技大學楊安 石教授等人[23-24]，以體外胸主動脈實驗和電腦數值模擬將血管內流場穩態以及暫態的結 果相互對照。以上結果皆證實了數值模擬可合理精確的應用在預測血管內部流場特性及血 流動力因素，並將觀察到的結果以及病變做結合。 3.2 脈波傳導速度
脈波傳導速度(Pulse Wave Velocity，PWV)為量測動脈硬度時，廣泛被使用的指標，PWV 定義如下： C 1 ρ A PWV (1) 式(1)中 A 為主動脈截面面積，ρ 為血液密度，C 則為血管壁的順應性(Compliance)。而血 管區域順應性(Local Area Compliance) 的定義如下:
ΔP
ΔA
C
(2)
式(2)中 ΔA 是收縮前後主動脈截面面積變化，ΔP 則是收縮前後血管內壓力差(Artial Pulse Pressure)。C 越小則 PWV 值越高，代表動脈血管壁越硬，血管彈性越差。2006 年，Yu, et al.[25] 以一套全新的軸向速度分佈方法來計算 PWV。藉由 MRI 影像找出降胸主動脈的軸中心位 置，並計算軸中心上的血流速度分佈。將軸心血流速度(Z-axis)沿軸心位置(X-axis)和時間 (Y-axis)的相對關係繪製成三維脈波波形行進圖(如圖 2 所示)，利用 MRI 所獲得數據所歸 納出來的方程式，進而計算出降胸主動脈的脈波傳導速度。 由 X 和 Y 平面投影方向可以清楚的發現在各個時間點下速度的等高線圖。從圖中可明 顯的觀察出波峰波谷的效果即為造成 PWV 值差異的關鍵。因此血流速度與軸心位置和時 間的關係可以一線性式近似。 z=ax+by+c (3) PWV=b/a (4) 利用 MRI 所獲得數據所歸納出來的方程式(3)、(4)，即可計算出主動脈波速。這種使用軸 向速度分佈圖的方法不僅可直接評估脈波的傳導速度，並且只需較少的假定即能以較直覺 與客觀的方式決定 PWV。 (4) 研究方法: 首先藉由PC-MRI影像分析並獲得一般正常人的胸主動脈外型，以Mimic®軟體將此胸主
動脈外型轉成iges檔，再利用電腦輔助繪圖(Computer Aided Drawing,CAD)讀入轉成快速成 型(Rapid Prototyping, RP)技術所需之CAM檔，製作成人體胸主動脈的透明模型，用以進行 體外(In-Vitro)胸主動脈實驗。圖3所示為體外胸主動脈流場實驗設備的管路系統，實驗初期 驗證階段係使用水為工作流體，環境溫度約27℃條件下，其密度與黏滯係數為997 kg/m3 與 0.855 cp。圖中之沉水馬達M201為三星公司所製造，最大流量為22 L/min，用以推動液體流 動。根據Francis[26]所引述Martini 的方法得知，主動脈內各部位的流量比例分別為；上行 主動脈佔100%，下行主動脈佔77.4%，頭臂動脈佔11.3%，左頸動脈佔6.8%，左鎖骨下動脈 佔4.5%，為了量測各部位的流量比例，所以於上行主動脈入口前端加裝了電磁式流量計 (ISOIL ML-190)。其餘部位皆裝設四只浮子式流量計(Dwyer)，依序放置於降胸主動脈、頭 臂動脈、左頸動脈和左鎖骨下動脈之出口段，另使用球閥控制各部位流量。另外，電磁式 流量計之內部管路裝有整流蜂巢段(Honeycombs)，將入口剖面流速調整為均勻速度之流場 流入。並於升胸主動脈入口、降胸主動脈、頭臂動脈、左頸動脈和左鎖骨下動脈之下游處 亦放置5組壓力感測器(Druck PTX-1400)以量測穩態條件下之靜壓力值，並作為數值模擬之 邊界條件設定。 體外與活體胸主動脈之流場量測與追蹤血管壁面邊界位置方面，係藉由台中榮民總醫 院提供之1.5Tesla PC-MRI灰階影像分析系統(Sonata, Siemens, Erlangen, Germany)，以心電 訊號圖(ECG)的Q波做為觸發訊號，針對穩態與暫態條件下之體外胸主動脈作一系列之掃 描，而流場影像分析之參數設定如表1所示。若以一次心臟週期1T=825ms為例，各軸影像 皆取76張，每張影像間隔則為11ms。在PC-MRI掃描胸主動脈過程中，針對上型主動脈至下 行主動脈以矢狀切面(Sagittal Plane)進行一連串的掃描、編碼，最後得到XYZ三軸的血流速
度數據。在PC-MRI掃描過程中之Y軸切面；即矢狀切面，所表現的色階變化最為明顯，容 易以視覺分辨出血管壁的大略位置，所以在各個時間點下皆選用Y軸速度的PC-MRI胸主動 脈灰階影像，作為B-Spline曲線近似主動脈邊界的依據。在體外胸主動脈壁面邊界搜尋方法 中，圖4所示為搜尋主動脈邊界的程式流程圖。本計畫係利用數學軟體Matlab®作為搜尋胸 主動脈邊界的工具，其過程先以B-Spline曲線去近似胸主動脈邊界，再以血流速度及速度分 布之曲率變化作為判斷血管邊界依據，即可定義出胸主動脈於一個心血管週期下的壁面位 置，接著計算出胸主動脈壁面剪應力值與震盪剪應 力指標。 取PC-MRI掃描後的Y軸灰階速度影像作為B-Spline曲線的依據，於血管內靠近管壁的 地方選取適當的控制點，並將其連成一條圓滑的B-Spline曲線，在B-Spline曲線上各點的法 線方向，於血管內部以間距0.05mm之距離，依序沿著血管內腔壁面的上下表面向外各取100 點，並由內差計算出法線方向上各點的XYZ三軸的血流速度，並將各軸速度投影至主軸方 向，以取得流體速度分布圖及其曲率變化圖，如圖5所示為體外胸主動脈模型之PC-MRI影 像與B-Spline曲線。然而於追蹤血管壁面邊界的部份，依據流體力學的No-Slip boundary condition的精神，得知血管壁面上的血液流速與血管壁之相對速度零，所以接觸血管壁面 上的血流速度分布與壁面之間會有一個明顯的曲率(血流速度分佈對法線方向的二次微分) 變化，而且在血管壁面的位置上，理想狀況下沿壁面方向(切線方向)速度應為零，如圖6(a)(d) 所示，圖6(d)代表有迴流現象時之速度分佈情形。圖6(b)(c)(e)(f)則代表血管外面速度有一近 乎均勻背景雜訊之情形，其中圖6(e)(f)代表有迴流現象時之速度分佈情形。因此，透過血流 速度分布及其曲率變化的探討即可定義出血管壁面正確的位置。又因曲率為斜率的變化 率，及速度的二次微分，故可利用內差求出向外延伸每100點上的個別速度，再以公式(5) 計算出曲率值[2-3]。有了速度分布與曲率變化可正確的判別出體外主動脈模型的邊界位 置，再利用公式(6)求出三維體外主動脈的壁面剪應力。 (5) (6) 其中ν為各點的速度值，μ為血液的黏滯係數，u為流體切線方向的速度，s與n各為切線 方向與法線方向的距離。 本研究使用計算流體力學軟體 ACE+®來探討人體胸主動脈內的流場情形，先求解穩態 的流場，而後暫態流場，最後延伸至考慮血液以及血管壁流固耦合(Fliud-Structure Interaction) 的現象，在此表示的為考慮至流固耦合的完整研究方法，其穩態的為忽略固體的部分以及 假定入口速度固定為一定值，暫態的結果為僅忽略固體的部分。 在求解流體的部份，假定血液為脈動(Pulsatile)、不可壓縮以及牛頓流體，使用三維且 隨時間改變的連續方程式及動量方程式。由文獻中可知[27-29]，當血流剪切率(Shear Rate) 低於 100 s-1時，人體的血液傾向表現出非牛頓流體的行為，但在較大型的血管如胸主動脈 等其剪切率都大於 100 s-1，所以流體為牛頓流體在胸主動脈模型中為合理的假設;在求解固
體的部份，用位移的方程式來描述血管壁面的變形。 在暫態不可壓縮層流以及流體性質為常數的情況下，隨時間變化的連續方程式(7)以及 Navier-Stokes 方程式(8-10)如下: . 0 ) ( ) ( ) ( z u u y u u x u ux xm y ym z zm (7)
. . 2 ) ( ) ( ) ( x z x x y x x zm z x ym y x xm x x g x u z u z y u x u y x u x x p u u u z u u u y u u u x u t (8)
y z y y x y y zm z y ym y y xm x y g y u z u z y u y y u x u x y p u u u z u u u y u u u x u t
. 2 ) ( ) ( ) ( (9)
z z z y z x z zm z z ym y z xm x zg
z
u
z
y
u
z
u
y
x
u
z
u
x
z
p
u
u
u
z
u
u
u
y
u
u
u
x
u
t
.
2
)
(
)
(
)
(
(10) 其中 ux /uy /uz 以及 uxm /uym /uzm 分別代表在 x 軸、y 軸以及 z 軸上的流體以及網格速 度，而p、ρ、μ、g 分別為壓力、密度、黏滯係數、重力加速度等，其中流體密度為 1121 Kg/m3 。 用以求解偏微分方程式之數值方法大致上可分為有限差分法(Finite Difference)、有限體積法 (Finite Volume)以及有限元素法(Finite Element)。本研究中以有限體積法離散質量與動量守 衡方程式，進而將統御方程式以代數方程式表示，其中採用 Van Doormal et al.[30]所提出的 SIMPLEC (Semi-Implicit Method for Pressure-Linked Equations Consistent)演算法，以計算出 流體區域內部的速度與壓力分布。將 原 方 程 式 積 分 後 ， 其 中 對 流 項 (Convection Term) 採 用 二 階 中 央 差 分 法 (Central Difference)離散，擴散項(Diffusion Term)亦是用二階中央差分法，另外在暫態項的部份則是 用二階的科雷克-尼克生方法(Crank-Nicolson Scheme)以達到數值穩定的計算。在暫態流場 中，不同的時間間隔會影響解的穩定性，其需要滿足 Courant-Friedrichs- Lewy (CFL) Number 小於 1，在此 CFL Number =(區域速度值)×(最大時間間隔)/ (區域內網格長度)，在本研究的
計算中，流體變數如三軸速度以及壓力的標準殘值誤差在每個時間點都必頇小於 10-4，且
每個週期同一時間點之間質量的誤差必頇在 0.5%以內。
在固體的部份，我們假定胸主動脈壁面為無孔材料(Non-Porous)且僅會產生小形變 (Deformation)，壁面亦為速度不滑動邊界(No-Slip Boundary Condition Of Velocity)，且在邊 界上的壓力梯度為 0。由於假定壁面只產生小變形，我們假設固體模型為彈性(Elastic)、均 質(Homogeneous)和等向性(Isotropic)，以線彈性方法(Linear Elasticity Approach)來求解變
形，其統御方程式(11-19)如下:
.
2 2 x s zx yx x x sg
z
y
x
t
v
(11).
2 2 y s zy xy y y sg
z
x
y
t
v
(12).
2 2 z s yz xz z z sg
y
x
z
t
v
(13) s 代表壁面的密度， v、和分別代表位移、應力和應變向量(張量)。
v
xv
yv
z
v
,
,
(14)
T yz xz xy z y x
,
,
,
,
,
(15)
T z y z x y x z y x T yz xz xy z y xy
v
z
v
x
v
z
v
x
v
y
v
z
v
y
v
x
v
,
,
,
,
,
,
,
,
,
,
(16) 應力-應變的關係為線性，
D
。 G G G G G G D 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 2 0 0 0 2 (17)
2 1 1 E (18)
1 2 E G (19) D(Young’s Modulus)和包松比(Poisson’s Ratio)。Riley et al.[31] 利用非侵入式的超音波技術來 決定 3,321 位成人動脈的楊氏係數，其結果顯現出量測值介於 701 kpa 到 983 kpa 之間，本 研究的楊氏係數和包松比分別為 840kpa [32] 和 0.45 [33]，固體密度為 1161 Kg/m3。數值 模型中存在流體與固體的交介面，在流固耦合的邊界上必需要滿足兩個條件: (1) 流體移動 的量和固體移動的量必頇要相同; (2) 作用在流體上的力(壓力和剪應力) 和作用在固體上 的力相同。 詳細的計算步驟如下: 1. 將計算所得的穩態解 ( 整個脈動週期中，質量流率最大的時間點所獲得的數值解) 作 為暫態解 t = t0 時的初始條件。 2. 用方程式(7)-(10) 來求解流體的速度以及壓力，此後再解出作用在流體與固體交界面上 流體的力。 3. 用方程式(11)-(19) 來求解血管壁面變形的位移量。 4. 如果流體速度、壓力以及壁面變形位移量的標準殘值大於 10-4 (未收斂)，則壁面變形 位移量將加諸於流體部份的邊界上，使得流體部分和固體部分的網格有所變動。 5. 重複步驟(2)-(4) 直到流體速度、壓力以及壁面變形位移量的標準殘值小於 10-4 (收 斂)。本研究中，此步驟通常需要 500-700 個疊代以達到收斂標準。 6. 將步驟(5) 獲得的數值解當作新的時間點下的初始條件，且每個時間點的間距為 0.02 秒 (tnew= told + 0.02 秒)。 7. 重複步驟 (2)-(6) 直到計算區域的流場和血管壁面的變形量獲得週期性的解。 在本研究中，需要 5 個週期才能求得流場和固體變形週期性解，所以我們用第 6 個週 期的解做為最後呈現的結果。首先從入口流量在心臟收縮峰值(Peak Systole)開始計算脈動 的流場，接下來進入減速度相位(Decelerating Phase)以及舒張末期 (Later Diastole)，最後到 達加速度相位(Accelerating Phase)，經過 6 個週期總共有 300 個時間點 (心跳周期為 1 秒)。 此外，胸升主動脈的入口流量波形如圖7。 (5) 結果與討論: 建構完成的體外胸主動脈模型如圖8，在架設整個體外胸主動脈實驗設備與管路系統方 面，也已經完成如圖9。此外，用Matlab®建構程式協助PC-MRI分析體外胸主動脈穩態流場 速度分布，以及B-spline曲線，如圖10、11。 圖12-18為穩態數值分析的結果，係應用ACE+®計算流體力學分析軟體以模擬、驗證主 動脈在穩態條件下的流場特性。圖12展示數值計算之主動脈幾何尺寸、座標系統及網格分 布，計算區域的總格點數為244867。 圖13列出理論模擬與PC-MRI在平面上絕對速度等位線的分布，在(y*=y/D=0.8，其中D 為入口直徑)剖面上，理論模擬與MRI量測灰階影像的絕對速度分布的中，可得知流體沿著 主動脈弓到達下行主動脈時，靠近外側壁面之速度分布皆呈現較大的絕對速度，且最大之 絕對 速 度 相 對誤 差值 約為 1%。由圖14可發現粒子軌跡在頭臂動脈近側壁面 (Proximal Walls)、左頸動脈近側壁面與左鎖骨下動脈中間及近側壁面皆以渦旋流動，由圖中的流線可 以得知，在這三處形成數個逆時鐘方向的渦漩結構。 圖15為胸主動脈之全域流場截面速度分布曲線圖，圖中於截面a與b處可以發現，因上
升主動脈銜接主動脈弓的彎管幾何形狀因素，使得上升主動脈外側壁面附近的流體迴旋流 動後，速度較上升主動脈內側壁面大，直接作用於鄰近頭臂動脈近側壁面的底部與主動脈 弓外側壁面上的接何處。截面c、d、e可以發現，主動脈弓中心流體因歧管效應而偏離移動 並翻轉進入頭臂動脈、左頸動脈與左鎖骨動脈，所以主動脈弓外側壁面速度較大且直接作 用於各分岔歧管底部壁面上。 圖16為升主動脈、主動脈弓與下降主動脈各處橫截面流線圖。於截面a可發現，分別於 截面的右上方與左下方有二次流渦漩存在，又因為上升主動脈與主動脈弓之間的截面積逐 漸縮小，彎管曲率明顯改變且有扭轉與分岔歧管效應存在，所以於截面b、c、d皆可得知有 穿越流從主動脈弓內側逐漸向截面中心增加，再經由主動脈弓中心翻轉進入各分岔歧管 內。截面f上之右上方與左下方存在微小的二次流渦漩，且沿著出口方向逐漸減小，僅為 11.6%u0(其中u0為穩態入口條件下之入口流速) 。主動脈壁面壓力分布亦為重要的血流動力 參數之一。臨床發現，當血管壓力與血液流速突然升高，血管壁面容易於較薄處開始破裂 或是產生剝離的現象[34] 。 圖17為主動脈壁面動壓分布圖，由圖中發現，主動脈弓的三個分岔處之基底鄰近主動 脈弓外側壁面之交接處，和主動脈弓的左右兩側皆有相對的高壓存在，此最高之流體總壓 力與最大之壁面剪應力共同作用於此三個交界處，此流體之應變率與壁面剪應力成正比關 係。由圖18可發現此區域的應變率皆為最大值，故其可能對血管內皮細胞造成傷害，形成 主動脈剝離的現象。而於主動脈弓左右兩側與內側壁面上，和降主動脈前段之外側壁面上 顯現有較高之流體動壓力，且由圖18亦可得知這幾處有次高的應變率發生，合理推測得主 動脈弓左右兩側與內側壁面上，和降主動脈前段之外側壁面上也可能形成主動脈剝離的現 象，一旦血管壓力與血液流速突然升高時，皆有可能導致血管病變而危及到生命。而上述 幾處高危險區域皆與文獻[1][5]有高度的正相關性。 圖19-24為暫態數值分析的結果，圖19展示暫態模型的數值計算之主動脈幾何尺寸、座 標系統及網格分布，計算區域的總格點數為81603，整個模型分成6個部份，包含了升胸主 動脈(Ascending Aorta)、降胸主動脈(Descending Thoracic Aorta)、主動脈弓(Aortic Arch)、頭 臂動脈(Brachiocephalic Artery)、左頸動脈(Left Common Carotid Artery) 和左鎖骨下動脈 (Left Subclavian Artery)。
圖20以及圖21為流場矢狀截面的軸向速度在CFD模擬和PC-MRI實驗結果的比較，以及 四個相位流場的流線圖，由於模擬邊界條件將胸主動脈視為剛體，所以我們已知邊界的位 置，圖中實心黑點即為PC-MRI實驗邊界的位置，在Accelerating Phase (t/T= 1.82)以及Peak Systole (t/T= 1.0)這兩個相位中，Cross Section A處模擬和預測的結果顯示流體彎曲往Inferior 方向的主動脈弓壁面聚集。此外在Cross Section B和D處有複雜的幾何結構，軸向速度形狀 的改變受到二次流現象的影響，而造成二次流現象是由於主動脈弓的曲率以及其三個分 支，和文獻[36]有高度相關。在Decelerating Phase (t/T= 1.18)相位，Cross Section A-C處軸向 速度分佈呈現一個M型的趨勢，且流體在Cross Section D處彎曲往Posterior方向的主動脈弓 壁面聚集，此處呈現一個三度空間的複雜流場現象。
由圖22可發現在Accelerating Phase以及Peak Systole這兩個相位中，頭臂動脈和左頸動 脈流場為層流，而在左鎖骨下動脈有逆時鐘的渦漩的流場，且大部分的流體則是由左鎖骨 下動脈的Posterior方向壁面流過。在Decelerating Phase相位，逆時鐘渦漩出現在主動脈弓三 根分支的前端，且渦漩在各分支Posterior方向造成的影響延伸到下游。在Later Diastole相位 時，流體流出主動脈的三個分支並進入下行主動脈。實驗和模擬的結果皆清楚的顯示主動
脈弓三分支的逆流，以及順時針和逆時針的渦漩產生在頭臂動脈和左頸動脈的入口，而這 兩個渦漩的差異是在於大小以及位置，在頭臂動脈的順時針渦漩相對的較小，且位置較靠 近Anterior方向壁面，其餘順時針的渦漩延伸到左頸動脈Anterior方向的下游處。此外，頭 臂動脈和左頸動脈處的逆時針渦漩沿著Posterior方向的壁面延伸，而左鎖骨下動脈的逆時針 渦漩則是產生在Posterior方向的下游，模擬和實驗的結果皆為一致。 圖23為胸主動脈模型內表面的壁面剪應力分佈圖，在Accelerating Phase相位時，整個胸 主動脈模型呈現一個全面性的低剪應力分佈，然而在流體通過主動脈弓時，向Inferior方向 的主動脈弓壁面聚集，到了下行主動脈時，流體向Posterior方向的壁面聚集，所以此兩處有 較高的壁面剪應力。從圖中也顯示出在Accelerating Phase以及Decelerating Phase相位，壁面 剪應力的分佈以及極值的位置基本上沒有太大的變化，但在大小方面有不同的差異，在 Decelerating Phase相位時，由於主動脈弓的二次流現象還有流場分離的情形造成在主動脈弓 處的幾個部位有較高的壁面剪應力，在Peak Systole相位，整個流場呈現較高的壁面剪應力， 尤其是在頭臂動脈、左頸動脈和左鎖骨下動脈三根分岔Anterior方向的壁面有最高的壁面剪 應力(~22.8 N/m2 )，此一現象的成因為流體從上行主動脈通過主動脈弓造成強烈的二次流現 象。由於在Later Diastole相位時，流體的速度為整個週期中最緩慢的時刻，所以整個模型的 壁面剪應力比較前三個相位，相對的低很多(最大值~2.43 N/m2 )，目前得到的這些結果和文 獻[10][30]的結果有高度相關。 圖24為胸主動脈模型內表面的壓力分佈圖，且此動壓為無因次化後的結果(P =P*/ρuo2 )，符號P、P*和uo分別代表無因次壓力、壓力和平均進口速度。我們知道在管內快速的累 積動壓會導致管壁發生破壞或剝離。在Accelerating Phase時，最大的壓力發生在最大加速度 處，也就是上行主動脈的入口，最小壓力則是在下行主動脈，在Peak Systole相位時，在主 動脈弓的三個分支Posterior方向的停滯點，在這個高壓區，高速的流體由於通過主動脈弓分 支後迅速的減速，另一個高壓區在上行主動脈的Anterior方向的壁面，成因是巨大的曲率變 化造成流體減速，然而最低的剪應力則是在下行主動脈，其餘較低的剪應力出現在主動脈 弓的Inferior方向壁面，從圖上的結果顯示，整個週期高壓的位置出現在主動脈弓的三個分 支Posterior方向的停滯點以及下行主動脈Anterior方向的壁面。 在邊界判讀的部份，我們首先設計一個圓直管實驗來驗證我們搜尋主動脈邊界的程式 正確性，實驗設計如圖25，以一內徑為22 mm，厚度1.5 mm的壓克力圓管，在一固定的流 速下進行PC-MRI拍攝。其中，為使流場為一穩定的完全發展流，根據南赫(Langhaar)之理 論，發展完全入口長度L為雷諾數Re之函數，亦即 L=0.058ReD (20) 式中，D=管的直徑，因層流之雷諾數必頇小於2000，故L=0.058ReD＜0.058×2000D=116D， 亦即流體從入口至拍攝的長度頇為大於100倍管徑的長度，才可獲得完全發展流。 以下的結果為一雷諾數為1500，圓直管入口速度為6.84 cm/s，圖26為圓直管經過 PC-MRI拍攝後的檔案，利用Matlab® 程式搜尋找出的邊界，其中綠色線為B-Spline曲線，紅 色線則為利用速度分佈之曲率變化所判別出的近似邊界。圖27則為壁面上剪應力的分布情 形，因其中心速度Vc=13.68 cm/s，故可以式(21)求出理論上之壁面剪應力，其大小為0.02 kpa。由圖中可知，以法線距離(dr)為2 Pixel 時，所求出的壁面剪應力值與理論值最接近(如 表2所示)。即受到MRI雜訊的影響較小；可獲的較精確之壁面剪應力。 (21)
接下來的部份為模擬暫態流場以及流固耦合等邊界條件皆加入之模型，圖 28 展示數 值計算之主動脈幾何尺寸、座標系統以及網格分佈。計算區域的總格點數為 290750，其中 流體的部份佔 234358，固體的部份佔 56372，最密的格點分布於流體與固體的交界面上以 便精確的求得壁面的剪應力值，其間隔為 5 µm。整個模型分為六個部份(A-F)，分別為升胸 主動脈、頭臂動脈、左頸動脈、左鎖骨下動脈、主動脈弓以及降胸主動脈，其壁面厚度分 別為 1.63, 0.86, 0.63, 0.67, 1.63 and 1.63 mm [35]，由於模型建構為多區塊結構性網格，為了 提高格點的品質，在升胸主動脈入口端、頭臂動脈、左頸動脈、左鎖骨下動脈以及降胸主 動脈出口端截面，及各個歧管與主動脈弓的交界面上，均採用橢圓平滑法 ESM (Elliptical Smoothing Method)以維持複雜幾何形狀的格點正交化。 圖 29 為考慮流固耦合情況下的胸主動脈矢狀截面的軸向速度分布以及流線圖，此一結 果與先前的一些文獻有許多相似之處[16, 37]，在比對了考慮流固耦合的模型，以及未考慮 流固耦合的模型的軸向速度以及流線，我們發現兩者之間的差異在大部分的時間點都不 大。另外在收縮峰值(Peak Systole, t/T=1)、減速度相位(Decelerating Phase, t/T=1.18)以及加 速度相位(Accelerating Phase, t/T=1.82)時，由升胸主動脈入口端流入的流體會流向主動脈弓 的 Inferior 方向的壁面，形成這樣的現象是因為升胸主動脈到主動脈弓間為一彎管的構造， 流體撞擊到 Anterior 方向的壁面後會由角向往 Inferior 方向的壁面聚集，造成此處有較高的 軸向速度以及較多的流體，除此之外，我們發現有考慮流固耦合的模型在這幾個相位中， 在 Inferior 方向主動脈弓處的流體有較長的分離區域，這樣的現象對此位置的壁面剪應力 (WSS)亦會有所影響，關於這個部份在後面的章節會再提及。另外，由於在減速度相位 (Dccelerating Phase, t/T=1.18)時，主動脈弓三根分支的逆流造成了明顯的渦流出現在分支與 主動脈弓交界處的 Anterior 方向壁面附近。當流體進入降胸主動脈時，會偏向 Posterior 方 向的降胸主動脈壁面聚集。比較考慮流固耦合以及剛體模型的結果，考慮流固耦合的模型 計算出來的結果在升胸主動脈入口處有較明顯的旋轉流，另外在頭臂動脈、左頸動脈以及 左鎖骨下動脈與主動脈弓交界處的渦流也較明顯，以上的結果跟文獻 Kilner et al. [38] 以及 Shahcheraghi et al. [14] 的結果相似。 近年來電腦模擬已經可以被運用在探討主動脈的血流動力因素上，主動脈的曲率、扭 轉、三根分支的存在以及截面面積的改變都會影響胸主動脈內流體流動的情形。圖 30 為考 慮流固耦合情況下的胸主動脈特定截面軸向速度分布曲線隨時間的變化，圖中的α 定義為 順時鐘方向的旋轉角，當其為方向法向量為水平時 α = 0°，在截面 a、b、c、d、e 和 f 的 α 值分別為 0°、0°、0°、0°、90°和 180°，模擬的結果顯示在截面 a、e 和 f 的流體會逐漸的 發展成 M 型的高流速區域，這樣的情形在收縮峰值(Peak Systole, t/T=1)、減速度相位 (Decelerating Phase, t/T=1.18)以及舒張末期(Later Diastole, t/T=1.70)較為明顯，除此之外我們 觀察到截面 b、c 和 d 靠近分支的入口的 Posterior 方向處的小區域在減速度相位(Decelerating Phase, t/T=1.18)會出現逆流的現象，在舒張末期(Later Diastole, t/T=1.70)時，升胸主動脈入 口附近的軸向速度分佈相當均勻，另外在主動脈的三根分支處在此相位下有均勻的負流量 出現，此主動脈三根分支的負流將在在主動脈弓處與由升胸主動脈來的主血流結合後流向 降胸主動脈。
圖 29 和圖 30 所呈現的結果和文獻中其他的數值模擬的結果有高度的相似[14, 16] ， 另外和文獻中以雷射都普勒測速儀所量出的實驗結果也有高度相似 [39]，Shahcheraghi et al.[14] 和 Liepsch et al.[39] 等人指出在減速度相位時，在主動脈弓處流體會發展成為 M 型 的高流速區域，Shahcheraghi et al.[14] 等人在設定主動脈弓三根分支的邊界條件時，假定
每個分支佔總流量的 5%，其結果證明了流過主動脈弓處的流體的扭曲以及因為主動脈弓曲 率而使得流體有部份區域流速變大等現象，Jin et al.[16] 等人則是利用 MRI 和速度繪圖結 合，製出電腦數值模擬的模型，用來檢視主動脈弓的曲率以及壁面的順應性和移動量對流 體的影響，在此模擬的設定中由 MRI 所找出的幾何形狀以及速度分佈做為模擬的邊界條 件，其中在升胸主動脈處亦考慮了徑向的主動脹縮以及由於心臟跳動所帶動的主動脈壁面 移動量，除了升胸主動脈以外的其他部分(接近主動脈弓的入口處開始一直延伸到三根分支 以及降胸主動脈)則是假設為剛體且固定在空間中，其模擬的結果顯現出在接近舒張末期時 升胸主動脈的入口處流體有旋轉的現象，我們的結果中亦有這樣的現象。由 Jin et al.[16] 等 人的研究假設可以推測，在一些分支以及剛體和彈性體結合處，應該會產生反射波，此反 射波應會對流場造成相當的影響。 壁面剪應力為一個重要的血流動力因素，且其與主動脈剝離以及粥樣動脈硬化斑塊的 形成息息相關，在此壁面剪應力的定義為合成剪應力(即以三方向的速度向量合成的總速率 來計算剪應力值)，圖 31 為考慮流固耦合之胸主動脈流體-固體交界壁面剪應力分佈圖，我 們發現在收縮峰值(Peak Systole, t/T=1)時，存在於主動脈弓處高度扭曲的流體在三根分支 Posterior 方向的壁面產生了最大的壁面剪應力，其形成的原因可能為流體從升胸主動脈進 入主動脈弓時，流體的軸向速度直接衝擊到主動脈弓三根分支 Posterior 方向的壁面，其值 大約為 9.53 N/m2，為整個週期中的最大值。壁面剪應力次高的區域產生於降胸主動脈 Posterior 方向的壁面以及主動脈弓 Inferior 方向的壁面上，其形成的原因為流體經過主動脈 弓後流向降胸主動脈時，會偏向 Posterior 方向的壁面聚集，造成此處的壁面剪應力較高， 而主動脈弓 Inferior 方向則是有較明顯的二次流現象，造成較高的壁面剪應力，上述壁面剪 應力的分布的趨勢在收縮峰值(Peak Systole, t/T=1)、減速度相位(Decelerating Phase, t/T=1.18) 和加速度相位(Accelerating Phase, t/T=1.82)等時間點的差異性不大，而在舒張末期(Later Diastole, t/T=1.70)時，主動脈弓三根分支的逆流造成在此時間點下的最高壁面剪應力存在於
主動脈弓 Inferior 方向的壁面，但由於流量很小，所以此處的壁面剪應力(~0.71 N/m2
)相較 於其他時間點下的結果低了許多。整體而言，考慮流固耦合所預測出的壁面剪應力分佈與 剛體的結果[24] 在趨勢上是相似的，不過考慮流固耦合的結果顯示出高剪應力區域的分布 在降胸主動脈的 Left 和 Right side 方向比剛體的分布結果廣，形成此現象的原因為考慮流 固耦合的結果在 Inferior 方向主動脈弓處的流體有較長的分離區域所造成。 圖 32 為彈性體與剛體胸主動脈特定位置壁面剪應力隨時間變化之比較圖，此圖共針對 胸主動脈內四個特定位置 a、b、c 和 d 分別做壁面剪應力隨時間變化的探討，圖中虛線為 剛體的結果，實線為考慮流固耦合之彈性體的結果，其取樣時間為 0.02 秒，總取樣的時間 為 1 秒，所選擇的 4 個點皆為主動脈剝離以及粥樣動脈硬化的好發位置，從此圖我們觀察 到 4 個點的壁面剪應力隨時間分布趨勢大致上與升胸主動脈入口的流量波形相似，從收縮 峰值(Peak Systole, t/T=1)開始，壁面剪應力從最大值開始下降，到了 1.3-1.7 之間壁面剪應 力降低到了較低的值，而在 1.7-2 中因為進入了加速度相位(Accelerating Phase, t/T=1.82)， 所以壁面剪應力在此急速的上昇，此現象亦與升胸主動脈入口的流量波形相似。比較彈性 體與剛體之差異性，我們觀察到彈性體的壁面剪應力值在時間 1.0-1.3 以及 1.7-2 之間大致 上會略低於剛體的壁面剪應力值，相反地剛體的壁面剪應力值在時間 1.3-1.7 之間大致上會 略低於彈性體的壁面剪應力值，上述的結果與文獻中的數值模擬的結果有高度的相似[16]。 本研究中對人體胸主動脈壁面剪應力所做的探討如圖 31 以及圖 32 所示，數值模擬的結 果和文獻中的結果相似[14, 17]，由先前模擬的經驗得知出口(頭臂動脈、左頸動脈、左鎖骨
下動以及降胸主動脈)的質量流率的分佈會影響壁面剪應力的大小，所以胸主動脈出口的質 量流率設定相當重要，Shahcheraghi et al.[14] 的研究中將整個週期內主動脈弓三根分支的 質量流率都分別設為總進口流量的 5%，本研究將 Shahcheraghi et al.[14] 對邊界條件的設定 再進一步的改善，主動脈弓三根分支的質量流率是由體外胸主動脈實驗進行量測後的結 果，來作為模擬的邊界條件，而在 Nakamura et al.[17] 的研究中，利用在升胸主動脈前段 加上心室的部份，藉由模擬心室的變形將血液從心臟打出，進而探討舒張和收縮血流以及 心室流對胸主動脈的影響，不過其模型並未包含主動脈弓的三根分支，其研究最主要的部 份為比較在不同入口邊界條件下，模型升胸主動脈的壁面剪應力分佈，而本研究將 Nakamura et al.[17] 的研究改善再加上了主動脈弓的三根分支，來探討其對血流動力因素壁 面剪應力的影響。 圖 33 為考慮流固耦合情況下的之胸主動脈流體-固體交界壁面壓力分佈圖，圖中的壓力 為經過無因次化後的結果，其無因次化的方程式為 P = P* / Pmax，符號 P、P*以及 Pmax 分別 代表無因次化壓力、原始壓力以及原始最大壓力值，在整個週期中壓力最大值發生在胸主 動脈弓處三根分支 Posterior 方向的壁面，與壁面剪應力最大值的位置幾乎一致，推測原因 亦是流體從升胸主動脈進入主動脈弓時，軸向的流體直接衝擊到主動脈弓三根分支 Posterior 方向的壁面，次高的壓力則是出現在於降胸主動脈 Posterior 方向的壁面以及主動 脈弓 Inferior 方向的壁面，與壁面剪應力次高處的位置也幾乎一致，除此之外，最小的壓力 則是發生在減速度相位 Decelerating Phase (t/T=1.18)時的升胸主動脈入口。
為了預估隨時間變化的壁面剪應力震盪值(Oscillatory WSS Index, OSI)，本研究利用計
算過後的剪應力值，利用繪圖軟體 Tecplot®針對每個時間點下的剪應力進行後處理，獲得 壁面剪應力震盪值，其定義為 OSI=
T T dt dt 0 0 1 2 1
，其中 T 為整個心跳週期 (Taylor et al.[40]). ， 而 其 中 的 剪 應 力 τ 是 沿 著 胸 主 動 脈 壁 面 的 值 ， 其 表 示 式 為 τ=
y x v y n n x v ，符號 n 代表胸主動脈內垂直曲面的流體速度方向，μ 代表模擬流 體的動黏滯係數。圖 34 為考慮流固耦合之胸主動脈流體-固體交界壁面剪應力震盪值分佈 圖，結果顯示高壁面剪應力震盪值出現於主動脈弓，而其最大值出現在主動脈弓三根分支 的交界處的 Anterior 方向壁面，一般而言存在低壁面剪應力而有高震盪壁面剪應力的區域 較易產生壁面內膜的增厚以及斑塊增生。 如同之前所提及的，血流動力因素和許多與血管相關的疾病息息相關，雖然這些疾病 和血流動力因素間的生化學機制有許多部分尚未釐清，Malek et al.[41] 假設胸主動脈瘤的 發生是由於血管舒張以及動脈壁面的組織在承受高剪應力下，造成內皮細胞氧化氮生成累 積的影響而產生改變，先前的病理學研究如 Roberts[42] 和 Svensson et al.[43] 指出壁面剪 應力和主動脈剝離有相當大的關連性，本研究的結果顯示壁面剪應力的分佈 (圖 31)和 Debakey et al.[1] 主動脈剝離好發位置的結果有非常高的相關性，其研究中的數種主動脈剝 離的分類，極具有醫學上的代表性。另一方面，血管截面存在低壁面剪應力以及高壁面剪 應力震盪值為粥樣動脈硬化的高危險群[44]，粥樣動脈硬化好發的區域大致上在一些相當好 辨識的位置，如彎曲處、分岔處以及血管的分支，而本研究的壁面剪應力震盪值分佈(圖 34) 在彎曲處、分岔處以及血管的分支處都較高，和粥樣動脈好發的位置吻合。胸主動脈的脈波傳導速度為一個重要的臨床指標，通常用 PC-MRI 技術來擷取所需的 資料，應用在判斷主動脈順應性以及一些心血管方面的疾病，為一個有效率的指標，本研 究率先嘗詴用數值計算流固耦合之胸主動脈的結果與計算脈波傳導速度的方法做結合，將 計算臨床人體脈波傳導速度的方法應用在計算本主動脈數值研究的脈波傳導速度，如此預 估我們模型的主動脈順應性，並且可以將之與文獻中之臨床人體的結果做比對，進一步驗 證我們數值模型以及計算法則的正確性。 計算脈波傳導速度的方法有許多種，在先前的研究[25]中提到一個較新的方法，以 PC-MRI 量得的降胸主動脈軸心的軸向速度來預估降胸主動脈的脈波傳導速度，此方法將降 胸主動脈的脈波傳遞隨時間可視化，並且加上一些假設獲得脈波傳導速度。圖 35 顯示胸主 動脈矢狀截面的 B-spline 曲線(紅色邊線)和中心軸(藍色中心線)，左圖是實際人體的 PC-MRI 影像圖，右圖則是數值計算模型示意圖，邊線代表所選擇的降胸主動脈壁面邊界，而圖中 的藍線代表的是選擇的降胸主動脈的中心線，當我們在計算脈波傳導速度時，軸向速度皆 沿著降胸主動脈的中心線擷取，在本研究所選擇的降胸主動脈中心線長度約在 25 公分左 右，而影像的時間解析為 0.02 秒。 圖 36 代表的是三維脈波傳導速度分佈示意圖，圖中的 x、y 以及 z 軸分別代表在脈動 週期經過的時間、降胸主動脈的距離(圖 35 的藍線)，以及軸向速度，其中降胸主動脈距離 的起點為主動脈弓接降胸主動脈的位置，接近腹動脈的那一側為尾端，圖 36 中將三維脈波 傳導的過程分為三個階段，第一個為心臟預收縮階段(圖中灰色區範圍)，其定義為小於最大 軸向速度 10%的範圍，其平面方程式為 S1；而心臟收縮早期(圖中黑色範圍)的定義則為軸 向速度上升區段內，最大軸向速度的 10~50%，其平面方程式為 S2，平面方程式 S1以及 S2 如下式(22)以及(23)所示: S1 = a1 × x(space) ＋ b1 × y(time) ＋ c1. (22) S2 = a2 × x(space) ＋ b2 × y(time) ＋ c2. (23) 而此兩平面的相交線的斜率即為 PWV。故將此兩平面的法向量 N1(a1,b1,-1)與 N2(a2,b2,-1)做 外積，可得到此相交線的方向向量 I(b2-b1,a1-a2,a1b2-b1a2)，所以脈波傳導速度即為 PWV = (b2 - b1) / (a1 - a2). (24) 在本研究中以上計算脈波傳導速度的過程皆擷取計算流體力學軟體 ACE+® 所算出的軸向 速度結果，再運用數學軟體 Matlab®編寫的算則程式計算出結果。 圖 37 為模擬之降胸主動脈內軸向速度(z 軸)隨空間(x 軸)和時間(y 軸)變化圖(楊氏係數 = 840kpa、包松比=0.45)，利用方程式(24)所算出的脈波傳導速度為 552 公分/秒，Yu et al.[25] 利用 PC-MRI 量測到臨床人體胸主動脈的脈波傳導速度，其中一個測詴群組是 25 歲的健康 人，此族群計算出來的脈波傳導速度介於 500-600 公分/秒，由上述的結果顯示，本研究考 慮流固耦合之胸主動脈模擬的結果和 Yu et al.[25] 使用 PC-MRI 量測再計算出的結果有良 好的一致性，這也可以間接證明本研究模型的胸主動脈壁面順應性和臨床人體的胸主動脈 壁面順應性相似，這個創新的方法用來計算數值模型的脈波傳導速度，並且可以用來估計 臨床個別人體胸主動脈壁面的結構力學性質，可以用來減低以往在計算人體胸主動脈壁面 的難度以及不確定性，並可以更簡單的針對個別人體的脈波傳導速度來製出對應楊氏係數 的數值模型，更增加模擬的正確性。
計畫成果自評
本研究執行至此共發表了四篇學術期刊:
1.―In Vitro Characterization of Aortic Flow Using Numerical Simulation, PC-MRI and Particle Tracking Images,‖ Proceedings of the Institution of Mechanical Engineers, Part C, Journal of
Mechanical Engineering Science, Vol. 222, No.12, pp 2455-2462, 2008. (SCI Journal).
2.―Investigation of Pulsatile Flowfield in Healthy Thoracic Aorta Models,‖ Annals of Biomedical
Engineering, Vol. 38 (2), pp. 391–402, Feb., 2010. (SCI Journal).
3.―Clamping Analysis of Aorta for Open Heart Surgery,‖ Journal of Mechanics, accepted (2010). (SCI Journal).
4.―An Axial Velocity Profile Method to Derive the Aortic Pulse-Wave Velocity in Diabetic Patients,‖ Journal of Mechanics, accepted, Aug., (2010). (SCI Journal).
還有另一篇學術期刊已投稿，正在審核中
1.―Investigation of Pulsatile Flowfield and Pulse Wave Velocity in a Healthy Thoracic Aorta Model with Considering Fluid-structure Interaction,‖ Submitted to Annals of Biomedical
圖、表附錄
圖 1 Debakey et al. [1]之主動脈剝離分類
圖 2 三維軸向速度分布圖[25]。(I)下行主動脈軸向速度(Z 座標)在空間(X 座標)
與時間(Y 座標)的分布圖。(II) X-Z 平面：軸向速度在不同時間點的空間分佈
圖。(III)X-Y 平面：相同速度輪廓的傳導圖。(IV)Y-Z 平面：軸向速度在不同
空間的時間分佈圖。
圖3 主動脈實驗設備與管路系統示意圖
圖5 體外胸主動脈之MRI影像與B-Spline曲線
圖 7 心臟流速脈動波形
圖8 主動脈管路系統實體圖
圖10 B-spline曲線(穩態) 圖11 速度分布分布圖(穩態)
圖13 理論模擬與MRI量測結果之X-Z平面上絕對速度分布圖
圖15 主動脈之全域瞬時流場 圖16 主動脈各處之橫截
圖17 主動脈壁面總壓力分布圖(a)矢狀切面之右側壁面(b)矢狀切面之左側壁
面
圖18 主動脈壁面應變率分布圖(a)矢狀切面之右側壁面(b)矢狀切面之左側壁
面
圖19 暫態主動脈網格圖
圖21 流場矢狀截面的流線圖
圖23 主動脈壁面剪應力分布圖(暫態)
圖25 圓直管實驗平台
(a)
(b)
圖26 圓直管邊界判讀結果(a)左側(b)右側
(a)
(b)
圖 28 胸主動脈網格圖
圖 30 考慮流固耦合情況下的胸主動脈特定截面軸向速度分布曲線隨時間的
變化
圖 34 考慮流固耦合情況下的胸主動脈流體-固體交界壁面剪應力震盪值分佈
圖
圖 35 胸主動脈矢狀截面的 B-spline 曲線(紅色邊線)和中心軸(藍色中心線)。左
圖是實際人體的 PC-MRI 影像圖，右圖則是數值計算模型示意圖。
圖 36 三維脈波傳導速度分佈示意圖
表1 MRI影像設定參數
表2 壁面剪應力與理論值比較
P Piixxeel l WWSSSSMMeeaann((PPaa))((左左側側) ) SSttdd((左左側側) ) WWSSSSMMeeaann((PPaa))((右右側側) ) SSttdd((右右側側) ) 0 0..1 1 00..22887 7 00..2288556 6 00..1188887 7 00..11881122 0 0..5 5 00..00665511 00..0055449 9 00..0066000 0 00..00334400 1 1 00..00337744 00..0011990 0 00..0044008 8 00..00114444 1 1..5 5 00..00228888 00..0011119 9 00..0033112 2 00..00009 9 2 2 00..0022447 7 00..0000774 4 00..0022661 1 00..00110 0 理 理論論值 值 00..0022 00..002 2參考文獻
[1]. Debakey, M E., McCollum, C.H., Crawford, E.S., Howell, J., Noon G. P., Lawrie G. L., ―Dissection and Dissecting Aneurysms of the Aorta: Twenty-one Year Follow-up of 527 Patients
Treated Surgically.‖ Surg 92:1118-1134, 1982.
[2]. 劉清吉，人體主動脈流場之特性研究，私立大葉大學機械與自動化工程學系碩士論文， 民國九十三年。
[3]. 蔡坤晃，主動脈剝離現象之流場分析，私立大葉大學機械與自動化工程學系碩士論文， 民國九十四年。
[4]. 行政院衛生署網站：http://www.doh.gov.tw/statistic/index.htm ,民國九十七年。
[5]. Sahs, A L.,―Observations on the Pathology of Saccular Aneur- ysms.‖ J.Neurosurg. 24:79-806, 1996.
[6]. Erbel R., Delert H., Meyer J., ―Effect of Medical and Surgical Therapy on Aortic Dissection Evaluated by Transesophageal Echocardiography : Implications for Progonosis and Therapy.
Circulation.‖ 87:1604-1615, 1993.
[7]. Barakat, A.I., Karino, T.,Colton, C. K., ―Microcinematographic Studies of the Flow Field in the Excised Rabbit Aorta and its Major Branches, ‖ Journal of Biomechanics.31:217-228, 1997.
[8]. Keren, A., Kim, C B., Hu, B S., ―Accuracy of Biplane and Multiplane Transesophageal Echocardiography in Diagnosis of Typical Acute Aortic Dissection Intramural Hematoma.‖ J.
AmColl Cardiol. 28:627-636, 1996.
[9]. Perktold, K., Hofer, M., Rappitsch, G., Loew, M., Kuban, B. D., Friedman, M. H., ―Validated Computation of Physiologic Flow in a Realistic Coronary Artery Branch,‖ Journal of Biomechanics. 31: 217– 228, 1998.
[10]. Gijsen, F. J. H., Allanic, E., Vosse, F. N. van de., Janssen, J. D., ―The Influence of the Non-Newtonian Properties of Blood on the Flow in Large Arteries: Unsteady Flow in a 90° Curved Tube,‖ Journal of Biomechanics. 32:705-713, 1999.
[11]. Hugo G., Bogren, MD, PhD, and Michael H. Buonocore, MD, PhD, ―4D Magnetic Resonance Velocity Mapping of Blood Flow Patterns in the Aorta in Young vs. Elderly Normal Subjects.‖ Journal of Magnetic Resonance Imaging. 10:861-869, 1999.
[12]. Zabielski, L. and Mestel, A. J., ―Helical Flow Around Arterial Bendes for Varying Bond Mass.‖ Journal of Biomechanical Engineering. Vol. 122, April, 2000.
[13]. Zhao, S. Z., Xu, X. Y., Hughes, A. D., Thom, S. A., Stanton, A. V., Ariff, B., Long, Q., ―Blood flow and Vessel Mechanics in a Physiologically Realistic Model of a Human Carotid Arterial Bifurcation,‖ Journal of Biomechanics. 33:975-984, 2000.
[14]. Shahcheraghi, N., Dwyer, H. A., Cheer, A. Y., Barakat, A. I., Rutaganira, T., ―Unsteady and Three-Dimensional Simulation of Blood Flow in the Human Aortic Arch,‖ Journal of Biomechanical Engineering. Vol. 124, 378-387, 2002.
[15]. Buchanana, J. R., Kleinstreuera, C., Hyunb, S., Truskey, G. A., ―Hemodynamics Simulation and Identification of Susceptible Sites of Atherosclerotic Lesion Formation in a Model
Abdominal Aorta.‖ Journal of Biomechanics. 36:1185–1196, 2003.
[16]. Jin, S., Oshinski, J., Giddens, D. P., ―Effects of wall motion and compliance on flow patterns in the ascending aorta,‖ ASME J. Biomech. Eng 125:347-354, 2003.
[17]. Nakamura, M., Wada, S., Yamaguchi, T., ―Computational analysis of blood flow in an integrated model of the left ventricle and the aorta,‖ ASME J. Biomech. Eng 128:837-843, 2006. [18]. 丁大為，吳秉勳。植入枝架造成血管內壁呈皺摺狀變形後之血液動力及非牛頓流體效 應分析。第十一屆全國計算流體力學學術研討會。中華民國九十三年八月。 [19]. 丁大為，林柏宏。於人體生理條件下左冠狀動脈之幾何形狀變化對其壁面剪應力分布 之影響。行政院國家科學委員會專題研究計畫成果報告書。NSC91-2213-E-014-008。執 行期限 91/08/01~92/07/31。 [20]. 劉通敏，丁大為，陳禹銘。顱內彎形母管與其上側向動脈瘤之脈動流場特性數值模擬。 第十一屆全國計算流體力學學術研討會。中華民國九十三年八月。 [21]. 牛仰堯，伍邦銓，虞希禹，曾文毅，彭旭霞，李隆政，鄭守成。人體主動脈之磁振造 影良策與數值模擬。第十一屆全國計算流體力學學術研討會。中華民國九十三年八月。 [22]. 李明龍，周朝宜，李隆政，沈澄宇，施仁傑，林錫慶。冠狀動脈繞道血管之三維數值 模擬。第十一屆全國計算流體力學學術研討會。中華民國九十三年八月。
[23]. Yang, A. S., Wen, C. Y., Tseng, L. Y., ―In vitro characterization of aortic flow using numerical simulation, phase-contrast magnetic resonance imaging, and particle tracking images,‖ Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science 222(12):2455-2462, 2008.
[24]. Wen, C. Y., Yang, A. S., Tseng, L. Y., Chai, J. W., ―Investigation of pulsatile flowfield in healthy thoracic aorta models,‖ Annals of Biomedical Engineering 38:391-402, 2010.
[25]. Yu, H. Y., Peng, H. H., Wang, J. L., Wen, C. Y., Tseng, W. Y. I., ―Quantification of the Pulse Wave Velocity of the Descending Aorta Using Axial Velocity Profiles From Phase-Contrast Magnetic Resonance Imaging,‖ Magnetic Resonance in Medicine 56:876-883, 2006.
[26]. Francis, J. G. F., ―The QR Transformation I,‖ Comput. J., vol. 4:265-271, 1961. [27]. Whitmore, R. L., ―Rheology of the Circulation,‖ Oxford: Pergamon Press. 1968.
[28]. Yearwood, T. L., ―Steady and pulsatile flow analysis in a model of the human aortic arch,‖ Ph.D. thesis, Tulane University, New Orleans, 1979.
[29]. Fung, Y. C., ―Biomechanics. Mechanical Properties of Living Tissues,‖ New York: Springer, 1981.
[30]. Van-Doormaal, J. P., Raithby, G. D., ―Enhancements of the SIMPLE method for predicting incompressible fluid flows,‖ Numer. Heat Transfer 7:147-163, 1984.
[31]. Riley, W. A., Barnes, R. W., Evans, G. W., Burke, G. L., ―Ultrasonic measurements of the elastic modulus of the common carotid artery: the Atherosclerosis Risk in Communities (ARIC) Study,‖ Stroke 23:952-956, 1992.
[32]. Gao, F., Watanabe, M., Matsuzawa, T., ―Stress analysis in a layered aortic arch model under pulsatile blood flow,‖ Biomed. Eng. Online 5:25, 2006.
A., ―Fluid-structure interaction within realistic three-dimensional models of the aneurysmatic aorta as a guidance to assess the risk of rupture of the aneurysm,‖ Med. Eng. Phys 23:647–655, 2001.
[34]. Ferguson, G.G. Physical Factors in the Initiation, Growth, and Rupture of human Intracranial saccular aneurysms, J. Neurosurg. Vol.37,pp.666-677,1972.
[35]. Avolio, A. P., ―Multi-branched model of the human arterial system,‖ Med. & Biol. Eng. & Comput 18:709-718, 1980.
[36]. Nichols W., O’Rourke M., McDonald’s, ―Flow in Arterier, Lea & Febiger, Philadelphia,‖ London,1990.
[37]. Kim, Y. H., Kim, J. E., Ito, Y., Shih, A. M., Brott, B., Anayiotos., ―A Hemodynamic Analysis of a Compliant Femoral Artery Bifurcation Model using a Fluid Structure Interaction Framework,‖ Annals of Biomedical Engineering 36:1753-1763, 2008.
[38]. Kilner, P. J., Yang, G. Z., Mohiaddin, R. H., Firmin, D. N., Longmore, D. B., ―Helical and retrograde secondary flow Investigation of Pulsatile Flowfield in Healthy Thoracic Aorta Models patterns in the aortic arch studied by three-dimensional magnetic resonance velocity mapping,‖ Circulation 88:2235-2247, 1993.
[39]. Liepsch, D., Moravec, S. T., Baumgart, R., ―Some flow visualization and laser-Doppler velocity measurements in a tube-to-scale elastic model of a human arotic arch—a new model technique,‖ Biorheology 29:563-580, 1992.
[40]. Taylor, C. A., Cheng, C. P., Espinosa, L. A., Tang, B.T., Parker, D., Herfkens, R. J., ―In vivo quantification of blood flow and wall shear stress in the human abdominal aorta during lower limb exercise,‖ Ann. Biomed. Eng. 30:402-408, 2002.
[41]. Malek, A. M., Alper, S. L., Izumo, S., ―Hemodynamic shear stress and its role in atherosclerosis,‖ J. Am. Med. Assoc. 282:2035-2042, 1999.
[42]. Roberts, W. C., ―Aortic dissection: anatomy consequences and causes,‖ Am. Heart J. 101:195-214, 1981.
[43]. Svensson, L. G., Grawford, E. S., ―Aortic dissection and aortic aneurysm surgery: clinical observations, experimental investigations, and statistical analyses,‖ Part II. Curr. Probl. Surg. 29:913-1057, 1992.
[44]. Shaaban, A. M., Duerinckx, A. J., ―Wall shear stress and early atherosclerosis: a review,‖ American Journal of Roentgenology 174:1657-1665, 2000.
