孫凱順，張洪明，張桂敏，夏健明，孫毅

1.昆明理工大學(xué)工程力學(xué)系，云南昆明650500；2.云南省阜外心血管病醫(yī)院，云南昆明650032

前言

ABAQUS是一款具有強(qiáng)大非線性分析能力和模擬復(fù)雜系統(tǒng)高可靠性的有限元軟件。在計(jì)算流體軟件中，ABAQUS/CFD模塊有著強(qiáng)大的模擬計(jì)算能力，能模擬層流與湍流，穩(wěn)態(tài)與瞬態(tài)的內(nèi)部和外部流動。雷諾數(shù)范圍很寬，對復(fù)雜幾何有著很強(qiáng)的適應(yīng)性。有限元法和混合有限體積法的求解方法，使其在不可壓縮流的層流與湍流問題上有著較高的求解精度［1］。主動脈夾層（Aortic Dissection,AD）是一種嚴(yán)重的心血管疾病，主動脈內(nèi)壁的內(nèi)膜撕裂口發(fā)生后，血液流入主動脈中膜外層或中外膜交界處，可引起急性主動脈破裂而致患者死亡?；蛞驃A層真腔受到假腔的擠壓而致腔體狹窄和閉塞，腸管、腎臟、下肢等由真腔供血的重要臟器發(fā)生缺血性改變，并伴隨嚴(yán)重并發(fā)癥［2］。Stanford分類法將主動脈夾層分為A型和B型。Stanford A型：夾層破口累及升主動脈，無論第一破口源自何處；Stanford B型：夾層第一破口不累及升主動脈。當(dāng)今Stanford A型發(fā)病率顯著升高，嚴(yán)重威脅人民生命健康。主動脈外科手術(shù)的術(shù)式種類繁多，目前推行按照內(nèi)膜破口位置、數(shù)量和夾層涉及范圍制定個(gè)性化手術(shù)方案［3］，人工血管置換術(shù)就是其中一種有效治療主動脈夾層的術(shù)式［4］。深入研究主動脈血流動力學(xué)規(guī)律的意義重大，通過體內(nèi)實(shí)測獲得一系列主動脈血流動力學(xué)指標(biāo)的分布較為困難?，F(xiàn)如今醫(yī)學(xué)影像學(xué)、計(jì)算流體動力學(xué)的發(fā)展以及高性能計(jì)算機(jī)并行計(jì)算的應(yīng)用，使數(shù)值模擬仿真的研究方法實(shí)現(xiàn)體外模擬主動脈血流成為高效手段，進(jìn)而通過后處理呈現(xiàn)相關(guān)流域的血流動力學(xué)指標(biāo)［5］。主動脈血流數(shù)值模擬的研究已有較多進(jìn)展。邊界條件的差異直接影響模擬計(jì)算的收斂性［6］。關(guān)于主動脈入口流速，較多的研究采用理想化流速［7］。研究表明理想化的流速分布會使模型輸出與體內(nèi)真實(shí)血液動力學(xué)參數(shù)產(chǎn)生一定偏差［8］。流型假設(shè)上，層流被認(rèn)為適用于人類動脈中的血流［9-10］，且被廣泛運(yùn)用到主動脈血流模擬［11-12］，實(shí)際上主動脈血流存在湍流［13］。本文研究對象為一Stanford A型主動脈夾層病例，破口位于升主動脈，治療方式采用人工升主動脈置換術(shù)。獲取患者術(shù)后胸部CT圖像中的AD數(shù)據(jù)以構(gòu)建流域幾何模型［14］，主動脈入口流速使用術(shù)后實(shí)測流速，湍流模型采用Spalart-Allmaras，該湍流模型收斂性較好，可呈現(xiàn)二次流等基本流動特征［15］。經(jīng)模型網(wǎng)格劃分后進(jìn)行數(shù)值模擬。

1 計(jì)算模型和數(shù)值方法

1.1 控制方程

假定血液流動為不可壓縮牛頓流體的粘性流動，其連續(xù)性和動量方程的無量綱形式可表示如下［16］：

其中，u為流速矢量，p為流體壓力，Re為流體雷諾數(shù)。

1.2 數(shù)值模擬

