• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

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

    2016-02-17 06:26:40阿都建華尹立雪謝盛華
    中國生物醫(yī)學工程學報 2016年5期
    關鍵詞:心室流場左心室

    阿都建華 尹立雪 謝盛華

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

    6)量化三維心室?guī)缀文P椭亟ㄟ^程中復雜變化產(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

    猜你喜歡
    心室流場左心室
    心電向量圖診斷高血壓病左心室異常的臨床應用
    大型空冷汽輪發(fā)電機轉子三維流場計算
    大電機技術(2021年2期)2021-07-21 07:28:24
    轉杯紡排雜區(qū)流場與排雜性能
    基于HYCOM的斯里蘭卡南部海域溫、鹽、流場統(tǒng)計分析
    房阻伴特長心室停搏1例
    基于瞬態(tài)流場計算的滑動軸承靜平衡位置求解
    冠心病伴心力衰竭合并心房顫動的心室率控制研究
    初診狼瘡腎炎患者左心室肥厚的相關因素
    『孫爺爺』談心室為孩子上網(wǎng)撐起『保護傘』
    中國火炬(2014年6期)2014-07-24 14:16:34
    基于改進多尺度ASM和非剛性配準的4D-CT左心室分割
    日韩有码中文字幕| 精品国产国语对白av| 久久天躁狠狠躁夜夜2o2o| 欧美国产精品va在线观看不卡| 国产成人av教育| 国产成人欧美在线观看| 亚洲狠狠婷婷综合久久图片| 亚洲精品国产区一区二| 亚洲av片天天在线观看| svipshipincom国产片| 黄片大片在线免费观看| 亚洲精品美女久久久久99蜜臀| 在线十欧美十亚洲十日本专区| 熟女少妇亚洲综合色aaa.| 91国产中文字幕| 亚洲精品国产精品久久久不卡| 国产又爽黄色视频| 国产主播在线观看一区二区| 大香蕉久久成人网| 欧美性猛交黑人性爽| 久久精品夜夜夜夜夜久久蜜豆 | 色在线成人网| 淫妇啪啪啪对白视频| av福利片在线| 亚洲性夜色夜夜综合| 亚洲专区国产一区二区| 亚洲欧洲精品一区二区精品久久久| 欧美成人一区二区免费高清观看 | 美女午夜性视频免费| 日日爽夜夜爽网站| 一进一出抽搐动态| 90打野战视频偷拍视频| 午夜影院日韩av| 国产1区2区3区精品| 久久久久久亚洲精品国产蜜桃av| 香蕉国产在线看| 啦啦啦免费观看视频1| 久久精品人妻少妇| 一级a爱视频在线免费观看| 天天一区二区日本电影三级| 日日摸夜夜添夜夜添小说| 色尼玛亚洲综合影院| 两性夫妻黄色片| 中文字幕另类日韩欧美亚洲嫩草| 天堂动漫精品| 在线视频色国产色| 夜夜看夜夜爽夜夜摸| 免费无遮挡裸体视频| 亚洲av成人一区二区三| 欧美另类亚洲清纯唯美| 国产视频内射| 国产av又大| 欧美激情极品国产一区二区三区| 欧美zozozo另类| 熟女少妇亚洲综合色aaa.| 无遮挡黄片免费观看| 成人午夜高清在线视频 | 一级作爱视频免费观看| 在线观看日韩欧美| 真人一进一出gif抽搐免费| 久久精品成人免费网站| 色播亚洲综合网| 国产亚洲欧美在线一区二区| 日韩三级视频一区二区三区| 国产欧美日韩一区二区三| 国产成人系列免费观看| 亚洲成人久久性| 麻豆久久精品国产亚洲av| 精品少妇一区二区三区视频日本电影| 在线国产一区二区在线| 正在播放国产对白刺激| 啪啪无遮挡十八禁网站| 黄色女人牲交| 国产三级在线视频| 一区二区三区高清视频在线| 一边摸一边做爽爽视频免费| 国产黄片美女视频| 欧美乱码精品一区二区三区| 黄色毛片三级朝国网站| 免费在线观看日本一区| 国产高清videossex| tocl精华| 在线观看免费日韩欧美大片| 欧美日韩中文字幕国产精品一区二区三区| 国产欧美日韩一区二区三| 日本免费一区二区三区高清不卡| 欧美日韩瑟瑟在线播放| 亚洲午夜理论影院| 一进一出抽搐动态| 午夜福利视频1000在线观看| av超薄肉色丝袜交足视频| 最近最新中文字幕大全免费视频| 十八禁网站免费在线| 日韩 欧美 亚洲 中文字幕| 人人妻人人澡欧美一区二区| 自线自在国产av| 午夜福利高清视频| 国产精品综合久久久久久久免费| 啦啦啦观看免费观看视频高清| 午夜激情福利司机影院| 欧美成人一区二区免费高清观看 | 欧美最黄视频在线播放免费| 黄色a级毛片大全视频| 日日摸夜夜添夜夜添小说| 久久青草综合色| 9191精品国产免费久久| 9191精品国产免费久久| 精品国产国语对白av| 美国免费a级毛片| 中文资源天堂在线| 国产亚洲精品综合一区在线观看 | 成熟少妇高潮喷水视频| 欧美一级毛片孕妇| 成人三级做爰电影| 人妻久久中文字幕网| 特大巨黑吊av在线直播 | 老熟妇仑乱视频hdxx| 91字幕亚洲| 亚洲精品中文字幕在线视频| 欧美亚洲日本最大视频资源| 国内少妇人妻偷人精品xxx网站 | 午夜亚洲福利在线播放| 香蕉丝袜av| 亚洲欧美一区二区三区黑人| 日韩视频一区二区在线观看| 久久久国产精品麻豆| 久久人妻福利社区极品人妻图片| av视频在线观看入口| 国产精品久久久久久人妻精品电影| 老鸭窝网址在线观看| 日韩av在线大香蕉| 国产成人精品无人区| 免费看十八禁软件| 久久久久九九精品影院| 18禁裸乳无遮挡免费网站照片 | 在线观看日韩欧美| 热99re8久久精品国产| 在线观看免费午夜福利视频| 精品一区二区三区四区五区乱码| 午夜激情福利司机影院| 欧美一级毛片孕妇| 观看免费一级毛片| 男男h啪啪无遮挡| 日韩精品中文字幕看吧| 欧美成人午夜精品| 亚洲精品一卡2卡三卡4卡5卡| 757午夜福利合集在线观看| 又黄又爽又免费观看的视频| 国产精品一区二区精品视频观看| 亚洲国产精品合色在线| 国产精品一区二区精品视频观看| 一边摸一边做爽爽视频免费| 91字幕亚洲| 老熟妇仑乱视频hdxx| 国产男靠女视频免费网站| 亚洲国产精品999在线| 久久精品人妻少妇| 国产男靠女视频免费网站| 亚洲美女黄片视频| 男人舔女人下体高潮全视频| 日韩三级视频一区二区三区| 十八禁网站免费在线| 90打野战视频偷拍视频| 日日干狠狠操夜夜爽| 免费搜索国产男女视频| 在线观看免费日韩欧美大片| 午夜精品久久久久久毛片777| 日本撒尿小便嘘嘘汇集6| 国内久久婷婷六月综合欲色啪| 国内精品久久久久精免费| 91大片在线观看| 成人永久免费在线观看视频| 曰老女人黄片| 熟女电影av网| 亚洲自拍偷在线| 国产欧美日韩精品亚洲av| 国产亚洲精品久久久久5区| 久久久精品国产亚洲av高清涩受| 精品久久蜜臀av无| 精品午夜福利视频在线观看一区| 精品第一国产精品| 一级毛片女人18水好多| 久久久久久久久中文| 亚洲av中文字字幕乱码综合 | 国产蜜桃级精品一区二区三区| 久久天堂一区二区三区四区| 国产亚洲欧美在线一区二区| 一夜夜www| 日韩高清综合在线| 久99久视频精品免费| www国产在线视频色| 欧美大码av| 啦啦啦韩国在线观看视频| 国产免费男女视频| 亚洲国产精品合色在线| 啦啦啦 在线观看视频| 日本 av在线| 亚洲三区欧美一区| 国产激情欧美一区二区| 久久精品国产综合久久久| 91成年电影在线观看| 叶爱在线成人免费视频播放| 变态另类丝袜制服| 国产又色又爽无遮挡免费看| 色av中文字幕| 欧美一级毛片孕妇| 老司机午夜福利在线观看视频| 99精品久久久久人妻精品| 欧美成人性av电影在线观看| 日本精品一区二区三区蜜桃| 18禁国产床啪视频网站| 哪里可以看免费的av片| 国产精品免费一区二区三区在线| 一区二区日韩欧美中文字幕| 成人三级黄色视频| 在线视频色国产色| 亚洲美女黄片视频| 国产在线观看jvid| 91成人精品电影| 亚洲 国产 在线| 亚洲熟女毛片儿| 99re在线观看精品视频| 天堂动漫精品| 亚洲七黄色美女视频| 男人操女人黄网站| 久久精品91无色码中文字幕| 18美女黄网站色大片免费观看| 一边摸一边抽搐一进一小说| 亚洲精品中文字幕一二三四区| 色综合站精品国产| 国产精品日韩av在线免费观看| 日韩av在线大香蕉| 18美女黄网站色大片免费观看| 99riav亚洲国产免费| e午夜精品久久久久久久| 99久久国产精品久久久| 国产成人av教育| 日本一区二区免费在线视频| 亚洲欧美激情综合另类| 国产区一区二久久| 成人三级做爰电影| 亚洲 欧美一区二区三区| 日本一本二区三区精品| 校园春色视频在线观看| 午夜免费激情av| 午夜福利在线观看吧| 亚洲五月天丁香| 黄频高清免费视频| 亚洲真实伦在线观看| 久久精品国产清高在天天线| 欧美激情 高清一区二区三区| 欧美黑人欧美精品刺激| 国产精品一区二区三区四区久久 | 精品国产乱子伦一区二区三区| 99国产精品一区二区蜜桃av| 精品国产乱码久久久久久男人| 日韩精品青青久久久久久| 亚洲真实伦在线观看| www.自偷自拍.com| 色老头精品视频在线观看| 手机成人av网站| 又黄又爽又免费观看的视频| 又黄又粗又硬又大视频| 成人精品一区二区免费| 人妻丰满熟妇av一区二区三区| 亚洲欧美日韩高清在线视频| 国产1区2区3区精品| 亚洲 欧美一区二区三区| 免费观看精品视频网站| 久久 成人 亚洲| 久久久久久久久中文| 色综合婷婷激情| 色综合欧美亚洲国产小说| 高清在线国产一区| 日本成人三级电影网站| 在线十欧美十亚洲十日本专区| 国产一级毛片七仙女欲春2 | 亚洲av片天天在线观看| 国产高清videossex| 夜夜躁狠狠躁天天躁| 欧美一级a爱片免费观看看 | 亚洲一码二码三码区别大吗| av在线天堂中文字幕| 欧美激情高清一区二区三区| 啪啪无遮挡十八禁网站| 俺也久久电影网| √禁漫天堂资源中文www| 国产精品久久电影中文字幕| 日本一本二区三区精品| 亚洲av片天天在线观看| 97碰自拍视频| 天堂√8在线中文| АⅤ资源中文在线天堂| 久久精品国产亚洲av高清一级| 午夜福利成人在线免费观看| 啦啦啦观看免费观看视频高清| 91在线观看av| 美女国产高潮福利片在线看| 成人av一区二区三区在线看| 国产av不卡久久| 色精品久久人妻99蜜桃| 欧美一级毛片孕妇| 观看免费一级毛片| 国产精品乱码一区二三区的特点| 免费看a级黄色片| 国产免费av片在线观看野外av| 淫秽高清视频在线观看| 国产精品久久视频播放| 国产精品亚洲av一区麻豆| 热re99久久国产66热| 一级毛片女人18水好多| 日韩欧美在线二视频| 可以在线观看毛片的网站| 成在线人永久免费视频| 黄频高清免费视频| 亚洲七黄色美女视频| 天天添夜夜摸| 给我免费播放毛片高清在线观看| 亚洲av成人一区二区三| 日韩中文字幕欧美一区二区| 国产91精品成人一区二区三区| 精品久久久久久久毛片微露脸| 国产久久久一区二区三区| 久久久久国产精品人妻aⅴ院| 亚洲第一av免费看| 日韩高清综合在线| 欧美另类亚洲清纯唯美| 精品国产国语对白av| 久久国产精品影院| 久久香蕉国产精品| 国产精华一区二区三区| 欧美国产日韩亚洲一区| 亚洲成人久久性| 国产精品一区二区三区四区久久 | 国产精品av久久久久免费| 国产色视频综合| 精品高清国产在线一区| 日韩免费av在线播放| АⅤ资源中文在线天堂| 成人国语在线视频| 热re99久久国产66热| 天天躁夜夜躁狠狠躁躁| 亚洲av熟女| 性欧美人与动物交配| 欧美绝顶高潮抽搐喷水| 黄色 视频免费看| a级毛片a级免费在线| 日本一区二区免费在线视频| av片东京热男人的天堂| 最新在线观看一区二区三区| 国产精品久久久久久人妻精品电影| av超薄肉色丝袜交足视频| av在线播放免费不卡| 亚洲成av片中文字幕在线观看| 色精品久久人妻99蜜桃| 国产片内射在线| 91九色精品人成在线观看| 哪里可以看免费的av片| 亚洲av片天天在线观看| 精品电影一区二区在线| 久久国产精品男人的天堂亚洲| 19禁男女啪啪无遮挡网站| 日韩高清综合在线| 精品一区二区三区四区五区乱码| 国产精品 欧美亚洲| 亚洲av五月六月丁香网| 成人特级黄色片久久久久久久| 亚洲全国av大片| 在线观看一区二区三区| 国产精品免费一区二区三区在线| 国产精品精品国产色婷婷| bbb黄色大片| 精品久久久久久久人妻蜜臀av| 黄色 视频免费看| 国产伦人伦偷精品视频| 国产极品粉嫩免费观看在线| 亚洲人成网站高清观看| 自线自在国产av| 亚洲精品美女久久久久99蜜臀| 久久午夜亚洲精品久久| www.熟女人妻精品国产| 97碰自拍视频| 丝袜美腿诱惑在线| 中文字幕精品亚洲无线码一区 | 69av精品久久久久久| 国产成人欧美在线观看| 一级片免费观看大全| 高潮久久久久久久久久久不卡| 欧美av亚洲av综合av国产av| 欧美成人性av电影在线观看| 中出人妻视频一区二区| 18禁观看日本| 琪琪午夜伦伦电影理论片6080| 在线看三级毛片| 精品高清国产在线一区| 欧美激情 高清一区二区三区| 又大又爽又粗| 脱女人内裤的视频| 日韩免费av在线播放| 亚洲欧美一区二区三区黑人| 国产成人啪精品午夜网站| 国产免费男女视频| 亚洲精品国产精品久久久不卡| 无人区码免费观看不卡| 亚洲国产欧美网| 韩国精品一区二区三区| 俄罗斯特黄特色一大片| 欧美成人免费av一区二区三区| 国产高清视频在线播放一区| 90打野战视频偷拍视频| 日韩有码中文字幕| 女同久久另类99精品国产91| 色综合亚洲欧美另类图片| 91国产中文字幕| 成人18禁在线播放| 美女国产高潮福利片在线看| 777久久人妻少妇嫩草av网站| videosex国产| 女人被狂操c到高潮| www日本在线高清视频| 日韩欧美国产一区二区入口| 午夜免费激情av| 日韩欧美国产在线观看| 好男人电影高清在线观看| 两性夫妻黄色片| 男男h啪啪无遮挡| 国产不卡一卡二| 日韩一卡2卡3卡4卡2021年| 搞女人的毛片| 国产一区二区激情短视频| bbb黄色大片| 真人做人爱边吃奶动态| 国产三级在线视频| 亚洲av熟女| 成人精品一区二区免费| 精品久久久久久,| 国产麻豆成人av免费视频| 亚洲自偷自拍图片 自拍| 亚洲人成伊人成综合网2020| 国产成人精品久久二区二区91| 亚洲成国产人片在线观看| 亚洲专区国产一区二区| 亚洲第一青青草原| 一边摸一边抽搐一进一小说| 脱女人内裤的视频| 每晚都被弄得嗷嗷叫到高潮| 日本精品一区二区三区蜜桃| 亚洲精品久久国产高清桃花| 精品一区二区三区av网在线观看| 黄网站色视频无遮挡免费观看| 亚洲精品色激情综合| 国产精品久久视频播放| 午夜免费鲁丝| 国产私拍福利视频在线观看| 亚洲av电影不卡..在线观看| 亚洲国产毛片av蜜桃av| 18禁黄网站禁片午夜丰满| 在线观看www视频免费| 婷婷六月久久综合丁香| 桃色一区二区三区在线观看| 日韩精品中文字幕看吧| 久久狼人影院| 亚洲一码二码三码区别大吗| 久久国产精品人妻蜜桃| 免费在线观看影片大全网站| 不卡一级毛片| 精品乱码久久久久久99久播| 中文资源天堂在线| 精品人妻1区二区| 母亲3免费完整高清在线观看| 观看免费一级毛片| av免费在线观看网站| 日韩欧美三级三区| 亚洲五月婷婷丁香| 国产99久久九九免费精品| 久久久国产欧美日韩av| 1024视频免费在线观看| 亚洲 国产 在线| 亚洲av成人av| 一个人免费在线观看的高清视频| 亚洲国产精品999在线| 免费在线观看视频国产中文字幕亚洲| 特大巨黑吊av在线直播 | 欧美三级亚洲精品| 国产精品一区二区免费欧美| 国产精品久久久av美女十八| 男女之事视频高清在线观看| 国产成人系列免费观看| 国产亚洲精品久久久久5区| 操出白浆在线播放| 91九色精品人成在线观看| 国产av在哪里看| 欧美日韩乱码在线| 欧美zozozo另类| 久久国产精品人妻蜜桃| 自线自在国产av| 午夜亚洲福利在线播放| 国产高清videossex| 日本黄色视频三级网站网址| 90打野战视频偷拍视频| 亚洲成国产人片在线观看| 超碰成人久久| 日韩欧美一区视频在线观看| 欧美在线黄色| 制服丝袜大香蕉在线| 国产成人啪精品午夜网站| 亚洲专区国产一区二区| 国产成年人精品一区二区| 久久久国产成人免费| 欧美 亚洲 国产 日韩一| 国产不卡一卡二| 久久久久久久久久黄片| 久久久精品国产亚洲av高清涩受| www.自偷自拍.com| 免费人成视频x8x8入口观看| 级片在线观看| 亚洲一区二区三区色噜噜| 黄片小视频在线播放| 国产精品久久视频播放| 看免费av毛片| 日韩 欧美 亚洲 中文字幕| 成年人黄色毛片网站| 国产免费男女视频| 亚洲国产欧美日韩在线播放| а√天堂www在线а√下载| 久久午夜综合久久蜜桃| 国产又色又爽无遮挡免费看| 真人一进一出gif抽搐免费| 亚洲午夜理论影院| 在线天堂中文资源库| 成年免费大片在线观看| 国产亚洲精品第一综合不卡| 久久中文字幕一级| 非洲黑人性xxxx精品又粗又长| 男人舔奶头视频| 日日干狠狠操夜夜爽| 黄色a级毛片大全视频| 亚洲久久久国产精品| 可以免费在线观看a视频的电影网站| 日本成人三级电影网站| 最近最新免费中文字幕在线| 国产av一区在线观看免费| 午夜福利高清视频| 美女 人体艺术 gogo| 色哟哟哟哟哟哟| 久9热在线精品视频| 欧美精品亚洲一区二区| 听说在线观看完整版免费高清| 亚洲人成伊人成综合网2020| 国产高清视频在线播放一区| 国产精品精品国产色婷婷| 国产精品美女特级片免费视频播放器 | 欧美午夜高清在线| 99久久精品国产亚洲精品| 国产视频一区二区在线看| 在线观看午夜福利视频| 正在播放国产对白刺激| 精品人妻1区二区| 日韩一卡2卡3卡4卡2021年| 久久国产精品影院| 国产91精品成人一区二区三区| 亚洲国产精品久久男人天堂| 男人舔女人的私密视频| 黄片小视频在线播放| 国产高清videossex| 亚洲成人国产一区在线观看| 亚洲av第一区精品v没综合| 在线十欧美十亚洲十日本专区| 日本精品一区二区三区蜜桃| 精品第一国产精品| 国产午夜福利久久久久久| 一二三四社区在线视频社区8| 人妻丰满熟妇av一区二区三区| 男人舔女人下体高潮全视频| 十八禁人妻一区二区| 一边摸一边做爽爽视频免费| 成人手机av| 女性被躁到高潮视频| 丝袜美腿诱惑在线| 制服人妻中文乱码| 国产成年人精品一区二区| 久久精品国产99精品国产亚洲性色| 国产三级在线视频| 1024手机看黄色片| ponron亚洲| 最近最新免费中文字幕在线| 99国产精品99久久久久| 亚洲av日韩精品久久久久久密| 国产精品久久电影中文字幕| 亚洲电影在线观看av| 一本一本综合久久| 黄色成人免费大全| 俺也久久电影网| 日本熟妇午夜| 久久久久久久午夜电影| 韩国av一区二区三区四区| 亚洲精品美女久久av网站| 亚洲欧洲精品一区二区精品久久久| 亚洲av五月六月丁香网| 午夜激情av网站| 午夜精品在线福利| 午夜福利在线在线| 欧美一级毛片孕妇| 免费电影在线观看免费观看| 伦理电影免费视频| 男女床上黄色一级片免费看| 黑人欧美特级aaaaaa片| 90打野战视频偷拍视频| 欧美激情久久久久久爽电影| cao死你这个sao货| 国产精品亚洲一级av第二区| 日本一区二区免费在线视频| 美女大奶头视频| www国产在线视频色| 国产精品免费一区二区三区在线|