左 濤，彭紅梅，張東威，劉寶治，趙洪明

（1.內(nèi)蒙古民族大學 數(shù)理學院，內(nèi)蒙古 通遼 028043；2.內(nèi)蒙古民族大學 附屬醫(yī)院，內(nèi)蒙古 通遼 028007）

近年來，人口老齡化水平加劇，根據(jù)《中國心血管健康與疾病報告2020概要》［1］，中國居民患腦血管疾病人群的數(shù)目和死亡率在增加，由此可見，腦血管疾病對人類健康危害越來越大。各發(fā)達國家從上世紀末開始極力重視從血流動力學的角度研究腦血管疾病的成因。目前我國也在大力支持和發(fā)展基礎醫(yī)學研究，利用計算機模擬血流動力學狀況代替血流在體實驗也取得了很大的成效，但是多數(shù)研究為腦動脈血管血流狀況，對于腦靜脈血管血流狀況的研究較少，因此有必要對其進行研究。

大腦靜脈包括：上矢狀竇、下矢狀竇、直竇、竇匯、橫竇、乙狀竇、海綿竇、巖上竇及巖下竇。下矢狀竇向下連接直竇，直竇與上矢狀竇向下匯入竇匯，海綿竇連接巖下竇和巖上竇（圖1）。筆者所研究的橫竇是顱內(nèi)最大靜脈竇之一，也是顱內(nèi)靜脈竇構(gòu)造連接的重要樞紐。

靜脈橫竇又名側(cè)竇或水平竇，位于枕骨內(nèi)面的橫竇溝內(nèi)。橫竇（圖1）起于竇匯，延續(xù)為乙狀竇，接受巖上竇和上矢狀竇回流的血液［2］。因其部位存在弧狀和乙狀2種形狀彎曲，其血流環(huán)境復雜，臨床發(fā)現(xiàn)，在此處易發(fā)生顱內(nèi)靜脈竇血栓、靜脈竇狹窄等腦血管疾?。?］。因此對靜脈竇從血流動力學角度觀察橫竇處的血流狀況，以分析出現(xiàn)血管堵塞、狹窄等疾病的原因，從而更好地為臨床診治和醫(yī)務人員提供有效信息。

圖1 腦靜脈Fig.1 Cerebral vein

應用計算流體力學的方法，采用內(nèi)蒙古民族大學附屬醫(yī)院神經(jīng)內(nèi)科提供的二維CT數(shù)據(jù)，利用醫(yī)學建模軟件MIMICS21.0進行血管重建，將得到橫竇三維幾何模型導入Geomagic Studio中進行切割分離和光滑處理，利用ICEM CFD進行網(wǎng)格劃分，將模型結(jié)果導入FLUENT2020 R1中進行數(shù)值模擬，觀察模擬結(jié)果顯示在不同血液入口速度下血液流線分布、血液壓力分布、血液壁面切應力分布［3-4］，通過模擬結(jié)果顯示的血流動力學特性可以分析患者橫竇的血流動力學原因。同時對橫竇靜脈血流狀況進行計算機仿真，可以數(shù)值模擬分析出在顱內(nèi)靜脈竇處發(fā)生血栓、狹窄的血流動力學原因，為顱內(nèi)腦血管疾?。ㄈ缒X卒中）的診斷和治療提供幫助。

1 材料和方法

1.1 三維模型幾何構(gòu)建

數(shù)據(jù)獲取的圖像采集方法為計算機斷層掃描技術(shù)（CTA），圖像數(shù)據(jù)參數(shù)為：平面分辨率為512×512，像素大小為0.488 3 mm，層間距為0.900 mm，共256張切片。醫(yī)學影像三維重建技術(shù)是指將二維影像數(shù)據(jù)轉(zhuǎn)換成三維可視化圖像的技術(shù)，方便研究人員和醫(yī)務人員對圖像進行分析［4］。將CT數(shù)據(jù)（DICOM格式）導入Mimics后，首先閾值提取輪廓：選取合適的閾值范圍使得輪廓最清晰，閾值范圍為267~2 378，形成蒙面后進行3D計算得到骨骼、血管及組織模型（圖2（a）），通過分離、分割去除骨骼和組織，從模型中提取出橫竇靜脈（圖2（b）），對模型進行邊緣切割、去除冗余結(jié)構(gòu)，在3-matic中對其進行后處理，獲得較為光滑的3D血管模型（圖2（c））。