采集本研究病例的術(shù)后胸部CT圖像AD數(shù)據(jù)，在高性能計(jì)算機(jī)上應(yīng)用MIMICS軟件進(jìn)行三維重構(gòu)，經(jīng)過自動化閾值分割和手動分離，刪除所需模型以外的組織，得到包括主動脈根部、升主動脈、主動脈弓、降主動脈以及主動脈弓上部的3大分支（頭臂干、左頸總動脈、左鎖骨下動脈），再進(jìn)行模型表面光順化，獲得術(shù)后流場CAD模型。

劃分網(wǎng)格是流體數(shù)值模擬的重要一步，精細(xì)的網(wǎng)格使計(jì)算準(zhǔn)確性較高，但會增加計(jì)算機(jī)硬件成本和時(shí)間成本。將CAD模型文件導(dǎo)入網(wǎng)格劃分軟件，如圖1所示，劃分高度復(fù)雜幾何體網(wǎng)格有如下3步主要步驟：（1）幾何預(yù)處理，將幾何體抽殼后留下壁面塊狀曲面，將壁面塊狀曲面之間的共享邊轉(zhuǎn)化為抑制邊，可實(shí)現(xiàn)網(wǎng)格劃分時(shí)壁面塊狀曲面被視為整體連續(xù)曲面。（2）在連續(xù)曲面的基礎(chǔ)上進(jìn)行2D網(wǎng)格劃分和網(wǎng)格質(zhì)量控制。（3）在2D網(wǎng)格的基礎(chǔ)上進(jìn)行CFD網(wǎng)格劃分，并通過網(wǎng)格質(zhì)量檢查。邊界層網(wǎng)格為三菱柱單元，非邊界層網(wǎng)格為四面體單元［17］。為了在保證良好收斂性的同時(shí)降低計(jì)算量，選用如圖2所示的網(wǎng)格模型，共計(jì)單元數(shù)為292 012。

圖1 主動脈流域幾何模型Fig.1 Geometric model of aortic drainage area

圖2 主動脈流域網(wǎng)格模型Fig.2 Mesh model of aortic drainage area

將上述網(wǎng)格文件導(dǎo)入高性能計(jì)算機(jī)的ABAQUS。主要前處理設(shè)置如下：（1）應(yīng)用幅值曲線命令（Amplitude）的Smooth類型，將患者術(shù)后主動脈進(jìn)口實(shí)測流速進(jìn)行瞬態(tài)加載，如圖3所示。（2）根據(jù)文獻(xiàn)［18］，血液密度設(shè)定為1 060 kg/m3，血液運(yùn)動粘度設(shè)定為3.71 mPa·s。（3）所有壁面假定為無滑移邊界條件［19］。（4）本文的流體域設(shè)定為不可壓縮牛頓流體［20］。（5）根據(jù)流動充分發(fā)展條件將主動脈弓3大分支出口和主動脈出口設(shè)為自由出流。（6）定義雷諾數(shù)Re=ρvd/μ，其中ρ為血液密度，v為主動脈進(jìn)口流速，d為進(jìn)口截面的特征長度，μ為血液運(yùn)動粘度，得平均雷諾數(shù)為9 900，可知該流動存在湍流。使用上述湍流模型，進(jìn)行10個(gè)心動周期的計(jì)算，共計(jì)6.6 s，采用自動分析步長，得到穩(wěn)定收斂解，取第10個(gè)周期進(jìn)行分析。

圖3 主動脈進(jìn)口流速曲線Fig.3 Aortic inlet velocity

2 結(jié)果與討論

經(jīng)過10個(gè)周期的模擬計(jì)算，獲得主動脈壁面壓力云圖和主動脈縱向剖面流速矢量云圖。取第10個(gè)周期6個(gè)關(guān)鍵時(shí)刻的計(jì)算結(jié)果進(jìn)行分析，即t1=0 s、t2=0.105 6 s、t3=0.217 8 s、t4=0.231 0 s、t5=0.257 4 s、t6=0.369 6 s，這些時(shí)刻分別是開始射血時(shí)刻、射血早期時(shí)刻、射血加速最快時(shí)刻、射血峰值時(shí)刻、射血速度平穩(wěn)時(shí)刻、射血降速最快時(shí)刻。

2.1 血管壁面壓力分布

