基于CFD的左心室流場數(shù)值模擬研究現(xiàn)狀與趨勢

阿都建華 尹立雪 謝盛華

1(成都信息工程大學軟件工程學院，成都 610225)2(四川省醫(yī)學科學院·四川省人民醫(yī)院心血管超聲及心功能科超聲心臟電生理學與生物力學四川省重點實驗室，成都 610072)

基于CFD的左心室流場數(shù)值模擬研究現(xiàn)狀與趨勢

阿都建華1,2尹立雪2*謝盛華2

1(成都信息工程大學軟件工程學院，成都 610225)2(四川省醫(yī)學科學院·四川省人民醫(yī)院心血管超聲及心功能科超聲心臟電生理學與生物力學四川省重點實驗室，成都 610072)

隨著計算技術的快速發(fā)展，計算流體力學在心臟血流流場模擬中得到廣泛應用，取得了有價值的相關成果，在心臟功能的研究中扮演著越來越重要的角色。目前，運用計算流體力學實現(xiàn)心臟流體動力學特性及其流場數(shù)值模擬分析，已成為心血管基礎和臨床研究領域的一個重要研究方向，對心血管系統(tǒng)疾病預防、臨床診斷和治療具有重要意義。簡述基于CFD的左心室流場數(shù)值模擬方法的發(fā)展過程、研究現(xiàn)狀和方法分類，分析基于CFD的左心室流場數(shù)值模擬計算存在的主要問題，提出該研究領域未來的發(fā)展方向。

計算流體力學；左心室流場；數(shù)值模擬；流固耦合；

引言

心臟病學本質上是一門關于流的科學，心血管系統(tǒng)的主要功能是力學驅動、控制和維持整個人體系統(tǒng)中各個組織器官的血流灌注[1]。心室內(nèi)的血液流場更為復雜多變，心動周期內(nèi)不同時相血液流場主要表現(xiàn)為動態(tài)連續(xù)變化的層流和渦流，但在每個心動周期特定時相的流體總體特征通常又保持一定的相對穩(wěn)定性和可重復性，特定時間點的血流渦流參數(shù)、血流速度分布、壓力分布等流體特征量具有一定的不變性。不同的心臟病理解剖和病理生理異常，均會產(chǎn)生不同的異常血流流場。為了檢測心臟疾病導致的血流流場異常，掌握心臟內(nèi)血流變化過程及其時空分布就顯得非常重要[2]。

近年來，研究人員已將先進的醫(yī)學測量技術與計算機仿真技術、計算流體力學緊密結合起來，用于研究人體心臟的血液動力學狀態(tài)，通過計算機數(shù)值計算和圖像顯示的方法，在時間和空間上定量描述流場的數(shù)值解。常規(guī)心臟血流動力學實驗研究不僅風險高、花費大，而且也僅能獲取速度、壓力等有限的流動參數(shù)數(shù)據(jù)，而速度矢量、流線、壁面剪切應力等復雜的流場參數(shù)數(shù)據(jù)是不可能得到的。數(shù)值計算仿真技術成本低、效率高，不僅能有針對性地研究各種心臟幾何變形、力學及流固耦合因素對血流動力學的影響，而且可以完成常規(guī)實驗無法實現(xiàn)的生物力學多參數(shù)時空描述和可視化表達。當在體或離體實驗研究不可行時，數(shù)值仿真就成為唯一可行的血流動力學觀測手段。雖然僅僅依靠血流動力學的數(shù)值計算結果，還不能完全解釋心血管疾病相關的病理生理機制，觀測結果需要與臨床醫(yī)學實驗結果相結合，但是這些計算結果將彌補實驗研究的不足，有助于更好地揭示心血管疾病發(fā)生、發(fā)展的病理生理和病理解剖原因，從而可以對此類疾病的預防和治療提供有價值的參考信息[3]。各種心臟疾病往往會在心動周期的全部或某些特定時相呈現(xiàn)出心腔內(nèi)血液流場上的異常變化，因此深入研究心臟血液流場特征有可能為各種心血管疾病的早期精確診斷提供重要線索和依據(jù)。醫(yī)學臨床研究也已經(jīng)證實，精確的心臟血流流場及其血流動力學信息(如血流速度、剪切率、剪切應力等)對于心力衰竭、心肌病、動脈粥樣硬化、中風等眾多心血管疾病的早期預警和診斷具有重要意義，同時也能為心臟外科手術的風險評估、手術時機指導及術后治療提供重要的參考信息[4-9]。目前，運用計算流體力學實現(xiàn)心臟的流體動力學特性及其流場數(shù)值模擬分析已成為心血管研究領域的一個重要研究方向，基于CFD的左心室血流流場數(shù)值模擬研究在心臟功能的研究中扮演著越來越重要的角色，將在心臟功能狀態(tài)的評價及預測、治療決策、治療方案優(yōu)化、手術輔助等方面發(fā)揮重要作用。

1 基于CFD的左心室流場數(shù)值模擬分析方法類型

在已有研究中，心臟流場數(shù)值模擬方法的研究大多都以左心室內(nèi)流場的變化為研究對象，基于CFD的左心室流場數(shù)值模擬方法大致可以分為3種[10]：規(guī)定幾何結構法[11-23]、虛構流固耦合法(fluid structure interaction，F(xiàn)SI)[24-29]、真實流固耦合法[30-33]。這三種方法都具有良好的研究價值和臨床應用價值，其中真實FSI方法具有更為廣闊的應用前景，但同時也最為復雜。