圖2 血管模型Fig.2 Vascular models

1.2 計算模型與方法

1.2.1 計算模型網(wǎng)格劃分

網(wǎng)格劃分是計算流體力學的基礎，網(wǎng)格的質(zhì)量直接影響后面在fluent軟件的數(shù)值模擬。選用在ICEM CFD軟件中設定邊界層網(wǎng)格，生成面網(wǎng)格。圖（3）為橫竇血管的網(wǎng)格劃分模型圖，在模型中共生成了553 975個網(wǎng)格和205 830個節(jié)點。

1.2.2 控制方程與邊界條件

將橫竇靜脈中的血液假設為不可壓縮牛頓流體，并且血液滿足絕熱、均勻等條件，血液流動滿足的（Navier-Stokes）偏微分方程如下［5］：

圖3 網(wǎng)格劃分模型Fig.3 Grid division model

其中，t為時間，U為流體的速度矢量，P為流場壓強，血液密度ρ=1 050 kg·m-3，血液黏性系數(shù)μ=0.003 5 Pa·s。

設定計算模型為穩(wěn)態(tài)流動，顱內(nèi)靜脈竇的血流入口速度設為0.12 m·s-1和0.22 m·s-1［6］。血管的入口直徑d=8.56 mm，根據(jù)雷諾數(shù)公式：，在血液入口速度為0.12 m·s-1時，血液的雷諾數(shù)為313，在血液入口速度為0.22 m·s-1時，血液的雷諾數(shù)為565，雷諾數(shù)Re<2 300，判斷血液流動為層流［4］。在計算過程中，設置出口壓強為0 Pa，假設血管壁剛性無滑移，不考慮重力的影響。

1.2.3 計算數(shù)值模擬

使用計算流體力學軟件Fluent采用有限體積法進行數(shù)值模擬，有限體積法是將模型區(qū)域劃分為一系列控制區(qū)域，每個控制區(qū)域選一個節(jié)點作為代表，通過守恒型的控制方程，對控制體積做積分導出離散方程。

計算收斂條件時，將殘差設置為10-5，最大迭代次數(shù)設置為200。分別設置血流入口速度為0.12 m·s-1和0.22 m·s-1，計算顱內(nèi)橫竇靜脈的血液流場分布、壁面切應力及壁面壓力分布等血流動力學特性。同時可以將Fluent軟件數(shù)值模擬的結(jié)果代入CFD-post軟件進行可視化處理，得到橫竇血管的血液流場分布、壁面壓力及壁面切應力等云圖。

2 計算結(jié)果

2.1 血液流場分布

本研究選用健康成年人腦CT數(shù)據(jù)，得到的結(jié)果為正常血液流動狀態(tài)分布，流場的空間分布見圖4。圖4（a）是入口速度為0.12 m·s-1血液在靜脈橫竇中血液流場分布圖，除血管彎曲區(qū)域血液流動模式為直線流動，在橫竇靜脈乙狀彎曲處（箭頭1所指）出現(xiàn)明顯渦旋流動現(xiàn)象，并且血液流速很低，血流方向錯綜復雜，在橫竇靜脈弧狀彎曲處（箭頭2所指）其流動模式為螺旋流動，血液流速緩慢且流動不規(guī)律，箭頭（1、2）所指處均出現(xiàn)低流速區(qū)，血液流速在此區(qū)域存在最小值?；顝澢毫魅攵耍^3所指）血液流速存在最大值。圖4（b）是入口速度為0.22 m·s-1血液在靜脈橫竇中血液流場分布圖，在靜脈橫竇中血液流動模式大部分為層流流動，在橫竇靜脈乙狀彎曲處（箭頭1所指）血液流動不規(guī)律并伴有紊流現(xiàn)象，但血液流速明顯加快，在橫竇靜脈弧狀彎曲處（箭頭2所指）血液流速增加且流動不規(guī)律，血管彎曲處（箭頭1、箭頭2處）低流速區(qū)域減小，但仍為低流速區(qū)域，在此區(qū)域存在最小速度，弧狀彎曲血液流入端（箭頭3所指）血液流速存在最大值。

