顱內橫竇靜脈血流動力學計算機仿真

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

（1.內蒙古民族大學 數理學院，內蒙古 通遼 028043；2.內蒙古民族大學 附屬醫院，內蒙古 通遼 028007）

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

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

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

圖1 腦靜脈Fig.1 Cerebral vein

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

1 材料和方法

1.1 三維模型幾何構建

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

圖2 血管模型Fig.2 Vascular models

1.2 計算模型與方法

1.2.1 計算模型網格劃分

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

1.2.2 控制方程與邊界條件

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

圖3 網格劃分模型Fig.3 Grid division model

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

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

1.2.3 計算數值模擬

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

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

2 計算結果

2.1 血液流場分布

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

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

2.2 血管壁面壓力分布

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

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

2.3 壁面切應力分布

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

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

3 討論

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

通過對比圖5（a）入口速度為0.12 m·s-1和圖5（b）入口速度為0.22 m·s-1的血管壁面壓力數值模擬結果可得，從模型入口到出口血管壁面壓力逐漸遞減，在圖5（a）中弧狀彎曲和乙狀彎曲處壓力值較小，屬于低壓力區域，圖5（b）相對于圖5（a）血液入口速度增大，弧狀彎曲和乙狀彎曲雖然血管壁面壓力明顯增大，但乙狀彎曲（箭頭1所指）屬于低壓力區域，弧狀彎曲（箭頭2所指）屬于高壓力區域。高壓力環境（圖5（b）中箭頭2所指）會導致血管結構改變，損傷壁面［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的血管壁面切應力數值模擬結果可得，血液入口速度為0.12 m·s-1和0.22 m·s-1的分布規律相似，從切應力空間分布來看，切應力分布與血液流速和血管管徑有關。圖6（a）弧狀彎曲（箭頭2所指）和乙狀彎曲（箭頭1所指）內外側的壁面切應力值都很小，均為低切應力區域，圖6（b）相對于圖6（a）血液入口速度增大，乙狀彎曲（箭頭1所指）內側切應力明顯大于外側切應力，而弧狀彎曲（箭頭2所指）外側切應力明顯大于內側切應力，隨著血液入口速度的增大，弧狀彎曲和乙狀彎曲區域產生高切應力區域，橫竇靜脈持續存在高壁面切應力區域會增大血液與血管的摩擦，損傷血管內膜。橫竇靜脈持續存在低壁面切應力區域，血管內脂類顆粒及其他物質容易滯留、堆積，造成血管局部狹窄，易發生病變。

4 結論

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