引用本文： 韓露, 王偉, 劉金龍, 俞曉青, 蔡及明. 碟形血泵內部流場的數值仿真與優化. 中國胸心血管外科臨床雜志, 2016, 23(4): 380-385. doi: 10.7507/1007-4848.20160088 復制
目前心血管疾病數量約占所有疾病的30%，預計到2020年，該數值將增長到40%。伴隨復雜先天性心臟病外科技術的迅速發展，術后嚴重低心排血量綜合征的發生已成為影響手術預后的關鍵。部分患者雖經使用正性肌力藥物，仍不能成功脫離體外循環，需要機械輔助裝置進行心臟輔助循環支持，以使心肌功能得到恢復[1-2]。
我們應用計算機流體仿真技術（CFD）對一種兒童用碟形血泵內的流場進行數值模擬，通過分析數值計算結果對其產生的血栓和溶血現象進行系統的分析研究[3-5]，最后根據分析結果對血泵結構方面的不足進行了針對性優化。
1 材料與方法
1.1 計算機輔助建模
CFD計算的研究對象是血泵內部的流體，因此采用UG NX 8.0軟件對血泵內部的流體區域進行了幾何建模。血泵流體域的CFD模型如圖 1所示。
圖1
血泵CFD計算幾何模型
由于血泵是旋轉機械，其內部流體存在相對運動的區域。因此，需要對幾何模型進行合理的劃分，選擇合適的動靜交界面，將流體域劃分為4個區域：葉輪上下蓋板所包圍的流體區域，impeller域；血泵進口流入旋轉葉輪中心和葉輪上蓋板上方間隙的流體區域，top域；從葉輪中心流出，流經葉輪地盤下方間隙流體域，bottom域；經葉輪甩出后，進入壓水室，然后從血泵出口流出的流體，duct域。動靜區域的交界面采用frozen rotor邊界條件[6-8]。為了與試驗結果進行更準確的對比，對模型的進、出口分別延長了400 mm（圖 1中未顯示延長段）。
1.2 網格劃分
采用ICEM v14.5對各個區域進行網格劃分。由于top和bottom域的幾何尺度較小，為了提高網格劃分質量和計算精度，采用結構化網格。impeller和duct域采用非結構化網格。各個區域的網格如圖 2所示。進出口延長段也采用結構化網格。
圖2
生成網格
注：A為top域結構化網格；B為bottom域結構化網格；C為top域結構化網格；D為duct域非結構化網格
1.3 邊界條件
根據試驗工況，選擇參考壓力為1個大氣壓。試驗中，進口為儲血瓶自由液面，因此給定進口邊界條件為總壓邊界，值為0 atm。出口設為流量出口，根據試驗測得的流量給定不同轉速工況下的邊界值。計算工質為血液，忽略血液中的懸浮物質，假設其為不可壓的牛頓流體[9]，其密度為ρ=1056 kg/m3，動力粘度為μ=3.5×10-3 Pa·s。旋轉域impeller與其它3個區域均有接觸面，將這些面設置為interface邊界，采用frozen rotor模型。進出口延長段和血泵原進出口接觸面也設置為interface，因為沒有轉靜變換，因此不采用frozen rotor模型[10-12]。
1.4 計算
采用商業CFD軟件ANSYS?-FLUENT 12.0 計算泵內穩流情況，選用SIMPLE算?法的壓力-速度耦合方式，隱式分離式求解器(implicit segregated solver) [13]，以欠松弛方式求解迭代過程，欠松弛因子設置為：壓力0.3，體積力1.0，密度1.0，動量0.7 [14-15]。求解器對殘差收斂控制在10-5次方，運算中監控流體的x、y、z殘差曲線進行以及連續性參數，最終計算均收斂[16-17]。
2 結果
2.1 流量特性曲線
選取2 500 rpm、3 000 rpm和4 000 rpm三個不同轉速對血泵流場進行了CFD模擬，每個轉速下選取5個不同的流量出口邊界條件，計算血泵出口的靜壓，獲得了不同流量出口邊界時的原型機流場特征和流量特性曲線。將CFD模擬得到的出口靜壓與試驗測量結果進行了比較，見圖 3。結果表明，CFD仿真結果和試驗結果非常吻合。
圖3
流量特性
2.2 流場分析
圖 4給出了2 500 rpm轉速，不同流量工況下的血泵葉輪內部流線圖（y=-0.008 m）。葉輪中心是從血泵進口進來的流體，尚未得到葉輪的加速，所以速度很低。特別是在小流量（Q=2 lpm）的時候，葉輪中心的速度明顯低于其他地方。隨著流體向外圈流動，逐漸被加速，在葉片表面存在流動分離。從圖 4中可以看出，在不同流量下，葉輪內部均出現明顯的流動分離。隨著流量的增加，分離區域的范圍有減小趨勢，表明分離隨著流量的增加而有所改善。顯然，流動分離應該盡量避免，這一問題將在后面的優化中得到處理。隨著轉速的增加3 000 rpm和4 000 rpm，Q=4，10 lpm時的葉輪內部流場進行了比較。從圖 5、6上看出，葉輪內部還是存在流動分離，10 lpm流量下流場分離較4 lpm時略有改善。
圖4
通道內流線圖
注：A為Q=2 lpm；B為Q=4 lpm；C為Q=10 lpm
圖5
葉輪通道流線圖（3 000 rpm)
注：A為Q=4 lpm；B為Q=10 lpm
圖6
葉輪通道流線圖（4 000 rpm)
注：A為Q=4 lpm；B為Q=10 lpm
獲得不同流量工況下的血泵葉輪內部剪切應變率云圖。圖 7給出了2 500 rpm轉速下的切應力云圖，從圖中可以看出，剪切較強的區域集中在葉輪和壓水室交界的區域和葉片的頭部區域。在Q=2 lpm工況，壓水室離血泵出口較近的地方出現大面積高剪切應變率區域。這一現象隨著流量的增加很快消失，在Q=4 lpm及以上工況均未發現大面積的高剪應變區域，且隨著流量的增加，葉輪內部的剪切應變未發現明顯增加。
圖7
通道剪切應變率圖
注：A為Q=2 lpm；B為Q=4 lpm；C為Q=10 lpm
圖 8血泵壁面剪切力分布云圖。從數值上看，隨著流量的增加，剪切力的最大值并無明顯變化，基本維持在100 Pa左右。剪切力強弱分布和剪切應變率一致，高剪切力出現在壓水室四周和血泵進口進入到葉片腔室的頸部位置。3 000 rpm和4 000 rpm，Q=4，10 lpm時血泵壁面剪切力分布。血泵壁面剪切力在相同轉速時隨流量變化不大，而隨著轉速的提高，剪切力急劇增大。
圖8
壁面剪切力
注：A為Q=2 500 rpm，10 lpm；B為Q=Q=3 000 rpm，10 lpm；C為Q=Q=4 000 rpm，4 lpm；D為Q=Q=4 000 rpm，10 lpm
選擇轉速2 500 rpm，流量10 lpm工況下的血泵流場計算結果進行詳細分析。圖 9所示是血泵進口截面（z=0）的流線圖。流體從血泵進口進入，經葉輪加速后進入壓水室，然后從出口流出。從圖中可以很明顯地看出，旋轉的葉輪上下兩端，即top和bottom域內，存在著大尺度流動分離區，形成一個流動死區。流體在這兩個區域內的通流能力很弱，因此絕大部分的流體將繞過這一區域，從葉輪之間的通道流過。這一發現與我們之前在網格敏感性分析時的預測一致。流體經過葉輪通道加速后，被葉片甩出，進入壓水室，由于流體具有很強的旋轉速度，因此在壓水室內的流體形成很強的漩渦。
圖9
進口截面流線圖
流體通過葉片通道被逐漸加速，因此流體的總壓逐漸增加，葉輪旋轉的機械能轉化為流體的壓力勢能。圖 10所示為進口截面的湍動能分布云圖。從圖中可以看出，在旋轉部件的進口位置，特別是在易形成臺階流動的區域，這些地方流動分離嚴重，湍動能大，是能量高損失區域，在后期應予以關注。葉輪通道內的湍動能云圖表明，葉輪前緣的湍動能較高，表明流體在這一區域的損失較大。應力較高的區域集中在葉片的前緣和葉輪的四周，這些區域都是轉靜區域交界的地方，流體存在很大的速度梯度，因此剪切應變率較大。
圖10
進口截面湍動能云圖
3 討論
根據原型機CFD數值模擬的結果，從流線、壓力分布和湍動能分布等云圖看出，原型機流場的葉輪通道間存在較大的流動分離，局部地方流動損失較大，對血泵的性能造成不利的影響。將根據流場分析的結果，對原型機結構進行局部修改和優化，以改善血泵流場品質、提高醫學性能。
3.1 葉片數
由于葉片數少，葉柵稠度低，導致葉輪進出口面積變化較大，葉輪通道間擴壓劇烈，導致流動分離。進一步觀察葉輪通道間的速度矢量圖，可以看到，分離主要發生在葉背前緣。通過選擇合理的葉片數，可獲得合理的葉柵稠度，能夠有效控制葉輪通道間的流動分離。同時，由于分離嚴重導致流場紊亂，流動剪切增加，血泵內部的剪切與流動分離關系密切，分離越嚴重，剪切越強，反之越小。
經過計算，隨著葉片數的增加，血泵出口的靜壓值不斷提高，表明血泵的增壓能力隨著葉片數的增加而不斷增強。考慮到葉片數的增加，也會帶來葉型損失等的增加，因此從血泵的效率上看，當葉片數為6片時，血泵的效率最佳，所以優化機型采用6葉片方案。
3.2 攻角
葉片前緣的速度矢量圖，反應了葉片的攻角特性。由于來流速度矢量并未正對葉片前緣，因此在葉背形成了一個分離區。通過調整葉片的安裝角，可以改變葉片的攻角，從而改善葉片的攻角特性，使葉背的流動得到極大改善，流動分離區域消失。
3.3 曲率
血泵出口的流線和速度分布。出口流場有很強的旋轉，同時局部出現速度大的區域，造成局部剪切增強。血泵出口管與壓水室腔通過圓弧連接，改變圓弧的曲率半徑，有助于改善出口流體的預旋。優化將朝對壓力損失最有利、出口靜壓最大的方向進行，同時出口管道的速度梯度減小，有利于減小剪切。
3.4 總結
根據血泵原型機的幾何建模，對蝶形血泵流場進行了CFD模擬和分析。通過對碟形血泵的特性和流場分析，并與試驗結果進行對比，表明了測繪和建模過程的準確性，根據分析結果，尋找優化血泵的多個影響因素進行分析，為最終的優化血泵模型確定計算方法的合理性。小兒碟形離心式血泵具有較好的表現，葉輪壓力分布較為均勻，葉輪剪切力均在溶血性能可控范圍內，表明本設計方案合理可行，同時為本設計的優化提供了可預期的效果。
目前心血管疾病數量約占所有疾病的30%，預計到2020年，該數值將增長到40%。伴隨復雜先天性心臟病外科技術的迅速發展，術后嚴重低心排血量綜合征的發生已成為影響手術預后的關鍵。部分患者雖經使用正性肌力藥物，仍不能成功脫離體外循環，需要機械輔助裝置進行心臟輔助循環支持，以使心肌功能得到恢復[1-2]。
我們應用計算機流體仿真技術（CFD）對一種兒童用碟形血泵內的流場進行數值模擬，通過分析數值計算結果對其產生的血栓和溶血現象進行系統的分析研究[3-5]，最后根據分析結果對血泵結構方面的不足進行了針對性優化。
1 材料與方法
1.1 計算機輔助建模
CFD計算的研究對象是血泵內部的流體，因此采用UG NX 8.0軟件對血泵內部的流體區域進行了幾何建模。血泵流體域的CFD模型如圖 1所示。
圖1
血泵CFD計算幾何模型
由于血泵是旋轉機械，其內部流體存在相對運動的區域。因此，需要對幾何模型進行合理的劃分，選擇合適的動靜交界面，將流體域劃分為4個區域：葉輪上下蓋板所包圍的流體區域，impeller域；血泵進口流入旋轉葉輪中心和葉輪上蓋板上方間隙的流體區域，top域；從葉輪中心流出，流經葉輪地盤下方間隙流體域，bottom域；經葉輪甩出后，進入壓水室，然后從血泵出口流出的流體，duct域。動靜區域的交界面采用frozen rotor邊界條件[6-8]。為了與試驗結果進行更準確的對比，對模型的進、出口分別延長了400 mm（圖 1中未顯示延長段）。
1.2 網格劃分
采用ICEM v14.5對各個區域進行網格劃分。由于top和bottom域的幾何尺度較小，為了提高網格劃分質量和計算精度，采用結構化網格。impeller和duct域采用非結構化網格。各個區域的網格如圖 2所示。進出口延長段也采用結構化網格。
圖2
生成網格
注：A為top域結構化網格；B為bottom域結構化網格；C為top域結構化網格；D為duct域非結構化網格
1.3 邊界條件
根據試驗工況，選擇參考壓力為1個大氣壓。試驗中，進口為儲血瓶自由液面，因此給定進口邊界條件為總壓邊界，值為0 atm。出口設為流量出口，根據試驗測得的流量給定不同轉速工況下的邊界值。計算工質為血液，忽略血液中的懸浮物質，假設其為不可壓的牛頓流體[9]，其密度為ρ=1056 kg/m3，動力粘度為μ=3.5×10-3 Pa·s。旋轉域impeller與其它3個區域均有接觸面，將這些面設置為interface邊界，采用frozen rotor模型。進出口延長段和血泵原進出口接觸面也設置為interface，因為沒有轉靜變換，因此不采用frozen rotor模型[10-12]。
1.4 計算
采用商業CFD軟件ANSYS?-FLUENT 12.0 計算泵內穩流情況，選用SIMPLE算?法的壓力-速度耦合方式，隱式分離式求解器(implicit segregated solver) [13]，以欠松弛方式求解迭代過程，欠松弛因子設置為：壓力0.3，體積力1.0，密度1.0，動量0.7 [14-15]。求解器對殘差收斂控制在10-5次方，運算中監控流體的x、y、z殘差曲線進行以及連續性參數，最終計算均收斂[16-17]。
2 結果
2.1 流量特性曲線
選取2 500 rpm、3 000 rpm和4 000 rpm三個不同轉速對血泵流場進行了CFD模擬，每個轉速下選取5個不同的流量出口邊界條件，計算血泵出口的靜壓，獲得了不同流量出口邊界時的原型機流場特征和流量特性曲線。將CFD模擬得到的出口靜壓與試驗測量結果進行了比較，見圖 3。結果表明，CFD仿真結果和試驗結果非常吻合。
圖3
流量特性
2.2 流場分析
圖 4給出了2 500 rpm轉速，不同流量工況下的血泵葉輪內部流線圖（y=-0.008 m）。葉輪中心是從血泵進口進來的流體，尚未得到葉輪的加速，所以速度很低。特別是在小流量（Q=2 lpm）的時候，葉輪中心的速度明顯低于其他地方。隨著流體向外圈流動，逐漸被加速，在葉片表面存在流動分離。從圖 4中可以看出，在不同流量下，葉輪內部均出現明顯的流動分離。隨著流量的增加，分離區域的范圍有減小趨勢，表明分離隨著流量的增加而有所改善。顯然，流動分離應該盡量避免，這一問題將在后面的優化中得到處理。隨著轉速的增加3 000 rpm和4 000 rpm，Q=4，10 lpm時的葉輪內部流場進行了比較。從圖 5、6上看出，葉輪內部還是存在流動分離，10 lpm流量下流場分離較4 lpm時略有改善。
圖4
通道內流線圖
注：A為Q=2 lpm；B為Q=4 lpm；C為Q=10 lpm
圖5
葉輪通道流線圖（3 000 rpm)
注：A為Q=4 lpm；B為Q=10 lpm
圖6
葉輪通道流線圖（4 000 rpm)
注：A為Q=4 lpm；B為Q=10 lpm
獲得不同流量工況下的血泵葉輪內部剪切應變率云圖。圖 7給出了2 500 rpm轉速下的切應力云圖，從圖中可以看出，剪切較強的區域集中在葉輪和壓水室交界的區域和葉片的頭部區域。在Q=2 lpm工況，壓水室離血泵出口較近的地方出現大面積高剪切應變率區域。這一現象隨著流量的增加很快消失，在Q=4 lpm及以上工況均未發現大面積的高剪應變區域，且隨著流量的增加，葉輪內部的剪切應變未發現明顯增加。
圖7
通道剪切應變率圖
注：A為Q=2 lpm；B為Q=4 lpm；C為Q=10 lpm
圖 8血泵壁面剪切力分布云圖。從數值上看，隨著流量的增加，剪切力的最大值并無明顯變化，基本維持在100 Pa左右。剪切力強弱分布和剪切應變率一致，高剪切力出現在壓水室四周和血泵進口進入到葉片腔室的頸部位置。3 000 rpm和4 000 rpm，Q=4，10 lpm時血泵壁面剪切力分布。血泵壁面剪切力在相同轉速時隨流量變化不大，而隨著轉速的提高，剪切力急劇增大。
圖8
壁面剪切力
注：A為Q=2 500 rpm，10 lpm；B為Q=Q=3 000 rpm，10 lpm；C為Q=Q=4 000 rpm，4 lpm；D為Q=Q=4 000 rpm，10 lpm
選擇轉速2 500 rpm，流量10 lpm工況下的血泵流場計算結果進行詳細分析。圖 9所示是血泵進口截面（z=0）的流線圖。流體從血泵進口進入，經葉輪加速后進入壓水室，然后從出口流出。從圖中可以很明顯地看出，旋轉的葉輪上下兩端，即top和bottom域內，存在著大尺度流動分離區，形成一個流動死區。流體在這兩個區域內的通流能力很弱，因此絕大部分的流體將繞過這一區域，從葉輪之間的通道流過。這一發現與我們之前在網格敏感性分析時的預測一致。流體經過葉輪通道加速后，被葉片甩出，進入壓水室，由于流體具有很強的旋轉速度，因此在壓水室內的流體形成很強的漩渦。
圖9
進口截面流線圖
流體通過葉片通道被逐漸加速，因此流體的總壓逐漸增加，葉輪旋轉的機械能轉化為流體的壓力勢能。圖 10所示為進口截面的湍動能分布云圖。從圖中可以看出，在旋轉部件的進口位置，特別是在易形成臺階流動的區域，這些地方流動分離嚴重，湍動能大，是能量高損失區域，在后期應予以關注。葉輪通道內的湍動能云圖表明，葉輪前緣的湍動能較高，表明流體在這一區域的損失較大。應力較高的區域集中在葉片的前緣和葉輪的四周，這些區域都是轉靜區域交界的地方，流體存在很大的速度梯度，因此剪切應變率較大。
圖10
進口截面湍動能云圖
3 討論
根據原型機CFD數值模擬的結果，從流線、壓力分布和湍動能分布等云圖看出，原型機流場的葉輪通道間存在較大的流動分離，局部地方流動損失較大，對血泵的性能造成不利的影響。將根據流場分析的結果，對原型機結構進行局部修改和優化，以改善血泵流場品質、提高醫學性能。
3.1 葉片數
由于葉片數少，葉柵稠度低，導致葉輪進出口面積變化較大，葉輪通道間擴壓劇烈，導致流動分離。進一步觀察葉輪通道間的速度矢量圖，可以看到，分離主要發生在葉背前緣。通過選擇合理的葉片數，可獲得合理的葉柵稠度，能夠有效控制葉輪通道間的流動分離。同時，由于分離嚴重導致流場紊亂，流動剪切增加，血泵內部的剪切與流動分離關系密切，分離越嚴重，剪切越強，反之越小。
經過計算，隨著葉片數的增加，血泵出口的靜壓值不斷提高，表明血泵的增壓能力隨著葉片數的增加而不斷增強。考慮到葉片數的增加，也會帶來葉型損失等的增加，因此從血泵的效率上看，當葉片數為6片時，血泵的效率最佳，所以優化機型采用6葉片方案。
3.2 攻角
葉片前緣的速度矢量圖，反應了葉片的攻角特性。由于來流速度矢量并未正對葉片前緣，因此在葉背形成了一個分離區。通過調整葉片的安裝角，可以改變葉片的攻角，從而改善葉片的攻角特性，使葉背的流動得到極大改善，流動分離區域消失。
3.3 曲率
血泵出口的流線和速度分布。出口流場有很強的旋轉，同時局部出現速度大的區域，造成局部剪切增強。血泵出口管與壓水室腔通過圓弧連接，改變圓弧的曲率半徑，有助于改善出口流體的預旋。優化將朝對壓力損失最有利、出口靜壓最大的方向進行，同時出口管道的速度梯度減小，有利于減小剪切。
3.4 總結
根據血泵原型機的幾何建模，對蝶形血泵流場進行了CFD模擬和分析。通過對碟形血泵的特性和流場分析，并與試驗結果進行對比，表明了測繪和建模過程的準確性，根據分析結果，尋找優化血泵的多個影響因素進行分析，為最終的優化血泵模型確定計算方法的合理性。小兒碟形離心式血泵具有較好的表現，葉輪壓力分布較為均勻，葉輪剪切力均在溶血性能可控范圍內，表明本設計方案合理可行，同時為本設計的優化提供了可預期的效果。