圖4 血液流場分布Fig.4 Blood flow field distribution

2.2 血管壁面壓力分布

本例CT數(shù)據(jù)模擬顯示的顱內(nèi)橫竇靜脈血管壁面壓力分布圖見圖5。圖5（a）是入口速度為0.12 m·s-1時血管壁面壓力分布圖，血液從模型入口到出口壁面壓力逐漸遞減，在壁面壓力云圖中各個壓力值范圍分布規(guī)則，壓力差（圖中最大壓力與最小壓力的差值）范圍為0~81.250 Pa，在乙狀彎曲處（箭頭1所指）壓力值較小，屬于低壓力區(qū)，在弧狀彎曲處（箭頭2所指）雖然血管壁面壓力值相對于乙狀彎曲處較大，但依然屬于低壓力區(qū)。圖5（b）是入口速度為0.22 m·s-1時的血管壁面壓力分布圖，在云圖中各壓力值分布范圍規(guī)則，壓力差范圍為0~130.000 Pa，從模型中可以看出血管壁面壓力從入口到出口逐級遞減，在乙狀彎曲處（箭頭1所指）增大血液入口速度，血管壁面壓力明顯增大，但仍屬于低壓力區(qū)，在弧狀彎曲處（箭頭2所指）增大血流入口速度，血管壁面壓力顯著增大，屬于高壓力區(qū)。如圖5（b）所示，由于血管彎曲處血液流動不規(guī)律并伴有紊流現(xiàn)象，血液呈三維螺旋流動，相對于附近直線流動血管壁面壓力較大。

圖5 血管壁面壓力分布圖Fig.5 Distribution diagram of vascular wall pressure

2.3 壁面切應力分布

血管壁面切應力又稱內(nèi)皮剪切應力，是腔內(nèi)流動的黏性血流對其流經(jīng)的血管壁表面產(chǎn)生的摩擦力［7］。本例CT數(shù)據(jù)模擬顱內(nèi)橫竇靜脈血管壁面切應力分布見圖6。圖6（a）為入口速度為0.12 m·s-1時血管壁面切應力分布圖，血管乙狀彎曲（箭頭1所指）壁面切應力很小，屬于低切應力區(qū)域，弧狀彎曲（箭頭2所指）的內(nèi)側(cè)和外側(cè)壁面切應力都很小，都屬于低切應力區(qū)域。圖6（b）為入口速度為0.22 m·s-1時血管壁面切應力分布圖，橫竇靜脈乙狀彎曲（箭頭1所指處）內(nèi)外側(cè)壁面切應力明顯增大，屬于高壁面切應力區(qū)域?；顝澢^2所指處）外側(cè)壁面切應力變化明顯，存在高切應力區(qū)域，內(nèi)側(cè)壁面切應力變化很小，屬于低切應力區(qū)域，外側(cè)壁面切應力明顯大于內(nèi)側(cè)。

圖6 血管壁面切應力分布Fig.6 Distribution of vascular wall shear stress

3 討論

通過對比圖4（a）入口速度為0.12 m·s-1和圖4（b）入口速度為0.22 m·s-1的血液流場數(shù)值模擬結(jié)果可得，弧狀彎曲和乙狀彎曲區(qū)域血液流動環(huán)境復雜，圖4（a）中弧狀彎曲和乙狀彎曲位置（箭頭1、箭頭2所指）血液流速很低，均為低流速區(qū)域，圖4（b）相對于圖4（a）血液入口速度增大，血液在此位置（箭頭1、箭頭2所指）血液流速雖然有增大的現(xiàn)象，但是仍然屬于低流速區(qū)域。流速持續(xù)很低，血液中血小板和脂類等成分附壁堆積，易造成血管狹窄，影響血液流動［8］。當增大血液入口速度時，在弧狀彎曲和乙狀彎曲位置（箭頭1、箭頭2所指）血液流速明顯增大，但仍有低流速區(qū)域。

