血流数值模拟方法及装置

1.本技术属于计算流体力学技术领域，尤其涉及一种血流数值模拟方法及装置。


背景技术：

2.由于血管的几何形态、生理条件和生物力学等因素会对血管内部的血液流动产生影响，而复杂的血液流动会导致一些血管的生理变化，这样与病变的发生和发展有一定关联。因此，血流动力学模拟可以给脑血管疾病的预防和治疗提供重要参考依据。
3.目前，计算流体力学(computational fluid dynamics,cfd)方法作为一种常用的血流动力学定量分析工具，已广泛应用于血流动力学模拟，来评估血管内部的血流状态。在使用cfd方法对血管的流体力学性质研究时，由于远端血管的血流量无法一一测量，而现有的根据血管截面积确定血管中血流量方法准确性不足，从而在以血流量为依据进行血流动力学模拟时也会导致血流模拟结果不够准确。


技术实现要素：

4.有鉴于此，本技术提供了一种血流数值模拟方法及装置，可以提高血管的血流数值模拟结果的准确性。
5.第一方面，本技术实施例提供了一种血流数值模拟方法，包括：
6.获取待评估血管的三维模型和血液流速测量数据；所述血液流速测量数据包括所述三维模型中多个目标点的血液流速；所述三维模型包括多个血管区域；
7.根据所述多个目标点的血液流速和所述三维模型的入口处的压力差，确定每个所述血管区域中出口的阻力；
8.根据所述血液流速测量数据和多个所述血管区域中出口的阻力，确定所述待评估血管的流体控制方程的边界条件；
9.根据所述边界条件对所述流体控制方程求解，获取所述待评估血管的血流数值模拟结果。
10.本技术实施例根据目标点的血液流速和三维模型的入口处的压力差，确定血管区域中出口的阻力，提高了确定出口的阻力时的准确性，并根据出口的阻力确定流体控制方程的边界条件，提高了求解流体控制方程得出的血流数值模拟结果的准确性。
11.在一种可能的实现方式中，所述根据所述多个目标点的血液流速和所述三维模型的入口处的压力差，确定每个所述血管区域中出口的阻力，包括：
12.根据各个所述目标点的血液流速分别获取各个所述目标点的血流量；
13.根据多个所述目标点的血流量，确定每个所述血管区域相对于所述三维模型的血流量占比；
14.根据每个所述血管区域相对于所述三维模型的血流量占比和所述三维模型的入口处的压力差，确定每个所述血管区域中出口的阻力。
15.本技术实施例通过目标点的血液流速确定出血管区域的血流量占比，并根据血流
量占比和三维模型的入口处的压力差，确定出口的阻力，提高了确定血管的出口阻力的准确性。
16.在一种可能的实现方式中，所述根据每个所述血管区域相对于所述三维模型的血流量占比和所述三维模型的入口处的压力差，确定每个所述血管区域中出口的阻力，包括：
17.根据所述三维模型的总血流量和所述三维模型的入口处的压力差，获取所述三维模型的总阻力；
18.针对每个所述血管区域，根据所述血管区域相对于所述三维模型的血流量占比和所述三维模型的总阻力，获取所述血管区域的总阻力；并根据所述血管区域的总阻力和所述血管区域中出口的横截面尺寸，获取所述血管区域中出口的阻力。
19.本技术实施例先确定出三维模型的总阻力，再根据血管区域的血流量占比分配血管区域的总阻力，然后根据血管区域中出口的截面尺寸获取血管区域中的出口阻力，提高了得出的血管区域中的出口阻力的准确性。
20.在一种可能的实现方式中，所述流体控制方程为纳维-斯托克斯方程，所述边界条件包括出口边界条件，所述根据所述血液流速测量数据和多个所述血管区域中出口的阻力，确定所述待评估血管的流体控制方程的边界条件，包括：
21.根据所述血液流速测量数据和每个所述血管区域中出口的阻力，确定每个所述血管区域中每个出口对应的弹性腔模型的参数；
22.根据多个所述血管区域中每个出口对应的弹性腔模型的参数，确定所述流体控制方程的出口边界条件。
23.本技术实施例通过弹性腔模型来对出口进行模拟，确定出弹性腔模型的参数后得出流体控制方程的出口边界条件，提高了出口边界条件的准确性。
24.在一种可能的实现方式中，所述根据所述边界条件对所述流体控制方程求解，包括：
25.将所述三维模型划分为多个网格；
26.根据所述多个网格，在空间域上和时间域上分别离散所述流体控制方程，得到稀疏非线性系统；
27.根据所述边界条件对所述稀疏非线性系统求解。
28.本技术实施例通过将三维模型划分为多个网格，根据网格离散流体控制方程，提高了流体控制方程的可解性。
29.在一种可能的实现方式中，所述根据所述边界条件对所述流体控制方程求解，包括：
30.采用牛顿-克雷洛夫-施瓦兹算法对所述流体控制方程求解。
31.本技术实施例通过牛顿-克雷洛夫-施瓦兹算法求解流体控制方程，提高了用计算机求解流体控制方程时的计算效率。
32.在一种可能的实现方式中，所述获取待评估血管的三维模型，包括：
33.获取所述待评估血管的断层扫描影像数据；
34.根据所述断层扫描影像数据，获取所述待评估血管的三维模型。
35.本技术实施例通过断层扫描影像数据获取待评估血管的三维模型，提高了待评估血管的三维模型的准确性。
36.第二方面，本技术实施例提供了一种血流数值模拟装置，包括：
37.第一获取模块，用于获取待评估血管的三维模型和血液流速测量数据；所述血液流速测量数据包括所述三维模型中多个目标点的血液流速；所述三维模型包括多个血管区域；
38.第一确定模块，用于根据所述多个目标点的血液流速和所述三维模型的入口处的压力差，确定每个所述血管区域中出口的阻力；
39.第二确定模块，根据所述血液流速测量数据和多个所述血管区域中出口的阻力，确定所述待评估血管的流体控制方程的边界条件；
40.第二获取模块，用于根据所述边界条件对所述流体控制方程求解，获取所述待评估血管的血流数值模拟结果。
41.第三方面，本技术实施例提供了一种血流数值模拟装置，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，用于执行上述第一方面的任意可能的实现方式中的方法。
42.第四方面，本技术实施例提供了一种计算机可读存储介质，所述可读存储介质用于保存计算机程序，所述计算机程序被处理器执行时，能够实现上述第一方面的任意可能的实现方式中的方法。
43.第五方面，本技术实施例提供了一种计算机程序产品，当计算机程序在血流数值模拟装置上运行时，使得血流数值模拟装置执行上述第一方面的任意可能的实现方式中的方法。
44.可以理解的是，上述第二方面至第五方面的有益效果可以参见上述第一方面中的相关描述，在此不再赘述。
附图说明
45.图1是本技术实施例提供的血流数值模拟方法的流程示意图；
46.图2为本技术实施例提供的三维模型区域划分及目标点所在血管截面的示意图；
47.图3是本技术实施例提供的经颅多普勒超声图；
48.图4是本技术实施例提供的确定血管区域中出口的阻力的方法的流程示意图；
49.图5为本技术实施例提供的根据每个血管区域相对于三维模型的血流量占比确定出口阻力的流程示意图；
50.图6是本技术实施例提供的三元弹性腔模型的示意图；
51.图7是本技术实施例提供的网格划分的示意图；
52.图8是本技术实施例提供的血流数值模拟装置的示意性框图；
53.图9是本技术实施例提供的血流数值模拟装置的结构示意图。
具体实施方式
54.以下描述中，为了说明而不是为了限定，提出了诸如特定系统结构、技术之类的具体细节，以便透彻理解本技术实施例。然而，本领域的技术人员应当清楚，在没有这些具体细节的其它实施例中也可以实现本技术。在其它情况中，省略对众所周知的系统、装置、电路以及方法的详细说明，以免不必要的细节妨碍本技术的描述。
55.应当理解，当在本技术说明书和所附权利要求书中使用时，术语“包括”指示所描述特征、整体、步骤、操作、元素和/或组件的存在，但并不排除一个或多个其它特征、整体、步骤、操作、元素、组件和/或其集合的存在或添加。
56.还应当理解，在本技术说明书和所附权利要求书中使用的术语“和/或”是指相关联列出的项中的一个或多个的任何组合以及所有可能组合，并且包括这些组合。
57.如在本技术说明书和所附权利要求书中所使用的那样，术语“如果”可以依据上下文被解释为“当...时”或“一旦”或“响应于确定”或“响应于检测到”。类似地，短语“如果确定”或“如果检测到[所描述条件或事件]”可以依据上下文被解释为意指“一旦确定”或“响应于确定”或“一旦检测到[所描述条件或事件]”或“响应于检测到[所描述条件或事件]”。
[0058]
另外，在本技术说明书和所附权利要求书的描述中，术语“第一”、“第二”、“第三”等仅用于区分描述，而不能理解为指示或暗示相对重要性。
[0059]
在本技术说明书中描述的参考“一个实施例”或“一些实施例”等意味着在本技术的一个或多个实施例中包括结合该实施例描述的特定特征、结构或特点。由此，在本说明书中的不同之处出现的语句“在一个实施例中”、“在一些实施例中”、“在其他一些实施例中”、“在另外一些实施例中”等不是必然都参考相同的实施例，而是意味着“一个或多个但不是所有的实施例”，除非是以其他方式另外特别强调。术语“包括”、“包含”、“具有”及它们的变形都意味着“包括但不限于”，除非是以其他方式另外特别强调。
[0060]
发明构思简介
[0061]
本技术实施例提供的血流数值模拟方法适用于对血管内的血流情况进行分析，本技术实施例不对具体的应用场景做出限定，能通过本技术实施例提供的方法确定出血管内血液的流体控制方程的边界条件，从而求解得出血流的数值模拟结果的应用场景都属于本技术的保护范围。在一种可能的应用场景中，本技术实施例提供的血流数值模拟方法可以应用于人体器官中的血管，示例性的，人体器官中的血管可以为脑血管。本技术实施例将以脑血管为例对血流数值模拟方法进行具体说明。
[0062]
图1为本技术实施例提供的血流数值模拟方法的流程示意图，下面结合图1对本技术实施例提供的方法进行说明，这种血流数值模拟方法包括如下步骤：
[0063]
步骤101：获取待评估血管的三维模型和血液流速测量数据，血液流速测量数据包括三维模型中多个目标点的血液流速，三维模型包括多个血管区域。
[0064]
本技术实施例对获取三维模型的方法不做限定，在一种可能的实现方式中，获取待评估血管的断层扫描影像数据，然后根据断层扫描影像数据，获取待评估血管的三维模型。示例性的，断层扫描影像数据为电子计算机断层扫描(computed tomography,ct)影像数据。
[0065]
获取到断层扫描影像数据后，为了通过断层扫描影像数据构建出三维模型，还需要对影像数据进行进一步的处理。以脑血管为例，由于脑部血管分支繁多且医学断层扫描影像数据的显影质量的精度有限，无法实现对于直径过小的远端血管的重建，因此重建出的三维模型主要包括直径较大的血管。
[0066]
值得注意的是，重建出的血管的三维模型应该包括获取血液流速测量数据的目标点。其中，本技术实施例所述的血液流速测量数据，指的是通过仪器、测量方法等，对待评估血管进行真实测量后得到的有效数据。目标点则是待评估血管中，通过血液流速测量数据
得到有效数据的点。示例性的，图2为本技术实施例提供的三维模型区域划分及目标点所在血管截面的示意图，下面结合图2对三维模型进行示例性的描述。图2中所示的三维模型为脑动脉血管的三维模型，其中包括了左右侧大脑中动脉和前动脉，s1、s2、s3和s4为四个目标点。
[0067]
在血液流动的过程中，血液由三维模型底部较粗的血管流向三维模型顶部较细的血管，因此，图2中底部为三维模型的入口，顶部为出口。示例性的，图2中包括入口a和出口b，血液由入口a流入三维模型，通过出口b流出三维模型。类似的，按照血液的流动方向，可以看出图2的三维模型中包括2个入口和多个出口。
[0068]
本技术也不对根据断层扫描影像数据构建三维模型的技术手段做出限定，在一种可能的实现方式中，通过mimics或geomagic软件进行三维模型的构建。
[0069]
值得注意的是，通过三维模型的构建，可以计算得出三维模型中血管的半径、直径或横截面积中的至少一项。
[0070]
从图2中可以看出，本技术实施例在三维模型中划分出了多个血管区域，不同血管区域的确定将在后续步骤中用于根据流量分配出口阻力。因此，血管区域的确定以目标点的位置为基础。具体的，血管区域内的总血液流量可以通过目标点的血液流速测量数据进行确定。
[0071]
以图2为示例，其中包括1，2和3三个血管区域。
[0072]
目标点的血液流速和截面积已知，因此可以计算得出目标点的血流量。此外，血管内存在流量守恒定律，即流入血管分叉的血流量之和等于流出的血流量之和。如图2所示，血液由s1流入血管区域1，血管区域1内的总血流量为s1的血流量，血液由s2和s3流入血管区域2，血管区域2内的总血流量为s2和s3的血流量之和，血液由s4流入血管区域3，血管区域3内的总血流量为s4的血流量。图2只是进行示例性的描述，当三维模型中存在一定的目标点时，可以通过类似的方法确定三维模型内的血管区域。本技术实施例不对血液流速测量数据的获取方法做出限定，以脑血管为例，在一种可能的实现方式中，血液流速测量数据可以包括经颅多普勒超声(transcranial doppler,tcd)数据。tcd通过人类颅骨自然薄弱部位作为超声波入射点，利用多普勒效应获取目标点的血液流速，因此，tcd方法能够获取血液流速的目标点有限。
[0073]
值得注意的是，由于本技术实施例所使用的方法为cfd方法，因此，本技术实施例需要的血液流速为能够体现出血液流速随时间变化的具体数据，以及周期时长。周期时长可以理解为血液流速变化的时间周期，本技术实施例需要测量的血液流速为一个或多个周期时长内的血液流速。示例性的，图3为本技术实施例提供的经颅多普勒超声图，其中包括s1，s2，s3和s4四个目标点的血液流速，其中横轴为时间，纵轴为血液流速，根据图3的波形可以看出其中包括四个周期时长的血流速度。在本技术实施例中，也就是获取图2的三维模型中所示的四个目标点对应的血液流速。
[0074]
步骤102：根据多个目标点的血液流速和三维模型入口处的压力差，确定每个血管区域中出口的阻力。
[0075]
在一种可能的实现方式中，步骤102包括以下步骤，图4为本技术实施例提供的确定血管区域中出口的阻力的方法的流程示意图，下面结合图4对这些步骤进行说明：
[0076]
步骤201：根据各个目标点的血液流速分别获取各个目标点的血流量。
[0077]
如前文所述，在步骤101中能够得出目标点的血液流速，且血流量指的是单位时间流经血管某一截面的血量，因此可以计算出目标点的血流量。示例性的，在图2所示的三维模型中，可以根据目标点s1、s2、s3和s4的血液流速和截面积分别计算出对应的血流量。
[0078]
步骤202：根据多个目标点的血流量，确定每个血管区域相对于三维模型的血流量占比。
[0079]
在确定每个血管区域的流量占比时，需要首先确认待评估血管的总血流量。在一种可能的实现方式中，通过tcd方法确认三维模型入口处对应的血液流速，计算出流入待评估血管的总血流量。
[0080]
随后，通过目标点的血流量确定每个血管区域的血流量，示例性的，图2中血管区域1的血流量为s1的血流量。
[0081]
在得出总血流量和每个血管区域的血流量后，可以得出每个血管区域的血流量占比。
[0082]
步骤203：根据每个血管区域相对于三维模型的血流量占比和三维模型入口处的压力差，确定每个血管区域中出口的阻力。
[0083]
在一种可能的实现方式中，步骤203包括图5中的步骤，图5为本技术实施例提供的根据每个血管区域相对于三维模型的血流量占比确定出口阻力的流程示意图，下面结合图5对这些步骤进行说明：
[0084]
步骤301：根据三维模型的总血流量和三维模型的入口处的压力差，获取三维模型的总阻力。
[0085]
在前文中提到了总血流量的确认方式，因此可以获取入口处的压力差，从而通过欧姆定律确认三维模型的总阻力。其中，欧姆定律即血流量与血管两端的压力差成正比，与血流阻力成反比。
[0086]
以脑血管为例，由于出口处的血管的压力在一个周期时长内较为稳定，在使用欧姆定律时可以将出口处的血管压力近似为常数。因此，获取两个不同时刻之间的入口处压力差，以及三维模型中对应的血流量，就可以计算出三维模型的总阻力。具体的，三维模型的总阻力指的是三维模型中所有出口的总的阻力。
[0087]
在一种可能的实现方式中，入口处压力差使用收缩期和舒张期的左臂处或颈动脉处的测量压力差，这一数值较为容易测量且准确性较高。
[0088]
步骤302：针对每个血管区域，根据血管区域相对于三维模型的血流量占比和三维模型的总阻力，获取血管区域的总阻力。
[0089]
根据步骤202中得出的血流量占比，以及步骤301中得出的三维模型的总阻力，可以计算得出血管区域的总阻力。其中，血管区域的总阻力指的是一个血管区域内所有出口总的阻力。
[0090]
以图2为例，可以看出不同血管区域彼此之间为并联关系，因此可以根据阻力的并联关系计算出每个血管区域的总阻力。
[0091]
步骤303：根据血管区域的总阻力和血管区域中出口的横截面尺寸，获取血管区域中出口的阻力。
[0092]
具体的，这里的血管区域中出口的阻力，指的是血管区域中每一个出口各自的阻力。
[0093]
示例性的，假设一个三维模型中存在1，2和3三个血管区域，一共存在9个出口，血管区域1中包括出口1-3，血管区域2中包括出口4-6，血管区域3中包括出口7-9，接下来基于这一假设对步骤301-步骤303进行说明。
[0094]
步骤301中根据三维模型的总血流量和三维模型的入口处的压力差，获取的三维模型的总阻力，也就是出口1-9并联形成的阻力。步骤302根据血管区域相对于三维模型的血流量占比和三维模型的总阻力，获取的血管区域的总阻力，也就是血管区域1-3各自的总阻力，具体的，血管区域1的总阻力为出口1-3并联形成的阻力。而在步骤303要确定的血管区域中出口的阻力，以血管区域1为例，也就是要确定出口1-3各自的阻力，因此在步骤303中要确定出9个阻力，对应三维模型的9个出口。
[0095]
应当理解，血管区域中出口的横截面尺寸越大，出口的阻力越小，血管区域中出口的横截面尺寸越小，出口的阻力越大。
[0096]
可选的，根据血管区域中出口的半径、直径或横截面积中的一项或多项获取出口的阻力。
[0097]
在一种可能的实现方式中，本技术实施例通过下式获取血管区域中出口的阻力：
[0098][0099]
其中，i表示血管区域中的第i个出口，ri是第i个出口的阻力，r
t
是血管区域的总阻力，ri是第i个出口的半径，m是血管区域中的出口数目。
[0100]
可以看出，本技术实施例通过目标点的血液流速确定出血管区域的血流量占比，根据血流量和入口处的压力差计算出三维模型的总阻力，再根据血流量占比分配血管区域的总阻力，最后根据出口半径确定出口的阻力，提高了确定血管的出口阻力的准确性。
[0101]
步骤103：根据血液流速测量数据和多个血管区域中出口的阻力，确定待评估血管的流体控制方程的边界条件。
[0102]
流体控制方程在数学上大多是由非线性偏微分方程藕合而成的方程组，具体的，流体的流动由三个基本的物理原理控制，即质量守恒定律，牛顿第二定律和能量守恒定律，反映对应定律的方程可以被理解为流体控制方程。
[0103]
本技术实施例不对流体控制方程的选择进行限定，在一种可能的实现方式中，流体控制方程为纳维-斯托克斯(navier-stokes,ns)方程。
[0104]
可选的，将血管中的血液近似为不可压缩牛顿流体。其中，不可压缩流体指的是在流动过程中密度不发生变化的流体，牛顿流体指的是流体相邻的两层平行流动的液体间产生的剪切应力与垂直于流动方向的速度梯度成正比的流体。在这种近似下，可以采用非稳态不可压缩ns方程：
[0105][0106]
其中，u是血流速度，t是时间，是血流速度u相对于时间t的偏导数，ρ是血液的密度，是哈密顿算子，σ是柯西应力张量。ω是ns方程的空间域，即三维模型所处的空间，(0,t]是ns方程的时间域，t是周期时长，σ的计算满足下式：
[0107]
σ＝-pi 2με(u)
[0108]
其中，p是血管的压力，i是3
×
3单位矩阵，μ是动态粘度，ε(u)是变形张量，ε(u)的计算满足下式：
[0109][0110]
其中，u
t
表示u的转置。
[0111]
应当理解，上述非稳态不可压缩ns方程所涉及的各项物理量中，血流速度u和压力p为方程中的未知量，也就是本技术实施例要进行求解的物理量。
[0112]
本技术实施例中确定出的血管区域中的出口阻力，即ri，主要用于流体控制方程的边界条件的确定。在一种可能的实现方式中，通过弹性腔模型来模拟血管区域中的每一个出口，根据血液流速测量数据和每个血管区域中出口的阻力，确定每个血管区域中每个出口对应的弹性腔模型的参数。
[0113]
弹性腔模型是一种通过电路模型模拟人体中血管的方法，这种方法将局部血管系统简化为阻力元件、顺应性元件和惯性元件。其中，阻力元件用于模拟血管引起的阻力，通常为电阻，用r表示；顺应性元件，也可以称为弹性元件，用于模拟血管的收缩舒张，通常为电容，用c表示；惯性元件用于模拟血液流动的惯性，通常为电感，用l表示。值得注意的是，弹性腔模型中通常都包括阻力元件和顺应性元件，但可以不包括惯性元件。
[0114]
本技术实施例不对弹性腔模型的选取进行限定，在一种可能的实现方式中，通过三元弹性腔模型对血管区域的出口进行模拟。图6为三元弹性腔模型的示意图，下面结合图6对三元弹性腔模型进行说明。
[0115]
如图6所示，可以看出其中包括两个阻力元件和一个顺应性元件，为了便于描述，将表述为第一阻力，表述为第二阻力，ci是动脉顺应性，pi(t)是压力，具体在本技术实施例中是第i个出口的出口边界压力。
[0116]
在一种可能的实现方式中，通过下式计算弹性腔模型的参数：
[0117][0118][0119][0120]
其中，c
t
是血管区域的总顺应性，在一种可能的实现方式中，将c
t
假设为一个常数。
[0121]
具体的，和ci也就是本技术实施例需要确定的弹性腔模型的参数。
[0122]
在得出弹性腔模型的参数后，根据血管区域中每个出口对应的弹性腔模型的参数，确定流体控制方程的出口边界条件。
[0123]
在本技术实施例采取的弹性腔模型为三元弹性腔模型的情况下，出口边界条件满足下式：
[0124][0125]
其中，pi(t)是第i个出口的出口边界压力，p
ib
(t)是第i个出口处的远端压力，qi(t)是第i个出口处的血流量。在一种可能的实现方式中，p
ib
(t)可以被假设为零或一个常
数。在每个出口的阻力可以通过步骤102确定的情况下，结合前文所述的欧姆定律可以得出qi(t)。因此，pi(t)为上述出口边界条件中的未知物理量。
[0126]
本技术实施例通过弹性腔模型来模拟出口，并通过出口的阻力确定弹性腔模型的参数，从而确定流体控制方程的出口边界条件，提高了得出的出口边界条件的准确性。
[0127]
值得注意的是，如果要求解控制方程，还需要方程具有已知的入口边界条件，壁面边界条件和初始边界条件。
[0128]
本技术实施例不对入口边界条件，壁面边界条件和初始边界条件的确定方式进行限定。在一种可能的实现方式中，入口边界条件和壁面边界条件通过下式确定：
[0129][0130]
其中，vi是入口处的流入速度，γi表示入口，γw表示壁面，应当理解，壁面即血管壁。在一种可能的实现方式中，通过经颅多普勒超声获取入口处的流入速度，在另一种可能的实现方式中，通过流量守恒定律得出入口的总血流量，再根据总流量和入口处的截面积计算流入速度，流入速度采用下式计算：
[0131]
vi＝n
·
q/s
[0132]
其中，n是入口的内向法线，s是入口的横截面积，q是入口处的血流量。在一种可能的实现方式中，初始条件通过下式确定：
[0133]
u|
t＝0
＝u0,ω
[0134]
其中，u0是所述血液在初始时间的速度，在一种可能的实现方式中，将u0假设为0，在另一种可能的实现方式中，也可以将u0假设为其它常数。
[0135]
步骤104：根据边界条件对流体控制方程求解，获得待评估血管的血流数值模拟结果。
[0136]
在一种可能的实现方式中，本技术实施例对流体控制方法求解使用的是基于有限元分析的数值求解方法，因此，在进行求解之前需要先将三维模型划分为多个网格。
[0137]
其中，有限元分析方法将方程的求解域看成是由许多个小的互连子域组成，对每一单元假定一个合适的近似解，然后推导求解这个域总的满足条件，从而得到问题的解。在求解本技术实施例的ns方程时，也可以理解为整个三维模型为求解域，三维模型中的每一个网格为一个子域。
[0138]
本技术实施例不对网格划分的参数进行限定，在一种可能的实现方式中，对三维模型做非结构化四面体网格划分。其中，非结构化指的是网格区域内的内部点不具有相同的毗邻单元，四面体指每个网格的形状为四面体。图7是本技术实施例提供的网格划分的示意图，可以看出，图7采用了非结构化四面体网格划分。
[0139]
在一种可能的实现方式中，将三维模型划分为4.73
×
106或2.05
×
107个网格。
[0140]
具体的，划分网格的过程中需要先设置划分所需的参数，在完成划分之后还要对网格进行细化和光滑处理，并在软件内检查网格质量。
[0141]
在一种可能的实现方式中，使用ansys中的icem模块对三维模型进行网格划分。
[0142]
完成网格划分后，根据多个网格，在时间域上和空间域上分别离散流体控制方程，得到稀疏非线性系统，然后根据边界条件对稀疏非线性系统求解。
[0143]
在一种可能的实现方式中，在空间域上用有限元方法离散流体控制方程，在时间
域上用欧拉方法离散流体控制方程。具体的，在空间域上将流体控制方程离散到前文所述的三维模型中的网格，而欧拉方法则是一种逐次迭代进行求解的方法。
[0144]
可选的，在空间域上通过p1-p1有限原方法离散流体控制方程，在时间域上通过隐式向后欧拉方法离散流体控制方程。其中，两个p1分别指对血流速度u和压力p都使用p1有限元方法，具体的，p1有限元方法指的是一阶线性的离散。
[0145]
具体的，离散后的流体控制方程为一个大型稀疏非线性系统。
[0146]
在完成方程的离散后，流体控制方程已经可以通过计算机进行求解。在一种可能的实现方式中，使用天河二号a超算(tianhe-2a)进行求解。
[0147]
在一种可能的实现方式中，采用牛顿-克雷洛夫-施瓦兹(newton-krylov-schwarz,nks)算法来求解流体控制方程。
[0148]
具体的，在得出稀疏非线性系统后，也可以先通过nks算法对稀疏非线性系统进行进一步的处理，再进行求解。
[0149]
具体的，nks算法指的是将用于求解非线性方程系统的牛顿(newton)算法，与克雷洛夫(krylov)子空间技术结合，得到用于求解非线性系统的牛顿-克雷洛夫(newton-krylov)子空间迭代方法，再在迭代的过程中与施瓦兹(schwarz)预条件子技术相结合，得到的一种提高方程系统可解性的方法。这里的子空间，可以理解为上述有限元方法中子域的另一种表述。
[0150]
本技术实施例通过nks算法对稀疏非线性系统进行处理，提高了计算机求解流体控制方程的计算效率。
[0151]
本技术实施例对ns方程进行求解后，由于血流速度u和压力p是ns方程中的未知参数，因此可以直接求解出血流速度u和压力p。
[0152]
在求解出血流速度u和压力p的基础上，可以基于血液的流体力学性质通过血流速度u和压力p计算出血液的其它流体力学参数，得出血流数值模拟结果。在一种可能的实现方式中，可以计算血液的剪应力，剪应力指的是物体由于外因而变形时，在它内部任一截面的两方出现的相互作用力。
[0153]
可选的，为了使求解流体控制方程得出的数值模拟结果更为形象，可以对数值模拟结构进行后处理，用图像的形式表现数值模拟结果。在一种可能的实现方式中，使用paraview软件进行后处理。
[0154]
本技术实施例根据目标点的血液流速和三维模型的入口处的压力差，确定血管区域中出口的阻力，提高了确定出口的阻力时的准确性，并根据出口的阻力确定流体控制方程的边界条件，提高了求解流体控制方程得出的血流数值模拟结果的准确性。
[0155]
图8示出了根据本技术实施例的血流数值模拟装置的示意性框图。如图8所示，所述装置800包括：
[0156]
第一获取模块810，用于获取待评估血管的三维模型和血液流速测量数据；血液流速测量数据包括三维模型中多个目标点的血液流速；三维模型包括多个血管区域；
[0157]
第一确定模块820，用于根据多个目标点的血液流速和三维模型的入口处的压力差，确定每个血管区域中出口的阻力；
[0158]
第二确定模块830，根据血液流速测量数据和多个血管区域中出口的阻力，确定待评估血管的流体控制方程的边界条件；
[0159]
第二获取模块840，用于根据边界条件对流体控制方程求解，获取待评估血管的血流数值模拟结果。
[0160]
在一种可能的实现方式中，第一确定模块820用于：
[0161]
根据各个目标点的血液流速分别获取各个目标点的血流量；
[0162]
根据多个目标点的血流量，确定每个血管区域相对于三维模型的血流量占比；
[0163]
根据每个血管区域相对于三维模型的血流量占比和三维模型的入口处的压力差，确定每个血管区域中出口的阻力。
[0164]
在一种可能的实现方式中，第一确定模块820用于：
[0165]
根据三维模型的总血流量和三维模型的入口处的压力差，获取三维模型的总阻力；
[0166]
针对每个所述血管区域，根据血管区域相对于三维模型的血流量占比和三维模型的总阻力，获取血管区域的总阻力；
[0167]
根据血管区域的总阻力和血管区域中出口的横截面尺寸，获取血管区域中出口的阻力。
[0168]
在一种可能的实现方式中，流体控制方程为纳维-斯托克斯方程，边界条件包括出口边界条件，第二确定模块830用于：
[0169]
根据血液流速测量数据和每个血管区域中出口的阻力，确定每个血管区域中每个出口对应的弹性腔模型的参数；
[0170]
根据多个血管区域中每个出口对应的弹性腔模型的参数，确定流体控制方程的出口边界条件。
[0171]
在一种可能的实现方式中，第二获取模块840用于：
[0172]
将三维模型划分为多个网格；
[0173]
根据多个网格，在空间域上和时间域上分别离散流体控制方程，得到稀疏非线性系统；
[0174]
根据边界条件对稀疏非线性系统求解。
[0175]
在一种可能的实现方式中，第二获取模块840用于：
[0176]
采用牛顿-克雷洛夫-施瓦兹算法对流体控制方程求解。
[0177]
在一种可能的实现方式中，第一获取模块810用于：
[0178]
获取待评估血管的断层扫描影像数据；
[0179]
根据断层扫描影像数据，获取待评估血管的三维模型。
[0180]
图9是本技术实施例提供的血流数值模拟装置的结构示意图。如图9所示，血流数值模拟装置900包括：至少一个处理器90(图9中仅示出一个)处理器、存储器91以及存储在所述存储器91中并可在所述至少一个处理器90上运行的计算机程序92，所述处理器90执行所述计算机程序92时用于实现上述任意各个血流数值模拟方法实施例(比如图1中的方法)中的步骤。
[0181]
所称处理器90可以是中央处理单元(central processing unit，cpu)，该处理器90还可以是其他通用处理器、数字信号处理器(digital signal processor，dsp)、专用集成电路(application specific integrated circuit，asic)、现成可编程门阵列(field-programmable gate array，fpga)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、
分立硬件组件等。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。
[0182]
所述存储器91在一些实施例中可以是所述装置900的内部存储单元，例如血流数值模拟装置900的硬盘或内存。所述存储器91在另一些实施例中也可以是所述装置900的外部存储设备，例如所述装置900上配备的插接式硬盘，智能存储卡(smart media card,smc)，安全数字(secure digital,sd)卡，闪存卡(flash card)等。进一步地，所述存储器91还可以既包括所述装置900的内部存储单元也包括外部存储设备。所述存储器91用于存储操作系统、应用程序、引导装载程序(bootloader)、数据以及其他程序等，例如所述计算机程序的程序代码等。所述存储器91还可以用于暂时地存储已经输出或者将要输出的数据。
[0183]
所属领域的技术人员可以清楚地了解到，为了描述的方便和简洁，仅以上述各功能单元、模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分配由不同的功能单元、模块完成，即将所述装置的内部结构划分成不同的功能单元或模块，以完成以上描述的全部或者部分功能。实施例中的各功能单元、模块可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中，上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。另外，各功能单元、模块的具体名称也只是为了便于相互区分，并不用于限制本技术的保护范围。上述系统中单元、模块的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。
[0184]
本技术实施例还提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现可实现上述各个方法实施例中的步骤。
[0185]
本技术实施例提供了一种计算机程序产品，当计算机程序产品在血流数值模拟装置上运行时，使得血流数值模拟装置执行时实现可实现上述各个方法实施例中的步骤。
[0186]
所述集成的单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本技术实现上述实施例方法中的全部或部分流程，可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一计算机可读存储介质中，该计算机程序在被处理器执行时，可实现上述各个方法实施例的步骤。其中，所述计算机程序包括计算机程序代码，所述计算机程序代码可以为源代码形式、对象代码形式、可执行文件或某些中间形式等。所述计算机可读介质至少可以包括：能够将计算机程序代码携带到拍照装置/终端设备的任何实体或装置、记录介质、计算机存储器、只读存储器(rom，read-only memory)、随机存取存储器(ram，random access memory)、电载波信号、电信信号以及软件分发介质。例如u盘、移动硬盘、磁碟或者光盘等。在某些司法管辖区，根据立法和专利实践，计算机可读介质不可以是电载波信号和电信信号。