1.1 規(guī)定幾何結構法

規(guī)定幾何結構法通常利用臨床采集的圖像進行幾何結構建模，以此模型作為規(guī)定幾何結構構建移動網(wǎng)格，并以心臟內(nèi)壁為邊界，將心臟內(nèi)壁的移動作為邊界條件，同時不考慮心臟的幾何結構對血流產(chǎn)生的反饋作用力，最后進行心室流場的數(shù)值模擬。

圖1 流場模式和左心室充盈期的等壓線[23]Fig.1 Flow patterns and isobars in the LV during filling[23]

Vierendeels等以二維對稱幾何結構對左心室舒張期的血流流場進行建模，通過任意拉格朗日歐拉公式求解移動網(wǎng)格，實現(xiàn)流場域的模擬[23]。在心臟等容擴張開始后，開始進行模擬計算。在心臟的等容擴張階段，假定血流流體在心室內(nèi)部是靜止的，并且速度為零，不存在壓力差。在這個階段，所有節(jié)點的壓力都將根據(jù)心室壁模型內(nèi)的楊氏模量變化進行計算。當心室壓力下降到低于心房壓力時，二尖瓣膜就會打開，此時心室內(nèi)壓力由心室壁的舒張和血流的流體動力學共同決定。在研究中，作為邊界條件的二尖瓣血流模式以及二尖瓣打開時的心房壓力都來自集中參數(shù)模型的獨立計算。模擬結果如圖1所示，分別展示了左心室充盈期4次時間步長的速度矢量和等壓線，分別是充盈早期速度峰值(E峰)、舒張期后期、心房收縮期速度峰值(A峰)和A峰減速期。其中，第一個早期充盈E峰由心室舒張產(chǎn)生，第二個峰值是心房壓縮產(chǎn)生的A峰，從圖中可以充分觀察到充盈期4個階段的漩渦形成過程和壓力分布。

Saber等基于CFD方法利用MRI數(shù)據(jù)實現(xiàn)左心室血流的數(shù)值模擬，通過采集一個心動周期內(nèi)的MRI數(shù)據(jù)建立一個3D時變心臟內(nèi)膜圖像模型,再使用CFD方法模擬左心室的流場[17]。在建模過程中，使用了CMRTOOLS工具實現(xiàn)縱軸和橫軸MRI切面的手工圖像描繪處理；在血流流體的模擬過程中，則使用有限體積法進行計算。該模型的模擬結果展示了三維左心室內(nèi)壁的收縮和舒張過程，并能很好地模擬流場的渦流等主要特征。心室血流流場結果都基本與活體臨床實驗一致。由于沒有考慮瓣膜的影響，其流場速度模擬的準確性有一定局限性，圖2展示了該模型模擬的橫切面速度矢量分布，其中1～4是收縮期階段速度矢量分布，5～6是舒張期階段速度矢量分布。

圖2 橫切面速度向量分布[17]Fig.2 Velocity in short-axis plane[17]

2003年，Saber等又完善了以前的研究內(nèi)容，不但改善了縱向和橫向切面圖像的分割方法，而且在CFD模型中增加了近端左心房和升主動脈[18]。由于MRI無法提供瓣膜的清晰信息，因此研究中僅對二尖瓣和主動脈瓣做了近似模擬，無法對瓣膜口的精確位置、瓣膜的具體形態(tài)等細節(jié)問題做更為準確的建模和模擬，圖3為該模型在4個舒張期階段的橫切面速度矢量分布，能明顯看出形成了一個主要的大型前漩渦和一個后部小漩渦。

圖3 橫切面速度向量分布[18]Fig.3 Vectors of simulated diastolic velocities in a short axis plane[18]

Domenichini等建立了一個三維左心室模型，對舒張期的血流流體動力學進行研究[34]。該三維左心室模型是一個長橢圓形，流體動力來源于簡化的波形輸入血流，物理和數(shù)學模型如圖4所示。該模型的流動方程采用貼體坐標系，并在方位角方向用傅里葉函數(shù)展開，最終方程使用混合有限譜差分技術進行數(shù)值化求解。在研究中，為了理解心室擴展期的主要血流變化現(xiàn)象，通過變化控制參數(shù)來實現(xiàn)血流動力學的分析。模擬結果顯示，具有典型渦流結構特征的血流流場在控制參數(shù)發(fā)生變化時大致保持一致，直到用較低的斯特魯哈爾數(shù)，血流流場才變化為弱湍流的形式。圖5中展示了在特定常量參數(shù)時渦流隨著時間的變化在一個對稱平面中的演變過程，其中的(a)～(c)分別為3個不同時間步長時的y平面渦流矢量分布。

圖4 物理和數(shù)學模型輪廓[34]Fig.4 Sketch of the physical problem and of the mathematical model[34]

圖5 對稱平面的wy分布[34](渦流矢量的y分量, 斯特魯哈爾數(shù)=0.072，斯托克斯數(shù)=144，ε=0.125)。(a) t=20/128, (b) t=10/128,(c)t=50/128Fig.5 Distribution of wy[34] (y component of vorticity vector) on the symmetry plane, Strouhal number=0.072, Stokes number=144, ε=0.125). (a) t=20/128, (b) t=10/128,(c)t=50/128

TORSTEN SCHENKEL等提出的方法[22]和Saber等的方法[17-18]類似，但其建立了更為準確的瓣膜模型。該方法使用了Karlsruhe心臟模型，并按標準的臨床MRI圖像采集流程采集數(shù)據(jù)，通過一個二維時變注入孔口描述心臟瓣膜，同時根據(jù)采集的縱軸圖像估計瓣膜的位置和大小。模型中通過一個半自動的網(wǎng)格生成器創(chuàng)建計算網(wǎng)格，使用實現(xiàn)了任意拉格朗日歐拉公式求解的納維-斯托克斯(Navier-Stokes)解算器完成計算過程。由于核磁共振掃描的時間分辨率無法滿足CFD計算的需要，因此模型中采用了一個三階貝賽爾曲線近似方案來實現(xiàn)平滑的邊界和網(wǎng)格移動。模擬結果顯示，在心臟舒張期間，通過二尖瓣的血流形成一個不對稱的噴射，導致初始渦環(huán)不對稱發(fā)展。這些流場特征不僅和活體實驗結果完全一致，而且對注入血流邊界條件及心房幾何形狀有極高的敏感性，如圖6所示。其中，A部分表示速度概況、矢量及流量的測量，B部分表示曲線狀心房時用不同彩色線條描繪的血流流線速度，C部分是直線狀心房時用不同彩色線條描繪的血流流線速度。

圖6 注入條件變化血流模擬比較[22]Fig.6 Sensitivity to inflow conditions[22]

SIMON J. SONNTAG等將計算流體力學模擬研究方法和實驗研究方法結合起來，量化研究心臟二尖瓣回流的影響[35]。該研究開發(fā)了一套模擬左心室的循環(huán)流動腔室實驗研究系統(tǒng)，系統(tǒng)中采用3個不同的孔口模擬二尖瓣回流的變化，并同時利用二維多普勒超聲和PIV測量記錄流場域。同時，采用與實驗系統(tǒng)完全相同的幾何結構和邊界條件利用CFD進行模擬。PIV實驗結果和CFD計算模擬結果展現(xiàn)了高度的一致性，模擬的CDI也表現(xiàn)出與彩色多普勒圖像一樣的特性。在該研究中，數(shù)值模擬方法結果成功驗證了實驗研究方法的結果，并展示了實驗方法和計算模擬方法相結合帶來的靈活性。研究中使用的CFD模型如圖7所示，基于圓孔和矩形孔的時均PIV和CFD的模擬結果比較如圖8所示。

圖7 CFD的CAD模型、邊界條件和使用網(wǎng)格[35]Fig.7 CAD model, boundary conditions and used mesh (shown in one plane) of the CFD simulation[35]

圖8 基于圓孔和矩形孔的時均PIV和CFD的模擬結果比較[35]Fig.8 Comparison of time-averaged PIV and CFD for the circular orifice[35]

1.2 虛構FSI方法

虛構FSI法是一種模擬心臟血流流場的宏觀近似方法，利用浸入邊界法適用于模擬形狀復雜的彈性結構邊界的特點，用浸入邊界法將心臟結構簡化為一個彈性纖維體，并利用插值狄拉克函數(shù)來描述血流和心肌組織之間的相互作用。

圖9 模擬計算結果。(a)早期心室充盈階段；(b)心室射血階段；(c)通過二尖瓣膜的渦流運動；(d)右心室大渦流[27]Fig.9 Results from a computation.(a) A time early in ventricular filling, (b) A time during ventricular ejection, (c) Shows the vortex motion behind the mitral valve leaflets, (d) Shows a large vortex in the right ventricle[27]

浸入邊界模擬方法不僅能計算流體的運動，而且還能同時計算沉浸在流體之中又和流體發(fā)生相互作用的彈性邊界的運動[27]。在浸入邊界法中，血流被看做是牛頓液體，并通過歐拉速度表示，壓力分布在一個規(guī)則的三維計算格子中，血液動力則通過Navier-Stokes方程計算。在模擬計算過程中，不斷地進行求解計算并替換由外力導致變形的彈性邊界。彈性結構所描述的邊界在計算網(wǎng)格空間中能連續(xù)自由移動，可以通過狄拉克函數(shù)計算這個過程，同時這種移動變形產(chǎn)生的力量會通過彈性邊界傳遞給周邊的相鄰計算格子，這就表現(xiàn)為流體在外界動力的作用下運動，此時再次應用數(shù)值化狄拉克函數(shù)計算新的計算格子速度，然后通過計算格子的移動速度又來確定新的移動邊界，最后將邊界移動到新的位置。隨著這些計算過程的迭代執(zhí)行，新的網(wǎng)格邊界位置隨每個時間步長不斷重新定位，流體則表現(xiàn)為不斷的運動。

四大工程:給長城貼瓷磚、給赤道鑲金邊、給太平洋裝欄桿、給喜馬拉雅山安電梯間;四小工程:給蒼蠅戴手套、給蚊子戴口罩、給蟑螂戴避孕套、給老鼠戴腳鐐。

在文獻[27]中，僅要求提供血流流體的物理屬性、邊界的彈性屬性和邊界的初始幾何結構。在已經(jīng)建立的模型基礎上，采用GL(Graphics Library,GL)庫開發(fā)了一個可視化的自定義展示系統(tǒng)，并在SGI(Silicon Graphics,SGI)工作站上運行，用以顯示模擬計算的效果，如圖9所示。在這個系統(tǒng)中，通過計算流體標記的運動來表示心臟內(nèi)的血流運動，同時系統(tǒng)的控制管理非常靈活，可以自由控制計算過程，比如打開或關閉心臟模型中的任意一組纖維，并能隨意改變觀察的視角和放大興趣區(qū)域。

雷諾數(shù)是表征方程中遷移慣性項與黏性項比值大小的量，在它的不同范圍含有完全不同的兩種流動——層流和湍流。Lemmon等開發(fā)了一個薄壁心室計算模型，結合生理雷諾數(shù)模擬血流和心肌組織交互過程，該模型能有效減少使用的計算資源，使模擬心室運動的數(shù)值計算可以脫離超級計算機，在普通的用戶工作站上運行[24]。為了反映流體動力以及組織病理生理學上的變化，模擬過程中心肌組織和血流的相互作用使用浸入邊界法計算，并通過壓力耦合方程的半隱式方法來求解流體質量和動量守恒方程。在該模型中，并沒有包含一個實際的二尖瓣膜來模擬控制孔口的開閉，而是在二尖瓣膜的孔口位置直接假設一個閥門進行模擬。該研究將計算模擬結果和臨床數(shù)據(jù)進行了比較，顯示出二尖瓣流線、心室流場和心房流場均與臨床數(shù)據(jù)一致，心房流場的變化也和相關理論一致。此外，模型對心房和心室內(nèi)壓力場的跟蹤也顯示出了和活體測量一樣的結果。由此可以看出，該模型能檢查出由于心臟病變引起的舒張功能變化。

圖10 網(wǎng)格分區(qū)以及背景與重疊網(wǎng)格之間的通信[37]Fig.10 Mesh partitions and communications between background and OGs[37]

Lemmon等針對心臟舒張功能已減弱的病例進行建模[25]，研究中應用的3個舒張期功能障礙的研究病例和臨床觀察的舒張期流場變化一致。研究結果表明，心臟舒張延遲會減少心臟早期的充盈量，如果心室的僵硬程度不斷加強，充盈量就會持續(xù)減少。同時，在心房收縮期增強心房壓縮力度會導致更高的后期充盈速度和心房壓力。該研究結果顯示，心臟功能障礙狀態(tài)可以通過改變心肌纖維靜止狀態(tài)的長度和僵硬度的關聯(lián)關系來建模，這大大增強了研究人員未來對心臟疾病建模的信心。

由于血流和心肌組織的物質差異太大，瓣膜的數(shù)據(jù)信息又難以采集，致使生理條件下瓣膜運動的數(shù)值分析難度很大，因此在心臟數(shù)值建模中，主動脈瓣膜位置的流固耦合常常都被忽略。2003年，Hart等研究了心臟收縮過程中主動脈瓣膜位置的流固耦合效應[36]。拉格朗日乘子虛擬域方法作為一種數(shù)值模擬方法，常被用以描述計算流體領域的大瓣膜運動。在研究中，將該方法應用于一個主動脈瓣膜支架的三維有限元模型，然后為模型提供瓣膜及流經(jīng)血流的機械行為數(shù)據(jù)。模擬計算結果顯示，在心臟收縮期間，植入支架的瓣膜按流體運動方式隨血流移動。

2007年，Tai等提出了浸入對象法，實現(xiàn)包含流固耦合的三維非定常流模擬[37]，該方法中聯(lián)合應用了一個自由矩陣隱式雙時間步長和有限體積法的并行非結構多網(wǎng)格Navier-Stokes解算器，使用流固耦合方法來詳細研究三維非定常流的流動變化，并基于此研究人工心臟瓣膜的打開過程。在浸入對象方法在動量守恒方程的人工壓縮子循環(huán)的過程中，引入了一個表示人體動力能量的F來獲取速度分布。在流場域中，浸入了對象網(wǎng)格來定義對象的邊界，背景位置和重疊網(wǎng)格的通信處理過程如圖10所示。該方法的優(yōu)點是幾乎任意形狀都可以加入模型而無需重構網(wǎng)格，因為重構網(wǎng)格將消耗大量的時間和計算資源。這樣就能在生理條件下，利用人工心臟瓣膜完成復雜詳細的非定常血流和血流瓣膜流固耦合的三維模擬。實驗結果已經(jīng)驗證，這種并行隱式非結構多網(wǎng)格方法在計算三維非常定流的過程中具有高效性和精確性的特點。

1.3 真實FSI方法

真實FSI法分別利用有限元方法和計算流體力學方法對心臟系統(tǒng)的組織結構和血流進行建模，再通過流固耦合算法描述組織結構和血流之間的相互作用。一個成功的真實FSI心臟血流模擬方法最少需要合理地解決以下3個方面的問題。首先，表示心臟幾何結構的解算機必須能描述非線性、各向異性和不均勻的心臟壁組織結構特征；其次，CFD代碼必須能處理流場域的大變形，并且同時更新相應的計算網(wǎng)格；再次，流固耦合算法必須能確保固體結構和血流的正確耦合。因此，真實FSI方法的描述方法最為復雜，需要集成心肌力學、非線性有限元、計算流體力學、計算網(wǎng)格動力學、流固耦合算法來實現(xiàn)心臟流場的模擬。當前，這些方法的研究已經(jīng)取得了重要的進展，F(xiàn)ung等[38]、Lin等[39]提出了心臟壁肌肉結構的本購方程，Hunter等[40]、Smith NP等[41]提出了分析心室力學的有限元模型，Navier-Stokes方程的任意拉格朗日歐拉公式已經(jīng)成功地在計算流體力學移動網(wǎng)格和FSI中得到了應用[42-43]，Zhang Hou等[44]提出了同步交互的流固耦合算法。

Hunter等結合心室解剖、心肌組織的幾何結構和材料屬性，細胞膜離子通道，鈣調控和心臟肌細胞的肌絲力學來解釋心臟的集成功能，建立了分析心室力學的有限元模型[40]。他們在模型中專門設計了計算框架，用來建立心臟細胞和組織的結構及功能與整個心臟整體集成功能運轉之間的復雜關聯(lián)聯(lián)系。然而，如果希望該模型能夠洞察理解更多的心臟疾病過程，則還需要包含更多的生理功能，包括代謝和信號轉導通路。

Watanabe等將亞細胞分子和心臟器官生理水平關聯(lián)起來，開發(fā)了一個基于有限元的三維模擬程序[45]。該模擬程序結合了興奮收縮偶聯(lián)及其應用的細胞學機制，模擬包含人類左心室收縮功能的流固耦合。該模型在細胞元模型中采用了描述橫橋動力學的神經(jīng)元模型和四態(tài)模型，用有限元網(wǎng)格對心腔內(nèi)的心室墻和血流進行建模，應用任意拉格朗日-歐拉公式求解自動網(wǎng)格的更新，并采用了強耦合策略。在模型中，還使用電子模擬肺循環(huán)，左心房作為前負荷，彈性腔模型作為后負荷，模擬心室充盈的動力和噴射過程，如圖11所示。他們成功模擬再現(xiàn)早期快速充盈和心房壓縮的雙相填充流過程，近似于臨床的觀測結果。因此，該模擬器中的流固耦合模型能分析出充盈流的波傳播速度，同時對于在宏觀層面上建立分子異常和臨床疾病之間的聯(lián)系是一個有力的分析工具。

圖11 左心室收縮和心室內(nèi)血流的時延圖像[45]Fig.11 Time-lapse images of LV contraction and intraventricular blood flow[45]

程永光等為了驗證在左心室充盈過程中基于流固耦合架構的商業(yè)軟件的可行性，研究了心臟充盈階段的流固耦合模擬過程，顯示出其在心臟功能研究和臨床應用領域的巨大前景[10]。他們使用實現(xiàn)了任意拉格朗日-歐拉公式求解的商業(yè)軟件ADINA來模擬左心室的流場，具體解決方案是采用Navier-Stokes方程的任意拉格朗日-歐拉公式求解流體問題，采用有限元模型對心肌組織進行建模求解，并在每個時間步長里同時求解耦合方程。對于人類心臟的左心室充盈流，采用基于規(guī)定時變楊氏膜量的三維橢圓薄壁模型幾何結構來進行模擬，如圖12所示。雖然充盈過程中結構變形很大，但耦合能平穩(wěn)收斂。

圖12 用于流固耦合模擬的心室模型幾何結構。綠色表示心室薄壁，橙色表示流場、入口和二尖瓣口[10]Fig.12 Geometry of the model ventricle for the coupled simulation. ventricular thin wall; flow domain andin let, mitral orifice[10]

程永光等詳細分析了心室模型的壓力容積關系、壓力的空間時間分布、瞬態(tài)速度向量和渦流模式，認為分析結果和實際測量數(shù)據(jù)從定性和定量角度都一致。在圖12中，可以觀察到不同時間步長下的速度模式。在快速充盈階段，會產(chǎn)生一個強力噴射，將血流從二尖瓣噴射到心室腔內(nèi)。在這個強力噴射的驅動下，在二尖瓣口下方的區(qū)域會形成一個環(huán)流。隨著快速充盈的繼續(xù)，這個渦環(huán)平穩(wěn)地向心室中部移動。作為一個附加結果，在減速期會在主動脈流出口處形成一個小的渦環(huán)。隨后，在心臟舒張后期，渦環(huán)會到達膨脹的心室中部，同時渦環(huán)的前端向后移動，而后端仍然向前移動。同時，逐漸減弱的渦流開始擴大并移出通道，進而在心室舒張后期被輸入流摧毀。充盈階段幾個典型時刻的速度向量如圖13所示。

圖13 充盈階段幾個典型時刻的速度向量[10]Fig.13 Velocity vectors at several typical moments during the filling[10]

研究者對該模型的模擬過程做了許多簡化假設，如薄壁心室?guī)缀谓Y構、各向異性和同質材料，忽略二尖瓣以及統(tǒng)一輸入流速度，都明顯限制了模擬模型的實際臨床應用。但是，研究工作驗證了這個框架的結構變形耦合能力，為未來建立真實心臟充盈過程模型提供了有力的參考。

Sebastian等建立了一種心臟血流流場模擬的三維代碼耦合方法[33]。在該方法中，他們首次用心肌復合模型替換了原來規(guī)定幾何運動的Karlsruhe心臟模型，并且分別利用有限體積法和有限元法，實現(xiàn)任意拉格朗日-歐拉公式和固體力學有限元彈性方程的離散化，同時通過數(shù)據(jù)交換迭代策略，保證底層基礎控制方程的接口平穩(wěn)。在該模型使用的代碼耦合框架中，分別應用Abaqus 6.7-1(www.simulia.com)以及Fluent6.3.26 (www.fluent.com) 作為FEM和FVM的解算器，接口交互通信使用由弗勞恩霍夫研究所(www.scai.fraunhofer.de)開發(fā)的網(wǎng)格并行代碼耦合接口MpCCI 3.0.5 (www.mpcci.de)來實現(xiàn)，最后通過迭代計算完成心臟流場的模擬。圖14展示的是該研究模擬不同時間段的心室內(nèi)流場。

圖14 KaHMo FSI/MRT中不同時間段的心室內(nèi)流場[33](左上：充盈早期；右上：充盈后期；左下：心房收縮；右下：中期噴射)Fig.14 Inner-ventricular flow field in KaHMo FSI/MRT for different time-frames[33](Top left: early filling; top right: late filling;bottom left: atrium contraction; bottom right: mid ejection)

2 存在的問題及發(fā)展趨勢

綜上所述，基于CFD的左心室流場數(shù)值模擬研究工作已經(jīng)取得了巨大的進步，但仍然存在許多問題。

首先，基于采集圖像的建模方法大多都使用MRI數(shù)據(jù)來對左心室的形態(tài)和運動進行建模，MRI不僅檢測費用較高，而且通常臨床常規(guī)檢查所獲取數(shù)據(jù)的時間分辨率非常低(9～17幀/s)。超聲心動圖有更高的時間分辨率，在扇形B模式條件下通?？梢赃_到100～200幀/s，在實驗環(huán)境下甚至可以達到2 000幀/s，可以有效克服這個問題。此外，由于強磁場的原因,MRI對諸如體內(nèi)有磁金屬或起搏器的特殊患者不適用。因此，需要進一步研究基于超聲心動圖的左心室建模方法。另一個問題是如何通過實驗數(shù)據(jù)實現(xiàn)三維幾何結構的自動生成。通常，采集的圖像數(shù)據(jù)是二維切面圖，通過在這些圖像上插值來實現(xiàn)三維幾何結構的重建。當前，這個過程大多在計算機輔助設計軟件上手工完成或半自動生成，這將消耗大量的時間。建立自動化分析方法可以節(jié)省時間，排除或最小化人機交互過程產(chǎn)生的主觀性。最后，需要在通過成像數(shù)據(jù)重構三維幾何結構的過程中量化復雜的變化過程，比如心動周期內(nèi)心室的運動變化、成像過程中主體的運動以及成像體位位置的不確定，將來都需要進一步量化從而支持臨床決策。

總之，左心室數(shù)值模擬研究未來的發(fā)展方向將集中在以下幾個方面：

1)臨床圖像數(shù)據(jù)時間分辨率的提高。當前，核磁共振掃描圖像的空間分辨率較高，但時間分辨率非常低，需要采集多個心動周期的圖像才能實現(xiàn)左心室數(shù)值模擬建模，在此過程中容易遺漏一些細節(jié)信息。未來可以繼續(xù)通過提高磁共振成像設備的掃描速度來提高圖像的時間分辨率，但短期內(nèi)還難以實現(xiàn)。此外，可以結合計算機圖形圖像技術，充分利用超聲心動圖高時間分辨率的優(yōu)點，實現(xiàn)對左心室的數(shù)值模擬建模。

2)左心室流場流固耦合新方法的研究。雖然已有許多左心室流場數(shù)值模擬研究充分考慮了流固耦合因素，并對此做了深入研究，做出了較好的模擬結果，但這些研究在計算效率和模擬準確度上都需要進一步提升，因此建立一套在心肌力學上簡化但在血流流場上模擬效果好、計算效率高的左心室流固耦合計算模式非常重要。

3)左心室二尖瓣膜的精確模擬。由于現(xiàn)有的成像設備和圖像技術通常難以準確提取二尖瓣膜的位置和移動數(shù)據(jù)信息，現(xiàn)在的研究大多通過數(shù)據(jù)假設來完成二尖瓣膜的數(shù)值模擬。隨著科技的進一步發(fā)展，通過對各種成像設備圖像信息的研究，可以結合圖像處理技術來提取真實的二尖瓣膜數(shù)據(jù)信息，以實現(xiàn)更為準確的左心室流場模擬。

4)心臟扭曲運動的模擬。心臟的扭曲運動是導致升主動脈內(nèi)旋動流的主要原因，也是影響左心室內(nèi)部血流流場變化的重要因素，但目前的數(shù)值模擬研究基本都還未考慮到這一因素對流場變化的影響。

5)三維心室?guī)缀文Ｐ偷淖詣由?。計算圖像分割技術已經(jīng)日趨成熟，可以研究適合超聲圖像的圖像分割算法，自動分割一個心動周期內(nèi)的超聲圖像，實現(xiàn)自動構造三維心室?guī)缀文Ｐ汀?/p>

6)量化三維心室?guī)缀文Ｐ椭亟ㄟ^程中復雜變化產(chǎn)生的影響。在三維幾何模型重建過程中會產(chǎn)生各種復雜的細微變化，這些變化會造成心室內(nèi)模擬流場的變化，并最終影響模擬流場的準確性。這些復雜的細微變化可以通過反復實驗來確定，并建立模型量化其影響，通過建立參數(shù)調整模型來最小化這些變化產(chǎn)生的影響。

7)三維心臟流場數(shù)值模擬模型中臨床病人數(shù)據(jù)的注入配置方案。支持臨床應用是三維心臟流場數(shù)值模擬的主要目的，便捷高效地注入配置臨床病人數(shù)據(jù)，準確再現(xiàn)一個心動周期內(nèi)病人的三維心臟流場變化過程，從而滿足臨床診斷和治療的需要。

[1] Richter Y, Edelman ER. Cardiology is flow[J]. Circulation, 2006, 113(23):2679-2682.