通過對比圖5（a）入口速度為0.12 m·s-1和圖5（b）入口速度為0.22 m·s-1的血管壁面壓力數(shù)值模擬結(jié)果可得，從模型入口到出口血管壁面壓力逐漸遞減，在圖5（a）中弧狀彎曲和乙狀彎曲處壓力值較小，屬于低壓力區(qū)域，圖5（b）相對于圖5（a）血液入口速度增大，弧狀彎曲和乙狀彎曲雖然血管壁面壓力明顯增大，但乙狀彎曲（箭頭1所指）屬于低壓力區(qū)域，弧狀彎曲（箭頭2所指）屬于高壓力區(qū)域。高壓力環(huán)境（圖5（b）中箭頭2所指）會導致血管結(jié)構(gòu)改變，損傷壁面［8-9］。圖5（a）血管壁面壓力差范圍為0~81.250 Pa，圖5（b）血管壁面壓力差范圍為0~130.000 Pa，圖5（b）相對于圖5（a）血流入口速度增大，血管壁面壓力差范圍在變大。

通過對比圖6（a）入口速度為0.12 m·s-1和圖6（b）入口速度為0.22 m·s-1的血管壁面切應力數(shù)值模擬結(jié)果可得，血液入口速度為0.12 m·s-1和0.22 m·s-1的分布規(guī)律相似，從切應力空間分布來看，切應力分布與血液流速和血管管徑有關(guān)。圖6（a）弧狀彎曲（箭頭2所指）和乙狀彎曲（箭頭1所指）內(nèi)外側(cè)的壁面切應力值都很小，均為低切應力區(qū)域，圖6（b）相對于圖6（a）血液入口速度增大，乙狀彎曲（箭頭1所指）內(nèi)側(cè)切應力明顯大于外側(cè)切應力，而弧狀彎曲（箭頭2所指）外側(cè)切應力明顯大于內(nèi)側(cè)切應力，隨著血液入口速度的增大，弧狀彎曲和乙狀彎曲區(qū)域產(chǎn)生高切應力區(qū)域，橫竇靜脈持續(xù)存在高壁面切應力區(qū)域會增大血液與血管的摩擦，損傷血管內(nèi)膜。橫竇靜脈持續(xù)存在低壁面切應力區(qū)域，血管內(nèi)脂類顆粒及其他物質(zhì)容易滯留、堆積，造成血管局部狹窄，易發(fā)生病變。

4 結(jié)論

本研究通過醫(yī)學影像軟件MIMICS21.0基于CT數(shù)據(jù)成功建立了顱內(nèi)橫竇的三維血管模型，將得到的三維血管模型代入到計算流體力學軟件FLUENT中進行血流動力學數(shù)值模擬，同時計算出在入口速度為0.12 m·s-1和0.22 m·s-1時的橫竇靜脈的血液流場、壁面壓力、壁面切應力等血流動力學特征，數(shù)值模擬分析顯示在靜脈橫竇弧狀彎曲和乙狀彎曲處血流速度緩慢，出現(xiàn)渦旋流動等血流現(xiàn)象，在此區(qū)域壁面切應力值較小，為低切應力區(qū)，此現(xiàn)象說明了在靜脈橫竇彎曲處容易發(fā)生淤滯，血管內(nèi)流通區(qū)域減少，影響血液正常流動，增大血液入口速度后，弧狀彎曲處出現(xiàn)高壓力區(qū)域，高壓力環(huán)境會導致血管結(jié)構(gòu)改變，損傷血管壁面，在乙狀彎曲內(nèi)側(cè)和弧狀彎曲外側(cè)出現(xiàn)高壁面切應力區(qū)域，高壁面切應力容易使血管受損。橫竇由于其構(gòu)造的獨特性，在彎曲區(qū)域血流環(huán)境復雜，極易發(fā)生疾病，通過計算機仿真模擬真實地反映了橫竇靜脈血管內(nèi)血液流動，可以幫助醫(yī)務人員更清楚地認識橫竇靜脈的血流動力學特點，通過血流動力學參數(shù)對患病風險進行評估，為顱內(nèi)橫竇疾病的治療提供理論支持。