6個(gè)關(guān)鍵時(shí)刻的壁面壓力云圖如圖4所示。主動脈壁面壓力與邊界幾何形狀和血流速度密切相關(guān)。在心臟收縮期，研究對象的壁面壓力隨著入口血流速度增大而迅速升高。在心跳周期不同時(shí)間點(diǎn)，流域壁面壓力整體呈現(xiàn)從主動脈近心端到遠(yuǎn)心端的下降趨勢，局部壓力有差異分布。并且在流速較高時(shí)段，由于壁面幾何改變，血流流動使升主動脈流域和主動脈弓流域的外側(cè)壓力分布明顯高于內(nèi)側(cè)。當(dāng)入口血流速度到達(dá)0.231 s的最高流速2 m/s時(shí)，流域整體壁面壓力到達(dá)峰值。隨后進(jìn)入舒張期，整體壁面壓力伴隨入口血流速度的減少而降低。

圖4 主動脈壁面壓力分布（Pa）Fig.4 Wall pressure distribution at the aorta(Pa)

2.2 流域流速分布

6個(gè)關(guān)鍵時(shí)刻的主動脈縱向剖面流速矢量云圖如圖5所示。本例的入流流速存在病理性，從最低流速0.031 7 m/s到峰值流速2 m/s的循環(huán)中呈現(xiàn)非線性上升下降，入流平均流速為1.008 m/s。在心動收縮期初期，流域整體流速不高，升主動脈與主動脈弓流域流速較小，并伴隨著如圖所示的二次流，分別為升主動脈中部和降主動脈起始處內(nèi)側(cè)的二次回流和二次環(huán)流。心臟收縮期屬于快速射血期，隨著入口流速的提高，流域整體流速逐漸提高，升主動脈內(nèi)側(cè)流速與主動脈弓流速漲幅較小，且存在由于壁面復(fù)雜幾何和二次流狀態(tài)而表現(xiàn)出的局部性和波動性，升主動脈的二次流出現(xiàn)了向主動脈內(nèi)側(cè)偏移的趨勢。到心室最大收縮期0.231 s，入口流速到達(dá)峰值時(shí)，流域整體流速保持上升，不同位置點(diǎn)的流速高峰時(shí)刻不同，這說明流場流速與入口流速沒有同步性。當(dāng)進(jìn)入心室減慢射血期和舒張期，流域流速持續(xù)降低，升主動脈渦流和回流重新回靠。在幾何關(guān)系上，渦流較多出現(xiàn)在局部幾何復(fù)雜或腔體狹窄的動脈中，如心臟瓣膜、靜脈瓣膜、血管分支處以及血管急轉(zhuǎn)彎的內(nèi)側(cè)邊，可知本例存在兩處典型急轉(zhuǎn)彎，為主動脈根部與升主動脈交界處內(nèi)側(cè)和主動脈弓與降主動脈交界處內(nèi)側(cè)。流速而論，流速對于渦流形成的影響呈現(xiàn)雙向作用。當(dāng)血流流速較低時(shí)，血液粘滯性可致流動阻力較大，軸流不明顯而邊流作用增加，渦流形成可能性提高。當(dāng)流速較高時(shí)，血流剪切力增加且血流沖擊內(nèi)壁，也可致渦流出現(xiàn)的幾率上升。

圖5 主動脈縱向剖面流速矢量分布（m/s）Fig.5 Blood velocity vector distribution at aortic vertical section(m/s)

3 結(jié)論

血流力學(xué)行為與模型幾何形狀和進(jìn)出口邊界條件息息相關(guān)。本文介紹了ABAQUS/CFD在主動脈血流數(shù)值模擬上的應(yīng)用，得到主動脈壁面壓力云圖和剖面流速云圖。相比健康成年人體，本例的壁面壓力分布和流速分布整體上偏大，需要術(shù)后進(jìn)一步治療調(diào)理。在不同時(shí)刻，從近心端到遠(yuǎn)心端，壁面壓力整體呈現(xiàn)下降趨勢，局部壓力有差異分布，較大壓力位置為主動脈根部內(nèi)側(cè)和升主動脈起始處外側(cè)。二次流的主要分布位置為升主動脈內(nèi)側(cè)和降主動脈起始處內(nèi)側(cè)。心血管疾病的發(fā)病機(jī)理和病理過程與血流力學(xué)環(huán)境的異常密不可分。局部力學(xué)屬性的突變，包括血液的湍流、二次環(huán)流、二次回流、血管內(nèi)壁局部高壓等，會使血細(xì)胞、脂粒等附壁堆積（或?qū)е聞用}粥樣硬化類疾?。⒀軆?nèi)膜受損、血管局部擴(kuò)張等現(xiàn)象發(fā)生。本例的血流數(shù)值模擬可為人工血管置換術(shù)治療Stanford A型主動脈夾層后的病情發(fā)展研究提供參考。