CN109102568A - 基于区域分解的血管血流模拟方法及相关装置 - Google Patents
基于区域分解的血管血流模拟方法及相关装置 Download PDFInfo
- Publication number
- CN109102568A CN109102568A CN201810525629.6A CN201810525629A CN109102568A CN 109102568 A CN109102568 A CN 109102568A CN 201810525629 A CN201810525629 A CN 201810525629A CN 109102568 A CN109102568 A CN 109102568A
- Authority
- CN
- China
- Prior art keywords
- zoning
- grid
- sub
- blood
- zonings
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
- G16H50/00—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics
- G16H50/50—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics for simulation or modelling of medical disorders
Landscapes
- Engineering & Computer Science (AREA)
- Medical Informatics (AREA)
- Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Public Health (AREA)
- Epidemiology (AREA)
- Pathology (AREA)
- Databases & Information Systems (AREA)
- Data Mining & Analysis (AREA)
- General Health & Medical Sciences (AREA)
- Primary Health Care (AREA)
- Biomedical Technology (AREA)
- Computer Graphics (AREA)
- Geometry (AREA)
- Software Systems (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本申请公开了一种基于区域分解的血管血流模拟方法、血流模拟装置及计算机可读存储介质,其中方法包括获取血管的特征数据;根据特征数据构建血管的三维模型,三维模型定义有计算区域;对计算区域进行离散化处理,生成刻画计算区域的网格;将计算区域划分为多个子计算区域,每个子计算区域的网格点数量一致;基于网格,同时对至少两个子计算区域进行计算,以获得计算区域的血流参数。本申请可实现对血管中血流的精确高效的模拟。
Description
技术领域
本申请涉及血流数值模拟领域,特别是涉及一种基于区域分解的血管血流模拟方法、血流模拟装置以及计算机可读存储介质。
背景技术
血液流动的特征在一定程度能够反映出血管是否存在疾病以及患者是否存在血流动力学改变导致的疾病,例如血流储备分数(FFR)可反映缺血风险,血流速度可反映血管堵塞程度等。因而对血流进行模拟分析已经成为了当前血管疾病预防和诊断领域的研究热点
早期由于计算机能力的限制,在对血流进行模拟时均进行了很多简化,例如物理模型的简化,离散化网格的简化;这种简化计算虽然计算时效性高,但无法得出血流的一些重要特征,模拟的精度不高。而当前计算机硬件水平获得发展后,却没有与计算机能力匹配的模拟方法,无法同时保证血流模拟的效率和精度的提高,即当前对血流进行模拟的精度和效率依旧不高。
发明内容
本申请提供一种基于区域分解的血管血流模拟方法、血流模拟装置以及计算机可读存储介质,以解决现有技术中对血流进行模拟的精度和效率不高的问题。
为解决上述技术问题,本申请提出一种基于区域分解的血管血流模拟方法,包括获取血管的特征数据;根据特征数据构建血管的三维模型,三维模型定义有计算区域;对计算区域进行离散化处理,生成刻画计算区域的网格;将计算区域划分为多个子计算区域,每个子计算区域的网格点数量一致;基于网格,同时对至少两个子计算区域进行计算,以获得计算区域的血流参数。
为解决上述技术问题,本申请提出一种血流模拟装置,包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序,处理器执行计算机程序时实现上述方法的步骤。
为解决上述技术问题,本申请提出一种计算机可读存储介质,其上存储有计算机程序,计算机程序被处理器执行时实现上述方法的步骤。
本申请血管血流模拟方法包括:获取血管的特征数据,然后根据特征数据构建血管的三维模型,在三维模型上定义有计算区域;然后对计算区域进行离散化处理,生成刻画该计算区域的网格;将计算区域划分为多个子计算区域,且每个子计算区域的网格数量一致;使网格点数量一致,是为了使得在同时对多个子计算区域求解时,能够保证计算量的均衡;然后对多个子计算区域中的至少两个进行同时计算,最终获得计算区域的血流参数。本申请方法中将大规模的计算区域划分为多个平均的子计算区域,然后同时并行求解,使得大规模的计算问题变得分布式、小型化,同时求解使得计算更高效,而单个子计算区域的问题规模较小,计算精度会更高;因而本申请的方法提高了血流模拟的精度和效率。
附图说明
图1是血流模拟系统的结构示意图;
图2是血流模拟方法的流程示意图;
图3是本申请血管血流模拟方法一实施例的流程示意图;
图4是图3所示的实施例中计算区域的分区示意图;
图5是本申请血管血流模拟方法另一实施例的流程示意图;
图6是本申请血流模拟装置一实施例的结构示意图;
图7是本申请计算机可读存储介质一实施例的结构示意图。
具体实施方式
为使本领域的技术人员更好地理解本申请的技术方案,下面结合附图和具体实施方式对发明所提供的一种基于区域分解的血管血流模拟方法、血流模拟装置以及计算机可读存储介质做进一步详细描述。
对血管中血流进行模拟是利用流体力学的方法来模拟血液的流动,获得包括血流的力学特征血流参数。本申请中构建一实现血流模拟的计算机系统,请参阅图1,图1是血流模拟系统的结构示意图,该血流模拟系统100包括以下几大模块。
数据导入模块11。
用于获取血管的特征数据,例如当该血流模拟系统100应用于患者的疾病诊断时,可通过数据导入模块11导入患者血管的特征数据,包括影像数据、生理数据等。
其中影像数据可以是基于核磁造影MRA、计算机断层造影CTA、数字减影造影DSA或超声弹性成像等技术获得。即首先由其他设备或技术获得血管的影像数据或生理数据,然后血流模拟系统100再根据该影像数据或生理数据进行血流模拟。
三维建模模块12。
用于根据由数据导入模块11获取的血管特征数据构建血管的三维模型。具体即由血管的影像数据来构建血管的三维模型。在构建本三维模型的过程中,可以对三维模型进行平滑处理,使得三维模型更符合血管本身的形状;还可以对三维模型进行修剪处理,例如在对心脏的血管进行建模时,仅保留主动脉的部分,修剪掉三维模型中的其他血管。
在用户使用本血流模拟系统100时,三维建模模块12可将所构建的三维模型可视化向用户呈现,由用户判断该三维模型是否可行,若不可行,则重新构建三维模型,例如包括血管壁在内的血管模型。
在完成三维模型的构建后,三维模型上定义有计算区域,即需要进行血流模拟的目标区域,计算区域可以是整个三维模型,也可以是三维模型中某部分区域,例如可对患者心脏的整个血管建立三维模型,但仅对主动脉中的血流进行模拟,此时计算区域即主动脉。
网格生成模块13。
用于对计算区域进行离散化处理,生成刻画计算区域的网格。所生成的网格能够体现出计算区域的整个形状。
由于进行血流模拟的计算过程从物理数学意义上来看是对血流动力学控制方程求解,而求解血流动力学控制方程的求解即求解偏微分问题,在求解偏微分问题时一般采用离散化处理方法,例如有限元方法、有限体积方法、间断有限元方法等。因此,需要将计算区域划分出用于数值模拟的离散网格。
网格生成模块13可以将计算区域划分为结构化网格或非结构化网格,由于血管几何形态的复杂性,本申请中采用非结构网格。对于三维结构的计算区域,可采用Delaunay准则、前沿推进算法(Advancing Front Method)或Shephard-Yerry算法划分出四面体网格单元;也可采用映射法、子映射法、扫琼法、基于栅格法、中轴面法、Plastering法或WhiskerWeaving法划分出六面体网格单元。
边界条件模块14。
用于确定对计算区域进行模拟计算时的边界条件,其中边界条件可以人为设置,也可以根据血管的生理数据确定。
模型求解模块15。
基于所生成的网格以及所确定的边界条件,对计算区域进行求解,即求解施加了一定边界条件的血流动力学控制方程,在求解时,首先也需要对该血流动力学控制方程进行离散化处理,然后再通过算法对离散后的方程进行求解。最终获取该计算区域内的血流参数,包括血流速度、血压、剪切力、血管壁变形等信息。
上述模块之间的连接关系已经体现在对模块的功能描述中,在此不再赘述。通过血流模拟系统100能够实现对血流的模拟,最终获得计算区域内血流的参数,根据血流参数即可进行疾病诊断,分析病变进行病理研究,还可指导心脑血管的相关手术,例如血管搭桥手术,血管支架放入手术,还可对术后情况进行模拟评估,优化血管搭桥手术中使用桥梁和血管支架放入手术中支架结构的优化设计。
为了更为方便的实现用户交互,血流模拟系统100中还包括以下模块。
可视化模块16。
用于显示血流模拟后的血流参数,可结合三维模型进行模拟显示,例如对于血流参数FFR值或压力值,可基于三维模型呈现不同颜色分布图显示给用户,使得血流模拟结果更加直观。
报告生成模块17。
根据血流模拟结果生成报告,给出疾病诊断或治疗建议。在报告生成之前,可以先向用户呈现血流模拟结果,由用户判断该血流模拟结果是否正常,若不正常,例如医生认为该结果与其他方式的诊断结果差异较大,则可重新进行网格生成,模拟计算;若正常,则继续生成报告。
在血流模拟系统100进行实际应用时,以该血流模拟系统面向医生的应用为例。
医生可在该血流模拟系统100的网页界面或应用界面上进行操作,输入患者血管的特征信息,特征信息由数据导入模块11导入到系统中。
系统中三维建模模块12根据特征信息生成血管的三维模型,并在界面上进行可视化呈现。
此时医生判断该三维模型是否可行,并将判断结果反馈给系统;若可行,则系统网格生成模块13对计算区域进行离散化处理,生成刻画计算区域的网格。
模型求解模块15基于网格对计算区域进行求解,在求解过程中由边界条件模块14确定计算区域的边界条件及参数;最终获得计算区域的血流参数。
可视化模块16将血流参数进行可视化呈现。
此时医生判断该模拟结果,血流参数是否有异常,并将判断结果反馈给系统;若无异常,则系统的报告生成模块17根据模拟结果生成疾病诊断报告或治疗方案报告。
上述各个模块构成实现血流模拟的系统,从方法的角度来看,则主要通过以下步骤来实现血流模拟。请参阅图2,图2是血流模拟方法的流程示意图。
S11:获取血管的特征数据。
S12:根据特征数据构建血管的三维模型,在三维模型上定义计算区域。
S13:对计算区域进行离散化处理,生成刻画计算区域的网格。
S14:基于网格,对计算区域进行计算,从而获得计算区域的血流参数。
上述步骤均对应于血流模拟系统100中的各个模块,步骤中的具体过程不再赘述,其中,步骤S11对应数据导入模块11、步骤S12对应三维建模模块12、步骤S13对应网格生成模块13、步骤S14对应模型求解模块15和边界条件模块14。
上述步骤S11-S14为实现血流模拟的基本步骤,即本申请中血流模拟方法实施例均基于上述步骤S11-S14实现。本申请为了提高血流模拟的精度和效率,从多个方面提出对血流模拟过程进行了优化。例如图3和图5所示的实施例,图3所示的实施例在上步骤S13离散化处理中引入区域分解;并基于区域分解,在上述步骤S14中引入并行计算,以提高模拟效率。图5所示的实施例在上述步骤S14中引入流固全耦合计算,以提高模拟精度。
以下分别对本申请的两实施例进行具体介绍。首先参阅图3,图3是本申请血管血流模拟方法一实施例的流程示意图。本实施例中对计算区域进行分解,将大规模的计算区域分解为多个小规模的计算区域进行分别独立计算,从而提高计算效率,并且采用并行算法,保证了计算精度。本实施例是基于区域分解的血管血流模拟方法,其对血流进行模拟的方法包括以下步骤。
S21:获取血管的特征数据。
S22:根据特征数据构建血管的三维模型。其中,三维模型上定义有计算区域。
S23:对计算区域进行离散化处理,生成刻画计算区域的网格。
S24:将计算区域划分为多个子计算区域。
对于本步骤可结合图4理解,图4是图3所示的实施例中计算区域的分区示意图。将计算区域划分为多个子计算区域,即将大规模的计算区域划分为多个小规模计算区域,以减小计算规模,提高计算效率。
在本步骤S24中进行子计算区域划分时,需要保证每个子计算区域的网格点数量一致,由于计算规模是以网格点数量来评价的,为了保证整体的计算效率,本步骤中对计算区域的计算规模进行均等划分,即每个子计算区域的网格点数量一致。
当然在实际应用中,是无法保证每个子计算区域的网格点数量完全一样,因此本步骤中每个子计算区域之间的网格点数量差值在可接受范围内即可,即网格点数量一致并不限定为网格点数量完全一样,保证基本一致即可。
S25:基于网格,同时对至少两个子计算区域进行计算,以获得计算区域的血流参数。
在将计算区域划分为多个子计算区域后,子计算区域之间的计算相互独立,因而本步骤S25中,同时对至少两个子计算区域进行计算,从而获得计算区域的血流参数,即多个子计算区域并行计算,提高了计算效率。若计算能力足够强大,也可对所有的子计算区域同时并行计算。
本实施例中首先对计算区域进行区域划分,然后再对多个区域进行并行计算,提高了计算效率,且扩大了计算规模。本实施例充分利用高性能计算机提高血流动力学分析的精度和效率,满足了临床医学应用的准确性和实时性要求。
此外,本实施例中还进一步提出对上述步骤的优化,以提高计算效率和精度。
例如步骤S23可包括以下两步:
S231:对计算区域进行离散化处理,生成刻画计算区域的粗网格;
S232:对粗网格进行加密处理,生成刻画计算区域的细网格。
上述步骤S231-S232表示多层次的实现网格生成,包括以下多种情况。
可以是首先对计算区域进行小规模的离散化处理,生成一套粗网格,且保证该粗网格能够刻画计算区域的形状;然后再对该粗网格进行加密处理,继而生成刻画计算区域的细网格。
对应在实际应用中,本实施例的血流模拟计算规模较大,一般将其上传至超级计算中心进行计算,因而在本地计算机上应用血流模拟系统实现血流模拟方法时,可首先在本地计算机上生成一套粗网格,然后再将该粗网格上传到超级计算中心进行网格的加密等处理。该过程可实现在本地计算机上获得刻画计算区域较好的粗网格,基于该质量较好的粗网格进行加密,能够进一步保证后续网格质量。
还可以是对粗网格进行多次加密处理,例如对初始的粗网格进行第一次加密处理,得到第一细网格;然后对第一细网格进行第二次加密处理,得到第二细网格;第一细网格和第二细网格相对来说即分别为步骤S231中的粗网格和S232中的细网格。
采用这种方式获得多水平的粗网格和细网格,在计算时即可通过多水平的计算来提高计算的效率和精度。具体来说,即步骤S25中对计算区域进行计算可包括以下两步。
S251:基于粗网格,同时对至少两个子计算区域进行计算,获得至少两个子计算区域的粗网格血流参数。
S252:基于细网格,并根据至少两个子计算区域的粗网格血流参数,同时对至少两个子计算区域进行计算,以获得至少两个子计算区域的细网格血流参数。
该多水平计算可以是瀑布型,对应上述步骤S251和S252,首先对粗网格进行计算,然后根据粗网格计算结果对细网格进行计算,即粗网格计算—细网格计算;还可以是V循环型,首先对粗网格进行计算,然后根据粗网格计算结果对细网格进行计算,其次根据细网格计算结果对粗网格进行再次计算,即粗网格计算—细网格计算—粗网格计算;还可以是W循环型,即粗网格计算—细网格计算—粗网格计算—细网格计算—粗网格计算。
对于粗细网格多水平的计算过程,粗细网格之间有数值互换的步骤,即插值计算过程。为了保证边界处插值的精度,本实施例中在步骤S231-S234构造粗细网格时,保证所生成的细网格和粗网格在计算区域的边界具有相同的网格点数量。
以上对粗网格进行加密处理所采用的是一致加密算法,能够在不改变网格质量的情况下对网格快速加密,例如对于三角形的网格单元,连接其各条边的中点,从而将一个三角形网格单元划分为四个三角形网格单元;对于三维的四面体网格单元,同理可将其分给为八个四面体单元。
在对粗网格进行加密的过程中,还可对其进行粗化处理,使得所生成的细网格中保留有粗网格中部分网格点的几何信息,即保留粗网格中的几何信息。首先选择保留一些在几何上重要的网格点,例如曲面上所有的点,平面边上两端点和内部的等距点;然后将没有选择保留的网格点删除,具体过程采用Edge-contraction算法进行迭代筛选。在删除网格点后,对整个网格进行优化以保证网格质量。
对于步骤S24的分区,本申请还提出多水平分区技术,以获得规则的均衡的多个子计算区域,具体包括以下步骤。
S241:对计算区域进行第一次区域划分,获得多个第一子计算区域,每个第一子计算区域的网格点数量一致。
S242:同时对至少两个第一子计算区域进行第二次区域划分,获得多个第二子计算区域,每个第二子计算区域的网格点数量一致。
基于上述步骤S241和S242,还可对计算区域进行多次区域划分,直至获得所需要的分区数量。多次分区可使得每次所分区的数量不大,保证每次分区的子计算区域之间的网格点数量的一致性。因而可提高分区质量。
在步骤S25中采用并行算法对子计算区域进行求解,从而获得整个计算区域的计算结果。子计算区域由计算区域划分获得,子计算区域之间可重叠也可不重叠,对应采用不同的并行算法。
本实施例中,每相邻两个子计算区域之间具有重叠区域,对应采用重叠型区域分解算法。
本实施例中在对计算区域进行计算时,边界条件可以是人为设置,可以是根据邻近计算区域的计算所获得,也可以是根据所获取的生理数据确定。再此不做赘述。
本实施例血流模拟方法采用区域分解技术结合并行计算,提高了血流模拟的效率和精度。
请继续参阅图5,图5是本申请血管血流模拟方法另一实施例的流程示意图,本实施例中建立流固全耦合的物理数学模型,进行全耦合计算,即在进行血流模拟时考虑血流与血管之间的相互作用,因而本实施例提高了模拟计算精度;并且结合非线性系统求解算法,保证了计算的效率。本实施例是基于力学方程的血管血流模拟方法,本实施例对血流进行模拟的方法包括以下步骤。
S31:获取血管的特征数据。
S32:根据特征数据构建血管的三维模型。其中,三维模型中定义有计算区域。
S33:对计算区域进行离散化处理,生成刻画计算区域的网格。
本实施例中步骤S31-S33与上述对应步骤类似,具体不再赘述。
S34:构建计算区域的物理数学模型。
物理数学模型能够描述血液流动、血管壁变形及两者之间的相互作用力对血流参数的影响的物理现象。对该物理数学模型进行求解,所获得的血流参数,即代表了血流特征。物理数学模型即描述物理客观现象的控制方程,本实施例中所构建的物理数学模型包括全耦合的流体力学控制方程和固体力学控制方程。
根据流体力学控制方程可获得血液流速和压力等,流体力学控制方程包括:可压缩和不可压缩的Navier-Stockes方程及其对应的各种湍流模型,如雷诺时均、大涡模拟等。
根据固体力学控制方程可获得血管壁位移等,固体力学控制方程包括:线性弹性和非线性弹性的固体本构方程,以及黏弹性、弹塑性、多孔介质等模型。
所构建的物理数学模型中流固全耦合,因而物理数学模型中还包括:流固界面条件。
对应计算区域中的边界,物理数学模型中还包括:滑移或无滑移固壁边界条件,阻尼型出流边界条件、无压力边界条件等,不同边界条件对应不同的物理现象,且直接影响问题的计算复杂度和适应性。
S35:基于网格,对物理数学模型进行计算,获得计算区域的血流参数和血管参数。
上述步骤S34所构建的物理数学模型将流体力学控制方程和固体力学控制方程耦合到一个方程,本步骤S35在该方程中进行求解,计算过程无需两个方程之间的迭代,实现同步求解,保证求解精度。
具体来说,本步骤即基于网格,对流体力学控制方程和固体力学控制方程进行同时求解,且在求解时统一计算区域中流固交界面处血流网格的网格点信息和血管网格的网格点信息。
对于血流和血管壁来说,两者在流固交界面相互作用,且为相对作用力,例如在剪切力,位移等参数上具有一定的关系。因而在物理数学模型中采用流固界面条件来模拟交界面的情况;而在计算区域的网络中对流固交界面的网格点信息进行统一,具体来说,若交界面处血流网格和血管网格匹配,则通过点对点的信息转换来进行网格点信息的统一;若交界面处血流网格和血管网格不匹配,则利用插值法来进行网格点信息的统一,可以采用基于有限元基函数的线性和二次插值方法、径向基函数插值法、Mortar元方法等。
上述步骤建立流固耦合的物理数学模型,计算该物理数学模型以更为精确的对血管中的血流进行模拟。
在进行流固全耦合同步求解时,问题规模非常大,因此本实施例还提出利用非线性系统求解算法对物理数学模型进行计算,具体可采用牛顿-克雷洛夫-施瓦茨(Newton-Krylov-Schwarz)算法,包括以下步骤。
S351:对物理数学模型进行离散化处理,获得非线性方程组。
首先对物理数学模型进行离散化处理,即将偏微分方程离散化为非线性方程组。其中,对于流体力学控制方程,可采用稳定化的P1-P1元、经典Taylor-Hood P2-P1元等方法;对于固体力学控制方程,对给定弱对称应力张量的混合形式采用PEERS元,对位移采用非协调P1元;对于流固界面条件,离散格式采用mortar或杂交技术以及基于Lagrange乘子的新型方法。
S352:利用牛顿法对非线性方程组进行求解。
在本步骤的求解过程中,可采用线性搜索和可行域技术确定搜索方向和步长,在可行域中进行线性搜索;在牛顿法求解的迭代过程中,可采用网格序列法和非线性预处理技术,使得本步骤的非线性迭代过程具有网格无关的收敛性;而对于牛顿法中的Jacobian矩阵,则采用多色排序有限差分法、自动微分技术、Jacobian-free方法或显式生成等策略构造生成。
S353:利用克雷诺夫子空间迭代法对牛顿法中的线性方程组进行求解。
具体来说,在本步骤中采用GMRES,或短重现式(Short-Recurrence)的Lanczos双正交化方法求解牛顿法中具有非对称矩阵的线性方程组。
S354:利用区域分解法构造线性方程组中的预条件子。
构造预条件子对步骤S353中的线性求解进行加速,本实施例中采用重叠型施瓦茨(Schwarz)算法,具体可利用调节扩展加性Schwarz算法或限制型加性Schwarz算法构造预条件子。
若在本实施例中将计算区域划分为多个重叠子计算区域,此时需要引入对子计算区域的计算方法,本实施例中采用直接法或迭代法,包括LU分解算法、不完全LU分解算法、Gauss-Seidel迭代法等。子区域的矩阵是稀疏的,其非零元的存储和访问可按点块(point-block)方式进行,即直接法或迭代法可按网格节点的顺序同时存储和访问节点上多个变量。采用直接法求解子区域问题时,可采用不同的子区域矩阵排序方式,包括NestedDissection、One-way Dissection、Reverse Cuthill-McKee和Quotient Minimum Degree等方法。
利用非线性系统求解算法对物理数学模型进行求解,得到血流参数及血管参数。血管参数中包括血管壁的移动数值,血管壁的移动会对网格有一定的影响,在下次计算时,需要考虑网格的变化。因此本实施例中还进行以下步骤。
S36:根据血管参数计算网格的变化参数,从而更新计算区域。
当血液和血管壁相互作用使得计算区域发生变化,相应的计算区域的网格也需要对应调整,本步骤中具体采用动网格技术来捕获网格的移动,根据血管参数,利用可移动网格算法计算网格的变化参数,从而根据变化参数构造新网格。
本步骤中利用可移动网格算法的计算是基于可移动网格框架,例如ALE框架,ALE框架结合拉格朗日格式和欧拉格式,拉格朗日格式在固体表达中网格可随材料变形而变形,欧拉格式在流体表达中网格不移动不至于引起网格扭曲和交错。具体来说,当基于ALE框架来进行计算时,流固界面附近的流体区域采用ALE参考系建模,而远离流固界面的流体采用欧拉参考系建模。可移动网格算法则具体可采用弹簧体法、偏微分方程法或插值法。例如,弹簧体法假设计算网格域为弹性介质,选取恰当的弹性模量并利用位移边界条件,通过求解线性或非线性的弹性方程,得到计算网格随便捷的运动或变形。
为了提高网格更新质量,在计算过程中,构造杨氏模量和泊松比来获得高质量的新网格,即在利用可移动网格算法进行计算时,结合所构造的杨氏模量和泊松比。此外,对于移动量较大的问题,当移动后的网格质量较差时,通过网格点局部调整的方法来调整网格点的分布以提高网格质量。并且,在根据变化参数构造新网格后,还可基于网格优化算法,例如Delaunay算法对新网格进行优化。
本步骤S36中构建了一个动网格方程,以描述网格的移动,对动网格方程进行计算可获得网格的变化参数,由于计算过程涉及到血管参数,因而可将动网格方程集成在步骤S34中所构建的物理数学模型中,进行同时求解,即本步骤S36与步骤S35并没有严格的先后关系,可同时进行。
本实施例所构建的物理数学模型可以采用以下形式。
流体动力学方程:
流固界面条件:
σf·nf=-σs·ns on Γinterface,
d=x on Γinterface,
固体动力学方程:
动网格方程:
其中,为流场的Cauchy应力张量,u表示血流速度,pf为血流压力,ρf为血液密度,μ为血液的粘性系数(血液视为牛顿流体时,对应的μ为常数,视为非牛顿流体时,μ为一个复杂的函数)。
d表示血管壁的位移,σs=λtrace(ε)I+2με为血管壁的应力张量,其中λ和μs为Lame系数,
x表示网格移动的位移,σm为网格模型的应力张量,形式与σs一样,但对应的Lame系数取值不同。
为流体计算域的入流边界和出流边界,为固体(管壁)的除与流体交界外的边界,如血管的外壁,Γinterface为流体与固体的交界面(血液-血管壁交界面);α为稳定化常数,根据实验数据设定具体的值。Ωf为流体计算区域,Ωs为固体计算区域。
对于上述方程中的应力张量的选取、粘性系数的构造及边界条件的选取、流体和固体方程中应力张量的选取均需要根据血液及血管壁的具体性质而定,在实际临床应用中,针对每个病例,其数值各不相同。
本实施例的方法构建了流固全耦合的物理数学模型,并对计算区域的网格变化进行更新,使得血流模拟更为精确;在模拟计算中引入非线性系统求解算法,以保证计算效率。
图3所示实施例涉及网格生成后对计算区域的进行区域分解,并基于区域分解进行并行计算;图5所示实施例涉及建立流固全耦合的物理数学模型,以进行流固全耦合计算,并针对性设置求解算法;二者所涉及到的技术可结合应用。
例如血流模拟系统中网格生成模块采用图3所示实施例的网格生成技术及区域分解技术,模型求解模块则引入图5所示的流固全耦合模型的构建,及对应算法的采用。对于血流模拟过程来说,实现了模拟效率和精度的极大提高。
上述血流模拟方法均可基于图1所示的软件构架实现,在硬件结构方面,请参阅图6,图6是本申请血流模拟装置一实施例的结构示意图。
本实施例血流模拟装置200包括处理器21和存储器22,存储器22中存储有可在处理器21上运行的计算机程序,处理器21执行计算机程序时能够实现上述血流模拟方法。
本实施例中处理器21为广义的处理器,其可以包括多个处理器,可以是设置在不同设备中的处理器,例如可以包括本地计算机中的处理器和超级计算中心的处理器群。本实施例中在进行血流模拟的并行计算时,有多少个子计算区域,则对应设置多少个处理器,以分别对子计算区域进行计算。
本实施例中处理器21可采用异构体系结构,在使用本实施例处理器21进行大规模数据的并行计算时,采用一系列并行加速技术,即根据不同计算的特征,将计算过程放置到不同的处理器上,实现处理器能力的最大化利用。
例如在求解方面,将非线性求解和线性求解时耗时较多的非线性离散函数的计算、稀疏矩阵向量乘、稀疏矩阵快速分解等核心计算模块移植到GPU、MIC或其它众核处理器上。采用多色或多scheduling策略改进稀疏矩阵的分解和回溯求解,一方面提升算法并行度,另一方面保持求解器的收敛效率。将边界条件等涉及大量分支操作的部分从区域内部独立出来,由逻辑处理能力强的部件进行计算,以使任务更合理地得到分配。
在处理器计算指令方面,优化指令级并行、线程级并行和进程级并行,采用数据重用、计算与访存重叠、数据融合和对界访问、数据合并传输、向量化、科学运算函数优化等技术优化众核处理器上的操作,提高实际运行时的浮点效率。
处理器的指令实现编程语言方面,对于GPU,MIC或其它众核处理器,使用CUDA,OpenACC,openCL语言于GPU上执行;使用OpenMP语言在MIC上执行;使用pthread,athread或其它函数包在其它众核处理器上执行。
在多个处理器并行处理大规模数据时,不仅仅需要在计算方面提高效率,还需在数据的传输方面提高效率。例如本实施例中采用了一系列大规模离散数据并行处理技术。
分块并行I/O技术:建立分块的数据结构,使处理器之间负载平衡。即将表示血管和血液的物理量的离散数据按计算区域分块且并行读入自或输出于一个或多个数据文件,分块个数与所使用的处理器个数相同。本实施例中分块并行I/O技术基于MPI-2(及以上版本)函数库,采用指定显式偏移量、独立文件指针或共享文件指针的方式实现。数据文件包括HDF5、VTK等格式。
放粗或加密输出技术,即将表示血管和血液的网格进行放粗或加密,并将原网格上的物理量插值到新网格上,然后使用所述分块并行I/O技术输出。
采用向量压缩技术或基于MPI_pack数据打包技术降低数据通信量和I/O数据量,并设计I/O与计算的重叠机制以解决大规模离散数据的I/O瓶颈。
本实施例血流模拟装置实现了对配置数千上万的计算核心的高级计算机的充分利用,通过软件和算法的配合充分调用了计算资源,提高血流动力学分析的精度和效率,实现了60%以上的并行可扩展效率,提高了血流动力学模拟精度的同时也减少了计算时间。
上述血流模拟方法以软件形式实现并作为独立的产品销售或使用时,可存储在一个电子设备可读取存储介质中,即,本发明还提供一种计算机可读存储介质,请参阅图7,图7是本申请计算机可读存储介质一实施例的结构示意图,计算机可读存储介质300中存储有计算机程序,该计算机程序被处理器执行时实现上述方法的步骤。计算机可读存储介质可以为U盘、光盘、服务器等。
以上所述仅为本申请的实施方式,并非因此限制本申请的专利范围,凡是利用本申请说明书及附图内容所作的等效结构或等效流程变换,或直接或间接运用在其他相关的技术领域,均同理包括在本申请的专利保护范围内。
Claims (10)
1.一种基于区域分解的血管血流模拟方法,其特征在于,所述方法包括:
获取所述血管的特征数据;
根据所述特征数据构建所述血管的三维模型,所述三维模型定义有计算区域;
对所述计算区域进行离散化处理,生成刻画所述计算区域的网格;
将所述计算区域划分为多个子计算区域,每个子计算区域的网格点数量一致;
基于所述网格,同时对至少两个所述子计算区域进行计算,以获得所述计算区域的血流参数。
2.根据权利要求1所述的方法,其特征在于,所述对所述计算区域进行离散化处理,生成刻画所述计算区域的网格,包括:
对所述计算区域进行离散化处理,生成刻画所述计算区域的粗网格;
对所述粗网格进行加密处理,生成刻画所述计算区域的细网格。
3.根据权利要求2所述的方法,其特征在于,所述对所述粗网格进行加密处理,包括:
对所述粗网格进行加密及粗化处理,使得所生成的细网格中保留有所述粗网格中部分网格点的几何信息。
4.根据权利要求2所述的方法,其特征在于,所述对所述粗网格进行加密处理,包括:
对所述粗网格进行加密处理,使得所生成的细网格和所述粗网格在所述计算区域的边界具有相同的网格点数量。
5.根据权利要求2所述的方法,其特征在于,所述基于所述网格,同时对至少两个所述子计算区域进行计算,包括:
基于所述粗网格,同时对至少两个所述子计算区域进行计算,获得所述至少两个子计算区域的粗网格血流参数;
基于所述细网格,并根据所述至少两个子计算区域的粗网格血流参数,同时对所述至少两个子计算区域进行计算,以获得所述至少两个子计算区域的细网格血流参数。
6.根据权利要求1所述的方法,其特征在于,所述将所述计算区域划分为多个子计算区域,包括:
将所述计算区域划分为多个子计算区域,且每相邻两个子计算区域之间具有重叠区域。
7.根据权利要求1所述的方法,其特征在于,所述将所述计算区域划分为多个子计算区域,包括:
对所述计算区域进行第一次区域划分,获得多个第一子计算区域,每个所述第一子计算区域的网格点数量一致;
同时对至少两个所述第一子计算区域进行第二次区域划分,获得多个第二子计算区域,每个所述第二子计算区域的网格点数量一致。
8.根据权利要求1所述的方法,其特征在于,所述获取所述血管的特征数据,包括:获取所述血管的影像数据和生理数据;
所述根据所述特征数据构建所述血管的三维模型,包括:根据所述影像数据构建所述血管的三维模型;
所述方法进一步包括:根据所述生理数据确定所述计算区域的边界条件;
所述基于所述网格,同时对至少两个所述子计算区域进行计算,包括:基于所述网格,根据所述边界条件同时对至少两个所述子计算区域进行计算。
9.一种血流模拟装置,包括存储器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序,其特征在于,所述处理器执行所述计算机程序时实现权利要求1-8中任一项所述方法的步骤。
10.一种计算机可读存储介质,其上存储有计算机程序,其特征在于,所述计算机程序被处理器执行时实现权利要求1-8中任一项所述方法的步骤。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810525629.6A CN109102568A (zh) | 2018-05-28 | 2018-05-28 | 基于区域分解的血管血流模拟方法及相关装置 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810525629.6A CN109102568A (zh) | 2018-05-28 | 2018-05-28 | 基于区域分解的血管血流模拟方法及相关装置 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN109102568A true CN109102568A (zh) | 2018-12-28 |
Family
ID=64796543
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201810525629.6A Pending CN109102568A (zh) | 2018-05-28 | 2018-05-28 | 基于区域分解的血管血流模拟方法及相关装置 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN109102568A (zh) |
Cited By (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109992858A (zh) * | 2019-03-20 | 2019-07-09 | 五邑大学 | 一种基于sph的出入流边界计算方法、装置和存储介质 |
CN110634572A (zh) * | 2019-09-24 | 2019-12-31 | 杭州阿特瑞科技有限公司 | 基于力学方程的血管血流模拟方法及相关装置 |
CN110675957A (zh) * | 2019-09-24 | 2020-01-10 | 杭州阿特瑞科技有限公司 | 血管血流模拟方法及相关装置 |
CN110795808A (zh) * | 2019-10-31 | 2020-02-14 | 北京理工大学 | 一种脑环境参数确定装置、方法及电子设备 |
CN111223186A (zh) * | 2020-01-15 | 2020-06-02 | 中南大学 | 三维随机孔结构模型的物理建模方法、系统和设备 |
CN113139260A (zh) * | 2020-01-17 | 2021-07-20 | 中国石油化工股份有限公司 | 一种用于提高钻井仿真计算速度的系统及方法 |
Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP2002172668A (ja) * | 2000-12-08 | 2002-06-18 | Toray Ind Inc | キャビティ内流体の解析方法および装置 |
US20100298719A1 (en) * | 2007-10-31 | 2010-11-25 | Samuel Alberg Kock | Method for calculating pressures in a fluid stream through a tube section, especially a blood vessel with atherosclerotic plaque |
US20130156158A1 (en) * | 2010-08-27 | 2013-06-20 | Konica Minolta Medical & Graphic, Inc. | Thoracic diagnosis assistance system and computer readable storage medium |
WO2013171039A1 (en) * | 2012-05-16 | 2013-11-21 | Feops Bvba | Pre -operative simulation of trans - catheter valve implantation |
CN104933225A (zh) * | 2015-05-25 | 2015-09-23 | 中国科学院过程工程研究所 | 实现计算流体力学大规模实时模拟的方法 |
CN105095534A (zh) * | 2014-04-23 | 2015-11-25 | 北京冠生云医疗技术有限公司 | 对血管中血流进行模拟的方法和系统 |
CN106780477A (zh) * | 2016-12-30 | 2017-05-31 | 上海联影医疗科技有限公司 | 一种血流分析方法和系统 |
-
2018
- 2018-05-28 CN CN201810525629.6A patent/CN109102568A/zh active Pending
Patent Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP2002172668A (ja) * | 2000-12-08 | 2002-06-18 | Toray Ind Inc | キャビティ内流体の解析方法および装置 |
US20100298719A1 (en) * | 2007-10-31 | 2010-11-25 | Samuel Alberg Kock | Method for calculating pressures in a fluid stream through a tube section, especially a blood vessel with atherosclerotic plaque |
US20130156158A1 (en) * | 2010-08-27 | 2013-06-20 | Konica Minolta Medical & Graphic, Inc. | Thoracic diagnosis assistance system and computer readable storage medium |
WO2013171039A1 (en) * | 2012-05-16 | 2013-11-21 | Feops Bvba | Pre -operative simulation of trans - catheter valve implantation |
CN105095534A (zh) * | 2014-04-23 | 2015-11-25 | 北京冠生云医疗技术有限公司 | 对血管中血流进行模拟的方法和系统 |
CN104933225A (zh) * | 2015-05-25 | 2015-09-23 | 中国科学院过程工程研究所 | 实现计算流体力学大规模实时模拟的方法 |
CN106780477A (zh) * | 2016-12-30 | 2017-05-31 | 上海联影医疗科技有限公司 | 一种血流分析方法和系统 |
Non-Patent Citations (14)
Title |
---|
卢培培;许学军;: "高波数波动问题的多水平方法", 计算数学, no. 02 * |
张玮, 王元, 徐忠: "多重网格技术在SIMPLE内外迭代中的应用", 西安交通大学学报, no. 07 * |
朱文婕: "图像分割技术在医学图像分割中的应用", 《安徽科技学院学报》 * |
朱文婕: "图像分割技术在医学图像分割中的应用", 《安徽科技学院学报》, no. 03, 15 May 2011 (2011-05-15) * |
李薇等: "基于计算机仿真的血管搭桥术出口端血流分析", 《中国西部科技》 * |
李薇等: "基于计算机仿真的血管搭桥术出口端血流分析", 《中国西部科技》, no. 36, 25 December 2008 (2008-12-25) * |
杨媛媛: "网格"计算"流体力学", 《中国计算机用户》 * |
杨媛媛: "网格"计算"流体力学", 《中国计算机用户》, 8 May 2006 (2006-05-08) * |
王娜等: "基于改进SPH的皮肤表面血流模拟算法", 《山东大学学报(工学版)》 * |
王娜等: "基于改进SPH的皮肤表面血流模拟算法", 《山东大学学报(工学版)》, no. 01, 7 January 2016 (2016-01-07) * |
肖曼玉等: "基于区域分解的并行SIMPLER算法研究", 《应用基础与工程科学学报》 * |
肖曼玉等: "基于区域分解的并行SIMPLER算法研究", 《应用基础与工程科学学报》, no. 03, 30 September 2006 (2006-09-30) * |
衣同训, 姜勇, 索沂生: "多重网格法求解轴流压气机内部流场", 应用力学学报, no. 03, pages 105 - 111 * |
赵洪雷;王松涛;韩万金;冯国泰;: "多级涡轮三维气动优化设计的可行性分析与实现", 热能动力工程, no. 01 * |
Cited By (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109992858A (zh) * | 2019-03-20 | 2019-07-09 | 五邑大学 | 一种基于sph的出入流边界计算方法、装置和存储介质 |
CN110634572A (zh) * | 2019-09-24 | 2019-12-31 | 杭州阿特瑞科技有限公司 | 基于力学方程的血管血流模拟方法及相关装置 |
CN110675957A (zh) * | 2019-09-24 | 2020-01-10 | 杭州阿特瑞科技有限公司 | 血管血流模拟方法及相关装置 |
CN110795808A (zh) * | 2019-10-31 | 2020-02-14 | 北京理工大学 | 一种脑环境参数确定装置、方法及电子设备 |
CN111223186A (zh) * | 2020-01-15 | 2020-06-02 | 中南大学 | 三维随机孔结构模型的物理建模方法、系统和设备 |
CN111223186B (zh) * | 2020-01-15 | 2022-03-25 | 中南大学 | 三维随机孔结构模型的物理建模方法、系统和设备 |
CN113139260A (zh) * | 2020-01-17 | 2021-07-20 | 中国石油化工股份有限公司 | 一种用于提高钻井仿真计算速度的系统及方法 |
CN113139260B (zh) * | 2020-01-17 | 2024-02-09 | 中国石油化工股份有限公司 | 一种用于提高钻井仿真计算速度的系统及方法 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN109064559A (zh) | 基于力学方程的血管血流模拟方法及相关装置 | |
CN109102568A (zh) | 基于区域分解的血管血流模拟方法及相关装置 | |
Hennigh et al. | NVIDIA SimNet™: An AI-accelerated multi-physics simulation framework | |
Lacasta et al. | An optimized GPU implementation of a 2D free surface simulation model on unstructured meshes | |
Marsden et al. | A computational framework for derivative-free optimization of cardiovascular geometries | |
Li et al. | Adjoint sensitivity analysis for time-dependent partial differential equations with adaptive mesh refinement | |
Mayer et al. | 3D fluid–structure-contact interaction based on a combined XFEM FSI and dual mortar contact approach | |
Krafczyk et al. | Two-dimensional simulation of fluid–structure interaction using lattice-Boltzmann methods | |
Grinberg et al. | A new computational paradigm in multiscale simulations: Application to brain blood flow | |
CN110634572B (zh) | 基于力学方程的血管血流模拟方法及相关装置 | |
Biancolini et al. | Fast interactive CFD evaluation of hemodynamics assisted by RBF mesh morphing and reduced order models: The case of aTAA modelling | |
Safaei et al. | Roadmap for cardiovascular circulation model | |
Xu et al. | PyCAC: The concurrent atomistic-continuum simulation environment | |
Guzzetti et al. | Propagating uncertainties in large-scale hemodynamics models via network uncertainty quantification and reduced-order modeling | |
Du et al. | Deep learning-based surrogate model for three-dimensional patient-specific computational fluid dynamics | |
Nickerson et al. | Using CellML with OpenCMISS to simulate multi-scale physiology | |
Maas et al. | A plugin framework for extending the simulation capabilities of FEBio | |
Goswami et al. | Neural operator learning of heterogeneous mechanobiological insults contributing to aortic aneurysms | |
JP7221062B2 (ja) | 流体解析システム、流体解析方法、および流体解析プログラム | |
Dzwinel et al. | A concept of a prognostic system for personalized anti-tumor therapy based on supermodeling | |
Letov et al. | A geometric modelling framework to support the design of heterogeneous lattice structures with non-linearly varying geometry | |
Jenabidehkordi et al. | An open source peridynamics code for dynamic fracture in homogeneous and heterogeneous materials | |
Tam | Principal stress line computation for discrete topology design | |
Tai et al. | Numerical simulation of 3D fluid–structure interaction flow using an immersed object method with overlapping grids | |
Black et al. | Deep neural networks for parameterized homogenization in concurrent multiscale structural optimization |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination |