CN108170895A - 纹线可视化设备和方法 - Google Patents

纹线可视化设备和方法 Download PDF

Info

Publication number
CN108170895A
CN108170895A CN201711261432.8A CN201711261432A CN108170895A CN 108170895 A CN108170895 A CN 108170895A CN 201711261432 A CN201711261432 A CN 201711261432A CN 108170895 A CN108170895 A CN 108170895A
Authority
CN
China
Prior art keywords
streakline
point
time
fluid
information
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.)
Granted
Application number
CN201711261432.8A
Other languages
English (en)
Other versions
CN108170895B (zh
Inventor
栉部大介
渡边正宏
久田俊明
杉浦清了
鹫尾巧
冈田纯
冈田纯一
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Fujitsu Ltd
University of Tokyo NUC
Original Assignee
Fujitsu Ltd
University of Tokyo NUC
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Priority claimed from JP2016236734A external-priority patent/JP6741252B2/ja
Priority claimed from JP2017075671A external-priority patent/JP6760603B2/ja
Application filed by Fujitsu Ltd, University of Tokyo NUC filed Critical Fujitsu Ltd
Publication of CN108170895A publication Critical patent/CN108170895A/zh
Application granted granted Critical
Publication of CN108170895B publication Critical patent/CN108170895B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/02—Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
    • A61B6/03—Computed tomography [CT]
    • A61B6/032—Transmission computed tomography [CT]
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/50—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications
    • A61B6/504—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications for diagnosis of blood vessels, e.g. by angiography
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/50—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications
    • A61B6/507—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications for determination of haemodynamic parameters, e.g. perfusion CT
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5211—Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data
    • A61B6/5217—Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data extracting a diagnostic or physiological parameter from medical diagnostic data
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0012—Biomedical image inspection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/28—Fuselage, exterior or interior
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/06—Power analysis or power optimisation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G06T2207/30101—Blood vessel; Artery; Vein; Vascular
    • G06T2207/30104—Vascular flow; Blood flow; Perfusion
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0012—Biomedical image inspection
    • G06T7/0014—Biomedical image inspection using an image reference approach
    • G06T7/0016—Biomedical image inspection using an image reference approach involving temporal comparison
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Medical Informatics (AREA)
  • Theoretical Computer Science (AREA)
  • General Health & Medical Sciences (AREA)
  • Radiology & Medical Imaging (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Animal Behavior & Ethology (AREA)
  • Public Health (AREA)
  • Optics & Photonics (AREA)
  • Pathology (AREA)
  • Biophysics (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Veterinary Medicine (AREA)
  • General Physics & Mathematics (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Oral & Maxillofacial Surgery (AREA)
  • Dentistry (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Vascular Medicine (AREA)
  • Pulmonology (AREA)
  • Physiology (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Quality & Reliability (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种纹线可视化设备和方法。在纹线可视化设备中，处理单元通过使用表达式来计算多个网格点上的流体在多个第一时间点中的每一个处的速度的时间微分值，所述表达式包括用于校正由位置信息表示的多个网格点的加速运动引起的误差的校正值。处理单元基于多个网格点上的流体在多个第一时间点中的每一个处的速度的时间微分值和速度，计算随着分析时间的进行从粒子生成源依序输出的一系列多个粒子在多个第二时间点中的每一个处的位置，所述多个第二时间点具有短于第一时间间隔的第二时间间隔。处理单元生成关于指示一系列多个粒子的纹线的显示信息。

Description

纹线可视化设备和方法
技术领域
本文中讨论的实施方式涉及纹线可视化设备和方法。
背景技术
流体力学是力学中的学术领域之一，并且其描述流体的行为。流体力学已经被应用于各种工业领域，在所述工业领域中，不仅空气或水的流动而且诸如温度或浓度的物理量的转移都被作为问题来处理。例如，流体力学已经应用于风洞实验以评估机动车原型，并在实验结果的基础上已经对这些机动车的空气动力学特性进行了优化。然而，这些风洞实验是非常昂贵的。因此，替代风洞实验，已经通过使用计算流体力学进行了模拟风洞实验的计算机模拟(流体模拟)。
近来计算机性能的提高已经使流体模拟得到了迅速的进展。因此，流体模拟不仅被应用于飞行器、机动车、铁路车辆、船舶等的空气动力学特性的评估，而且还被应用于心脏、血管等的血流状态的分析。
当进行流体模拟时，分析结果被可视化，使得可以在视觉上容易地理解分析结果。对流体模拟的结果进行可视化的一种手段是显示纹线。纹线是通过连接已经通过空间中的特定点的流体粒子而形成的曲线。在风洞实验中，从预定位置喷出的烟雾的踪迹是纹线。即，通过在流体模拟中生成关于纹线的显示信息并显示所述纹线，可以将流体中的粒子的运动——如在风洞实验中的烟雾的踪迹中的粒子的运动——可视化，而无需执行任何风洞实验。
已经提出了与流体模拟相关的各种技术。例如，已经详细地提出了执行高速模拟并且快速且平滑地表示流体中的场景的技术。还提出了易于将结构-流体分析模拟的结果应用于对血管异常的诊断的技术。还提出了能够使不熟悉计算流体力学的使用者如医生进行适当的血流模拟的设备。此外，已经发表了与流体模拟相关的各种论文。例如参见以下文献：
日本特许专利公开第2003-6552号。
日本特许专利公开第2015-97759号。
国际公开小册子第WO2016/056642号。
Tino Weinkauf和Holger Theisel，“Streak Lines as Tangent Curves of aDerived Vector Field”，IEEE可视化与计算机图形学报(IEEE Transactions onVisualization and Computer Graphics)，2010年11月-12月，第16卷，第6期。
Erwin Fehlberg，“LOW-ORDER CLASSICAL RUNGE-KUTTA FORMULAS WITHSTEPSIZE CONTROL AND THEIR APPLICATION TO SOME HEAT TRANSFER PROBLEMS”，NASA技术报告(NASA TECHNICAL REPORT)，NASA TR R-315，1969年7月。
J.Donea，A.Huerta，J.-Ph.Ponthot和A.Rodríguez-Ferran，“ArbitraryLagrangian-Eulerian methods”，计算力学百科全书(Encyclopedia of ComputationalMechanics)，约翰威立股份有限公司(John Wiley&Sons Ltd.)，2004年11月，第413至437页。
Seiryo Sugiura，Takumi Washio，Asuka Hatano，Junichi Okada，HiroshiWatanabe，Toshiaki Hisada，“Multi-scale simulations of cardiacelectrophysiology and mechanics using the University of Tokyo heartsimulator”，生物物理学与分子生物学进展(Progress in Biophysics and MolecularBiology)，2012年10月-11月，第110卷，第380至389页。
William H.Press等，“Numerical Recipes in C：The Art of ScientificComputing”，剑桥大学出版社(Cambridge University Press)，1992年10月30日，第113-117页。
Joseph E.Flaherty，“Finite Element Analysis-Chapter 4Finite ElementApproximation”，2005年4月1日(2017年2月27日搜索的)，<URL:www.cs.rpi.edu/～flaherje/pdf/fea4.pdf>。
然而，这些常规的纹线分析技术——例如在分析机动车周围的空气流动的情况下——基于分析空间中的结构不会变形的假设。当结构变形时，难以精确地追踪纹线。
发明内容
在一个方面中，实施方式的目的是即使当结构变形时也实现对纹线的追踪。
根据一个方面，提供了一种纹线可视化设备，包括：存储装置，用于存储流体信息，流体信息包括位置信息，该位置信息指示多个网格点在具有第一时间间隔的多个第一时间点中的每个第一时间点处的位置，多个网格点随着流体模拟的分析时间的进行在分析空间中以加速运动移动，并且流体信息包括速度信息，该速度信息指示多个网格点上的流体在多个第一时间点中的每个第一时间点处的速度；以及处理装置，用于：基于流体信息并且通过使用表达式来计算多个网格点上的流体在多个第一时间点中的每个第一时间点处的速度的时间微分值，所述表达式包括用于校正由位置信息表示的多个网格点的加速运动引起的误差的校正值；基于多个网格点上的流体在多个第一时间点中的每个第一时间点处的速度的时间微分值和速度来计算随着分析时间的进行从粒子生成源依序输出的一系列多个粒子在多个第二时间点中的每个第二时间点处的位置，多个第二时间点具有短于第一时间间隔的第二时间间隔；以及生成关于指示一系列多个粒子的纹线的显示信息。
附图说明
图1示出根据第一实施方式的纹线可视化设备的示例；
图2示出根据第二实施方式的系统配置示例；
图3示出可视化设备的硬件配置示例；
图4示出纹线计算示例；
图5是示出可视化设备的功能的框图；
图6示出弹性体信息文件组的示例；
图7示出流体信息文件组的示例；
图8示出预分析文件的示例；
图9是示出纹线可视化处理的过程的示例的流程图；
图10是示出预分析处理的过程的示例的流程图；
图11示出当基于任意拉格朗日-欧拉(Arbitrary Lagrangian-Eulerian)法控制的网格点自身运动时出现的表观力的原理；
图12和图13是示出时间演化计算处理的过程的流程图；
图14示出龙格-库塔(Runge-Kutta)系数表的示例；
图15示出纹线的数据示例；
图16示出网格点的位置如何随时间不断变化；
图17是示出用于计算在时间点t＝tk时的场信息的处理的过程的示例的流程图；
图18是示出用于计算在时间点t＝tk时的弹性体场信息的处理的过程的示例的流程图；
图19示出四面体单元的示例；
图20是示出基于龙格-库塔法的时间演化处理的过程的示例的流程图；
图21示出能够作为计算结果而获得的位置的示例；
图22是指示状态变量的真值表；
图23是示出用于确定心脏内部的位置的处理的过程的示例的流程图；
图24是示出用于确定时间演化后位置是否落在流体内的处理的过程的示例的流程图；
图25是示出用于计数由移动向量形成的线与心肌的交点的处理的过程的示例的流程图；
图26是示出用于设置预测球体的半径的处理的过程的示例的流程图；
图27示出通过时间分割而实现计算量的减少的构思；
图28示出根据时间分割数目的总计算时间的变化的示例；
图29示出单元的边长的分布的示例；
图30示出基于预测的球体半径的计算成本的变化；
图31示出移动距离的概率分布；
图32示出在假设概率分布时的计算效率曲线；
图33示出纹线的显示的示例；
图34示出当时间步长改变时的轨迹的精确度的变化；
图35示出已经向其添加了对ALE网格点的运动的修正的速度场的示例；以及
图36示出防止纹线嵌入到心肌中的有利效果，该效果是作为计算精度提高的结果而获得的。
具体实施方式
下面将参照附图来描述实施方式，其中，类似的附图标记贯穿全文指代类似的要素。两个或更多个实施方式可以彼此组合而不会产生不一致。[第一实施方式]
首先，将描述第一实施方式。第一实施方式提供了即使当结构变形时也能够追踪纹线的纹线显示设备。
当结构变形时，难以精确地追踪纹线。其中一个原因是，当结构变形时，网格点在计算上处于关于移动边界问题的运动中，并且因此产生表观力。
科里奥利(Coriolis)力是表观力的简单示例。如果在地球上保持不动的观察者发射炮弹，那么在观察者看来就好像单靠离心力所不能解释的力正在弯曲炮弹的轨迹。这是因为，当从空间看时，由于观察者他或她自身也在与地球一起运动。即，观察者具有加速度，并且该加速度产生表观力。不是内在力而是看起来好像由于观察者的运动(坐标系和处理该坐标系的方式)而作用的这样的力被称为“表观力”。
即，当结构变形时，结构变形的影响需要被反映在流体运动的计算中。在该操作中，流体的运动与限定结构形状的单个网格点的移动相关联地被计算。如果网格点随加速运动而移动，则表观力存在于基于网格点的流体运动的计算中。除非通过反映该表观力来计算流体的运动，否则纹线被不精确地计算。例如，如果表示流体速度的速度场的计算精确度低，则呈现出不可能的行为，例如纹线嵌入到心肌中，并且计算失败。
因此，在第一实施方式中，通过校正由单个网格点的加速运动而引起的表观力，流体的速度场被精确地计算。因此，精确地得出了心肌附近的纹线的行为。因此，防止了纹线嵌入到心肌中的现象。因此，由于提高了计算稳定性，并且纹线没有嵌入到心肌中，因此避免了计算错误。
图1示出根据第一实施方式的纹线可视化设备的示例。该纹线可视化设备10包括存储单元11、处理单元12和显示单元13。例如，存储单元11是纹线可视化设备10的存储器或存储装置。例如，处理单元12是纹线可视化设备10的处理器或算术电路。例如，显示单元13是纹线可视化设备10的图形电路。
例如，存储单元11包括结构信息11a和流体信息11b。结构信息11a指示流体模拟分析空间中的结构的形状随时间的变化。流体信息11b包括关于网格点V1至V4的位置信息和网格点V1至V4上的流体速度信息。例如，位置信息指示多个网格点V1至V4在具有第一时间间隔的多个第一时间点(T0、T1、…)中的每个第一时间点处的位置。此外，速度信息指示多个网格点V1至V4上的流体在多个第一时间点(T0、T1、…)中的每个第一时间点处的速度。随着流体模拟分析时间的进行，多个网格点在分析空间中以加速运动的方式移动。虽然图1示出了四个网格点V1至V4，但是分析空间包括未示出的许多其他网格点。此外，在图1中，单个网格点处的流体速度由箭头指示。
处理单元12基于流体信息11b来计算纹线5a和5b。例如，通过使用包括用于校正由位置信息表示的多个网格点的加速运动引起的误差的校正值的表达式，处理单元12计算多个网格点V1至V4上的流体在多个第一时间点中的每个第一时间点处的速度的时间微分值。具体地，例如，处理单元12基于位置信息来计算多个网格点V1至V4中的每个网格点上的在多个第一时间点(T0、T1、…)中的每个第一时间点处的速度和加速度。接下来，通过使用包括作为变量的所计算的速度和加速度的表达式，处理单元12计算多个网格点V1至V4上的流体在多个第一时间点(T0、T1、…)中的每个第一时间点处的速度的时间微分值。以下表达式(6)和(7)是用于计算速度和加速度的表达式的示例。此外，以下表达式(4)和(5)是用于计算速度的时间微分值的表达式的示例，所述表达式包括作为变量的速度和加速度。
接下来，处理单元12定义随着分析时间的进行从粒子生成源4顺序输出的多个粒子。接着，基于多个网格点V1至V4上的流体在多个第一时间点(T0、T1、…)中的每个第一时间点处的速度和速度的时间微分值，处理单元12计算多个粒子在多个第二时间点(t0、t1、…)中的每个第二时间点处的位置。多个第二时间点(t0、t1、…)具有比第一时间间隔短的第二时间间隔。例如，处理单元12生成插值曲线，插值曲线中的每一个平滑地连接指示多个网格点V1至V4中的相应的一个网格点处的在多个第一时间点(T1、T2、…)中的每个第一时间点处的速度的点，并且插值曲线中的每一个表示流体速度随时间的变化。当插值曲线通过单个第一时间点处的点时，插值曲线具有基于第一时间点处的流体速度的时间微分值的斜率。接下来，基于多个网格点V1至V4的各个插值曲线，处理单元12计算多个粒子在多个第二时间点中的每个第二时间点处的位置。将利用下面的表达式(20)来描述这些插值曲线的示例。
此外，处理单元12基于结构信息11a来再现结构2的形状随时间的变化。处理单元12生成关于各自指示从粒子生成源4输出的一系列多个粒子的纹线5a和5b的显示信息。
显示单元13基于由处理单元12生成的显示信息来在显示屏1上显示纹线5a和5b。例如，显示单元13在流体区域3a或3b内叠加纹线5a或5b，所述流体区域3a或3b包括在显示屏1上的结构2的图像上的流体。在图1中，纹线5a是第二时间点“t0”处的纹线，而纹线5b是第二时间点“t1”处的纹线。
以这种方式，纹线可视化设备10通过反映由各个网格点的加速运动引起的表观力来精确地计算多个粒子的位置，并且通过生成关于精确纹线的显示信息来显示精确纹线。
此外，处理单元12可以被配置成确定纹线5a或5b是否进入结构，并且仅当纹线5a或5b未进入结构时才显示纹线5a或5b。在这种情况下，为了高效处理，处理单元12可以被配置成当基于第一分析时间点处的第一纹线来计算第二分析时间处的第二纹线时执行以下处理。
例如，处理单元12将包括分析空间中的第一纹线上的第一位置处的离散点的部分区域设置为离散点的分析目标区域。例如，分析目标区域是球形区域。接下来，处理单元12基于分析目标区域中的流体的速度来计算指示离散点上的粒子在第二分析时间点处的目的地的第二位置，所述速度由流体信息指示。接下来，处理单元12基于关于分析目标区域中的结构的信息来确定由分析目标区域的结构在第二分析时间点处占据的区域，所述信息由结构信息指示。接下来，处理单元12基于第一位置和第二位置来确定第二纹线进入或者未进入占据的区域。当确定第二纹线未进入占据的区域时，处理单元12生成关于第二纹线通过第二位置的显示信息并且显示第二纹线。
以这种方式，通过将分析范围限制于分析目标区域并且确定纹线是否已经进入结构的占据的区域，吞吐量减少，并且有效的处理被执行。
[第二实施方式]
接下来，将描述第二实施方式。第二实施方式提供了一种能够随着心脏的运动将心脏中的血流的纹线可视化的可视化设备。
例如，使用计算流体分析可以模拟流体的行为，甚至是技术上或伦理上难以测量的系统中的流体，例如心脏中的血流传输。因此，计算流体分析被用于讨论其中心脏中的血流传输存在故障的先天性心脏病等的治疗。即，计算流体分析是重要的技术。通过使用可视化设备来将这样的计算流体分析的结果可视化，卫生保健专业人员如医生能够容易地理解分析结果并制定治疗计划。
图2示出根据第二实施方式的系统配置示例。可视化设备100经由网络20连接到心脏模拟器200。心脏模拟器200是执行心肌运动和冠状循环的模拟的计算机。可视化设备100从心脏模拟器200获取模拟结果。接下来，可视化设备100基于模拟结果来生成关于纹线的显示信息并且显示所生成的纹线。例如，模拟结果包括关于指示以下项的三维(3D)模型的信息：每个时间点的心脏形状、血管中的血液的速度以及心肌或血液的物理性质值。
图3示出可视化设备100的硬件配置示例。可视化设备100由处理器101全面控制。处理器101经由总线109连接至存储器102和多个外围装置。处理器101可以是多处理器。处理器101是算术处理设备，例如中央处理单元(CPU)、微处理单元(MPU)或数字信号处理器(DSP)。通过使处理器101执行程序来实现的功能中的至少一部分可以通过使用诸如专用集成电路(ASIC)或可编程逻辑器件(PLD)的电子电路来实现。
存储器102用作可视化设备100的主存储装置。存储器102临时地存储由处理器101执行的操作系统(OS)程序或应用程序中的至少一部分。此外，存储器102存储处理器101执行的处理所需的各种类型的数据。例如，诸如随机存取存储器(RAM)的易失性半导体存储装置被用作存储器102。
连接至总线109的外围装置的示例包括存储装置103、图形处理装置104、输入接口105、光驱装置106、装置连接接口107以及网络接口108。
存储装置103电或磁地写入和读取其存储介质上的数据。存储装置103被用作可视化设备100的辅助存储装置。存储装置103存储OS程序、应用程序以及各种类型的数据。例如，可以将硬盘驱动器(HDD)或固态驱动器(SSD)用作为存储装置103。
图形处理装置104连接至监视器21。图形处理装置104根据来自处理器101的指令在监视器21的屏幕上显示图像。监视器21的示例包括阴极射线管(CRT)显示装置和液晶显示(LCD)装置。
输入接口105连接至键盘22和鼠标23。输入接口105将从键盘22或鼠标23发送的信号发送至处理器101。鼠标23是定点装置。也可以使用不同的定点装置，例如触控面板、平板计算机、触摸板或追踪球。
光驱装置106通过使用激光等来读取存储在光盘24上的数据。光盘24是存储通过光反射读取的数据的便携式存储介质。光盘24的示例包括数字通用光盘(DVD)、DVD-RAM、光盘只读存储器(CD-ROM)以及可记录(R)/可重写(RW)CD。
装置连接接口107是用于将外围装置连接至可视化设备100的通信接口。例如，存储器装置25或存储器读写器26可以连接至装置连接接口107。存储器装置25是能够与装置连接接口107通信的存储介质。存储器读写器26能够在存储器卡27上读写数据。存储器卡27是卡式存储介质。
网络接口108连接至网络20。网络接口108经由网络20与其他计算机或通信装置交换数据。
根据第二实施方式的处理功能可以由以上硬件配置来实现。第一实施方式中描述的设备还可以由与图3示出的可视化设备100的硬件配置等同的硬件配置来实现。
例如，可视化设备100通过执行存储在计算机可读存储介质中的程序来实现根据第二实施方式的处理功能。保存由可视化设备100执行的处理内容的程序可以存储在各种类型的存储介质中的任何一个中。例如，由可视化设备100执行的程序可以存储在存储装置103中。处理器101将存储装置103中的程序的至少一部分加载到存储器102上，并且执行所加载的程序。由可视化设备100执行的程序可以存储在便携式存储介质中，例如光盘24、存储器装置25或存储器卡27中。例如，在存储在便携式存储介质中的程序由处理器101安装在存储装置103中之后，程序由处理器101执行。处理器101可以直接从便携式存储介质读取程序并执行所读取的程序。
接下来，将描述纹线。
图4示出纹线计算示例。可视化设备100在分析目标空间中限定粒子生成源30。当开始分析时，从粒子生成源30不断发出粒子群33。当流场不随时间变化时，粒子群33形成固定的曲线(流线)。然而，当流场随时间变化时，由粒子群33形成的曲线随时改变。纹线31和纹线32是当流场随时间变化时形成的这样的曲线。在图4中，纹线31表示时间点t0处的一系列粒子群33，并且纹线32表示时间点t1处的一系列粒子群33。由于存在障碍物35，因此这些纹线31和纹线32变得非常弯曲。
这些纹线31和纹线32对于将粒子群33如何在时变流场中传输进行可视化是有用的。例如，将描述障碍物35是机动车的情况。为了使开发的机动车的空气阻力可视化，粒子生成源30被布置在机动车的前方，并且从被布置在布置有粒子生成源30的地方的风扇等朝向机动车供应空气。此外，在可视化设备100中的流体模拟中，连续地发射粒子群33，并且粒子群33的轨迹被测量作为纹线31和纹线32。纹线31和纹线32直接描述流体的传输并将其可视化。因此，纹线适用于各种领域。
对纹线的计算和可视化已经进行了大量的研究。此外，针对湍流、不稳定流动等也已经进行了大量研究。然而，在诸如心脏的弹性体(其被流体视为壁面)经历大的变形的模拟中，在纹线的可视化上先前没有进行太多的研究。由于心脏周期性搏动并且反复膨胀和收缩，因此它们是经历大变形的系统的典型示例。此外，由于该周期性运动在心脏的泵送作用中起重要作用，因此在考虑心脏病的治疗时，对弹性体周期性地经历大变形的系统中的血流的传输进行评估是重要的。
在生物模拟领域，已经在计算机上模拟了心脏的行为。通过在计算机上的模拟，评估由操作获得的治疗的有效性而无需实际执行操作。因此，使用生物模拟可以使医生在实际执行操作之前能够考虑最佳的治疗计划。特别地，心脏模拟针对具有复杂3D结构的心脏，并且心脏的行为动态地改变。如果将表示心脏中的血流的传输的纹线与心脏的行为协同地进行可视化，则医生就可以容易地在视觉上了解心脏的状态。在视觉上容易地显示心脏的状态对于防止判断中的错误是有效的。
以下几点是对心脏中的血流的纹线进行可视化要克服的障碍。
1.当心肌(弹性体)变形较大时，难以精确地追踪心肌周围的迹线和纹线的行为。
2.在迹线的情况下，只需要对单个点执行计算。但是，在纹线的情况下，需要对形成线的全部N个点执行计算，导致明显大的计算量。
3.有限元模型的一些低质量网格增加了总体计算量。
因此，通过使用以下功能，根据第二实施方式的可视化设备100以可行的计算量将精确的纹线可视化。
1-1：可视化设备100精确地确定纹线上的单个点是否已经进入运动区域外部的心肌或已经落在模拟目标系统的外部。
1-2：当可视化设备100确定纹线上的点已经落在运动区域的外部时，可视化设备100调整作为用于纹线的微分方程中的参数的时间步长，以防止点落在移动区域的外部。
1-3：为了估计关于场在任何时间点处的信息，可视化设备100通过使用插值方法对场进行插值。
1-4：为了提高纹线的运动的精确度，可视化设备100通过使用其中反映表观力的表达式来计算纹线。
2-1：由于1-3中的功能的应用增加了计算量，因此可视化设备100计算纹线上的点移动的最大距离，并且仅计算关于具有等于最大距离的半径的预测球体内部的场的信息。以这种方式，不管数据容量如何，可视化设备100均保持一定的计算量。
2-2：通过分割时间步长，可视化设备100减小预测球体的半径，需要比不使用预测球体时所需的计算量更小的计算量，并且同时提高精确度。
3-1：由于大多数计算是在高质量网格上执行的，因此可视化设备100通过假设所有网格都是高质量网格来执行推测计算。在计算失败的情况下，可视化设备100执行精确的计算。以这种方式，计算量减小。在该推测计算中，通过允许纹线上的点的目的地落在预测球体的外部的可能性，可视化设备100减小预测球体的半径。计算失败的情况是纹线上的点的目的地不存在于预测球体内的情况。
3-2：由于可视化设备100在3-1中执行推测计算，因此可视化设备100准备概率模型并且确定实现最小计算量的参数集，最小计算量包括计算失败时所需的处罚(penalty)。
4-1：为了有效地执行其中反映表观力的纹线计算，预先计算基于表观力的系数，并且通过参照每个时间步长的计算结果来计算纹线。
通过在可视化设备100上实现这些功能来获得以下有利效果。
1.即使当心肌(弹性体)变形较大时，可视化设备100也能够在将心肌(弹性体)的运动考虑在内时计算纹线。
2.可视化设备100能够通过使用预测球体和其中反映表观力的的表达式来快速且精确地计算在纹线上的单个点。
3.可视化设备100能够通过使用概率模型来设置使计算成本最小化的预测球体的半径。
在下文中，将详细描述可视化设备100的功能。
图5是示出可视化设备100的功能的框图。可视化设备100包括作为信息存储功能的模拟结果存储单元110和预分析结果存储单元120。模拟结果存储单元110存储从心脏模拟器200获取的模拟结果。例如，当心脏模拟器200执行计算流体动力学模拟时，关于动态变化的弹性体和流场在L个时间点t0、t1、...、tL(L是1或更大的整数)处的模拟结果被存储在文件中。例如，关于心肌的信息和关于血流的信息被作为分离的文件存储在模拟结果存储单元110中。在图5的示例中，每个时间点的关于心肌的信息被存储为弹性体信息文件组111，并且每个时间点的关于血流的信息被存储为流体信息文件组112。
预分析结果存储单元120存储通过表观力的预分析计算的基于表观力的系数。在下文中，通过预分析计算的一组系数将被称为预分析结果Ω。例如，预分析结果存储单元120存储包括预分析结果Ω的预分析文件121。
通过分析这些存储在模拟结果存储单元110中的模拟结果，可视化设备100计算描述关于血流的传输的信息的纹线。作为模拟结果输出的时间间隔Δti＝ti+1-ti不需要与心脏模拟器200求解微分方程时所使用的时间间隔相匹配。为了减少信息量，通常只输出一些模拟结果。因此，为了精确地获得纹线，可视化设备100通过使用多个时间点处的输出文件来使用插值方法等，并且估计目标时间点处的各种物理量。
接下来，将描述可视化设备100的处理功能。可视化设备100包括信息读取单元130、预分析单元140、纹线计算单元150、显示处理单元160以及分析单元170。
信息读取单元130从模拟结果存储单元110读取指示流体分析结果的文件。预分析单元140通过使用由信息读取单元130读取的信息来计算用于纹线计算的系数。预分析单元140将包括预分析结果Ω的预分析文件121存储在预分析结果存储单元120中。纹线计算单元150通过使用由信息读取单元130读取的信息和由预分析单元140计算的系数来计算纹线。显示处理单元160将获得的结果可视化。
分析单元170是通常由信息读取单元130、预分析单元140、纹线计算单元150以及显示处理单元160使用的一组功能。当执行特定分析处理时，信息读取单元130、预分析单元140、纹线计算单元150以及显示处理单元160请求分析单元170执行处理并获得结果。
分析单元170包括系数计算单元171、流体信息分析单元172、弹性体信息分析单元173、预测球体半径计算单元174、内部预测球体信息计算单元175、移动边界碰撞确定单元176以及显示格式确定单元177。系数计算单元171计算包括在其中反映表观力的插值多项式中的系数。流体信息分析单元172分析流体的速度场、离散点的位置以及边界表面。弹性体信息分析单元173分析非流体的诸如心肌的弹性体的离散点的位置和边界表面。预测球体半径计算单元174设置用于在计算纹线时提高计算速度和计算精确度的预测球体的半径。内部预测球体信息计算单元175例如计算预测球体内部的心肌位置和速度场。移动边界碰撞确定单元176确定纹线上的点是否已经进入心肌——作为计算错误的结果。显示格式确定单元177确定如何显示获得的纹线。
例如，图5所示的单个单元的功能可以通过使计算机执行与相应单元对应的程序模块来实现。
接下来，将详细描述作为模拟结果而获得的信息。
图6示出弹性体信息文件组111的示例。弹性体信息文件组111是每个模拟时间点处的一组弹性体信息文件111a、111b等。弹性体信息文件111a、111b等中的每一个被给予文件名如“stru(X).inp”。在这种情况下，单个文件名中的“X”表示编号，并且这些编号根据模拟时间点的时间顺序以升序给出。
弹性体信息文件111a、111b等包括指示心脏在各个时间点处的形状的心肌数据。心肌数据包括沿着单个网格(布置在3D空间中的顶点)的x轴、y轴和z轴的坐标值，单个网格ID指示包括在心肌中的四面体单元(TETRA)的四个顶点和施加于单个单元的力。
图7示出流体信息文件组112的示例。流体信息文件组112是每个模拟时间点处的一组流体信息文件112a、112b等。例如，流体信息文件112a、112b等中的每一个被给予文件名如“flui(Y).inp”。在这种情况下，单个文件名中的“Y”表示编号，并且这些编号根据模拟时间点的时间顺序以升序给出。
流体信息文件112a、112b等包括指示各个时间点处的血流的血流数据。血流数据包括沿着单个网格(布置在3D空间中的顶点)的x轴、y轴和z轴的坐标值，单个网格ID指示包括在血管中的四面体单元(TETRA)的四个顶点，并且单个速度场向量指示血液在单个网格上流动的方向和速度。
图8示出预分析文件的示例。在预分析文件121中，记录网格编号、时间点、位置坐标、速度场(仅血流)以及插值多项式系数，使得插值多项式(下面的表达式(2))在单个时间点和单个网格点处再现。例如，预分析文件121包括输出时间点指数表121a、心肌侧插值多项式系数表121b以及流体侧插值多项式系数表121c。
在输出时间点指数表121a中，流体分析模拟上的时间点(Time)与时间点指数(Time Index)相关联。
在心肌侧插值多项式系数表121b中，将心肌侧插值多项式中使用的轴向位置(Position)和系数(系数a至d)设置成与网格ID(GRID ID)、时间点指数(Time Index)以及轴指数(Direction Index)的集合相关联。
在流体侧插值多项式系数表121c中，将流体侧插值多项式中使用的轴向位置(Position)、轴向速度(Velocity)以及系数(系数a至d)设置成与网格ID(GRID ID)、时间点指数(Time Index)以及轴指数(Direction Index)的集合相关联。
流体侧插值多项式系数表121c可以包括单个处的流体的速度场的时间微分的计算结果。
可视化设备100基于图6和图7示出的模拟结果和图8示出的预分析结果Ω来计算纹线并将纹线可视化。在下文中，将详细描述纹线可视化处理。
图9是示出纹线可视化处理的过程的示例的流程图。在下文中，将逐步地描述图9所示的处理。
[步骤S101]纹线计算单元150确定预分析文件121是否被存储在预分析结果存储单元120中。如果纹线计算单元150确定预分析文件121被存储，则处理进行至步骤S102。如果否，则处理进行至步骤S103。
[步骤S102]纹线计算单元150从预分析结果存储单元120读取预分析文件121。接下来，处理进行至步骤S105。
[步骤S103]预分析单元140执行用于计算插值多项式系数的预分析处理。将参照图10在下面详细描述预分析处理。
[步骤S104]预分析单元140生成包括在预分析处理中计算出的系数的预分析文件121，并且将预分析文件121存储在预分析结果存储单元120中。
接下来，在步骤S105至S108中，纹线计算单元150执行单个纹线的初始设置。
[步骤S105]纹线计算单元150设置要被计算的纹线的数目M(M是1或更大的整数)。例如，纹线计算单元150将由用户输入的值设置为纹线的数目M。
[步骤S106]纹线计算单元150设置纹线计算的数目N(N是1或更大的整数)。在下文中，纹线计算的数目N将被称为“时间演化数目”。例如，纹线计算单元150将由用户输入的值设置为时间演化数目N。
由于纹线随时间变化，因此通过设置时间演化数目N来确定一系列时间点t0、t1、...、和tN，在这些时间点处纹线的结果被输出。当指定的时间演化数目N比文件的数目L大(L是1或更大的整数)时，纹线计算单元150可以将第(L+1)文件作为第二心动周期中的跳动，并且将时间点t0处的文件用于第(L+1)文件。例如，所述一系列时间点中的时间间隔被设置成0.01秒。然而，可替选地，所述一系列时间点可以具有不规则的时间间隔。
[步骤S107]纹线计算单元150设置纹线的粒子生成源的坐标。例如，纹线计算单元150将分析空间中由用户指定的点设置为粒子生成源的坐标。例如，用户在参考心肌信息和血流信息时指定空间中的点。纹线计算单元150读取指定点的坐标作为坐标向量X0。当纹线的数目是1时，纹线的粒子生成源被设置为具有坐标向量X0。当纹线的数目是多个时，纹线计算单元150随机地将纹线lj的粒子生成源设置在以下球体中，该球体具有坐标向量X0为其中心且具有半径r(r是正实数)。纹线的粒子生成源被从血流中的坐标中选择。纹线计算单元150将设置的粒子生成源的坐标向量Xj设置为纹线的粒子生成源。
接下来，纹线计算单元150如下执行第j(j＝1、2、...、M)纹线lj的初始设置。
第j纹线lj由与时间演化数目N匹配的离散点形成。因此，纹线计算单元150生成指示包括在纹线lj中的离散点的点Pij(i＝0、1、2、...、N)。纹线计算单元150将单个离散点的初始值的坐标设置为粒子生成源的坐标。
当计算在时间点t＝ti处的纹线lj时，对单个点Pij(i＝0、1、2、...、i)进行时间演化计算，作为从相应的生成源发出的纹线粒子的位置。由于没有从任何生成源发出与点Pij(i＝i+1、...、N)相对应的纹线粒子，因此当计算时间点t＝ti处的纹线lj时，不对这些纹线粒子进行计算。此外，纹线计算单元150以值i的升序来计算离散点。由于从相应的粒子生成源发出，因此较早计算的离散点具有较长的时间。
[步骤S108]纹线计算单元150针对以下情况执行设置：纹线lj上的点已经落在大动脉等中，即落在目标系统中的流体边界的外部。纹线lj上的点Pij可能落在大动脉等中，即通过流体边界落在系统的外部。在这种情况下，由于在系统的外部没有定义流体速度场，因此点Pij在下一时间点处的计算不能被执行。因此，纹线计算单元150为每个点Pij设置区域确定标记作为单个离散点的参数。当点Pij已经落在目标区域内时，区域确定标记指示“T”。相比之下，当点Pij已经通过大动脉等中的流体的流动而漂移并且落在目标区域的外部时，区域确定标记指示“F”。由于在初始设置中流体包括所有点Pij，因此纹线计算单元150将每个离散点的区域确定标记设置为“T”。
[步骤S109]纹线计算单元150从指数i＝1起以升序对指数i(i＝1、2、...、和N-1)中的每一个重复一组步骤S110至步骤S114。
[步骤S110]纹线计算单元150从指数j＝1起以升序对指数j(j＝1、2、...、和M)中的每一个重复一组步骤S111至步骤S113。
[步骤S111]纹线计算单元150设置作为时间演化的开始的时间点并将该时间点存储在存储器102中。在第i次计算中，计算开始时间点被设置为t＝ti。时间演化结束时间点被设置为t＝ti+1。
[步骤S112]纹线计算单元150在由ti≤t≤ti+1定义的时间点之间执行时间演化计算。基于时间演化计算，对每个时间点t＝ti处从纹线lj的粒子生成源发出的所有点Pij(i＝0、1、2、...和i)进行时间演化，并且线上的所有点被时刻地更新。坐标值被获取作为在时间点[ti，ti+1]处的纹线lj上的各个点Pij的时间演化计算的结果，其在下一时间点t＝ti+1处被设置为坐标Pi+1,j。
[步骤S113]纹线计算单元150将所获取的计算结果存储在存储器中。基于计算结果，显示处理单元160将纹线lj可视化。此外，纹线计算单元150能够将纹线lj的坐标值输出至文件。
[步骤S114]每当纹线计算单元150执行一组步骤S111至步骤S113时，纹线计算单元150向指数j加1。在对指数j＝M执行步骤S111至S113之后，纹线计算单元150执行步骤S115。
[步骤S115]每当纹线计算单元150执行一组步骤S110至步骤S114时，纹线计算单元150向指数i加1。在对指数i＝N-1执行步骤S110至S114之后，纹线计算单元150结束纹线可视化处理。
接下来，将详细描述预分析处理(步骤S103)。预分析处理是用于生成系数的处理，所述系数用于估计在任何时间点t处的血流部分的速度场和心肌有关的结构信息。
首先，将描述作为前提的输出数据。输出数据被给出模拟时间点t＝ti。在下文中，关于心肌的结构信息将被表示为向量M(向量r，ti)。在心肌上，给出有限数目的离散点向量rk。这些信息项已经被输出至如图6所示的关于心肌的弹性体信息文件111a、111b等。离散点由各个GRID ID唯一标识，并且相应的坐标被存储。每个轴的单个坐标值是指示与GRIDID对应的离散点的单个离散点向量rk的元素。
在有限元方法中，通过使用各自具有这些离散点作为顶点的图形(例如，四面体图形)作为单元来执行计算。单元中的每个单元由图6中的弹性体信息文件111a、111b等中的TETRA ID唯一地标识，并且构成单元的离散点的GRID ID被存储。虽然在图6的示例中使用了四面体单元，但是单元可以具有与四面体图形不同的图形。心肌的形状由离散点的位置和无限元的排列来确定。
关于流体部分的结构信息也以与心肌相同的方式确定。关于流体部分的结构信息将被表示为向量B(向量r，ti)。流体部分也由有限数目的离散点和有限元形成，并且这些信息项已经被输出至如图7所示的流体信息文件112a、112b等。流体与心肌之间的差异在于，指示单个离散网格点(向量rk)上的流体速度场的速度场信息(速度向量v(向量rk，ti))也已经被输出为流体信息。预分析单元140通过使用以下项来执行预分析处理：所有时间点t0、t1、…、tL处的速度场信息、血流的结构以及心肌的结构。
图10是示出预分析处理的过程的示例的流程图。在下文中，将逐步描述图10所示的处理。
[步骤S121]预分析单元140从i＝1至i＝N-1依次对指数i中的每一个执行步骤S122。
[步骤S122]预分析单元140经由信息读取单元130从模拟结果存储单元110中读取时间点ti处的向量B(向量r，ti)、向量v(向量r，ti)以及向量M(向量r，ti)。
[步骤S123]每当预分析单元140执行步骤S122时，预分析单元140向指数i加1。在预分析单元140对指数i＝N-1执行步骤S122之后，处理进行至步骤S124。
[步骤S124]预分析单元140计算关于网格坐标向量rl的时间微分向量drl/dt在构成向量M(向量r，ti)的时间点ti处的值，其是心肌结构信息。预分析单元140对所有时间点和所有网格点执行该计算。
[步骤S125]预分析单元140计算关于网格坐标向量rl的时间微分向量drl/dt在构成向量B(向量r，ti)的时间点ti处的值，其是血流结构信息。预分析单元140对所有时间点和所有网格点执行该计算。
[步骤S126]预分析单元140计算血流的网格坐标向量rl在时间点ti处的速度场向量v(向量rl，ti)的时间微分向量dv(向量rl，ti)/dt。预分析单元140对所有时间点和所有网格点执行该计算。
步骤S124至S126中的指数l指示网格编号。此外，本实施方式假设关于心肌总共存在NM,elem个网格坐标，并且关于血流总共存在NB,elem个网格坐标。
接下来，将详细描述用于计算在心肌和血流部分的网格坐标向量rl上的时间微分向量drl/dt在时间点ti处的值的方法。所有网格点处的时间微分向量drl/dt的值以相同的方式计算。
首先，预分析单元140指定网格编号l并获取关于位置向量序列rk(ti)(i＝0、1、2、…、n-1)的信息。接下来，预分析单元140从所获取的信息中提取向量的X分量、Y分量以及Z分量，并将所提取的分量存储为坐标序列Xl(ti)(i＝0、1、2、...、n-1)、Yl(ti)(i＝0、1、2、...、n-1)和Zl(ti)(i＝0、1、2、...、n-1)，其中，n表示文件的数目(n是1或更大的整数)。
接下来，关于X分量，预分析单元140使用插值方法来计算通过所有点Xl(ti)的平滑且连续的曲线Xl(t)。作为插值方法，例如可以使用三次样条。可替选地，可以使用除三次样条之外的任何插值方法。在三次样条中，曲线Xl(t)由表达式(1)使用每个区间ti≤t≤t(i+1)定义的插值多项式X(l，i)(t)定义。
Xl(t)＝Xl，i(t)(ti≤t≤ti+1)(i＝0，1，2，…，n-1) (1)
下面的表达式(2)是插值曲线Xl，i(t)的具体形式。
Xl，i(t)＝al，i+bl，i(t-ti)+cl，i(t-ti)2+dl，i(t-ti)3 (2)
然而，由于n表示一个心搏周期，因此预分析单元140设置循环边界条件，使得第(n+1)值将对应于第零值。例如，预分析单元140设置tn+m＝tm。此外，预分析单元140设置与插值曲线Xl，i(t)的函数值有关的Xl，n+m＝Sl，m。在表达式(2)中，虽然系数al，i、bl，i、cl，i以及dl，i是未知的，但是通过添加包括每个区间ti≤t≤ti+1的微分的用于平滑且连续的曲线的条件，建立线性联立方程。通过计算上求解该联立方程，计算出所有的系数al，i、bl，i、cl，i以及dl，i的集合(参见“Numerical Recipes in C：The Art of Scientific Computing”)。此外，通过以下表达式(3)直接计算dXl/dt。
dXl/dt＝bl，i+2cl，i(t-ti)+3dl，i(t-ti)2 (3)
以这种方式，计算dXl/dt在时间点ti处的值，使得一阶微分和二阶微分关于时间是连续的。
接下来，预分析单元140将系数al，i、bl，i、cl，i以及dl，i的序列存储在预分析结果Ω中。根据相同的过程，还可以针对Y分量和Z分量计算插值方程。预分析单元140针对所有心肌网格点和血流网格点计算这些系数。因此，计算出向量drl/dt的值。
接下来，将描述计算速度场向量v(向量rl，ti)的时间微分向量dv(向量rl，ti)/dt的过程。假设网格点固定在流体中并与流体一起移动，则该时间微分向量dv(向量rl，ti)/dt表示网格点的加速度。这在流体力学领域被称为拉格朗日坐标。
当分析目标是不可压缩流体时，如果加速度乘以流体的密度，则获得施加至网格点上的粒子的力。因此，计算网格点的加速度是计算施加至不可压缩流体中的网格点上的粒子的力的大小和方向。
可以将流体的速度场的分量显示为速度场向量v＝(vx，vy，vz)。尽管在下文中将描述速度场的x分量(vx)，但是速度场的y分量(vy)可以以相同的方式计算。当计算速度场的z分量(vz)时，预分析单元140计算复合函数的微分。因此，建立以下表达式(4)。
根据表达式(4)，获得单个网格点(网格编号k)处的速度场的x分量(vx)在时间点t＝ti处的时间微分。网格编号k的坐标向量rk处的向量drk/dt已经通过插值多项式被计算。▽vx可以通过参考“Finite Element Analysis Chapter 4Finite ElementApproximation”等来计算。通过使用以下表达式(5)来计算剩余速度场(remainingvelocity field)的关于时间的偏微分。
在该表达式(5)中，向量aave表示在拉格朗日坐标上的网格编号k的网格点处的平均加速度向量，并且向量vave表示平均速度向量。向量aave和向量vave通过以下表达式(6)和(7)被定义。
表达式(5)包括拉格朗日坐标上的单个网格点的平均加速度(向量aave)的项和平均速度(向量vave)的项。由于单个网格点移动，因此这是为了校正由于网格点运动引起的表观力。
作为单个网格点执行加速运动时使用的分析方法，例如存在ALE(任意拉格朗日-欧拉)法。
图11示出当基于ALE法控制的网格点自身处于运动中时出现的表观力的原理。在基于ALE法控制的网格点基于根据心脏收缩的跳动而移动之后，网格点返回至其原始位置并且根据心跳移动。在这种情况下，由于网格点的速度时刻改变，因此出现加速度。然而，如果观察装置设置在网格点处，则观察者解释观察者是静止的。看起来好像是由于观察者的加速运动而产生的效应作为表观力作用于流体。因此，需要做出基于观察者的加速运动的校正。校正表观力的过程包括在表达式(5)至(7)中。以这种方式，计算单个网格点处的速度场在时间点t＝ti处的时间微分。预分析单元140计算所有网格点处的血流的单个速度场的时间微分，并且将包括单个网格点处的速度场的时间微分的数据集存储在预分析结果存储单元中120中作为预分析结果Ω。
接下来，将详细描述时间演化计算处理(步骤S112)。
图12和图13是示出时间演化计算处理的过程的流程图。在下文中，将逐步地描述图12所示的处理。
[步骤S131]纹线计算单元150读取点Pij的区域确定标记。纹线计算单元150确定区域确定标记是否指示“T”。当区域确定标记指示“T”时，处理进行至步骤S132。当区域确定标记指示“F”时，纹线计算单元150确定该点已经落在区域的外部并且结束时间演化计算处理，而无需执行计算。因此，当区域确定标记指示“F”时，纹线计算单元150不更新坐标值。
[步骤S132]纹线计算单元150读取点Pij的坐标值。
[步骤S133]纹线计算单元150从存储器102读取计算开始时间点t＝ti和计算结束时间点t＝ti+1。
[步骤S134]纹线计算单元150经由信息读取单元130从文件flui(i).inp读取流体部分在计算开始时间点t＝ti处的网格信息向量B(向量r，ti)和速度场向量v(向量r，ti)。此外，纹线计算单元150经由信息读取单元130从文件stru(i).inp读取向量M(向量r，ti)，向量M是关于心肌部分的弹性体的网格信息(关于心肌的结构的信息)。
[步骤S135]纹线计算单元150经由信息读取单元130从文件flui(i+1).inp读取流体部分在计算结束时间点t＝ti+1处的网格信息向量B(向量r，ti+1)和速度场向量v(向量r，ti+1)。此外，纹线计算单元150经由信息读取单元130从文件stru(i+1).inp读取向量M(向量r，ti+1)。
[步骤S136]纹线计算单元150计算由文件flui(i).inp和文件flui(i+1).inp指示的网格点的速度场向量的范数(速度场向量的长度)，并且通过使用插值表达式来计算相应时间区间中的速度的最大值。纹线计算单元150将所计算的最大值作为Vmax存储在存储器102中。
[步骤S137]如果纹线计算单元150在单个时间演化中计算区间[ti，ti+1]，则精确度不够。因此，纹线计算单元150用Ndiv(Ndiv是1或更大的整数)将区间[ti，ti+1]等分，并且设置中间时间点(t＝tk(k＝1、2、...、和Ndiv))。以这种方式，该区间被分成[ti，ti+Δt]、[ti+Δt，ti+2Δt]、...、和[ti+(Ndiv-1)Δt，ti+1]。纹线计算单元150自身将最优值设置为分割数目Ndiv。分割数目Ndiv可以外部给出。接下来，处理进行至图13中的步骤S141。
在下文中，将逐步地描述图13所示的处理。
[步骤S141]纹线计算单元150通过从k＝1至k＝Ndiv在每个中间时间点(t＝tk(k＝1、2、...、和Ndiv))重复一组步骤S142至步骤S150来执行时间演化计算。因此，获得关于每个中间时间点t＝tk的坐标向量rk，并且获得时间点t＝ti+1处的点Pi+1,j的坐标向量rNdiv＝向量ri+1。
[步骤S142]纹线计算单元150计算关于场在时间点t＝tk处的信息(流体结构信息向量B(向量r，tk)和速度场向量v(向量r，tk))。下面将参照图17来详细描述该步骤。
[步骤S143]纹线计算单元150计算关于弹性体场在时间点t＝tk处的信息(心肌结构信息向量M(向量r，tk))。下面将参照图18来详细描述该步骤。
[步骤S144]纹线计算单元150仅通过时间dt＝tk+1-tk来执行时间演化，以计算向量rk+1。
纹线计算单元150可以如下在步骤S142至步骤S144中每个中间时间点执行时间演化。假设在时间点t＝tk时点Pi,j在向量rk(tk)处，则通过表达式(8)来表示纹线方程。
在表达式(8)中，向量v(向量r，t)是在时间点t和位置向量r处的速度(场)。因此，时间Δt之后的坐标通过数值求解作为常微分方程的表达式(8)来计算。通过用四阶龙格-库塔(Runge-Kutta)法求解表达式(8)得到以下计算表达式。该构思也适用于更高阶龙格-库塔型公式，并且实现了精确的计算(参见“LOW-ORDER CLASSICAL RUNGE-KUTTA FORMULASWITH STEPSIZE CONTROL AND THEIR APPLICATION TO SOME HEAT TRANSFER PROBLEMS”)。在本文中，将描述作为典型示例的四阶龙格-库塔公式。
利用四阶龙格-库塔公式，建立以下表达式(9)和(10)。
纹线计算单元150从龙格-库塔系数表中获取表达式(10)中的系数αI和系数βI。
图14示出龙格-库塔系数表的示例。在图14所示的龙格-库塔系数表中，系数αI和βI是按照I存储的(I是从1至4的整数)。例如，如图14所示的龙格-库塔系数表被预先存储在存储器102中。
可以根据以下项来计算时间点t＝tk时在任何位置向量r处的速度场向量v(向量r，tk)：向量v(向量r，ti)、向量v(向量r，ti+1)、向量B(向量r，ti)、向量B(向量r，ti+1)、向量M(向量r，ti)以及向量M(向量r，ti+1)。假设在时间点t＝tk处的移动边界表面由Sk表示，则还可以计算关于在步骤S142和S143中使用的移动边界表面Sk上的速度场的信息。在下文中，移动边界表面将被简单地称为“边界表面”。因此，可以通过使用表达式(9)和(10)以及图14所示的龙格-库塔系数表来计算时间演化所需的中间值。通过将这些结果代入表达式(8)，计算向量r(tk+1)。下面将参照图20来详细描述时间演化处理。以下描述返回至图13。
[步骤S145]纹线计算单元150执行向量rk+1在心脏内部的位置的确定。纹线计算单元150执行该处理，因为所计算的向量rk+1包括无限时间宽度误差并且可以进入心肌。在作为子例程的位置确定之后，获取指示确定的位置的结果的状态。例如，当向量没有穿过或进入心肌时，获取作为正常终止状态的“0”或“2”。当向量rk+1指示分析目标流体中的单元中的位置时(例如，在心房或心室中)，获取状态“0”。当向量rk+1指示在分析目标流体中的单元的外部的位置时(例如，动脉中)，获取状态“2”。下面将参照图23来详细描述该步骤。
[步骤S146]纹线计算单元150确定确定结果是否指示状态“0”，所述状态“0”指示正常终止。当确定结果指示状态“0”时，处理进行至步骤S150。当确定结果不指示状态“0”时，处理进行至步骤S147。
[步骤S147]纹线计算单元150确定确定结果是否指示状态“2”，所述状态“2”指示终止。获取状态“2”的情况对应于在计算期间点Pij已经通过流体边界移动至系统外部的外部单元(例如移动至大动脉)的情况。当确定结果指示状态“2”时，处理进行至步骤S148。否则，处理进行至步骤S149。
[步骤S148]当点Pij已经落在系统的外部时，纹线计算单元150将纹线上的点Pij(j＝1至j)的区域确定标记设置为“F”。单个区域确定标记“F”指示对应点Pij已经落在分析区域的外部。接下来，纹线计算单元150结束时间演化计算处理。以这种方式，当终止状态指示“2”时，纹线计算单元150确定点Pij已经落在系统的外部，并将区域确定标记设置成“F”。通过连接由点Pi1、Pi2、Pi3、...、和PiN依序连接而形成的曲线来绘制纹线。因此，当任何点Pij被确定为已经落在区域的外部时，没有理由绘制从粒子生成源发出的先前的点。因此，纹线计算单元150将点Pij(j＝1至j)的所有区域确定标记设置成“F”并且结束处理。
[步骤S149]当确定结果没有指示正常终止时(当状态既不是“0”也不是“2”时)，纹线计算单元150减小用作控制参数的时间步长。即，纹线计算单元150通过更多的时间点来进一步分割计算，并且重复一组步骤S142至步骤S147。纹线计算单元150继续充满(repletion)步骤S142至S147，直至计算根据可变时间步进方法收敛。在减小时间步长之后，当纹线计算单元150确定正常终止时，处理进行至步骤S150。
当向量rk+1已经落在预测球体的外部时，纹线计算单元150也没有确定正常终止。在这种情况下，例如，通过增加预测球体的半径并且执行重新计算，纹线计算单元150能够防止向量rk+1落在预测球体的外部。
[步骤S150]纹线计算单元150将向量rk+1存储在存储器中，并且将所存储的值设置为第(k+1)次重复计算的初始值。
[步骤S151]每当纹线计算单元150执行一组步骤S142至步骤S150时，纹线计算单元150向指数k加1并且重复该处理。当纹线计算单元150对所有中间时间点(k＝Ndiv)完成时间演化计算时，处理进行至步骤S152。
[步骤S152]纹线计算单元150将最终计算的向量rNdiv作为向量rk+1存储在存储器中。
通过图12和图13所示的处理来更新纹线上的点的坐标。
图15示出纹线的数据示例。如图15所示，在每个分析时间点设置纹线上的点的坐标值。粒子生成源的坐标值被设置为时间点t＝t0处的纹线上的点的初始值。当时间点t＝t1时，更新指示初始发出的粒子的位置的点的坐标值。接下来，随着时间点的更新，发出新的粒子，并且更新指示发出的新粒子和旧粒子的位置的点的坐标值。在图15中，位置已经从它们之前的时间点改变的点的坐标值被加下划线。
接下来，将描述获得时间点tk(没有针对其给出输出数据)处的场信息(向量B(向量r，tk)、向量v(向量r，tk)以及向量M(向量r，tk))的过程。这些场信息的项用于基于龙格-库塔法的时间演化中描述的上述表达式(9)和(10)的计算中。
仅将时间点t＝ti和t＝ti+1处的数据输出为模拟结果。即，由于没有定义通过龙格-库塔法计算的中间时间tk处的网格坐标和速度场，因此纹线计算单元150根据输出文件的速度场来计算场信息的近似值。该计算的关键考虑是在执行模拟时时刻地移动网格位置。尽管可以以任何方式确定网格位置，但是常常使用任意拉格朗日-欧拉(ALE)法来解决诸如心脏的对象的边界移动的问题。在ALE法中，模拟中使用的坐标被独立地确定，从而不会使所描述的偏微分方程的解的精确度劣化。在许多情况下，偏微分方程被用于该确定。然而，用于确定网格位置的控制方程对数据分析的位置中的一个是不可用的，而只有给出的网格点的输出值是可用的。在这种情况下，网格点的位置随时间不断变化。
图16示出网格点的位置如何随时间不断变化。在图16中，横轴表示时间，而纵轴表示网格点(顶点)的X坐标值。如图16所示，单个网格点的位置随时间不断变化。因此，鉴于该事实，纹线计算单元150通过使用插值方法来估计在任何时间点处的网格位置。
图17是示出用于在时间点t＝tk时计算场信息的处理的过程的示例的流程图。在下文中，将逐步描述图17所示的处理。
[步骤S161]纹线计算单元150在存储器中设置在时间点tk处的场信息。所设置的信息包括向量B(向量r，ti)、向量v(向量r，ti)、向量B(向量r，ti+1)以及向量v(向量r，ti+1)。
时间点tk满足ti≤tk≤ti+1。此外，向量B(向量r，ti)、向量v(向量r，ti)、向量M(向量r，ti)、向量B(向量r，ti+1)、向量v(向量r，ti+1)以及向量M(向量r，ti+1)是已知的。这些向量值已经在步骤S134和步骤S135中从文件读取。
首先，纹线计算单元150对限定流体部分的结构的网格点执行以下处理。然而，为了减少计算量，纹线计算单元150在理论上计算点Pi+1,j的最大移动距离，并且仅对具有等于最大移动距离的半径R的球体内部的网格点执行处理。在下文中，该球体将被称为“预测球体”。
[步骤S162]纹线计算单元150设置预测球体的半径R。以下将参照图26来详细描述该处理。
[步骤S163]纹线计算单元150搜索预测球体的半径R内部的流体网格点，并将预测球体内部的网格点的数目设置为NB,elem。
[步骤S164]纹线计算单元150对预测球体的半径R内部的NB,elem个网格点中的每个网格点执行一组步骤S165至步骤S168。
[步骤S165]纹线计算单元150确定网格点i是否在边界表面上。如果网格点i在边界表面上，则处理进行至步骤S167。如果否，则处理进行至步骤S166。
取决于在步骤S165中对网格点i是否在边界表面上的确定，处理进行至不同的步骤。这是通过插值来计算网格点i在任何时间点tk处的坐标。当网格点i在心肌与血流之间的边界表面S(t)上时，相对于心肌对血流设置无滑移(slip-free)边界条件。
[步骤S166]基于关于网格点i的插值表达式，纹线计算单元150计算网格位置。例如，插值表达式是上面的表达式(2)。预分析结果Ω保存有关网格点i的坐标关于时间的插值表达式(表达式(2))的结果。因此，纹线计算单元150能够根据预分析结果Ω确定任何时间点处的网格坐标。接下来，处理进行至步骤S169。
[步骤S167]如果网格点i在边界表面S(t)上，则纹线计算单元150如下计算网格点i的位置，使得边界条件被满足。
当网格点i在边界表面S(t)上时，由相应速度场形成的拓扑空间中的点被假定为(向量rID(t)，向量vID(t))。由于已经输出了时间点ti和ti+1处的值，因此可以计算通过两点(向量rID(ti)，向量vID(ti))和(向量rID(ti+1)，向量vID(ti+1))的曲线。关于网格位置向量rID(t)和速度场向量vID(t)建立由以下表达式(11)表示的关系。
因此，存在四个条件表达式。对于边界表面S(t)上的网格点i，通过利用关于时间的三阶方程来定义向量ID(t)，四个条件被满足。因此，根据以下表达式(12)确定向量rID(t)。
关于系数向量ai，假设向量rID(t)的分量由ξID(ξ＝x，y，z)来表示，并且相应的微分分量由vξ(ξ＝x，y，z)来表示，则获得以下表达式。
a1ξ＝vξ(ti) (15)
a0ξ＝ξID(ti) (16)
在这些表达式中，使用了以下关系。
Δti＝ti+1-ti (17)
ΔξID＝ξID(ti+1)-ξID(ti) (18)
Δvξ＝vξ(ti+1)-vξ(ti) (19)
以这种方式，获取在任何时间点处的边界表面S(t)上的一组位置向量rID(t)。因此，通过计算流体中所有网格点i的位置，纹线计算单元150能够确定作为关于流体的结构的信息的向量B(向量r，tk)。
[步骤S168]纹线计算单元150将通过使用插值表达式(表达式(12))计算的系数存储在存储器102中。利用表达式(12)，可以同时计算速度场。速度场可以在这里计算或者根据需要通过使用表达式(11)和(12)来计算。因此，纹线计算单元150仅存储通过使用插值表达式(表达式(12))计算的系数。
[步骤S169]每当纹线计算单元150执行一组步骤S165至S168时，纹线计算单元150向指数i加1，并且重复该组步骤S165至S168。当纹线计算单元150已经完成对指数i＝NB,elem的计算时，纹线计算单元150已经完成了在预测球体内部的所有网格点i的位置的计算。
如同场信息的情况，对于作为心肌结构信息的向量M(向量r，t)，纹线计算单元150也能够计算任何时间点处的网格点信息。
图18是示出用于获得在时间点t＝tk时的弹性体场信息的处理的过程的示例的流程图。在下文中，将逐步描述图18所示的处理。
[步骤S171]纹线计算单元150在存储器中设置读取的心肌结构信息。所设置的信息包括向量B(向量r，ti)、向量v(向量r，ti)、向量M(向量r，ti)、向量B(向量r，ti+1)、向量v(向量r，ti+1)以及向量M(向量r，ti+1)。
[步骤S172]纹线计算单元150基于预测球体的半径R来搜索预测球体内部的心肌的结构，并且搜索心肌的网格点i。预测球体的半径R与在流体部分上执行的处理中使用的半径R相同。纹线计算单元150将检测到的网格点i的数目设置为NM,elem。
[步骤S173]纹线计算单元150对预测球体的半径R内部的NM,elem个网格点中的每个网格点i执行一组步骤S174至步骤S177。
[步骤S174]纹线计算单元150确定网格点i是否在边界表面S(t)上。如果网格点i在边界表面S(t)上，则处理进行至步骤S176。如果否，则处理进行至步骤S175。
[步骤S175]纹线计算单元150通过使用预分析结果Ω根据插值表达式(表达式(2))来计算网格点i的位置。接下来，处理进行至步骤S178。
[步骤S176]纹线计算单元150根据表达式(12)来计算网格点i的坐标。通过计算流体中的所有网格点i的位置，确定作为流体结构信息的向量M(向量r，tk)。
[步骤S177]纹线计算单元150将通过使用插值表达式(表达式(12))计算出的系数存储在存储器102中。如果首先执行对流体的处理，则存在其位置与流体侧的网格点i的位置相同的网格点i。因此，替代执行步骤S176和S177，在步骤S167中，纹线计算单元150可以从计算结果中获取网格点i的坐标并存储该坐标。
[步骤S178]每当纹线计算单元150执行一组步骤S174至步骤S177时，纹线计算单元150向i加1，并且重复该组步骤S174至S177。当纹线计算单元150完成对指数i＝NM,elem的计算时，纹线计算单元150结束用于在时间点t＝tk时计算弹性体结构信息的处理。以这种方式，纹线计算单元150能够计算心肌结构信息向量M(向量r，tk)。
接下来，将描述用于计算在任何时间点t处的血流中的网格坐标向量ri(t)处的速度场向量vl(向量rl(t)，t)的方法。
该速度场向量vl(向量rl(t)，t)被用于表达式(9)和(10)的计算中，表达式(9)和(10)被用于使用龙格-库塔法来执行在纹线上的点的时间演化。当输出数据可用的时间点是ti时，获取满足ti≤t≤ti+1的整数i。预分析结果Ω保存单个顶点在时间点ti处的网格坐标向量rl(ti)的值、该点处的速度场向量vl(向量rl(ti)，ti)的值以及速度场的时间微分向量dv(向量rl(ti)，ti)/dt。因此，在时间点ti和ti+1处，已经给出了速度场和微分。因此，包括微分的连续曲线被假设为插值曲线。由于速度场和时间微分是在两个时间点处给出的，因此可以使用由表达式(20)表示的具有四个未知变量的三阶多项式作为插值曲线。
系数向量bk,i被确定，使得在两个时间点ti和ti+1处的速度场向量vi(t)平滑且连续地连接。当ξ＝x，y，z时，通过使用表达式(21)至(26)来确定系数。在下面的表达式中，通过使用“'”来表示各个轴向上的速度的时间微分，例如“v'x(t)”，“v'y(t)”和“v'z(t)”。
以这种方式，计算给定网格点i在任何时间点t处的速度场。
接下来，将描述用于计算在任何时间点t处的流体中的任何坐标向量r的速度场向量v(向量r，t)的方法。
由于给定血流中的任何坐标向量r属于血流，因此获得包括坐标向量r的血流中的无限元Il。假设属于无限元Il的顶点的坐标由向量rl(t)表示，并且速度场由向量vl(向量rl(ti)，t)表示，则坐标向量r处的速度场向量v(向量r，t)根据以下表达式(27)确定。
在表达式(27)中，Ni(向量r，t)被称为结构函数。接下来，用于计算一阶四面体单元的结构函数的具体方法。
图19示出四面体单元的示例。假设四面体单元Il的顶点的坐标是P1(x1，y1，z1)、P2(x2，y2，z2)、P3(x3，y3，z3)以及P4(x4，y4，z4)，则给出如表达式(28)的结构函数。
在表达式(28)中，V是四面体单元Il的体积，并且系数ai、bi和ci被获得为由除了i之外的三个点形成的平面的方程的法向量。例如，N1是根据点P1以外的三个点P2、P3、P4形成的平面的法向量确定的。假设法向量是向量n1，则计算可以被执行为向量n1＝(a1，b1，c1)＝向量(P2P3)×向量(P2P4)。向量(P2P3)是从点P2至点P3的向量。向量(P2P4)是从点P2至点P4的向量。“×”是向量的叉积。关于法向量的方向，当从底表面上的由P2、P3和P4形成的三角形看四面体单元时，顶点P1的方向是正的。此外，由于平面aix+biy+ciz+di＝0的方程通过点P2(或点P3或P4)，因此系数di被确定。因此，当任何点向量r是四面体单元Il中的点时，根据表达式(28)来确定结构函数的值。此外，根据表达式(27)计算位置处的速度场向量v(向量r，t)(参见“Finite Element AnalysisChapter 4Finite Element Approximation”)。
接下来，将详细描述基于龙格-库塔法的时间演化处理。
在图13中的步骤S142和S143中获得心肌和血流网格点i在任何时间点处的位置之后，在步骤S144中，基于龙格-库塔法来执行纹线上的点向量r(tk)的时间演化。
图20是示出基于龙格-库塔法的时间演化处理的过程的示例的流程图。接下来，将逐步描述图20所示的处理。
[步骤S181]纹线计算单元150从i＝1起依次对指数i(i＝1、2、3、4)中的每一个重复成对步骤S182至S183。
[步骤S182]纹线计算单元150获得向量r(tk)所属的血流单元数目Il，并且获取用于基于表达式(27)计算速度向量的信息要素。
[步骤S183]纹线计算单元150通过使用表达式(9)和(10)来计算龙格-库塔法中使用的中间速度向量vi(i＝1、2、3、4)。在向量vi的计算中，还计算了流体中的任何坐标在没有给出数据的任何时间点处的速度场。例如，当向量v1是在时间点tk的坐标向量r(tk)处的速度场的值时，时间点tk并不总是通过预分析给出数据的时间点。此外，坐标向量rk并不总是表示给出数据的流体中的网格点i的坐标。因此，计算出任何坐标向量r在任何时间点t处的速度场。
例如，由于给出坐标向量r(tk)和时间点tk，因此根据血流中的任何坐标向量r在任何时间点t处的速度场向量v(向量r，t)的以上计算方法(表达式(27))来计算速度场。向量v2包括关于向量v1的信息，并且已经计算出向量v1。因此，还可以根据血流中的网格坐标向量ri(t)在任何时间点t处的速度场向量vl(向量rl(t)，t)的计算方法(表达式(20)和(27)来计算速度场。
[步骤S184]每当纹线计算单元150执行成对步骤S182和S183时，纹线计算单元150向i加1，并重复成对步骤S182和S183。当纹线计算单元150完成对i＝4的计算时，纹线计算单元150完成所有向量vi(i＝1、2、3、4)的计算，并且处理进行至步骤S185。
[步骤S185]纹线计算单元150根据表达式(9)来计算纹线在下一个时间点tk+1处的坐标向量r(tk+1)。
接下来，将详细描述心肌内部的位置的确定。
在下文中，将描述确定作为心脏内部的时间演化的结果而获得的向量rk+1的位置的过程。由于龙格-库塔法等包括有限误差，因此作为计算结果而获得的位置向量rk+1可以落在不真实的位置上。
图21示出能够作为计算结果而获得的位置的示例。心脏40的右心系统41和左心系统42的内部部分是存在分析目标流体的区域。右心系统41包括右心房和右心室。在第二实施方式中，左心系统42包括左心房和左心室。连接至左心系统42的动脉43的内部的血流不是分析目标。然而，大动脉的一部分可以被包括模拟中。
当点Pkj在时间点t＝tk处存在于流体中时，在通过如图21所示的时间演化的移动之后存在五个可能的目的地。目的地中的每一个都被给出对应的状态变量(状态)。当点没有移动越过心肌壁并且已经落在分析目标流体的单元内时，状态变量表示“0”。当点已经移动越过心肌壁并且已经落在分析目标流体的单元的外部时，状态变量表示“1”。当点没有移动越过心肌壁并且已经落在分析目标流体的单元的外部时，状态变量表示“2”。当点已经移动越过心肌壁并且已经落在分析目标流体的单元的内部时，状态变量表示“3”。当点没有移动越过心肌壁并且已经落在心肌与分析目标流体的单元之间的边界上时，状态变量表示“4”。
图22是指示状态变量的真值表。以下描述假设“P”表示时间演化之前的点，其始终存在于流体中，而“Q”表示点P在时间演化之后的目的地点。此外，以下描述假设纹线计算单元150仅对点Q执行确定。纹线计算单元150能够通过执行两种类型的确定处理来设置点Q的状态变量。
第一确定处理是流体确定，其中，纹线计算单元150确定目的地点Q是否存在于流体中。如果点Q存在于流体中，则确定T。如果否，则确定F。
第二确定是线确定，其中，纹线计算单元150确定通过连接初始点P与目的地点Q而形成的线PQ是否穿过心肌或心肌表面。如果线PQ穿过心肌(表面)，则确定T。如果否，则确定F。在线确定中，交点的数目也被确定。
当点Q存在于流体中并且线PQ没有穿过心肌时，纹线计算单元150确定正常移动并将“0”设置为状态变量。
当点Q不存在于流体中并且线PQ穿过心肌(表面)时，以下两种情况是可能的：(1)点Q已经移过越过心肌并且落在系统的外部；(2)点Q已经被嵌入在心肌中。在任一情况下，由于需要执行重新计算，因此纹线计算单元150将“1”设置为状态变量。
当点Q不存在于流体中并且线PQ没有穿过心肌(表面)时，纹线计算单元150确定点Q已经经由大动脉等落在模拟系统的外部，并且将“2”设置为状态变量。
即使当点Q存在于流体中时，也存在以下情况，在该情况下，确定出不可能的移动，诸如从左心房至右心房的移动。在这种情况下，当T被确定为流体确定时，T也被确定为线确定，从而交点的数目总是多个。因此，当交点的数目是2或更大时，纹线计算单元150将“3”设置为状态变量。
当T被确定为流体确定并且被确定为线确定时，如果交点的数目是1，则纹线计算单元150确定点已经落在流体与心肌的边界表面上。因此，纹线计算单元150将“4”设置为状态变量。
接下来，将描述以上状态确定处理的过程。
图23是示出用于确定心脏内部的位置的处理的过程的示例的流程图。在下文中，将逐步描述图23所示的处理。
[步骤S191]纹线计算单元150确定坐标向量rk+1在时间演化之后是否存在于流体中。将参照图24来详细描述该处理。
[步骤S192]纹线计算单元150确定由移动向量dr＝向量rk+1-向量rk形成的具有无限长度的线是否穿过心肌，并且对该线与心肌(表面)的交点的数目进行计数。将参照图25来详细描述该处理。
[步骤S193]纹线计算单元150基于图22所示的真值表来确定状态。即，步骤S191和S192的每个结果被作为真值——即真(T)或假(F)——而获得。步骤S193的结果被作为0或更大的整数而获得。纹线计算单元150参考真值表，并且确定值“0”至“4”中的任何一个作为与作为步骤S191至步骤S193的结果而获得的返回值对应的心脏的状态变量。纹线计算单元150使用确定结果作为位置确定结果。
以这种方式，执行位置确定，并且确定状态变量。
接下来，将详细描述用于确定时间演化后位置是否落在流体内的处理(步骤S191)。
图24是示出用于确定时间演化后位置是否落在流体内的处理的过程的示例的流程图。在图24的示例中，纹线计算单元150确定坐标向量rk+1是否存在于流体中。在下文中，将逐步描述图24所示的处理。
[步骤S201]纹线计算单元150获取时间点tk处的坐标向量rk。
[步骤S202]纹线计算单元150获取时间点tk+1处的坐标向量rk+1。
[步骤S203]纹线计算单元150获取预测球体的半径R。
[步骤S204]纹线计算单元150获取具有坐标向量rk为其中心且具有半径R的预测球体的内部的单元的列表(流体单元列表Lf)。
[步骤S205]纹线计算单元150将sizef设置为流体单元列表Lf中的单元的数目。
[步骤S206]纹线计算单元150将“F(假)”设置为结果的初始值。
[步骤S207]纹线计算单元150对单元列表(i＝1、2、...、和sizef)中的单个单元Li重复步骤S208。
[步骤S208]纹线计算单元150确定单元列表中的单元Li是否包括坐标向量rk+1。如果是，则处理进行至步骤S210。如果否，则纹线计算单元150进行至步骤S209。
[步骤S209]每当纹线计算单元150执行步骤S208时，纹线计算单元150向指数i加1并且重复步骤S208。当纹线计算单元150完成对i＝sizef的处理时，纹线计算单元150结束用于确定时间演化后位置是否落在流体内的处理。
[步骤S210]纹线计算单元150将单元Li的单元编号ID存储在存储器中。
[步骤S211]纹线计算单元150将结果改变成“T(真)”。
以这种方式，当目的地点的坐标被包括在单元中的任何单元中时，纹线计算单元150确定坐标向量rk+1存在于流体中，并且将返回值设置成T(真)。然而，当不包括目的地点的坐标时，纹线计算单元150将返回值设置成F(假)。
接下来，将详细描述用于搜索预测球体中的移动向量dr通过的弹性体单元的处理(步骤S192)。
图25是示出用于对由移动向量和心肌形成的线的交点进行计数的处理的过程的示例的流程图。在下文中，将逐步描述图25所示的处理。
[步骤S221]纹线计算单元150获取时间点tk处的坐标向量rk。
[步骤S222]纹线计算单元150获取时间点tk+1处的坐标向量rk+1。
[步骤S223]纹线计算单元150根据所获取的信息来计算移动向量dr＝向量rk+1-向量rk。因此，定义了指示点的移动路径的线。
[步骤S224]纹线计算单元150获取预测球体的半径R的值。
[步骤S225]纹线计算单元150创建在预测球体的半径R内部的弹性体单元的列表(弹性体单元列表Le)。
[步骤S226]纹线计算单元150将sizee设置为弹性体单元列表Le中的单元数目。
[步骤S227]纹线计算单元150将“F(假)”设置为结果的初始值。
[步骤S228]纹线计算单元150对单元列表(i＝1、2、...、和sizee)中的单个单元Li重复步骤S229。
[步骤S229]纹线计算单元150确定第i单元Li是否与线向量dr相交。由于单元Li是多面体，因此纹线计算单元150获得线向量dr关于单元Li的所有表面的交点。如果纹线计算单元150确定在单元Li的任何表面上没有交点，则纹线计算单元150确定线向量dr没有与单元相交。在这种情况下，处理进行至步骤S230。如果纹线计算单元150确定线向量dr与单元相交，则处理进行至步骤S231。
[步骤S230]每当纹线计算单元150执行步骤S229时，纹线计算单元150向指数i加1并且重复步骤S229。当纹线计算单元150完成对i＝sizee的处理时，纹线计算单元150结束用于搜索移动向量dr通过的弹性体单元的处理。
[步骤S231]纹线计算单元150将单元Li的单元编号ID存储在存储器中。
[步骤S232]纹线计算单元150获得指示线的移动向量dr与心肌表面的交点的数目Z(Z是1或更大的整数)，并且将数目Z存储在存储器102中。
[步骤S233]纹线计算单元150将结果改变成真。
如上所述，当不存在交点时，纹线计算单元150确定线向量dr没有与单元相交并且将返回值设置成F。相比之下，当存在至少一个交点时，纹线计算单元150确定线向量dr与单元相交并且将返回值设置成T。在这种情况下，存储与线向量dr具有至少一个交点的单元的单元编号和与心肌表面的交点的数目Z。
通过执行以上处理，纹线计算单元150能够适当地确定点的目的地的状态。当纹线计算单元150在该处理中确定目的地时，通过检查作为目的地确定目标的预测球体中的单元，纹线计算单元150能够更高效地执行处理。接下来，将详细描述用于设置预测球体的半径R的方法。
纹线计算单元150基于时间演化之前的坐标向量rk能够在时间步长Δt内移动多久来设置预测球体。在四阶龙格-库塔法的情况下，根据表达式(9)建立以下不等式。
因此，通过如表达式(30)所指示地定义半径R，半径R表示在时间步长Δt内，纹线上的点P移动的最大距离。
此外，点P在时间步长Δt之后必然存在于球体的内部。此外，中间向量vi也是存在于具有半径R的球体内部的点，其将如下被指示。即，例如，假设在表达式(9)中I＝1来建立表达式(31)。
因此，由中间向量指示的坐标表示具有半径R的球体内部的点。通过当在表达式(9)中I＝1、2、3、4时执行同样的操作，可以看出，所有中间向量vI也是具有半径R的球体内部的点。该半径R被设置为预测球体的半径。如将在以下描述的，纹线计算单元150根据流体的速度场来计算最大值，预测球体半径针对所述流体的速度场而设置。
根据表达式(20)，通过三阶多项式对速度场进行插值。因此，在各个网格点i处的速度场区间[ti，ti+1]处数值地获得最大值。通过从所获得的各个网格点i处的速度场的最大值中选择最大值，将与这两个时间点之间的速度场的最大范数对应的点设置为最大速度|向量vmax|。由于仅需要对数据集执行一次该处理，因此如果最大速度|向量vmax|被存储为预分析结果Ω，则该处理不需要随后被执行。接下来，当设置龙格-库塔法中的时间步长宽度Δt时，最大移动距离被确定为R＝|Δt向量vmax|。
然而，如果时间步长Δt太小，则如上所述设置的预测球体半径R变得小于网格宽度(相邻网格点之间的距离)。因此，在网格点i的离散化中存在最小值Rmin。例如，该最小值Rmin的初始值被设置成作为经验最小值的0.001[m]。最小值可以根据网格的统计分析来计算。
当预测球体的半径R变得小于所计算的Rmin时，预测球体的半径R被设置成Rmin。以这种方式，避免了预测球体内部不存在单元的情况。图26示出如上所述的用于设置预测球体半径R的处理的过程。
图26是示出用于设置预测球体半径的处理的过程的示例的流程图。在下文中，将逐步描述图26所示的处理。
[步骤S241]纹线计算单元150根据速度场和单元的边长的统计分析来设置预测球体的半径的最小值Rmin。
[步骤S242]纹线计算单元150获取时间点ti处的速度向量vi的值。
[步骤S243]纹线计算单元150获取时间点ti处的速度的最大值|向量vi,max|。
[步骤S244]纹线计算单元150获取时间点ti+1处的速度的最大值|向量vi+1,max|。
[步骤S245]纹线计算单元150确定|向量vi,max|是否等于或大于|向量vi+1,max|。当|向量vi,max|等于或大于|向量vi+1,max|时，处理进行至步骤S246。当|向量vi,max|小于|向量vi+1,max|时，处理进行至步骤S247。
[步骤S246]纹线计算单元150将|向量vi,max|设置成|向量vmax|。接下来，处理进行至步骤S248。
[步骤S247]纹线计算单元150将|向量vi+1,max|设置成|向量vmax|。接下来，处理进行至步骤S248。
[步骤S248]纹线计算单元150获取龙格-库塔法中的时间步长dt。
[步骤S249]纹线计算单元150计算|向量vmax|dt，并且将计算结果设置为预测球体半径R。
[步骤S250]纹线计算单元150确定预测球体半径R是否小于预测球体半径的最小值Rmin。当预测球体半径R小于最小值Rmin时，处理进行至步骤S251。当预测球体半径R等于或大于最小值Rmin时，纹线计算单元150结束本处理。
[步骤S251]纹线计算单元150将最小值Rmin设置为预测球体半径R。接下来，纹线计算单元150结束本处理。
通过执行以上处理，纹线计算单元150能够针对预测球体设置适当的半径R。
在第二实施方式中，可以通过执行以下处理更快地执行处理。
[通过分割方法提高计算精确度和速度]
在图12中的步骤S137中，给出输出文件的时间区间ti≤t≤ti+1被分成Ndiv个时间点。将详细描述如何分割时间区间。
在所有单元的纹线计算中，只有少量单元与单个纹线上的点有关。因此，为了减小存储容量和计算时间，使用预测球体。纹线的计算成本随预测球体的半径R3成比例地增加。这将如下被解释。
假设空间单元数目的密度是作为坐标向量r的函数的ρ(向量r)，则作为计算目标的单元的数目Nelem由以下表达式给出。
假设密度ρ(向量r)近似均匀，则ρ(向量r)＝ρ0。因此，Nelem∝R3。
当时间区间ti≤t≤ti+1被分割成Ndiv个时间点时，基于四阶龙格-库塔法的时间演化的数目为Ndiv。同时，通过表达式(30)来确定预测球体半径。因此，当时间步长为1/Ndiv时，预测球体半径也变成1/Ndiv。由于单个计算量与预测球体的半径R3成比例，因此计算量由Ndiv -3表示。由于该计算被重复Ndiv次，因此总的计算量由Ndiv -2来表示。图27示意性地示出该计算的构思。
图27示出通过时间分割而实现计算量的减少的构思。图27假设计算作为点P的目的地的点Q的情况。当不执行时间分割并且使用长的时间步长时，分析具有点P为其中心的大的预测球体45。另一方面，当执行时间分割并且使用短的时间步长时，点P处的粒子的可移动范围减小。因此，优选地分析比预测球体45小的预测球体46。即，纹线计算单元150并非对大的预测球体45执行单个计算，而是通过将时间区间分割成Ndiv个时间点并且获得分割的路径来得到较小的预测球体46。因此，减少了作为计算目标的单元的数目，并且减少了计算量。
图28示出根据时间分割数目的总计算时间的变化的示例。在图28的曲线图中，横轴表示分割数目Ndiv，而纵轴是计算时间。在图28中，虚线表示理论计算时间51的转变，而实线表示实际计算时间52的转变。实际计算时间52是时间区间ti≤t≤ti+1中的总计算时间的测量结果。当Ndiv＝0时，不使用预测球体，并且所有的单元被用作计算目标。
理论上，分割数目Ndiv越大，计算时间将越短。然而，实际上，如果分割数目太大，则每分割后单位时间执行的处理的次数增加，例如创建预测球体中的单元的列表。从而，处理时间增加。因此，存在作为分割数目Ndiv的最优值。例如，在开始纹线的计算之前，纹线计算单元150提前针对目标系统来测量最优值作为分割数目Ndiv，同时改变分割数目Ndiv。
[通过统计分析实现的预测球体半径的最小值的确定方法]
通过使用有限数目的离散点信息项来执行流体模拟。因此，如果设置过小的预测球体半径，则没有离散点信息项会包括在预测球体中。这意味着存在预测球体半径的下限。同时，由于预测球体半径通过表达式(30)设置，因此存在时间步长的下限。当通过使用预测球体执行计算时，设置预测球体半径的下限Rmin。当预测球体半径等于或小于下限时，将下限Rmin设置为预测球体半径。当预测球体半径等于或小于下限时，由于使用比预测球体半径大的Rmin，因此纹线上的点没有落在预测球体的外部。因此，为了稳定地进行计算，设置下限Rmin是重要的。此外，下限Rmin的值与时间分割数目Ndiv的设置有关。由于R0＝|Δt向量vmax|是最大移动距离，当确定了下限Rmin的值时，通过使用上取整函数来根据以下表达式设置分割数目Ndiv。
由于分割数目Ndiv的值与计算速度和计算精确度有关，因此设置下限Rmin在计算速度和计算精度方面也是重要的。
然而，在设定下限Rmin时需要注意。具体地，引入概率模型，并且纹线计算单元150执行允许计算失败的推测性计算。当计算失败时，纹线计算单元150通过使用使计算必定成功的参数来执行计算。该重新计算所需的时间被认为是处罚。设置统计上使包含处罚的计算时间最小化的下限Rmin。接下来，将详细描述如何设置下限Rmin。
首先，纹线计算单元150将半径Rw设置为用作处罚的最差值。模拟中使用的单元的最长边长被用作半径Rw。
图29示出了单元的边长的分布的示例。在图29的曲线图中，横轴表示单元的边长，而纵轴表示指示边长中的每一个的概率密度。在曲线图中，由虚线来表示两种类型的心脏模拟#1和#2的概率密度函数61和62。
在图29的示例中，模拟中使用的一些单元非常粗糙，并且最长边是0.008[m]。因此，如果将这些单元直接用于计算，则分割数目Ndiv也减小，并且速度也降低。然而，大多数计算可以在不使用最长边的情况下执行。
在设置作为最差值的半径Rw之后，纹线计算单元150设置用于计算的下限Rmin。在该操作中，纹线计算单元150在假设计算失败的同时计算统计计算成本。此外，鉴于推测性计算的执行，纹线计算单元150选择使计算量最小化的半径作为下限Rmin，其中，所述计算量包括计算失败时出现的处罚。虽然这显著取决于假设的概率模型，但是将通过使用简单的示例来描述推测性计算。首先，使用纹线上的点的计算假设所有单元都可以是具有等概率的分析目标。在这种情况下，计算根据图29中的概率密度函数61和62而成功。
以下描述假设：当半径是R时，计算时间是T。在这种情况下，对应的计算量与R3成比例。此外，以下描述假设：通过将半径减小至βR(0<β≤1)计算成功的概率为p[％]。计算成功的情况是点的目的地落在具有半径βR的预测球体的内部的情况。当计算失败时，通过使用半径R来执行重新计算。在这种情况下，包括处罚的计算时间T'通过表达式(34)来表示。
图30示出基于预测球体半径的计算成本的变化。在图30中，纵轴表示计算成本(T'/T)，而横轴表示预测球体半径R。如图30所示，可以看出，当预测球体半径R大约为0.003[m]时的包括处罚的成本比当预测球体半径是最差值半径Rmax(0.008[m])时的包括处罚的成本小大约10％。也就是说，可以看出，存在实现包括处罚的在统计上的最小计算量的半径。在实践中，由于执行有限数目的计算次数，因此所有单元并非以等概率用作分析目标。因此，比以上值更小的值被用作半径。
在这种情况下，尽管可能指示存在最优预测球体半径R的事实，但是最优值本身倾向于被高估。因此，为了防止这样的高估，可以通过以下操作来获得如图31所示的移动距离的分布：获得单个离散点处的速度场的范数，并且将时间步长乘以范数。
图31示出移动距离的概率分布。在图31中，横轴表示移动距离，而纵轴表示出现的数目。即使当存在两个具有相同形状和尺寸的四面体时，如果一个四面体上的速度场的大小是另一个四面体上的速度场的大小的两倍，则点通过目标四面体所需的时间减半。图32基于移动距离的分布通过使用表达式(34)示出作为预测球体半径R的函数的包括处罚的计算效率。
图32示出在假设了概率分布时的计算效率曲线。在图32中，横轴表示预测球体半径R，而纵轴表示计算效率。较小的计算效率值指示较好的效率。根据图32，最佳的计算效率与近似从0.0012[m]至0.0016[m]的范围中的预测球体半径R相对应。这是支持经验最优值Rmin＝0.0015±0.0005[m]的模型。
以这种方式，获得预测球体半径R的最优值。即，虽然当离散点落在预测球体的外部时纹线计算单元150需要执行作为处罚的重新计算，但是纹线计算单元150能够鉴于该处罚来设置预测球体的最有效的半径R。因此，纹线计算的效率提高。
[具体计算示例]
在下文中，将具体描述在通过使心脏模拟器分析心肌和冠状循环来实际计算纹线时测量的计算速度。
首先，将描述计算中使用的数据。输出对应于单次心跳的模拟结果，并且每0.01[秒]输出100个心脏状态数据。由于心脏搏动，因此速度场随时间变化。因此，速度场是不稳定的流。此外，由于搏动，心肌重复舒张和收缩。因此，心肌表面也移动，导致移动边界问题。为了描述心脏中的血流在该系统中的传输，将粒子生成源布置在心脏中的多个点处，并且计算纹线。
当可视化设备100从外部读取关于心肌和流体的信息时，纹线的数目M和计算的次数N被输入至可视化设备100。此外，纹线生成源的位置被输入至可视化设备100。纹线生成源被布置在心脏中的左心室和右心房中。根据纹线生成源的布置的指示，例如，可视化设备100将20个纹线生成源随机地布置在具有0.05[m]的半径的球体的内部的各个流体部分中。此外，计算次数被设置成与0.3次心跳相对应的30。通过使用四阶龙格-库塔法来执行时间演化计算。
在设置以上初始条件之后，当可视化设备100开始纹线计算时，每个步长从粒子生成源生成粒子。根据图9所示的纹线计算流程图来计算单个粒子在时间步长Δt＝0.01[秒]之后的位置。所计算的纹线被显示在监视器21上。
图33示出纹线的显示的示例。在图33的示例中，多条纹线71被叠加在心脏70的截面图像上。
当计算纹线时，纹线计算单元150使用预测球体来减小计算量并提高计算精确度。在使用预测球体的计算中，通过分割数目Ndiv来进一步分割时间步长Δt＝0.01[秒]。将实现统计上的最小计算量的最优值自动设置为分割数目Ndiv。在该模拟中，分割数目Ndiv被设置成大约3至7。
通过使用分割数目Ndiv进行优化来定量地测量计算速度，并且图28示出测量的结果。图28中的示例还示出当没有使用预测球体时的纹线的计算时间，以便阐明通过使用预测球体而获得的有利效果。Ndiv＝0与当没有使用预测球体时——即当所有单元被用作计算目标时——所需的时间相对应。
当使用预测球体(Ndiv＝1)时，即，当基本上不执行时间分割时，需要时间来建立预测球体。因此，与不使用预测球体的情况相比，需要更多的时间(大约1.47倍)。然而，当Ndiv＝2时，与Ndiv -2成比例的计算成本的降低效应变得大于建立预测球体的固定成本。因此，与不使用预测球体的情况相比，计算时间减少大约22.9％。虽然计算时间随着分割数目Ndiv的增加而缩短，但是可以看出，由于建立预测球体的固定成本，所以存在性能改善的上限。在图28的示例中，最佳分割数目是4。如果分割数目Ndiv进一步增加，则由于固定成本的效应具有更大的影响，因此计算时间开始增加。当Ndiv＝10时，与最佳情况相比，性能恶化了大约53％。
如上所述，虽然分割数目Ndiv存在最佳值，但是在实际计算中，最佳分割数目在计算纹线之前设置。此外，在图28的计算中，由于网格点的数量不大于大约50000，因此在不使用预测球体的情况下的计算时间短于使用预测球体并且分割数目Ndiv为1的情况下的计算时间。然而，如果模拟系统增大，则不使用预测球体的情况下的计算量与网格点的数目成比例地增加。例如，如果网格数目增加10倍，则时间增加10倍。然而，当使用预测球体时，由于预测球体半径R仅基于时间步长和速度的最大值来设置并且始终恒定，因此计算量也是恒定的。因此，系统(对其执行模拟)越大，通过使用预测球体而获得的有利效果将越大。
此外，将描述通过对预测球体进行时间分割而获得的另一有利效果——精确度的提高。虽然表达式(9)至(10)是在四阶龙格-库塔法的情况下给出的，但是已知：误差阶数被给出为O(Δt4)。因此，当时间由分割数目Ndiv分割时，单个龙格-库塔运算中的误差的结果是Ndiv -4倍。即使当在时间区间ti≤t≤ti+1中基于龙格-库塔法的计算总计Ndiv次并且误差累积了Ndiv次时，总误差仍然是Ndiv -3。当分割数目Ndiv是4时，总误差是1/64倍。
图34示出当时间步长改变时的轨迹精确度的变化。在图34中，横轴表示模拟时间，而纵轴表示离散点的x坐标。图34中的示例示出当时间步长Δt被设置成“0.01”、“0.05”、“0.001”以及“0.0001”时某一离散点的x坐标的时间变化(轨迹)。
当积分后的坐标变化量小或者积分运算的数目少时，不管时间步长如何，轨迹都不改变。然而，在积分运算的数目较大的区域(时间t>2.5[秒])中，与当时间步长Δt＝0.01[秒]时对应的轨迹从与当时间步长Δt＝0.001[秒]和0.0001[秒]时分别对应的轨迹显著偏移。然而，与当时间步长Δt＝0.001[秒]和Δt＝0.0001[秒]时对应的轨迹彼此非常接近。因此，可以看出，通过减小时间步长防止了计算误差的累积。
如上所述，设置分割数目Ndiv对于计算速度和计算精确度是重要的。此外，分割数目Ndiv取决于预测球体半径的最小值Rmin来设置。在实际计算中，如图31所示，鉴于单元的边长和速度场二者来假设概率模型，并且计算概率分布。此外，根据概率分布，通过使用表达式(34)来计算计算效率，并且计算实现最佳计算效率的最小值Rmin。即，如图32所示，存在作为预测球体的半径的最小值，并且最小值Rmin大约为0.0015±0.0005[m]。因此，总是以最小的计算成本执行计算。
接下来，将描述计算稳定性。在第二实施方式中，通过考虑基于ALE网格坐标的运动的校正，根据表达式(4)来计算单个速度场。
图35示出添加了ALE网格点的运动的校正的速度场的示例。图35示出以下计算结果：当通过线性插值执行速度场的评估时的速度的改变(虚线)的计算结果，以及当由ALE网格坐标的运动引起的表观力的校正根据表达式(4)执行时的速度的改变(实线)的计算结果。
尽管存在线性插值近似达到非常好的结果的情况，但是如果线性近似用于评估速度场相对于时间快速变化的区域([0.03≤t≤0.1])中的速度场，则速度场的评估精确度恶化。在这样的速度场的评估精确度差的区域中，由于出现大的数值误差，因此在计算纹线时，纹线可能进入心肌。这在图36中示意性地示出。
图36示出防止纹线嵌入到心肌中的有利效果——由于计算精确度提高而获得的效果。当在心脏收缩期中，纹线存在于心肌附近时，纹线上的点靠近心肌的移动方向上的心肌运动。
如果速度场的计算精确度低，则速度场被低估，并且纹线上的点的速度可以被计算为低于精确值。如果心肌附近的纹线上的点的速度被低估并且低于心脏收缩速度，则心肌超过纹线上的点。因此，纹线嵌入到心肌中。一旦纹线上的点被嵌入到心肌中，由于在心肌上没有定义流体的速度场，因此计算失败。
为了避免这种情况，第二实施方式在计算速度场时根据表达式(4)来校正ALE网格坐标的运动。以这种方式，由于提高了速度场的计算精确度，因此心肌附近的纹线上的单个点的速度不太频繁地低于心脏收缩速度，并且心肌不太频繁地超过纹线上的任何点。因此，避免了纹线嵌入到心肌中的情况。
如上所述，通过校正由单个网格点自身的加速运动引起的表观力，提高了纹线可视化设备的计算稳定性。
[其他实施方式]
纹线上的点彼此独立，彼此之间没有相互作用。因此，由于并行化适当地适用，因此可以通过使用消息传递接口(MPI)或开放式多处理(OpenMP)来以并行方式执行纹线的计算。以这种方式，计算速度提高。
此外，虽然在第二实施方式中已经示出通过使用心脏模拟的结果来将血流的纹线可视化的示例，但是第二实施方式也可以应用于其他流体模拟的结果。例如，当将可变翼机构布置在机动车的后部时，可以模拟和分析当可变翼机构被激活时可变翼机构周围的空气的流动。通过使用根据第二实施方式的可视化设备100，通过纹线来将当可变翼机构被激活时可变翼机构周围的空气的流动可视化。此外，第二实施方式也适用于指示当飞机的翼摆动时在摆动翼周围的空气的流动的流体模拟结果。
根据一个方面，即使当结构变形时，也追踪纹线。

Claims (6)

1.一种纹线可视化设备，包括：
存储装置，用于存储流体信息，所述流体信息包括位置信息，所述位置信息指示多个网格点在具有第一时间间隔的多个第一时间点中的每一个处的位置，所述多个网格点随着流体模拟的分析时间的进行在分析空间中以加速运动移动，并且所述流体信息包括速度信息，所述速度信息指示所述多个网格点上的流体在所述多个第一时间点中的每一个处的速度；以及
处理装置，用于：
基于所述流体信息并且通过使用表达式来计算所述多个网格点上的流体在所述多个第一时间点中的每一个处的速度的时间微分值，所述表达式包括校正值，所述校正值用于校正由所述位置信息表示的所述多个网格点的加速运动引起的误差；
基于所述多个网格点上的流体在所述多个第一时间点中的每一个处的速度的时间微分值和速度，计算随着所述分析时间的进行从粒子生成源依序输出的一系列多个粒子在多个第二时间点中的每一个处的位置，所述多个第二时间点具有短于所述第一时间间隔的第二时间间隔；以及
生成关于指示所述一系列多个粒子的纹线的显示信息。
2.根据权利要求1所述的纹线可视化设备，其中，
所述处理装置基于所述位置信息来计算所述多个网格点中的每一个在所述多个第一时间点中的每一个处的网格点移动速度和网格点移动加速度，并且
所述处理装置通过使用包括所述网格点移动速度和所述网格点移动加速度作为变量的表达式来计算所述多个网格点上的流体在所述多个第一时间点中的每一个处的速度的时间微分值。
3.根据权利要求1或2所述的纹线可视化设备，其中，
所述处理装置生成插值曲线，所述插值曲线中的每一个连接指示所述多个网格点中的相应的一个网格点上的流体在所述多个第一时间点处的速度的点，所述插值曲线中的每一个具有基于由分别与所述多个第一时间点相对应的点指示的流体的速度的时间微分值的斜率，并且所述插值曲线中的每一个表示流体的速度随时间的变化，并且
所述处理装置基于所述多个网格点的各个插值曲线来计算所述一系列多个粒子在所述多个第二时间点中的每一个处的位置。
4.根据权利要求1所述的纹线可视化设备，
其中，所述存储装置还存储结构信息，所述结构信息指示结构的形状在所述分析空间中随时间的变化，并且
其中，当基于第一分析时间点处的第一纹线来计算第二分析时间点处的第二纹线时，所述处理装置还进行以下操作：
将包括所述分析空间中的所述第一纹线上的第一位置处的离散点的部分区域设置为所述离散点的分析目标区域；
基于所述分析目标区域中的流体的速度来计算第二位置，所述第二位置指示所述离散点上的粒子在所述第二分析时间点处的目的地，所述速度由所述流体信息指示；
基于关于所述分析目标区域中的结构的信息，确定由所述分析目标区域中的结构在所述第二分析时间点处占据的区域，所述信息由所述结构信息指示；
基于所述第一位置和所述第二位置来确定所述第二纹线进入或者未进入所述占据的区域；以及
当确定所述第二纹线未进入所述占据的区域时，生成关于通过所述第二位置的第二纹线的显示信息。
5.根据权利要求1所述的纹线可视化设备，还包括：
显示装置，用于基于所述显示信息来显示所述纹线。
6.一种由计算机执行的纹线可视化方法，所述纹线可视化方法包括：
基于流体信息并且通过使用表达式来计算多个网格点上的流体在具有第一时间间隔的多个第一时间点中的每一个处的速度的时间微分值，所述多个网格点随着流体模拟的分析时间的进行在分析空间中以加速运动移动，所述流体信息包括位置信息，所述位置信息指示所述多个网格点在所述多个第一时间点中的每一个处的位置，并且所述流体信息包括速度信息，所述速度信息指示所述多个网格点上的流体在所述多个第一时间点中的每一个处的速度，所述表达式包括校正值，所述校正值用于校正由所述位置信息表示的所述多个网格点的加速运动引起的误差；
基于所述多个网格点上的流体在所述多个第一时间点中的每一个处的速度的时间微分值和速度，计算随着所述分析时间的进行从粒子生成源依序输出的一系列多个粒子在多个第二时间点中的每一个处的位置，所述多个第二时间点具有短于所述第一时间间隔的第二时间间隔；以及
生成关于指示所述一系列多个粒子的纹线的显示信息。
CN201711261432.8A 2016-12-06 2017-12-04 纹线可视化设备和方法 Active CN108170895B (zh)

Applications Claiming Priority (4)

Application Number Priority Date Filing Date Title
JP2016-236734 2016-12-06
JP2016236734A JP6741252B2 (ja) 2016-12-06 2016-12-06 流脈線可視化装置、流脈線可視化方法、および流脈線可視化プログラム
JP2017-075671 2017-04-06
JP2017075671A JP6760603B2 (ja) 2017-04-06 2017-04-06 流脈線可視化装置、流脈線可視化方法、および流脈線可視化プログラム

Publications (2)

Publication Number Publication Date
CN108170895A true CN108170895A (zh) 2018-06-15
CN108170895B CN108170895B (zh) 2021-12-14

Family

ID=60569597

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711261432.8A Active CN108170895B (zh) 2016-12-06 2017-12-04 纹线可视化设备和方法

Country Status (3)

Country Link
US (1) US10799191B2 (zh)
EP (1) EP3333738A1 (zh)
CN (1) CN108170895B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111210522A (zh) * 2020-01-14 2020-05-29 西南石油大学 利用fem在三维非结构网格流场内追踪流线分布的方法

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP3166032B1 (en) * 2014-07-03 2018-10-03 Fujitsu Limited Biometric simulation device, method for controlling biometric simulation device, and program for controlling biometric simulation device
CN116611176B (zh) * 2023-07-21 2023-09-15 中国空气动力研究与发展中心空天技术研究所 一种三维前缘可定制的流线追踪类乘波体设计方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101667137A (zh) * 2009-10-12 2010-03-10 哈尔滨工程大学 使用方向滤波技术的手指静脉纹路提取方法
CN104318601A (zh) * 2014-10-22 2015-01-28 电子科技大学 一种流体环境下人体运动仿真方法
CN104392105A (zh) * 2014-10-27 2015-03-04 西北工业大学 一种ib方法
US20150242545A1 (en) * 2014-02-21 2015-08-27 Junghyun Cho Method of Simulation of Moving Interfaces using Geometry-Aware Volume of Fluid Method
WO2016002054A1 (ja) * 2014-07-03 2016-01-07 富士通株式会社 生体シミュレーション装置、生体シミュレーション装置の制御方法、および生体シミュレーション装置の制御プログラム

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2003006552A (ja) 2001-06-19 2003-01-10 Toshiba Corp シミュレーションシステム、シミュレーション方法及びシミュレーションプログラム
US7472046B2 (en) * 2003-06-27 2008-12-30 Lucasfilm Entertainment Company Ltd. Apparatus and method of simulating the movement of elements through a region of 3D space
JP5240132B2 (ja) * 2009-09-04 2013-07-17 富士通株式会社 熱流体シミュレーション解析装置
CN105580019B (zh) * 2013-07-30 2018-09-14 哈特弗罗公司 为优化诊断性能利用边界条件模型化血流的方法和系统
US9790770B2 (en) * 2013-10-30 2017-10-17 The Texas A&M University System Determining performance data for hydrocarbon reservoirs using diffusive time of flight as the spatial coordinate
JP6362853B2 (ja) 2013-11-20 2018-07-25 キヤノンメディカルシステムズ株式会社 血管解析装置、および血管解析装置の作動方法
US20170311916A1 (en) 2014-10-08 2017-11-02 EBM Corporation Blood-flow analysis device for blood-flow simulation, method therefor, and computer software program

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101667137A (zh) * 2009-10-12 2010-03-10 哈尔滨工程大学 使用方向滤波技术的手指静脉纹路提取方法
US20150242545A1 (en) * 2014-02-21 2015-08-27 Junghyun Cho Method of Simulation of Moving Interfaces using Geometry-Aware Volume of Fluid Method
WO2016002054A1 (ja) * 2014-07-03 2016-01-07 富士通株式会社 生体シミュレーション装置、生体シミュレーション装置の制御方法、および生体シミュレーション装置の制御プログラム
CN104318601A (zh) * 2014-10-22 2015-01-28 电子科技大学 一种流体环境下人体运动仿真方法
CN104392105A (zh) * 2014-10-27 2015-03-04 西北工业大学 一种ib方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
CHIH-YUNG WEN等: "investigation of pulsatile flowfield in healthy thoracic aorta models", 《ANNALS OF BIOMEDICAL ENGINEERING》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111210522A (zh) * 2020-01-14 2020-05-29 西南石油大学 利用fem在三维非结构网格流场内追踪流线分布的方法

Also Published As

Publication number Publication date
US10799191B2 (en) 2020-10-13
US20180153483A1 (en) 2018-06-07
EP3333738A1 (en) 2018-06-13
CN108170895B (zh) 2021-12-14

Similar Documents

Publication Publication Date Title
Praveen et al. Low cost PSO using metamodels and inexact pre-evaluation: Application to aerodynamic shape design
Krone et al. Fast Visualization of Gaussian Density Surfaces for Molecular Dynamics and Particle System Trajectories.
CN102592133B (zh) 使用概率推进树进行评估的方法和系统
Giannakoglou et al. Aerodynamic shape design using evolutionary algorithms and new gradient-assisted metamodels
CN108170895A (zh) 纹线可视化设备和方法
CN108153934A (zh) 纹线可视化设备和方法
JP6760603B2 (ja) 流脈線可視化装置、流脈線可視化方法、および流脈線可視化プログラム
JP2022158542A (ja) 推論プログラム、推論方法および情報処理装置
Estep et al. Fast and reliable methods for determining the evolution of uncertain parameters in differential equations
Tutmez et al. Spatial estimation of some mechanical properties of rocks by fuzzy modelling
US9001117B2 (en) Selective flow visualization of traced particles
Dagliati et al. Inferring temporal phenotypes with topological data analysis and pseudo time-series
Khrissi et al. A performant clustering approach based on an improved sine cosine algorithm
CN116597294A (zh) Slam地图拓扑评估方法、装置、电子设备及存储介质
KR101068937B1 (ko) 기능성 게임의 개인화된 학습경로 생성방법
Nieto Ramos et al. Bayesian inference for fitting cardiac models to experiments: estimating parameter distributions using Hamiltonian Monte Carlo and approximate Bayesian computation
Liu et al. Parameter estimation in computational biology by approximate bayesian computation coupled with sensitivity analysis
Li et al. An improved algorithm optimization algorithm based on RungeKutta and golden sine strategy
Kent et al. Bayesian quality diversity search with interactive illumination
Michaelides et al. Statistical abstraction for multi-scale spatio-temporal systems
Gaier et al. Feature space modeling through surrogate illumination
Krispel et al. A benchmark for distance measurements
Mohsen et al. MULTI-POPULATION CUCKOO SEARCH ALGORITHM FOR SOLVING GLOBAL PROBLEMS OF OPTIMIZATION
Mavhemwa et al. Uniform spatial subdivision to improve boids algorithm in a gaming environment
Eero Optimal design for disc golf by computational fluid dynamics and machine learning

Legal Events

Date Code Title Description
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant