一种基于几何多尺度模型的一个半心室模型中奇静脉分流比的计算方法与流程

本发明提供了一种基于几何多尺度模型的一个半心室模型中奇静脉分流比的计算方法，其定义为公式(1)，

属于血流动力学计算领域。

背景技术：

研究一个半心室模型中奇静脉的分流比对于认识右心无效循环十分重要。而纵观大多数研究均是基于超声、磁共振成像方式得到奇静脉分流比。这种方法通过成像设备虽然能快速得到各血管的流速，但仪器本身有一定误差；其次对患者来说医疗花费较大，有些病人也不适于进行磁共振操作。现采用几何多尺度模型，利用有限元分析，既可以改变3d模型中奇静脉的内径，又能方便地通过调节0d集中参数模型的各部分参数(如血流阻力r、血管顺应性c、血液惯性l等)，来模拟不同肺阻力的实际生理，从而直观地观察3d模型内部的血液流动状况。在ansys-cfx后处理中可以方便地提取奇静脉流量、上腔静脉流量和肺主动脉流量，依据公式(1)计算出奇静脉分流比；此外还可以提取奇静脉的壁面切应力，从生物力学角度评价计算所得的奇静脉分流比。

几何多尺度建模是一种特殊的仿真血液循环系统的策略，本发明采用的是0d/3d耦合模型，其中3d模型是由ct图像重建的一个半心室模型，用来模拟感兴趣的局部流场细节，而0d集中参数模型是由r、c、l三元件组成的电路结构，用来仿真外周的血液循环系统。这种模型可以避免固定边界所带来的不良影响。

0d/3d耦合算法，它的重点在于0d集中参数模型与3d模型之间的数据交换。具体来说，0d集中参数模型为3d模型入口处提供流量波形边界，出口处提供压力波形边界，此时3d模型便可以在计算流体力学(cfd)数值模拟软件中进行流场计算，3d模型计算得到的入口的压力值与出口的流量值返回到0d集中参数模型。耦合算法中每个3d模型计算时间步进行一次数据交换，同时进行残差检测。通常定义3d模型出口压力和入口流量在不同时间步之间的误差作为残差检测项。当残差小于预先设定值时认为计算结果收敛，仿真结束，具体流程如图1所示。

ansys公司的ansys-cfx是cfd的仿真软件，它的前、后处理和流体计算较为完善。但想要实现0d/3d耦合模型流体计算，不仅需要掌握0d/3d耦合算法，还需了解ansys-cfx的二次开发、内存管理及多进程计算。以usercelfunction函数和userjunctionboxroutine程序块两种形式应用于cfx的3d模型计算之中。

技术实现要素：

本发明提出的一种基于几何多尺度模型的一个半心室模型中奇静脉分流比的计算方法，使用计算机仿真一个半心室模型的血流动力学状况。相较于常用的超声、磁共振成像方法，本发明的计算方法可以方便地改变奇静脉的内径，使用0d/3d耦合的算法，实现不同肺阻力时循环系统的仿真；还可以通过奇静脉的壁面切应力来评价奇静脉分流比的合理性。

一种基于几何多尺度模型的一个半心室模型中奇静脉分流比的计算方法，包括如下步骤：

1.构建一个半心室的3d模型；

2.构建血液循环系统0d集中参数模型并建立0d/3d耦合模型的几何多尺度模型；

3.对几何多尺度模型进行有限元分析，仿真不同肺阻力时一个半心室模型的血流动力学状况。提取奇静脉流量、上腔静脉流量和肺主动脉流量，由公式(1)计算得到奇静脉的分流比；提取不同肺阻力时奇静脉的壁面切应力，并与静脉系统承受的壁面切应力范围0.1-0.6pa进行比较，从生物力学角度评价奇静脉分流比的合理性，有关合理性的描述详情见下文步骤3中的内容。

步骤1具体为：

将胸部ct图像导入mimics中，经过灰度范围选择、阈值分割以及区域生长的图像处理后，初步得到格式为“.stl”的3d模型(包括奇静脉、上下腔静脉、肺主动脉、左右肺动脉)；再在freeform中对模型不连续和明显存在突变的部分进行修饰，得到顺滑的3d模型；之后在geomagic中对模型进行曲面光滑操作，得到格式为“.igs”的几何模型；最后用solidworks裁剪3d模型边界出入口，使其平整，并将模型导出格式为“.x-t”的数据，用于有限元计算。

步骤2具体为：

本发明构建的0d集中参数模型具体包括心脏模块、肺动脉模块、主动脉模块、上半身阻力与上腔静脉构成上半身循环模块、下半身阻力与下腔静脉构成下半身循环模块以及奇静脉模块。本发明中任意一段不含分叉且结构简单的血管用“倒l型”电路结构代替，这种结构如图2所示。心脏部分结构如图3所示，左心与右心的结构相同，使用前人研究中表达式e(t)来表示可变电容c(t)上的电压，以模拟心脏的收缩和舒张，e(t)是电容的倒数，单位为mmhg/ml；用二极管代替心脏瓣膜，以模拟生理上该位置的血液单方向流动的特点；电阻rla、rlv分别表示左心房、左心室的血流阻力。emax和emin分别是心脏收缩舒张功能的最大值与最小值。en(tn)是一个随时间变化的弹性模量，tc为心动周期。将0d集中参数模型中上半身阻力连接至入口上腔静脉，下半身阻力连接至入口下腔静脉，入口肺主动脉连接至心脏模块，出口左右肺动脉连接至心脏模块，与3d模型相连的外周循环部分仍以0d集中参数模型表示，从而组成0d/3d耦合的几何多尺度模型。

e(t)＝(emax-emin)·en(tn)+emin

en(tn)＝1.55[(tn/0.7)^1.9/(1+(tn/0.7)^1.9)][1/(1+(tn/1.17)^21.9)]

tn＝t/tmax，tmax＝0.15tc+0.2

步骤3具体为：

参考文献设置0d集中参数模型中各模块的参数：血流阻力r、血管顺应性c、血液惯性l及心脏收缩舒张功能最大值emax和最小值emin；再调试各模块的参数使左右肺动脉压力均值分别为10mmhg、20mmhg、30mmhg，波峰处出现重搏波，而在左右肺动脉压力为不同值时，上腔静脉流量波形、下腔静脉流量波形、肺主动脉流量波形均稳定收敛，分别对应肺阻力为2wood、3wood、4wood时的生理实际。

调试参数的具体过程：首先调节emax、emin使得左、右心室压力波形的波峰达到正常值120mmhg和25mmhg，若左心室压力波形呈下降趋势，调大上半身与下半身循环的阻力，反之，调小各阻力值；若右心室压力波形呈下降趋势，则调大肺动脉阻力，反之，调小其值，使得心室压力波形稳定收敛；其次增大主动脉电容，使主动脉压力波形中波峰与波谷之间的压差增大，直至压差为40mmhg，增大肺动脉电容，使得肺动脉压力波形中波峰与波谷之间的压差增大，直至压差为20mmhg；然后调节各模块的电阻使各模块流量波形稳定收敛；最后增大肺动脉电感，使得肺动脉波形出现重搏波；在每一次调试中以增大或减少20％的比例改变上述参数。

依据此时的参数将格式为“.x-t”的3d模型与0d集中参数模型出入口边界条件数据进行0d/3d耦合：首先在ansys-cfx中使用userjunctionboxroutine程序块完成0d/3d模型初始化；其次在ansys-cfx前处理中对3d模型进行网格划分，并设置血液密度、血液粘度、流体类型、时间步长、计算时间，将从0d集中参数模型得到不同肺阻力时的上腔静脉流量波形、下腔静脉流量波形、肺主动脉流量波形以及左右肺动脉压力波形作为3d模型的边界条件；然后在ansys-cfx中对3d模型进行流体计算；与此同时进行流体计算的收敛性判定，即定义流体计算相邻时间步之间3d模型中上腔静脉流量、下腔静脉流量、肺主动脉流量以及左右肺动脉压力的误差为残差检测项，当残差小于0.0001时，计算结果收敛，继续进行下一时间步的迭代计算直至指定时间，即三个心动周期2.4s，则仿真结束；如果残差不小于预先设定值0.0001，则使用ansys-cfx中的usercelfunction函数将从3d模型计算中得到的上腔静脉压力、下腔静脉压力、肺主动脉压力以及左右肺动脉流量作为0d集中参数模型计算的所缺项，使用userjunctionboxroutine程序块完成0d集中参数模型计算，继续进行此时间步的运算，直至收敛。

待计算完成后，在ansys-cfx后处理中，在奇静脉中段垂直于血管纵轴方向上建立横截面，分别提取奇静脉流量、上腔静脉流量和肺主动脉流量，按照公式(1)计算可得到一个半心室模型中的奇静脉分流比；提取奇静脉的壁面切应力，从生物力学角度评价奇静脉分流比的合理性，其合理性指标是奇静脉的壁面切应力大于静脉系统承受的壁面切应力范围0.1-0.6pa，且随肺阻力增大，奇静脉的壁面切应力与奇静脉分流比正相关，二者均逐渐增大；对比奇静脉的壁面切应力与计算所得的奇静脉分流比，若符合前述指标，则合理；若肺阻力为4wood时提取的奇静脉壁面切应力小于肺阻力为2wood时的壁面切应力，则肺阻力为4wood时计算出的奇静脉分流比是不合理的，这时将一个半心室模型中奇静脉内径增大5％后重复计算步骤1.3所述的内容；若肺阻力为2wood时提取的奇静脉壁面切应力大于肺阻力为4wood时的壁面切应力，则肺阻力为2wood时计算出的奇静脉壁面切应力是不合理的，这时将一个半心室模型中奇静脉内径减小5％后重复计算步骤1.3所述的内容。

附图说明

图10d/3d耦合算法流程图；

图2“倒l型”电路结构示意图；

图3左心电路结构示意图；

图4真实3d模型示意图；

图5构建的集中参数模型示意图；

图6构建的几何多尺度模型示意图；

图7集中参数模型调试参数后边界条件示意图；

inlet----入口

outlet----出口

l/ra(left/rightatrium)----左/右心房

l/rv(left/rightventricle)----左/右心室

r1----上半身阻力

r2----下半身阻力

o(aorta)----主动脉

svc(superiorvenacaval)----上腔静脉

ivc(inferiorvenacaval)----下腔静脉

azg(azygosvein)----奇静脉

l/rp(leftpulmonaryartery)----左/右肺动脉

mpa(themainpulmonaryartery)----肺主动脉

具体实施方式

下面结合具体实施方式解释本发明，实施例1

1.构建一个半心室模型

1.1在mimics中利用灰度阈值分割、区域生长的方法对ct图像数据进行处理，初步得到奇静脉内径为5mm的一个半心室模型，保存格式为“.stl”。

1.2将1.1步得到的模型导入freeform中，利用力反馈器，对模型不连续和明显存在突变的部分进行修饰，导出顺滑的3d模型。

1.3将1.2步得到的模型导入geomagic中，进行曲面光滑操作，以达到简化模型的目的，导出“.igs”格式的3d模型。

1.4将1.3步得到的模型导入solidworks中，在出/入口处建立基准面，裁剪边界口，使3d模型边界口平整；将“.igs”格式的模型转为可用于有限元计算的“.x-t”数据。如图4。

2.建立循环系统0d集中参数模型与0d/3d耦合模型

实例1建立的0d集中参数模型如图5所示。将一个半心室模型以3d模型表示，与它相连的外周循环部分以0d集中参数模型表示，从而组成0d/3d耦合模型，实例1建立的0d/3d耦合模型如图6所示。

3.几何多尺度有限元分析

调试0d集中参数模型的参数emax、emin、r、c和l，使得出/入边界口的波形符合肺阻力为2wood时的生理实际，如图7所示。

依据此时的参数将格式为“.x-t”的3d模型与0d集中参数模型出入口边界条件数据进行0d/3d耦合：首先在ansys-cfx中使用userjunctionboxroutine程序块完成0d/3d模型初始化；在ansys-cfx前处理中对3d模型划分网格，设置血液密度为1050kg/m³，血液动力学粘度为0.0035pa·s，流体类型设置为不可压缩的牛顿流体，设置时间步长为0.0025s，计算时长为2.4s，将得到的入口流量和出口压力波形作为3d模型的边界条件；然后在ansys中对3d模型进行流体计算；与此同时进行收敛性判定；定义相邻时间步之间3d模型中上腔静脉流量、下腔静脉流量、肺主动脉流量以及左右肺动脉压力的误差为残差检测项，当残差小于0.0001时，计算结果收敛，继续下一时间步的迭代计算直至2.4s，则仿真结束；若残差不小于预先设定值0.0001，则使用ansys-cfx中的usercelfunction函数将从3d模型计算得到的上腔静脉压力、下腔静脉压力、肺主动脉压力以及左右肺动脉流量作为0d集中参数模型计算的所缺项，使用userjunctionboxroutine程序块完成0d集中参数模型计算，继续进行此步运算，直至收敛。

计算结束后在ansys-cfx后处理中在奇静脉中段垂直于血管纵轴方向上建立横截面plan1，之后新建四个表达式expression，分别表示提取的奇静脉流量、上腔静脉流量、肺主动脉流量以及奇静脉分流比，分别定义为q3＝massflow()@plane1，qq1＝massflow()@inletsvc，qq2＝massflow()@inletmpa，additionpercent＝q3/(qq1+qq2)。计算得到奇静脉分流比为0.237。其中inletsvc和inletmpa为上腔静脉和肺主动脉的入口。定义表达式wss_azy＝ave(wallshear)@azy，提取奇静脉的壁面切应力为3.5pa，超过静脉系统承载的壁面切应力范围0.1-0.6pa，表明奇静脉受到应力后适应性生长，使得上腔静脉有部分血液经奇静脉弓流向右心室，当肺阻力超过2wood时，提取的奇静脉壁面切应力大于3.5pa，计算的分流比才是合理的，否则不合理，此时将一个半心室模型中奇静脉内径增大5％后重复计算步骤3中的内容。