[2] Carlh?ll CJ, Bolger A. Passing strange flow in the failing ventricle[J]. Circulation: Heart Failure, 2010, 3(2):326-331.

[3] 劉有軍, 刁越, 高松. 血流動力學與動脈粥樣硬化 (一)[J]. 北京生物醫(yī)學工程, 2004, 23(4):293-295.

[4] Foster GE, Deng Z, Boulet LM et al. Changes in left ventricular function and coronary blood flow velocity during isocapnic hypoxia: A cardiac magnetic resonance imaging study[J]. Journal of Cardiovascular Magnetic Resonance, 2016, 18(S1)： 1-2

[5] Zhu Y, Qian M, Niu L et al. Two-dimensional measurement of blood flow velocity in rat arteries based on ultrasonic particle image velocimetry[J]. Journal of Southern Medical University, 2014, 34(9):1305-1309.

[6] Zhou F, Cui YY, Wu LL et al. Analysis of flow field in mechanical aortic bileaflet heart valves using finite volume method[J]. Journal of Medical & Biological Engineering, 2016, 36(1): 110-120

[7] Shang EK, Nathan DP, Sprinkle SR et al. Peak wall stress predicts expansion rate in descending thoracic aortic aneurysms[J]. Annals of Thoracic Surgery, 2013, 95(2):593-598.

[8] Altnji HE, Bou-Sa?d B, Berre WL. Morphological and stent design risk factors to prevent migration phenomena for a thoracic aneurysm: A numerical analysis[J]. Medical Engineering & Physics, 2014, 37(1):23-33.

[9] 隋成龍. 左心室舒張早期功能模型評價方法研究[D]. 西安：西南交通大學， 2013.

[10] Cheng Y, Oertel H, Schenkel T. Fluid-structure coupled CFD simulation of the left ventricular flow during filling phase[J]. Annals of Biomedical Engineering, 2005, 33(5):567-576.

[11] Baccani B, Domenichini F, Pedrizzetti G. Vortex dynamics in a model left ventricle during filling[J]. European Journal of Mechanics-B/Fluids, 2002, 21(5):527-543.

[12] Baccani B, Domenichini F, Pedrizzetti G. Model and influence of mitral valve opening during the left ventricular filling[J]. Journal of Biomechanics, 2003, 36(3):355-361.

[13] Baccani B, Domenichini F, Pedrizzetti G et al. Fluid dynamics of the left ventricular filling in dilated cardiomyopathy[J]. Journal of Biomechanics, 2002, 35(5):665-671.

[14] Long Q, Merrifield R, Yang G et al. The influence of inflow boundary conditions on intra left ventricle flow predictions[J]. Journal of Biomechanical Engineering, 2003, 125(6):922-927.

[15] Nakamura M, Wada S, Mikami T et al. A Computational Fluid Mechanical Study on the Effects of Opening and Closing of the Mitral Orifice on a Transmitral Flow Velocity Profile and an Early Diastolic Intraventricular Flow[J]. JSME International Journal Series C Mechanical Systems, Machine Elements and Manufacturing, 2002, 45(4):913-922.

[16] Nakamura M, Wada S, Mikami T et al. Computational study on the evolution of an intraventricular vortical flow during early diastole for the interpretation of color M-mode Doppler echocardiograms[J]. Biomechanics and Modeling in Mechanobiology, 2003, 2(2):59-72.

[17] Saber NR, Gosman A, Wood NB et al. Computational flow modeling of the left ventricle based on in vivo MRI data: initial experience[J]. Annals of Biomedical Engineering, 2001, 29(4):275-283.

[18] Saber NR, Wood NB, Gosman A et al. Progress towards patient-specific computational flow modeling of the left heart via combination of magnetic resonance imaging with computational fluid dynamics[J]. Annals of Biomedical Engineering, 2003, 31(1):42-52.

[19] Schoephoerster RT, Silva CL, Ray G. Evaluation of left ventricular function based on simulated systolic flow dynamics computed from regional wall motion[J]. Journal of Biomechanics, 1994, 27(2):125-136.

[20] Taylor TW, Okino H, Yamaguchi T. Three-dimensional analysis of left ventricular ejection using computational fluid dynamics[J]. Journal of Biomechanical Engineering, 1994, 116(1):127-130.

[21] Long Q, Merrifield R, Xu X et al. Subject-specific computational simulation of left ventricular flow based on magnetic resonance imaging[J]. Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine, 2008, 222(4):475-485.

[22] Schenkel T, Malve M, Reik M et al. MRI-based CFD analysis of flow in a human left ventricle: methodology and application to a healthy heart[J]. Annals of Biomedical Engineering, 2009, 37(3):503-515.

[23] Vierendeels J, Riemslagh K, Dick E et al. Computer simulation of intraventricular flow and pressure gradients during diastole[J]. Journal of Biomechanical Engineering, 2000, 122(6):667-674.

[24] Lemmon JD, Yoganathan AP. Three-dimensional computational model of left heart diastolic function with fluid-structure interaction[J]. Journal of Biomechanical Engineering, 2000, 122(2):109-117.

[25] Lemmon JD, Yoganathan AP. Computational modeling of left heart diastolic function: examination of ventricular dysfunction[J]. Journal of Biomechanical Engineering, 2000, 122(4):297-303.

[26] McQueen D, Peskin C. Shared-memory parallel vector implementation of the immersed boundary method for the computation of blood flow in the beating mammalian heart[J]. The Journal of Supercomputing, 1997, 11(3):213-236.

