文摘
降低复杂性的计算方法结合血管血流动力学的模拟为研究对象,建立了一维(1 d)波传播模型与零维(0 d)集总模型血管微循环。尽管降低维度,当前算法用于解决模型方程,模拟压力和流量是相当复杂的,从而限制接受医学领域。这种复杂性主要来自方法结合1 d和0 d模型方程。本文提出了一种数值方法,不再需要额外的耦合方法,使随机组合1 d和0 d模型只使用压力作为状态变量。方法应用于60的主要动脉血管树组成的身体和头部。仿真结果是现实的。数值方法稳定,显示了良好的收敛性。
1。介绍
血液流动包括压力和流波,通过血管系统传播。作为妥协之间的计算需求和物理细节,一维(1 d)网络模型被用来研究压力和流量波形在正常和病理条件下(1]。这种建模方法已应用于系统性动脉系统(2- - - - - -6),冠状动脉树(7,8),和大脑血管树(9,10]。
这些一维网络模型由本地的元素描述压力和流量之间的关系。压力之间的关系、区域和血管中流动的一维波传播方程,即一维偏微分方程派生通过集成的质量和动量navier - stokes方程在血管的横截面积11]。血管口径减少和对周边血管数量的增加,达到一个点,它不再是实现单独模型血管。此时的血管是截断和贡献远端血管压力和流量是由0维集总模型,描述如windkessel模型(5,6,12)或结构树模型(3]。
解决方程组由0和1 d模型和模拟的压力和流传播波通过血管系统,各种各样的,相当复杂,算法存在。对于一维波传播方程,数值方法从同一关系的压力,区域和流或横截面平均速度。第一个状态变量的数值方法产生差异选择依然存在。面积和压力相关通过血管壁的本构关系,结果是一个压力(2,8],area-velocity [13),面积量(3),或压力速度公式9,11,14,15]。第二个来源的差异是选择的空间离散化方程。方法包括有限差分(3,8)和光谱/有限元方案(2,13]。结果是一组常微分方程的状态变量必须解决。第三,不同的方法被用来执行压力和流量的连续性在船分岔和血管之间的接口和外围。方法包括0和1 d的弱耦合方程(15,16),计算黎曼不变量(9,14),或添加惩罚方程(2]。
本研究的贡献是开发一个简化的数值方法中,压力是唯一的状态变量。在这种方法中,一维波传播和0维集总模型方程转换为相同的形式。因此,1 d和0 d模型结合而不需要指定附加耦合方程。这允许灵活的建模从0和1 d元素模拟血管网络的压力和流量。说明,应用数值方法模拟血管树中压力和流量波形由60的主要动脉的身体和头部。
2。方法
2.1。压力在大血管(1 d)的关系
在大型动脉血压(Pa),血液流动(m3·年代−1)、壁面切应力(Pa)和血管横截面积(m2)是由1 d质量方程和动量相关。当忽视通过血管壁渗漏以及引力,质量和动量的平衡是由(推导过程可以发现,例如,休斯和卢布林(11]和Van de Vosse Stergiopulos [1]) 在这个()制定,沿着船轴(m)表示的坐标,(m)表示容器半径,(m2·巴−1)表示容器area-compliance,(公斤·米−3)表示血液密度。
壁剪切应力和持续在(2)是由假设估计速度剖面。选择的速度剖面,几个选项是可能的(1]。在这项研究中近似假定速度概要文件(2]。使用近似速度资料,给出了壁面切应力 与(Pa·s)血液粘度,横截面积的比例与惯性控制流,和Womersley数量包括心动周期的持续时间[s],血管横截面积在参考压强。的近似速度资料,常数是由
完成质量和动量方程表达式区()和区域合规()作为压力的函数。在这项研究中,假设一个非线性弹性血管壁的压力依赖性的遵从性给出的 在哪里,,,指定的压力依赖性area-compliance(函数从Langewouters et al。17])。泊松比,扬斯模量和壁厚指定区域的半径依赖合规(函数从Bessems et al。2])。横截面积的压力依赖性的表达式是通过整合区域合规压力。
2.2。外围国家的压力关系(0 d)
的贡献各动脉末梢血管终点站是集中在三元素windkessel模型(6,18)(图1(c))。通常一个微分方程推导出相关的压力和流在入口处的Windkessel2,8,10]: 与特性阻抗,外周阻力,外围遵从性。然而,在采用这个方程,它是隐含假定静脉出口压力为零。因此,模型的应用范围是有限的具体情况。更普遍的方法是限制Windkessel方程那些与压力降穿过不同元素构成Windkessel模型。使用该离散化,如图2,
2.3。数值方法
确定压力和流量的血管网络,血管1 d压力关系(1)和(2)和0 d外围的压力关系(7)是在一个完全耦合的方式解决。完整的耦合是通过铸造方程到相同的离散形式。
windkessels的子元素,两个节点压力和两个节点流(图中定义2)。一个关键的选择是节点流都是导演向内。的组合(7)与 收益率(附录一个) 在列和分别包含节点压力和流动。矩阵包含了外围合规和矩阵包含特征阻抗和外围电阻的倒数。二阶向后差分格式用于一步在时间与时间步。结果是windkessel方程写成 与
以前,休伯特等。19)提出了一个方法,每个血管段的波传播方程也被扔在一个集总模型包括电阻、遵从性,和电感,即构成windkessel相同的块。这种方法与谱元法通过基准测试Bessems et al。2]。在一个类似的,但更直接的方法是遵循不需要推导一个集总参数模型,导致一个简化的实现。
首先,波传播方程线性化使用的估计区域合规、横截面积、壁面切应力和对流加速获得前一个时间步(象征): 随后,该船段分为更小的近似大小的两节点元素。实际的大小元素用于离散化可以不同和船段之间 与船的长度。随后,一个梯形规则用于空间积分方程沿船轴和表达方程的节点压力和流动。在这种离散化步骤中,节点为每个元素流动也选择直接向内(图2)。一次,二阶向后差分格式用于一步。因此,质量和动量方程写成(附录B)
通过定义都是直接向内流动,压力和流量的连续性0 d-1d接口和1 d-1d接口(例如,分岔)是自动满足方程的元素在组装的过程中到大型方程组(附录C)。结果,组装大型方程组是由: 在哪里为每个节点包含一个零值,除了节点外部流程的规定。这个方程组解决一次压力边界条件和外部流。
请注意任何非线性流使刚度矩阵流的依赖,因此要求流后重新计算每个压力计算。这可以很容易地通过使用(10)或(14在元素级别);也就是说, 因此,对于节点所有元素都流计算。
压力和流量计算后,模拟所得到下一个时间步。重复这个过程,直到心脏周期时间是达到了。在这一点上,它是检查是否模拟血流动力学稳定状态。血流动力学稳定状态被认为是如果节点最大相对均方根准则,表示压力和流量小于。为, 与心动周期和数量节点号。完整的概述图算法如图3。
2.4。模拟
2.4.1。模型参数的选择
应用数值方法模拟血流动力学在一个60的主要动脉动脉树组成的身体和头部(图1(a))。血液密度和血液粘度假定为1.05公斤·米吗−3和巴勒斯坦权力机构年代,分别。参考透壁的压力()被设置为13.3 kPa。压力依赖性的合规根据指定参数,kPa,kPa,,(8]。注意,这些参数意味着参考合规压力降低。假设不可压缩血管壁材料呈现。值船长度、壁厚、半径和杨氏模量以及windkessel参数取自穆德et al。10、表2和图3)。
主动脉流入被认为是唯一的外部流程和规定根据波形如图1(b)。将心动周期时间年代。在windkessels静脉和血管外的压力设置为0 kPa。
2.4.2。模拟执行
评估该数值方法的收敛性行为对时间和空间离散化,完成一系列的模拟组合的元素大小()约为2.5、10和40毫米和时间步骤()1、2、4、8女士。每个仿真开始零压力和流动。
为每个模拟,血流动力学收敛准则心动周期数的函数决定。在收敛,主动脉压力和流量波形可视化,左腿、左臂,在大脑中。空间和时间离散化的影响在模拟压力和流量波形是由相对均方根差量化,由 计算的波形,得到最密集的网格和最小时间步长(下标REF)作为参考。
3所示。结果
对于大多数模拟压力和流量都聚集在心脏周期(图124)。中间心脏周期,规范压力约一个数量级低于流动。元素的大小似乎没有效果的下降正常心动周期数的函数;也就是说,没有视觉元素大小之间的歧视是可能的。收敛略慢对于较大的时间步长,但只有在情况下的最大时间步(ms)额外的心动周期是必需的。
图5显示了模拟压力和流量波形。结果表明,压力波的振幅增加向外围。此外,相对高度的重搏切迹对压力脉冲减少向外围。动脉附近等外围的手臂给逆转流在一个心动周期的一部分。
如图6的影响,时间步长计算波形(1)通常是一个数量级的压力低于流动和(2)增加对外围;即最大流量差异发生在前沟通大脑中动脉、尺动脉的最大压差(图6 (b))。元素的大小只有轻微影响计算压力和流量波形相比时间步的效果(图6(一))。使用一个元素的大小40毫米,而不是在最小时间步长增加2.5毫米规范不到的流。在时间步2而不是1女士女士在最小的元素大小导致大约增长了。增加时间步长到8女士导致阻尼的压力和流量波形。
(一)
(b)
4所示。讨论
在这项研究中,开发了一个简化的数值方法对血压和流波形的时域模拟夫妇的血管系统非线性一维(1 d)波传播模型血管零维(0 d)集中(windkessel)模型的外围用压力自由度。
展示在生理环境中性能的方法,该方法应用于模拟血流动力学在包含60的主要动脉的血管网络的身体和大脑。船的具体选择行为,速度剖面,windkessel参数,和必要的边界条件是超出了本研究的范围。血管的压力区关系被认为非线性和对流加速包括评估行为的方法求解模型方程最非线性形式。
得到的压力和流量波形(图的方法5)是类似模拟(10)以及实验测量了被别人(20.];即计算波形演示的生理特征(1)压力波的振幅的增加和减少的相对高度与距离增加主动脉根重搏切迹,(2)流在一个心动周期的一部分的逆转动脉的手臂。
4.1。该方法收敛性的行为
通常需要12心脏周期达到收敛当从零压力和流量条件。这个融合相当独立使用的元素和时间步长离散化模型方程(图4)。当比较聚集的情况下,时间步长有显著影响的压力和流量波形比元素大小(图6(一)),时间步的影响通常是高一个数量级比压力流。的时间步8女士介绍了重要的阻尼波压力和流量相比获得的结果使用1毫秒的时间步(图6 (b))。增加对周边的影响;即最大均方根差异被发现在手臂和尺动脉的前沟通大脑中动脉。时间的增加效果一步外围最有可能造成的生理趋陡的压力和流量对边缘(图波形5)。2的时间步女士收益率的增加相对均方根差流的,这表明1女士是足够小的时间步(图6(一))。
4.2。该方法的好处
在引言部分上市,许多不同的数值方法存在两个一维波传播方程的大动脉血管的0 d集中windkessel方程外围血管树的一部分。通常,血管段的波传播方程都写在使用有限的离散形式/ spectral-element或有限差分方案。这些方法需要额外的缺点,轨耦合方程定义的黎曼不变量或惩罚函数。此外,外围模型的方程通常包含通过求解特征方程(如三元素的windkessel (6))的波传播方程。这种方法的缺点是,这种特征方程必须可用。例如,这是并非如此,当windkessels的终点站是连接到一个静脉循环。
在数值方法,windkessel和波传播方程转换为相同的形式强烈一些他们不需要额外的耦合方程或可用性的特征方程。事实上,windkessel的任意组合(或集中)和波传播元素是可能的,允许更大范围的应用程序结合动脉血管网络,微循环(外围),和静脉。
容易扩展的数值方法允许集总模型的心脏研究arterioventricular如由他人交往(21- - - - - -23]。心脏收缩可以考虑通过指定心室压力(时变)本质边界条件,而非处方主动脉流入。心室体积可以被更新在弹射使用主动脉流入的时间步。尽管这样,心室和主动脉流只是弱耦合,计算复杂性是有限的。
4.3。该方法的局限性
将方程的波传播模型转换为相同的形式作为windkessel这些模型,要求每个离散元素只包含两个节点的直接向内流动。因此,如用于高阶元素,例如,光谱元素离散化不再是可能的。然而,这种近似限制秩序被发现有小影响元素大小似乎是次要的收敛性以及压力和流量曲线。
如切图所示6(一)关于时间步,收敛阶小于二阶,虽然二阶向后差分格式用于集成的时间。这融合订单预计会减少由于非线性一维波传播方程,但有可能被放大使用估计前一个时间步的线性化的过程。线性化的手段,例如,牛顿迭代计划可能已经完成但不包括进一步简化数值方法(实现)。
血管网络的提出在这项研究中,我们把windkessel外围模型。然而,算法并不局限于这个特定的模式选择。等集总元件模型结构树模型(3)可以很容易地注册只要压力关系可以转换为相同的形式,对波传播模型。近似的方法也不是限制速度概要文件在本研究假设。使用的,例如,Womersley速度概要文件也是可能的,作为唯一要求的方法进行时间是这个区域合规,壁面切应力和对流加速在以前的时间步的。
5。结论
总之,开发的一种新型数值方法计算的压力和流量波形的血管系统。使用压力作为唯一的自由度,0 d集中(windkessel)元素和1 d可以随机波传播元素结合而不需要额外的耦合方程。这个属性促进灵活建模从0和1 d元素广泛应用在研究血管血流动力学。
附录
答:(0 d) Windkessel元素方程的推导
应用程序定义的节点压力和流如图2了(7)和(8)。矩阵形式,这些方程 与矩阵定义为 组装方程组的收益率 与 和,。因为应该没有净流或从节点3,节点的选择定向流动向内收益率。注意,剩下的流,,表示外部流。应用二阶向后差分格式的收益率 在这,,。
b .推导方程(1 d)血管的元素
描述的线性化版本1 d质量平衡和动量是由(12)的一个离散的点,血管分为许多重叠两节点元素。作为一个后果 与元素的数量。梯形规则用于空间整合元素方程从节点1到节点2。从传统的离散化与流直接从节点1到节点2如图2域的收益率、集成的元素 与元素的长度和。上标c是用来表明流采用的传统定义。以矩阵形式,(B.2)成为 在哪里和包含元素的节点压力和流量,分别。使用二阶向后差分格式和时间步长,(B.3)被编写为 与
接下来,提出了离散化的开关是由如图2,在第二节点直接向内流动,也就是说,和。注意,这个开关意味着改变的第二列的迹象。因此,(B.4)成为 与 分离后的矩阵从流列,(1 d)元素方程最后阅读
c .组装的大型方程组:耦合(0 d)和方程(1 d)元素
请注意,(本)和(B.8)都是形式 组装的大型方程组包括总和的1 d元素根据方程(责任),以及所有0 d元素方程。在这一过程中,元素流动将被添加。由于所有节点直接向内流动,内部节点,收益质量的平衡零进入组装流程列(表示)。事实上,只有非零流在外部节点在组装流将产生一个非零元素列。因此,方程组是由组装 因为外部流动,在主动脉根等规定,节点压力只剩下自由度,和系统可以解决一次剩余的压力容器或windkessel末端是规定。
承认
这项研究是由欧盟第七框架计划(ARCH ict - 224390)。