[27] McQueen DM, Peskin CS. A three-dimensional computer model of the human heart for studying cardiac fluid dynamics[J]. ACM SIGGRAPH Computer Graphics, 2000, 34(1):56-60.

[28] Peskin CS. Fluid dynamics of the heart and its valves[C]//Proceedings of APS meeting. New York: American Physical Society, 1997: 2.

[29] Vesier CC, Lemmon Jr JD, Levine RA et al. A three-dimensional computational model of a thin-walled left ventricle. In: Proceedings of the 1992 ACM/IEEE conference on Supercomputing, 1992: IEEE Computer Society Press, 1992: 73-82.

[30] Chahboune B, Crolet J. Numerical simulation of the blood-wall interaction in the human left ventricle[J]. The European Physical Journal Applied Physics, 1998, 2(3):291-297.

[31] Vierendeels J, Riemslagh K, Dick E, et al. Computer simulation of left ventricular filling flow: impact study on echocardiograms[J]. Computers in Cardiology, 1999, 26: 177-180.

[32] Watanabe H, Hisada T, Sugiura S et al. Computer simulation of blood flow, left ventricular wall motion and their interrelationship by fluid-structure interaction finite element method[J]. JSME International Journal Series C, 2002, 45(4):1003-1012.

[33] Krittian S, Janoske U, Oertel H, et al. Partitioned fluid-solid coupling for cardiovascular blood flow[J]. Annals of Biomedical Engineering, 2010, 38(4):1426-1441.

[34] Domenichini F, Pedrizzetti G, Baccani B. Three-dimensional filling flow into a model left ventricle[J]. Journal of Fluid Mechanics, 2005, 539:179-198.

[35] Sonntag SJ, Li W, Becker M et al. Combined computational and experimental approach to improve the assessment of mitral regurgitation by echocardiography[J]. Annals of Biomedical Engineering, 2014, 42(5):971-985.

[36] De Hart J, Peters G, Schreurs P et al. A three-dimensional computational analysis of fluid-structure interaction in the aortic valve[J]. Journal of Biomechanics, 2003, 36(1):103-112.

[37] Tai C, Liew K, Zhao Y. Numerical simulation of 3D fluid-structure interaction flow using an immersed object method with overlapping grids[J]. Computers & Structures, 2007, 85(11):749-762.

[38] Fung YC. Biomechanics: circulation[M]. Berlin: Springer Science & Business Media, 2013.

[39] Lin D, Yin F. A multiaxial constitutive law for mammalian left ventricular myocardium in steady-state barium contracture or tetanus[J]. Journal of Biomechanical Engineering, 1998, 120(4):504-517.

[40] Hunter PJ, Pullan AJ, Smaill BH. Modeling total heart function[J]. Annual Review of Biomedical Engineering, 2003, 5(1):147-177.

[41] Smith N, Nickerson D, Crampin E, et al. Computational mechanics of the heart. from tissue structure to ventricular function[J]. J Elasticity, 2000, 61(1):113-141.

[42] Souli M, Ouahsine A, Lewin L. ALE formulation for fluid-structure interaction problems[J]. Computer Methods in Applied Mechanics and Engineering, 2000, 190(5):659-675.

[43] Zhang Q, Hisada T. Analysis of fluid-structure interaction problems with structural buckling and large domain changes by ALE finite element method[J]. Computer Methods in Applied Mechanics and Engineering, 2001, 190(48):6341-6357.

[44] Zhang H, Bathe KJ. Direct and iterative computing of fluid flows fully coupled with structures[C]//First MIT Conference on Computational Fluid and Solid Mechanics. Boston: MIT, 2001: 1440-1443.

[45] Watanabe H, Sugiura S, Kafuku H, et al. Multiphysics simulation of left ventricular filling dynamics using fluid-structure interaction finite element method[J]. Biophysical Journal, 2004, 87(3):2074-2085.

The Recent Advances and Prospects on the Numerical Simulation of Flow Field within Left Ventricular Based on CFD Techniques

Adu Jianhua1,2Yin Lixue2*Xie Shenghua2

1(ChengduUniversityofInformationTechnology,Chengdu610225,China)2(SichuanAcademyofMedicalSciences&SichuanProvincialPeople′sHospital,SichuanProvincialKeyLaboratoryofUltrasoundinCardiacElectrophysiologyandofBiomechanics,Chengdu610072,China)

With the rapid development of computer technology, computational fluid dynamics has been widely used in the numerical simulation of flow field in heart and some valuable relevant outcomes has been obtained. CFD is playing an increasingly important role in the study of cardiac function. Currently, it has become an important research in basic and clinical cardiovascular research areas in which CFD is used to simulate and analyze fluid dynamics and flow field in heart, which is important for prevention, diagnosis and treatment of the cardiovascular disease. In this paper, we describe the development of numerical simulation of flow field within left ventricular based on CFD, the research status and classification method. We summarize the main issues of numerical simulation of flow field within left ventricular based on CFD and put forward the future direction of this research field.

computational fluid dynamics; flow field within left ventricular; numerical simulation; fluid-structure coupling

10.3969/j.issn.0258-8021. 2016. 05.011

2016-02-22， 錄用日期:2016-04-14

國家自然科學基金(30970698);四川省圖形圖像智能系統(tǒng)協(xié)同創(chuàng)新中心開放基金;成都信息工程大學引進人才項目(KYTZ201322)

R318

A

0258-8021(2016) 05-0587-011

*通信作者(Corresponding author)， E-mail: yinlixue@yahoo.com