纹线可视化设备和方法与流程

本文中讨论的实施方式涉及纹线可视化设备和方法。

背景技术：

流体力学是力学中的学术领域之一，并且描述流体的行为。流体力学已经被应用于各种工业领域，在所述工业领域中，不仅空气或水的流动而且诸如温度或浓度等的物理量的转移都被作为问题来处理。例如，流体力学已经应用于风洞实验以评估汽车原型，并在实验结果的基础上对这些汽车的空气动力学特性进行了优化。然而，这些风洞实验是非常昂贵的。因此，替代风洞实验，已经通过使用计算流体力学进行了模拟风洞实验的计算机模拟(流体模拟)。

近来计算机性能的提高已经使流体模拟得到了迅速的进展。因此，流体模拟不仅被应用于飞行器、汽车、铁路车辆、船舶等的空气动力学特性的评估，而且还被应用于心脏、血管等的血流状态的分析。

当进行流体模拟时，分析结果被可视化，使得可以在视觉上容易地理解分析结果。对流体模拟的结果进行可视化的一种手段是显示纹线。纹线是通过连接已经通过空间中的特定点的流体粒子而形成的曲线。在风洞实验中，从预定位置喷出的烟雾的痕迹是纹线。即，通过在流体模拟中计算纹线并显示所述纹线，可以将流体中的粒子的运动——如风洞实验中的烟雾的痕迹——可视化，而无需执行任何风洞实验。

已经提出了与流体模拟相关的各种技术。例如，已经详细地提出了执行高速模拟并快速且平滑地表示流体中的场景的技术。还提出了易于将结构-流体分析模拟的结果应用于对血管异常的诊断的技术。还提出了能够使不熟悉计算流体力学的使用者如医生进行适当的血流模拟的设备。此外，已经发表了与流体模拟相关的各种论文。例如参见以下文献：

日本特许专利公开第2003-6552号

日本特许专利公开第2015-97759号

国际公开小册子第WO2016/056642号

Tino Weinkauf and Holger Theisel，“Streak Lines as Tangent Curves of a Derived Vector Field”，IEEE Transactions on Visualization and Computer Graphics，2010年11月-12月，第16卷，第6期。

Erwin Fehlberg，“LOW-ORDER CLASSICAL RUNGE-KUTTAFORMULAS WITH STEPSIZE CONTROL AND THEIR APPLICATION TO SOME HEAT TRANSFER PROBLEMS”，NASATECHNICAL REPORT，NASATR R-315，1969年7月。

J.Donea，A.Huerta，J.-Ph.Ponthot and A.Rodríguez-Ferran，“Arbitrary Lagrangian-Eulerian methods”，Encyclopedia of Computational Mechanics，John Wiley&Sons Ltd.，2004年11月，第413-437页。

Seiryo Sugiura，Takumi Washio，Asuka Hatano，Junichi Okada，Hiroshi Watanabe，Toshiaki Hisada，“Multi-scale simulations of cardiac electrophysiology and mechanics using the University of Tokyo heart simulator”，Progress in Biophysics and Molecular Biology，2012年10月-11月，第110卷，第380-389页。

然而，这些传统纹线分析技术——例如在分析汽车周围的空气流动的情况下——是基于分析空间中的结构不会变形的假设。当结构变形时，难以精确地追踪纹线。

技术实现要素：

在一个方面中，实施方式的目的是即使当结构变形时也实现对纹线的追踪。

根据一个方面，提供了一种纹线可视化设备，其计算流体模拟时间中的多个分析时间点处的指示一系列多个粒子的纹线，并且显示纹线，纹线可视化设备包括存储装置，用于保持：结构信息，其指示分析空间中的结构的形状的时间变化；以及，流体信息，其指示分析空间中的存在流体的区域中的多个点处的流体的速度的空间变化和时间变化中的至少一个；以及，纹线可视化设备包括处理装置，用于：当基于第一分析时间点处的第一纹线来计算第二分析时间点处的第二纹线时，将包括分析空间中的第一纹线上的第一位置处的离散点的部分区域设置为离散点的分析目标区域；基于分析目标区域中的流体的速度来计算指示离散点上的粒子在第二分析时间点处的目的地的第二位置，所述速度由流体信息指示；基于关于分析目标区域中的结构的信息来确定由分析目标区域中的结构在第二分析时间点处占据的区域，所述信息由结构信息指示；基于第一位置和第二位置来确定第二纹线进入或者未进入占据区域；并且当确定第二纹线未进入占据区域时，显示通过第二位置的第二纹线。

附图说明

图1示出根据第一实施方式的纹线可视化设备的配置示例；

图2示出根据第二实施方式的系统配置示例；

图3示出可视化设备的硬件配置示例；

图4示出纹线计算示例；

图5是示出可视化设备的功能的框图；

图6示出弹性体信息文件组的示例；

图7示出流体信息文件组的示例；

图8是示出纹线可视化处理的过程的示例的流程图；

图9和图10示出示出时间演化计算处理的过程的流程图；

图11示出纹线的数据示例；

图12示出网格点的位置如何随时间不断变化；

图13是示出用于计算在时间点t＝tk时的场信息的处理的过程的示例的流程图；

图14是示出用于计算在时间点t＝tk时的弹性体场信息的处理的过程的示例的流程图；

图15示出能够作为计算结果而获得的位置的示例；

图16是指示状态变量的真值表；

图17是示出状态确定处理的过程的示例的流程图；

图18是示出用于确定时间演化后位置是否落在流体内的处理的过程的示例的流程图；

图19是示出用于搜索移动向量dr通过的弹性体单元的处理的过程的示例的流程图；

图20是示出用于设置预测球体的半径的处理的过程的示例的流程图；

图21示出通过时间分割而实现计算量的减少的构思；

图22示出根据时间分割数目的总计算时间的变化的示例；

图23示出单元的边长的分布的示例；

图24示出基于预测球体半径的计算成本的变化；

图25示出移动距离的概率分布；

图26示出在假设概率分布时的计算效率曲线；

图27示出纹线的显示的示例；以及

图28示出当时间步长改变时的轨迹的精确度的变化。

具体实施方式

下面将参照附图来描述实施方式，其中，相同的附图标记贯穿全文指代相同的要素。两个或更多个实施方式可以彼此组合而不会产生不一致。[第一实施方式]

首先，将描述第一实施方式。第一实施方式提供了即使当结构变形时也能够追踪纹线的纹线显示设备。

图1示出根据第一实施方式的纹线可视化设备的配置示例。该纹线可视化设备10计算在流体模拟时间中的多个分析时间点处的指示一系列多个粒子的纹线，并且显示纹线。例如，纹线可视化设备10是包括作为算术处理装置的处理器和作为主存储装置的存储器的计算机。纹线可视化设备10包括存储单元11和处理单元12。

存储单元11保持结构信息11a和流体信息11b。例如，存储单元11是纹线可视化设备10的存储器或存储装置。

结构信息11a指示分析空间1中的结构2的形状的时间变化。流体信息11b指示分析空间1中的存在流体的区域(流体区域3a、流体区域3b)中的多个点处的流体的速度的时间变化。虽然图1分开地示出了流体区域3a和3b，但是流体区域3a和流体区域3b可以在未示出的某处彼此相连接。存在右心和左心彼此相连接的情况，例如先天性心脏病。例如，结构2是心脏。在这种情况下，流体是心脏中的血液。当结构2是心脏时，例如通过对心脏的搏动和血液的冠状循环执行相互作用分析的模拟来生成结构信息11a和流体信息11b。

处理单元12基于存储在存储单元11中的结构信息11a和流体信息11b来计算纹线，并且显示所计算的纹线。例如，处理单元12是纹线可视化设备10的处理器。

当处理单元12在纹线可视化处理中基于第一分析时间点(t0)处的纹线9a(第一纹线)来计算第二分析时间点(t1)处的纹线9b(第二纹线)时，处理单元12执行以下处理。

处理单元12将包括分析空间1中的第一分析时间点(t0)处的纹线9a上的第一位置5处的离散点的部分区域设置为离散点的分析目标区域4(步骤S1)。分析目标区域4例如是在其中心具有第一位置5的球形区域。在这种情况下，例如，处理单元12计算第一分析时间点(t0)与第二分析时间点(t1)之间的时间段中的流体的最大速度，并且基于第一分析时间点(t0)与第二分析时间点(t1)之间的差和最大速度来计算分析目标区域4的半径。第一分析时间点(t0)与第二分析时间点(t1)之间的差值是用在纹线计算中的时间步长。

接下来，处理单元12基于分析目标区域4中的流体的速度来计算指示离散点上的粒子在第二分析时间点(t1)处的目的地的第二位置6(步骤S2)，其中，所述速度由流体信息指示。更具体地，第二位置6指示在第一分析时间点(t0)处已经存在于离散点处的虚拟粒子在被流体运送之后在第二分析时间点(t1)处存在于何处。

接下来，处理单元12基于与分析目标区域4中的结构2有关的信息来确定在第二分析时间点(t1)处由分析目标区域4中的结构2占据的区域(步骤S3)，其中，所述信息由结构信息指示。

接下来，处理单元12基于第一位置5和第二位置6来确定在第二分析时间点(t1)处纹线9b进入或者未进入占据区域(步骤S4)。例如，处理单元12执行以下确定：第二位置6是否已经落在流体的内部的第一确定；以及，连接第一位置5与第二位置6的线是否已经穿过结构2的表面的第二确定。接下来，处理单元12基于第一确定和第二确定的结果来确定纹线9b是否已经进入占据区域。例如，当第二位置6已经落在流体的内部并且当线未穿过结构2的表面时，处理单元12确定纹线9b未进入占据区域。

即使当第二位置6已经落在流体的内部时，如果线已经穿过结构2的表面，则处理单元12确定纹线9b已经进入占据区域。例如，在图1的示例中，假设第二位置6已经落在流体区域3a的内部，尽管第二位置6存在于流体的内部，然而纹线已经穿过结构2的存在于流体区域3b与流体区域3a之间的单元。因此，在这种情况下，处理单元12确定纹线9b已经进入结构2的占据区域。

此外，当第二位置6已经落在流体的外部并且当线没有穿过结构2的表面时，处理单元12确定纹线9b未进入占据区域。然而，在这种情况下，由于纹线9b已经落在流体的外部，因此处理单元12结束计算，而不显示纹线9b。

当处理单元12确定离散点未进入占据区域时，处理单元12显示通过第二位置6的纹线9b(步骤S5)。在图1的示例中，从粒子生成源8生成的第一分析时间点(t0)处的纹线9a和第二分析时间点(t1)处的纹线9b叠加在显示屏7上的结构2的截面上。

当处理单元12确定离散点已经进入占据区域时，处理单元12在第一分析时间点(t0)与第二分析时间点(t1)之间的时间段中设置至少一个第三分析时间点，并且计算指示离散点上的粒子在第三分析时间点处的目的地的第三位置。接下来，处理单元12基于第三位置来重新计算第二位置6。当处理单元12确定离散点未进入占据区域——作为重新计算的结果时，处理单元12显示纹线9a和纹线9b。

通过针对预定时间内的各个分析时间点如上所述地计算并显示纹线，形状随时间变化的结构2周围的流体的纹线9a和9b被精确地显示在显示屏7上。即，纹线可视化设备10不太频繁地显示错误纹线，例如已经进入随时间变化的结构2的纹线。换言之，纹线可视化设备10显示精确的纹线。

此外，计算纹线时要进行分析的结构信息11a和流体信息11b对应于分析目标区域4中的信息。因此，处理单元12能够高效地执行其处理。由于能够如此高效地计算纹线，所以纹线可视化设备10能够鉴于结构2的形状的变化来执行计算。此外，纹线可视化设备10能够容易地将其形状随时间变化的结构2周围的流体的纹线可视化。

当设置分析目标区域4时，例如，处理单元12可以基于存在流体的区域中的指示流体速度的多个点之间的间隔来设置分析目标区域4的半径的最小值。在这种情况下，当计算的半径小于最小值时，处理单元12将最小值设置为分析目标区域4的半径。例如，处理单元12可以将多个点中的两个相邻点之间的最大间隔值设置成作为分析目标区域4的半径的最小值。如果分析目标区域4的半径被设置得过小，则不能获得关于分析目标区域4中的流体的速度的信息和用于第二位置6的计算的信息。但是，以上述方式，这种情况不常发生。

处理单元12可以允许以下情况：在该情况中，分析目标区域4的半径小于通过将第一分析时间点(t0)与第二分析时间点(t1)之间的差乘以在第一分析时间点(t0)与第二分析时间点(t1)之间的时间段中的流体的最大速度而计算出的值。相乘结果是存在于与流体中的最大速度对应的位置处的离散点的移动距离(最大移动距离)。最大移动距离可以用作分析目标区域4的最大半径值。例如，处理单元12将分析目标区域4的最小半径值设置成小于最大移动距离。在这种情况下，离散点上的粒子的目的地可能落在分析目标区域4的外部。因此，当离散点上的粒子的目的地落在分析目标区域4的外部——作为第二位置6的计算的结果时，处理单元12通过扩展分析目标区域4来重新计算第二位置6。例如，处理单元12将分析目标区域4重新计算为具有与最大移动距离对应的半径的球形区域。

如上所述，由于处理单元12允许分析目标区域4的半径比离散点的最大移动距离小的情况，因此分析目标区域4的半径减小。由于要分析的信息量随着分析目标区域4的半径减小而减少，因此处理效率提高。然而，如果对在具有确定的半径的分析目标区域4中的第二位置6的计算失败，则处理单元12需要通过将分析目标区域4的半径设置成最大值(最大移动距离)来再次执行计算，这是额外的处理。因此，处理单元12可以计算使计算量最小化的半径，其中，所述计算量包括用于直到计算失败而重新计算第二位置6的计算量。当所计算的半径小于最小值时，处理单元12可以将最小值设置为分析目标区域4的半径。在这种情况下，例如，处理单元12基于离散点上的粒子的目的地落在分析目标区域4的外部的概率来设置达到最高处理效率的半径。即，处理单元12基于以下处理量来设置分析目标区域4的半径使得达到最小的处理量：当离散点上的粒子的目的地落在分析目标区域4的外部时增加的处理量；以及当对分析目标区域4设置较小的半径时减少的处理量。以这种方式，提高了处理效率。

此外，由于分析目标区域4的半径减小，因此分析目标区域4的最小半径值受到低质量网格——即有限元模型的质量——限制的情况不常发生。因此，由于存在低质量网格而导致的计算量的增加不常发生。

在以上描述中，执行了第一确定和第二确定二者，作为确定纹线9b是否已经进入占据区域的示例。然而，简单地，可以仅执行第二确定。即，通过仅执行确定连接第一位置5与第二位置6的线是否已经穿过结构2的表面的第二确定，处理单元12就能够确定纹线9b是否已经进入结构2的占据区域。在这种情况下，如果线已经穿过结构2的表面至少一次，则处理单元12确定纹线9b已经进入结构2的占据区域。例如，当流体区域3a和流体区域3b是封闭区域时并且当不存在流体通过其流出的出口例如心脏动脉时，处理单元12能够仅通过执行第二确定来精确地确定纹线9b是否已经进入占据区域。

[第二实施方式]

接下来，将描述第二实施方式。第二实施方式提供了一种能够随着心脏的运动将心脏中的血流的纹线可视化的可视化设备。

例如，使用计算流体分析使得可以模拟流体的行为，甚至是技术上或伦理上难以测量的系统中的流体，例如心脏中的血流传输。因此，计算流体分析被用于讨论其中心脏中的血流传输存在故障的先天性心脏病等的治疗。即，计算流体分析是重要的技术。通过使用可视化设备来将这样的计算流体分析的结果可视化，卫生保健专业人员如医生能够容易地理解分析结果并制定治疗计划。

图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)以及可记录CD(CD-R)/可重写CD(CD-RW)。

装置连接接口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通过使用插值方法对场进行插值。

2-1：由于1-3中的功能的应用增加了计算量，因此可视化设备100估计纹线上的点移动的最大距离，并且仅估计关于在具有等于最大距离的半径的预测球体内部的场的信息。以这种方式，不管数据容量如何，可视化设备100均保持一定的计算量。

2-2：通过分割时间步长，可视化设备100减小预测球体的半径，需要比不使用预测球体时所需的计算量更小的计算量，并且同时提高精确度。

3-1：由于大多数计算是在高质量网格上执行的，因此可视化设备100通过假设所有网格都是高质量网格来执行推测性计算。在计算失败的情况下，可视化设备100执行精确的计算。以这种方式，计算量减小。在该推测性计算中，通过允许纹线上的点的目的地落在预测球体的外部的可能性，可视化设备100减小预测球体的半径。计算失败的情况是纹线上的点的目的地不存在于预测球体内的情况。

3-2：由于可视化设备100在3-1中执行推测性计算，因此可视化设备100准备概率模型并且确定实现最小计算量的参数集，其中，所述最小计算量包括计算失败时所需的处罚。

通过在可视化设备100上实现这些功能获得了以下有利效果。

1.即使当心肌(弹性体)变形较大时，可视化设备100也能够在将心肌(弹性体)的运动考虑在内的同时来计算纹线。

2.可视化设备100能够通过使用预测球体来快速且精确地计算在纹线上的各个点。

3.可视化设备100能够通过使用概率模型来设置使计算成本最小化的预测球体的半径。

在下文中，将详细描述可视化设备100的功能。

图5是示出可视化设备100的功能的框图。可视化设备100包括：存储单元110，其保持从心脏模拟器200获取的模拟结果。例如，当心脏模拟器200执行计算流体动力学模拟时，将与在L个时间点t0、t1、...、tL(L是1或更大的整数)处的动态变化的弹性体和流场有关的模拟结果存储在文件中。例如，关于心肌的信息和关于血流的信息被作为不同的文件存储在存储单元110中。在图5的示例中，关于每个时间点的心肌的信息被存储为弹性体信息文件组111，并且关于每个时间点的血流的信息被存储为流体信息文件组112。

通过分析这些模拟结果，可视化设备100计算描述关于血流的传输的信息的纹线。作为模拟结果输出的时间间隔Δti＝ti+1-ti不需要与心脏模拟器200求解微分方程时所使用的时间间隔相匹配。为了减少信息量，通常只输出一些模拟结果。因此，为了精确地获得纹线，可视化设备100通过使用多个时间点处的输出文件来使用插值方法等，并且估计目标时间点处的各个物理量。

接下来，将描述可视化设备100的处理功能。可视化设备100包括信息读取单元120、纹线计算单元130、显示处理单元140以及分析单元150。

信息读取单元120从存储单元110读取指示流体分析结果的文件。纹线计算单元130通过使用所读取的信息来计算纹线。显示处理单元140将获得的结果可视化。

分析单元150是由信息读取单元120、纹线计算单元130以及显示处理单元140通用的一组功能。当执行特定分析处理时，信息读取单元120、纹线计算单元130以及显示处理单元140请求分析单元150执行处理并获得结果。

分析单元150包括流体信息分析单元151、弹性体信息分析单元152、预测球体半径计算单元153、预测球体内部信息计算单元154、移动边界碰撞确定单元155以及显示格式确定单元156。流体信息分析单元151分析流体的速度场、离散点的位置以及边界表面。弹性体信息分析单元152分析非流体的诸如心肌的弹性体的离散点的位置和边界表面。预测球体半径计算单元153设置用于在计算纹线时提高计算速度和计算精确度的预测球体的半径。预测球体内部信息计算单元154例如计算预测球体内部的速度场和心肌位置。移动边界碰撞确定单元155确定纹线上的点是否已经进入心肌——作为计算错误的结果。显示格式确定单元156确定如何显示获得的纹线。

例如，图5所示的各个单元的功能可以通过使计算机执行与相应单元对应的程序模块来实现。

接下来，将详细描述作为模拟结果而获得的信息。

图6示出弹性体信息文件组111的示例。弹性体信息文件组111是一组每个模拟时间点处的弹性体信息文件111a、111b等。弹性体信息文件111a、111b等中的每一个被给予文件名如“stru(X).inp”。在这种情况下，各个文件名中的“X”表示数字，并且这些数字根据模拟时间点的时间顺序以升序给出。

弹性体信息文件111a、111b等包括指示心脏在相应时间点处的形状的心肌数据。心肌数据包括沿着各个网格(布置在3D空间中的顶点)的x轴、y轴和z轴的坐标值、指示包括在心肌中的四面体单元(TETRA)的四个顶点的各个网格ID和应用于各个单元的力。

图7示出流体信息文件组112的示例。流体信息文件组112是一组每个模拟时间点处的流体信息文件112a、112b等。例如，流体信息文件112a、112b等中的每一个被给予文件名如“flui(Y).inp”。在这种情况下，各个文件名中的“Y”表示数字，并且这些数字根据模拟时间点的时间顺序以升序给出。

流体信息文件112a、112b等包括指示相应时间点处的血流的血流数据。血流数据包括沿着各个网格(布置在3D空间中的顶点)的x轴、y轴和z轴的坐标值，指示包括在血管中的四面体单元(TETRA)的四个顶点的各个网格ID以及指示在各个单元中流动的血液的方向和速度的各个速度场向量。

可视化设备100基于图6和图7所示的模拟结果来计算纹线并将纹线可视化。在下文中，将详细描述纹线可视化处理。

图8是示出纹线可视化处理的过程的示例的流程图。在下文中，将逐步地描述图8所示的处理。

首先，在步骤S101至S104中，纹线计算单元130执行各个纹线的初始设置。

[步骤S101]纹线计算单元130设置要被计算的纹线的数目M(M是1或更大的整数)。例如，纹线计算单元130将由用户输入的值设置为纹线的数目M。

[步骤S102]纹线计算单元130设置纹线计算的数目N(N是1或更大的整数)。在下文中，纹线计算的数目N将被称为“时间演化数目”。例如，纹线计算单元130将由用户输入的值设置为时间演化数目N。

由于纹线随时间变化，因此通过设置时间演化数目N来确定一系列时间点t0、t1、...、和tN，在这些时间点处纹线的结果被输出。当指定的时间演化数目N比文件的数目L大(L是1或更大的整数)时，纹线计算单元130可以将第(L+1)文件作为第二心动周期中的跳动，并且将时间点t0处的文件用于第(L+1)文件。例如，所述一系列时间点中的时间间隔被设置成0.01秒。然而，可替选地，所述一系列时间点可以具有不规则的时间间隔。

[步骤S103]纹线计算单元130设置纹线生成源的坐标。例如，纹线计算单元130将分析空间中由用户指定的点设置为纹线生成源的坐标。例如，用户在参考心肌信息和血流信息的同时指定空间中的点。纹线计算单元130读取指定点的坐标作为坐标向量X0。当纹线的数目是1时，纹线生成源被设置为具有坐标向量X0。当纹线的数目是多个时，纹线计算单元130随机地将纹线lj的生成源设置在具有以坐标向量X0为其中心且具有半径r(r是正实数)的球体中。生成源是从血流中的坐标中选择的。纹线计算单元130将设置的生成源的坐标向量Xj设置为纹线生成源。

接下来，纹线计算单元130如下来执行第j(j＝1、2、...、M)纹线lj的初始设置。

第j纹线lj由与时间演化数目N匹配的离散点形成。因此，纹线计算单元130生成指示包括在纹线lj中的离散点的点Pij(i＝0、1、2、...、N)。纹线计算单元130将各个离散点的初始值的坐标设置为粒子生成源的坐标。

当计算在时间点t＝ti处的纹线时，各个点Pij(i＝0、1、2、...、N)经历时间演化计算，作为从相应的生成源发出的纹线粒子的位置。由于没有从任何生成源发出与点Pij(i＝i+1、...、N)相对应的纹线粒子，因此当计算时间点t＝ti处的纹线时，这些纹线粒子不会经历计算。此外，纹线计算单元130以值i的升序来计算离散点。由于从相应的粒子生成源发出，因此较早计算的离散点具有较长的时间。

[步骤S104]纹线计算单元130针对以下情况执行设置：纹线上的点已经落在大动脉等中，即落在目标系统中的流体边界的外部。纹线上的点Pij可能落在大动脉等中，即穿过流体边界落在系统的外部。在这种情况下，由于在系统的外部没有定义流体速度场，因此点Pij在下一时间点处的计算不能被执行。因此，纹线计算单元130为每个点Pij设置区域确定标记作为各个离散点的参数。当点Pij落在目标区域内时，区域确定标记指示“T”。相比之下，当点Pij通过大动脉等中的流体的流动而漂移并且落在目标区域的外部时，区域确定标记指示“F”。由于在初始设置中流体包括所有点Pij，因此纹线计算单元130将每个离散点的区域确定标记设置为“T”。

[步骤S105]纹线计算单元130从指数i＝1起以升序对每个指数i(i＝1、2、...、和N-1)重复步骤S106至步骤S110的组。

[步骤S106]纹线计算单元130从指数j＝1起以升序对每个指数j(j＝1、2、...、和M)重复步骤S107至步骤S109的组。

[步骤S107]纹线计算单元130设置作为时间演化的开始的时间点并将该时间点存储在存储器102中。在第i次计算中，计算开始时间点被设置为t＝ti。时间演化结束时间点被设置为t＝ti+1。

[步骤S108]纹线计算单元130在由ti≤t≤ti+1限定的时间点之间执行时间演化计算。基于时间演化计算，每个时间点t＝ti处从单独的纹线生成源发出的所有点Pij(i＝0、1、2、...和i)经历时间演化，并且线上的所有点被立刻更新。作为在时间点[ti，ti+1]处纹线lj上的各个点Pij的时间演化计算的结果，获取了坐标值，所述坐标值被设置为在下一时间点t＝ti+1处的坐标Pi+1,j。

[步骤S109]纹线计算单元130将所获取的计算结果存储在存储器中。基于计算结果，显示处理单元140将纹线可视化。此外，纹线计算单元130能够将纹线的坐标值输出至文件。

[步骤S110]每当纹线计算单元130执行步骤S107至步骤S109的组时，纹线计算单元130向指数j加1。在对指数j＝M执行步骤S107至S109之后，纹线计算单元130执行步骤S111。

[步骤S111]每当纹线计算单元130执行步骤S106至步骤S110的组时，纹线计算单元130向指数i加1。在对指数j＝N-1执行步骤S106至S110之后，纹线计算单元130结束纹线可视化处理。

接下来，将详细描述时间演化计算处理(步骤S108)。

图9和图10示出示出时间演化计算处理的过程的流程图。在下文中，将逐步地描述图9所示的处理。

[步骤S121]纹线计算单元130读取点Pij的区域确定标记。纹线计算单元130确定区域确定标记是否指示“T”。当区域确定标记指示“T”时，处理进行至步骤S122。当区域确定标记指示“F”时，纹线计算单元130确定该点已经落在区域的外部并且结束时间演化计算处理而无需执行计算。因此，当区域确定标记指示“F”时，纹线计算单元130不更新坐标值。

[步骤S122]纹线计算单元130读取点Pij的坐标值。

[步骤S123]纹线计算单元130从存储器102读取计算开始时间点t＝ti和计算结束时间点t＝ti+1。

[步骤S124]纹线计算单元130经由信息读取单元120从文件flui(i).inp读取流体部分在计算开始时间点t＝ti处的网格信息向量B(向量r，ti)和速度场向量v(向量r，ti)。此外，纹线计算单元130经由信息读取单元120从文件stru(i).inp读取向量M(向量r，ti)，向量M是关于心肌部分的弹性体的网格信息(关于心肌的结构的信息)。

[步骤S125]纹线计算单元130经由信息读取单元120从文件flui(i+1).inp读取流体部分在计算开始时间点t＝ti+1处的网格信息向量B(向量r，ti+1)和速度场向量v(向量r，ti+1)。此外，纹线计算单元130经由信息读取单元120从文件stru(i+1).inp读取向量M(向量r，ti+1)。

[步骤S126]纹线计算单元130计算由文件flui(i).inp和文件flui(i+1).inp指示的网格点的速度场向量的范数(速度场向量的长度)，并且通过使用插值表达式来计算相应时间区间中的速度的最大值。纹线计算单元130将所计算的最大值作为Vmax存储在存储器102中。

[步骤S127]如果纹线计算单元130在单个时间演化中计算区间[ti，ti+1]，则精确度不够。因此，纹线计算单元130由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]。纹线计算单元130自身将最优值设置为分割数目Ndiv。分割数目Ndiv可以从外部给出。接下来，处理进行至图10的步骤S131。

在下文中，将逐步地描述图10所示的处理。

[步骤S131]纹线计算单元130通过对从k＝1至k＝Ndiv的每个中间时间点(t＝tk(k＝1、2、...、和Ndiv))重复步骤S132至步骤S140的组来执行时间演化计算。因此，获得关于每个中间时间点t＝tk的坐标向量rk，并且获得时间点t＝ti+1处的点Pi+1,j的坐标向量rNdiv＝向量ri+1。

[步骤S132]纹线计算单元130计算关于在时间点t＝tk处的场的信息(流体结构信息向量B(向量r，tk)和速度场向量v(向量r，tk))。下面将参照图13来详细描述该步骤。

[步骤S133]纹线计算单元130计算关于在时间点t＝tk处的弹性体场的信息(心肌结构信息向量M(向量r，tk))。下面将参照图14来详细描述该步骤。

[步骤S134]纹线计算单元130仅通过时间dt＝tk+1-tk来执行时间演化，以计算向量rk+1。

纹线计算单元130可以对每个中间时间点如下来执行步骤S132至步骤S134中的时间演化。假设在时间点t＝tk时点Pi+1,j在向量rk处，则通过表达式(1)来表示纹线方程。

在表达式(1)中，向量v(向量r，t)是在时间点t和位置向量r处的速度(场)。因此，时间Δt之后的坐标通过数值求解作为常微分方程的表达式(1)来计算。通过用四阶龙格-库塔(Runge-Kutta)法求解表达式(1)得到以下计算表达式。

可以根据向量v(向量r，ti)、向量v(向量r，ti+1)、向量B(向量r，ti)、向量B(向量r，ti+1)、向量M(向量r，ti)以及关于在时间点t＝tk处的移动边界表面Sk的信息来计算时间点t＝tk时在任何位置向量r处的速度场向量v(向量r，tk)。在下文中，移动边界表面将被简称为“边界表面”。可以根据向量M(向量r，ti)来计算关于在时间点t＝tk处的边界表面Sk的信息。因此，计算表达式(3)至表达式(6)。通过将这些结果代入表达式(2)，纹线计算单元130计算向量rk+1。

[步骤S135]纹线计算单元130执行向量rk+1的位置是否已经落在心脏的内部的确定。纹线计算单元130执行该处理，因为所计算的向量rk+1包括无限时间宽度误差并且可以进入心肌。在作为子例程的位置确定之后，获取指示确定的位置的结果的状态。例如，当向量没有穿过或进入心肌时，获取作为正常终止状态的“0”或“2”。当向量rk+1指示分析目标流体中的单元中的位置(例如心房或心室中)时，获取状态“0”。当向量rk+1指示在分析目标流体中的单元的外部的位置(例如动脉中)时，获取状态“2”。下面将参照图17来详细描述该步骤。

[步骤S136]纹线计算单元130确定确定结果是否指示状态“0”，所述状态“0”指示正常终止。当确定结果指示状态“0”时，处理进行至步骤S140。当确定结果不指示状态“0”时，处理进行至步骤S137。

[步骤S137]纹线计算单元130确定确定结果是否指示状态“2”，所述状态“2”指示终止。获取状态“2”的情况对应于以下情况：在计算过程中点Pij已经穿过流体边界移动至系统外部的外部单元，例如移动至大动脉。当确定结果指示状态“2”时，处理进行至步骤S138。否则，处理进行至步骤S139。

[步骤S138]当点Pij已经落在系统的外部时，纹线计算单元130将纹线上的点Pij(j＝1至j)的区域确定标记设置为“F”。单独的区域确定标记“F”指示对应点Pij已经落在分析区域的外部。接下来，纹线计算单元130结束时间演化计算处理。以这种方式，当终止状态指示“2”时，纹线计算单元130确定点Pij已经落在系统的外部并将区域确定标记设置成“F”。通过连接由点Pi1、Pi2、Pi3、...、和PiN顺序连接而形成的曲线来绘制纹线。因此，当任何点Pij被确定为已经落在区域的外部时，没有理由绘制从粒子生成源发出的先前的点。因此，纹线计算单元130将点Pij(j＝1至J)的所有区域确定标记设置成“F”并且结束处理。

[步骤S139]当确定结果没有指示正常终止时(当状态既不是“0”也不是“2”时)，纹线计算单元130减小用作控制参数的时间步长。即，纹线计算单元130通过更多的时间点来进一步分割计算，并且重复步骤S132至步骤S137的组。在减小时间步长之后，当纹线计算单元130确定正常终止时，处理进行至步骤S140。

当向量rk+1已经落在预测球体的外部时，纹线计算单元130也确定没有正常终止。在这种情况下，例如，通过增加预测球体的半径并且执行重新计算，纹线计算单元130能够防止向量rk+1落在预测球体的外部。

[步骤S140]纹线计算单元130将向量rk+1存储在存储器中，并且将所存储的值设置为第(k+1)次重复计算的初始值。

[步骤S141]每当纹线计算单元130执行步骤S132至步骤S140的组时，纹线计算单元130向指数k加1并且重复该处理。当纹线计算单元130对所有中间时间点(k＝Ndiv)完成时间演化计算时，处理进行至步骤S142。

[步骤S142]纹线计算单元130将最终计算的向量rNdiv作为向量rk+1存储在存储器中。

通过图9和图10所示的处理来更新纹线上的点的坐标。

图11示出纹线的数据示例。如图11所示，在每个分析时间点设置纹线上的点的坐标值。将粒子生成源的坐标值设置为时间点t＝t0处的纹线上的点的初始值。当时间点t＝t1时，更新指示初始发出的粒子的位置的点的坐标值。接下来，随着时间点的更新，发出新的粒子，并且更新指示发出的新粒子和旧粒子的位置的点的坐标值。在图11中，位置已经从它们之前的时间点改变的点的坐标值被加下划线。

接下来，将详细描述用于计算在时间点t＝tk时的以下信息的处理(步骤S132和步骤S133)：关于流体部分的结构(网格坐标)和速度场的信息；以及关于弹性体(心肌)部分的结构(网格坐标)的信息。以下信息被用于表达式(3)至表达式(6)的计算：关于结构(网格坐标)的信息——即向量B(向量r，tk)；以及关于流体部分的速度场的信息——即向量v(向量r，tk)。关于弹性体的结构的信息被用于确定纹线是否已经穿过了心肌并且由向量M(向量r，tk)表示。

仅在时间点t＝ti和t＝ti+1处输出数据。因此，由于通过龙格-库塔法计算的网格坐标和速度场在中间时间tk处没有被定义，因此纹线计算单元130根据输出文件的速度场来计算场的近似值。该计算的关键考虑是在执行模拟时时刻地移动网格位置。尽管可以以任何方式确定网格位置，但是常常使用任意拉格朗日-欧拉(ALE)法来解决诸如心脏的对象的边界移动的问题。在ALE方法中，模拟中使用的坐标被独立地确定，从而不会使所描述的偏微分方程的解的精确度劣化。在许多情况下，偏微分方程被用于该确定。然而，用于确定网格位置的控制方程对数据分析的位置中的一个是不可用的，而只有给出的网格点的输出值是可用的。在这种情况下，网格点的位置随时间不断变化。

图12示出网格点的位置如何随时间不断变化。在图12中，横轴表示时间，并且纵轴表示网格点(顶点)的X坐标值。如图12所示，各个网格点的位置随时间不断变化。因此，鉴于该事实，纹线计算单元130通过使用插值方法来估计在任何时间点处的网格位置。

图13是示出用于计算在时间点t＝tk时的场信息的处理的过程的示例的流程图。在下文中，将逐步描述图13所示的处理。

[步骤S151]纹线计算单元130将在时间点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)、向量B(向量r，ti+1)以及向量v(向量r，ti+1)是已知的。这些向量值已经在步骤S124和步骤S125中从文件读取。

首先，纹线计算单元130对限定流体部分的结构的网格点执行以下处理。然而，为了减少计算量，纹线计算单元130在理论上计算点Pi+1,j的最大移动距离，并且仅对具有等于最大移动距离的半径R的球体内部的网格点执行处理。在下文中，该球体将被称为“预测球体”。

[步骤S152]纹线计算单元130设置预测球体的半径R。以下将参照图20来详细描述该处理。

[步骤S153]纹线计算单元130搜索在预测球体的半径R内部的流体网格点，并将预测球体内部的网格点的数目设置为NB,elem。

[步骤S154]纹线计算单元130对预测球体的半径R内部的NB,elem个网格点中的每个网格点执行步骤S155至步骤S160的组。

[步骤S155]纹线计算单元130确定网格点i是否在边界表面上。如果网格点i在边界表面上，则处理进行至步骤S159。如果网格点i不在边界表面上，则处理进行至步骤S156。

[步骤S156]当网格点不在边界表面上时，纹线计算单元130计算网格点i的平均移动速度向量vave,i。纹线计算单元130通过对时间点tk处的各个网格点的坐标进行插值来执行该计算。在这种情况下，由于输出各个网格点的时间短，因此计算由一阶方程来近似。因此，由以下表达式来计算具有特定ID的网格点的在时间点t处的坐标向量rID(t)的平均移动速度向量vave。

[步骤S157]通过使用平均移动速度向量vave，纹线计算单元130计算网格点i的位置向量rID(t)。位置向量rID(t)通过以下表达式来计算。

[步骤S158]纹线计算单元130计算非边界表面上的网格点的速度场向量V(向量rID(t)，t)。纹线计算单元130通过使用以下表达式(9)来执行该计算。

[步骤S159]纹线计算单元130通过使用以下表达式(10)来计算指示边界表面上的网格位置的向量rID(t)。

[步骤S160]纹线计算单元130通过使用以下表达式(11)来计算指示边界表面上的网格位置的向量rID(t)的速度场向量V(向量rID(t)，t)。

[步骤S161]每当纹线计算单元130执行步骤S155至步骤S160的组时，纹线计算单元130向指数i加1，并且重复步骤S155至步骤S160的组。当纹线计算单元130已经完成对指数i＝NB,elem的计算时，纹线计算单元130已经完成预测球体内部的所有网格点的速度场的计算。

[步骤S162]纹线计算单元130计算包括要被计算的点向量r(t)的单元号码ID0。

[步骤S163]纹线计算单元130通过使用以下等式(12)根据包括在单元号码ID0中的网格点信息rID0(t)来计算在要被计算的点向量r(t)处的速度场向量V(向量r(t)，t)。

接下来，纹线计算单元130结束用于计算当时间点t＝tk时的场信息的处理。

在图13所示的处理中，除了在步骤S157中计算非边界表面上的位置向量rID(t)之外，在步骤S159中计算边界表面上的网格位置向量rID(t)。这是因为，对于边界表面S(t)，无滑移边界条件被在纳维叶-斯托克斯(Navier-Stokes)方程中定义，并用作约束条件。因此，纹线计算单元130确定要被计算的网格点是否在边界表面S(t)上。如果这样，则纹线计算单元130确定满足边界条件的网格位置向量。因此，纹线计算单元130执行以下计算。

假设由边界表面S(t)上的网格点位置和对应的速度场形成的拓扑空间中的点由(向量rID(t)，向量vID(t))表示，则在时间点ti和ti+1的值已经被输出。因此，(向量rID(ti)，向量vID(ti))和(向量rID(ti+1)，向量vID(ti+1))被确定。由于在关于粘性流体的纳维叶-斯托克斯方程中设置了无滑移边界条件，因此网格位置向量rID(t)与速度场向量vID(t)之间建立了以下关系。

因此，由于通过用基于时间的三阶方程定义向量rID(t)而总共存在四个关于边界表面S(t)上的网格点的数据，因此条件得到满足。因此，获得表达式(10)。

关于系数向量ai，假设向量rID(t)的分量由ξID(ξ＝x，y，z)来表示，并且相应的微分分量由vξ(ξ＝x，y，z)来表示，则获得以下表达式。

a1ξ＝vξ(ti) (16)

a0ξ＝ξID(ti) (17)

在这些表达式中，使用了以下关系。

Δti＝ti+1-ti (18)

ΔξID＝ξID(ti+1)-ξID(ti) (19)

Δvξ＝vξ(ti+1)-vξ(ti) (20)

以这种方式，估计在任何时间点处的边界表面S(t)上的位置向量rID(t)的组(图13的步骤S159)。因此，通过计算流体中所有网格点的位置，纹线计算单元130能够确定作为关于流体的结构的信息的向量B(向量r，tk)。

接下来，将描述纹线计算单元130如何计算时间点处的速度场。速度场可以被分割成X分量、Y分量和Z分量。这些分量中的每一个都是标量场。因此，通过以下操作得到下面的表达式：将速度场表示为Vj(向量r，t)(j＝X分量、Y分量、Z分量)；执行泰勒(Talyor)展开，直到点周围的一阶项(向量rki，ti)，其足够接近位置向量r并且其在空间时间中是作为输出文件的给出数据；并且忽略二阶项和后续项。

在表达式21中，指数k表示网格ID，并且指数i表示输出文件所存在的第i时间点。由于在有限元方法中，▽Vj是各个单元的顶点处的速度梯度，因此▽Vj通过下面的表达式来获得。

在表达式22中，Nk被称为结构函数并且被给出为多项式，并且Vjk是表示各个顶点处的速度场的值并且已经被输出至文件。因此，通过基于差来估计表达式(21)中的速度场中的时间导数项纹线计算单元130能够估计非边界表面上的速度场。即，纹线计算单元130能够通过使用上述表达式(9)来估计速度场(图13的步骤S160)。如同场信息的情况，对于作为心肌结构信息的向量M(向量r，t)，纹线计算单元130能够计算任何时间点处的网格点信息。

图14是示出用于获得在时间点t＝tk时的弹性体场信息的处理的过程的示例的流程图。在下文中，将逐步描述图14所示的处理。

[步骤S171]纹线计算单元130将读取的心肌结构信息设置在存储器中。所设置的信息包括向量B(向量r，ti)、向量v(向量r，ti)、向量M(向量r，ti)、向量B(向量r，ti+1)、向量v(向量r，ti+1)以及向量M(向量r，ti+1)。

[步骤S172]纹线计算单元130基于预测球体的半径R来搜索预测球体内部的心肌的结构，并且搜索心肌的网格点。预测球体的半径R是与在流体部分上执行的处理中使用的半径R相同的半径R。纹线计算单元130将检测到的网格点的数目设置为NM,elem。

[步骤S173]纹线计算单元130对预测球体的半径R内部的NM,elem个网格点中的每个网格点执行步骤S174至步骤S178的组。

[步骤S174]纹线计算单元130确定网格点i是否在边界表面S(t)上。如果网格点i在边界表面S(t)上，则处理进行至步骤S177。如果否，则处理进行至步骤S175。

[步骤S175]当网格点在非边界表面S(t)上时，如同流体部分的情况，纹线计算单元130通过使用表达式(7)来计算各个网格点i的平均移动速度向量vave,i。

[步骤S176]纹线计算单元130通过使用表达式(8)来计算网格点i的位置向量rID(t)。

[步骤S177]当网格点i在边界表面S(t)上时，由于网格位置与流体部分交叠，因此纹线计算单元130在流体部分搜索对应的网格。

[步骤S178]纹线计算单元130将发现的网格坐标设置为弹性体(心肌)在边界表面S(t)上的坐标。然而，该处理基于以下假设：流体部分在边界表面S(t)上的各个点首先被计算。

[步骤S179]每当纹线计算单元130执行步骤S174至步骤S178的组时，纹线计算单元130向t加1并且重复步骤S174至S178的组。当纹线计算单元130完成对指数i＝NM,elem的计算时，纹线计算单元130结束用于计算在时间点t＝tk时的弹性体结构信息的处理。以这种方式，纹线计算单元130能够计算心肌结构信息向量(向量r，t)。

接下来，将详细描述心肌内部的位置的确定。

在下文中，将描述确定心脏内部的作为时间演化的结果而获得的向量rk+1的位置的过程。由于龙格-库塔法等包括有限误差，因此作为计算结果而获得的位置向量rk+1可能落在不真实的位置上。

图15示出能够作为计算结果而获得的位置的示例。心脏40的右心系统41和左心系统42的内部部分是存在分析目标流体的区域。右心系统41包括右心房和右心室。在第二实施方式中，左心系统42包括左心房和左心室。连接至左心系统42的动脉43的内部的血流不是分析目标。然而，大动脉的一部分可以包括在模拟中。

当点Pkj在时间点t＝tk处存在于流体中时，如图15所示，在通过时间演化实现的移动之后存在五个可能的目的地。每个目的地都被给出对应的状态变量(状态)。当点没有移动越过心肌壁而且已经落在分析目标流体的单元内时，状态变量表示“0”。当点已经移动越过心肌壁而且已经落在分析目标流体的单元的外部时，状态变量表示“1”。当点没有移动越过心肌壁而且已经落在分析目标流体的单元的外部时，状态变量表示“2”。当点已经移动越过心肌壁而且已经落在分析目标流体的单元的内部时，状态变量表示“3”。当点没有移动越过心肌壁而且落在心肌与分析目标流体的单元之间的边界上时，状态变量表示“4”。

图16是指示状态变量的真值表。以下描述假设“P”表示时间演化之前的点，其始终存在于流体中，而“Q”表示点P在时间演化之后的目的点。此外，以下描述假设纹线计算单元130仅关于点Q执行确定。纹线计算单元130能够通过执行两种确定处理来设置点Q的状态变量。

第一确定处理是流体确定，其中，纹线计算单元130确定目的点Q是否存在于流体中。如果点Q存在于流体中，则确定T。如果否，则确定F。

第二确定是线确定，其中，纹线计算单元130确定通过连接初始点P与目的点Q而形成的线PQ是否穿过心肌或心肌表面。如果线PQ穿过心肌(表面)，则确定T。如果否，则确定F。在线确定中，交点的数目也被确定。

当点Q存在于流体中而且线PQ没有穿过心肌时，纹线计算单元130确定正常移动并将“0”设置为状态变量。

当点Q不存在于流体中而且线PQ穿过心肌(表面)时，以下两种情况是可能的：(1)点Q已经移过越过心肌而且落在系统的外部；以及(2)点Q已经嵌入在心肌中。在任一情况下，由于需要执行重新计算，因此纹线计算单元130将“1”设置为状态变量。

当点Q不存在于流体中而且线PQ没有穿过心肌(表面)时，纹线计算单元130确定点Q经由大动脉等已经落在模拟系统的外部，并且将“2”设置为状态变量。

即使当点Q存在于流体中时，也存在确定出诸如从左心房至右心房的移动的不可能的移动的情况。在这种情况下，当T被确定为流体确定时，T也被确定为线确定，从而交点的数目总是多个。因此，当交点的数目是2或更大时，纹线计算单元130将“3”设置为状态变量。

当T被确定为流体确定并且确定为线确定时，如果交点的数目是1，则纹线计算单元130确定点已经落在流体与心肌的边界表面上。因此，纹线计算单元130将“4”设置为状态变量。

接下来，将描述以上状态确定处理的过程。

图17是示出状态确定处理的过程的示例的流程图。在下文中，将逐步描述图17所示的处理。

[步骤S181]纹线计算单元130确定时间演化之后的坐标向量rk+1是否存在于流体中。将参照图18来详细描述该处理。

[步骤S182]纹线计算单元130确定由移动向量dr＝向量rk+1-向量rk形成的具有无限长度的线是否穿过心肌，并且对该线与心肌(表面)的交点的数目进行计数。将参照图19来详细描述该处理。

[步骤S183]纹线计算单元130基于图16所示的真值表来确定状态。即，步骤S181和S182的每个结果被获得作为真值，即真(T)或假(F)。步骤S183的结果被获得作为0或更大的整数。纹线计算单元130参考真值表，并且确定值“0”至“4”中的任何一个作为与步骤S181至步骤S183的结果而获得的返回值对应的心脏的状态变量。纹线计算单元130使用确定结果作为位置确定结果。

以这种方式，执行位置确定，并且确定状态变量。

接下来，将详细描述用于确定时间演化后位置是否落在流体内的处理(步骤S181)。

图18是示出用于确定时间演化后位置是否落在流体内的处理的过程的示例的流程图。在图18的示例中，纹线计算单元130确定坐标向量rk+1是否存在于流体中。在下文中，将逐步描述图18所示的处理。

[步骤S191]纹线计算单元130获取时间点tk处的坐标向量rk。

[步骤S192]纹线计算单元130获取时间点tk+1处的坐标向量rk+1。

[步骤S193]纹线计算单元130获取预测球体的半径R。

[步骤S194]纹线计算单元130获取具有以坐标向量rk为其中心且具有半径R的预测球体的内部的单元的列表(流体单元列表Lf)。

[步骤S195]纹线计算单元130将sizef设置为流体单元列表Lf中的单元的数目。

[步骤S196]纹线计算单元130将“F(假)”设置为结果的初始值。

[步骤S197]纹线计算单元130对单元列表(i＝1、2、...、和sizef)中的各个单元Li重复步骤S198。

[步骤S198]纹线计算单元130确定单元列表中的单元Li是否包括坐标向量rk+1。如果是，则处理进行至步骤S200。如果否，则纹线计算单元130进行至步骤S199。

[步骤S199]每当纹线计算单元130执行步骤S198时，纹线计算单元130向指数i加1并且重复步骤S198。当纹线计算单元130完成对i＝sizef的处理时，纹线计算单元130结束用于确定时间演化后位置是否落在流体内的处理。

[步骤S200]纹线计算单元130将单元Li的单元号码ID存储在存储器中。

[步骤S201]纹线计算单元130将结果改变成“T(真)”。

以这种方式，当目的点的坐标包括在单元中的任何单元中时，纹线计算单元130确定坐标向量rk+1存在于流体中，并且将返回值设置成T(真)。然而，当不包括目的点的坐标时，纹线计算单元130将返回值设置成F(假)。

接下来，将详细描述用于搜索预测球体中的移动向量dr通过的弹性体单元的处理(步骤S182)。

图19是示出用于搜索移动向量dr通过的弹性体单元的处理的过程的示例的流程图。在下文中，将逐步描述图19所示的处理。

[步骤S211]纹线计算单元130获取时间点tk处的坐标向量rk。

[步骤S212]纹线计算单元130获取时间点tk+1处的坐标向量rk+1。

[步骤S213]纹线计算单元130根据所获取的信息来计算移动向量dr＝向量rk+1-向量rk。以这种方式，定义了指示点的移动路径的线。

[步骤S214]纹线计算单元130获取预测球体的半径R。

[步骤S215]纹线计算单元130生成在预测球体的半径内的弹性体单元的列表Le。

[步骤S216]纹线计算单元130将sizee设置为弹性体单元列表Le中的单元数目。

[步骤S217]纹线计算单元130将“F(假)”设置为结果的初始值。

[步骤S218]纹线计算单元130对单元列表(i＝1、2、...、和sizee)中的各个单元Li重复步骤S219。

[步骤S219]纹线计算单元130确定第i单元Li是否与移动向量dr相交。由于单元Li是多面体，因此纹线计算单元130确定移动向量dr是否与单元Li的任何表面相交。如果纹线计算单元130确定在单元Li的任何表面上没有交点，则纹线计算单元130确定移动向量dr没有与单元相交。在这种情况下，处理进行至步骤S220。如果纹线计算单元130确定移动向量dr与单元相交，则处理进行至步骤S221。

[步骤S220]每当纹线计算单元130执行步骤S219时，纹线计算单元130向指数i加1并且重复步骤S219。当纹线计算单元130完成对i＝sizee的处理时，纹线计算单元130结束用于搜索移动向量dr通过的弹性体单元的处理。

[步骤S221]纹线计算单元130将单元Li的单元号码ID存储在存储器中。

[步骤S222]纹线计算单元130获得指示线的移动向量dr与心肌表面的交点的数目Z(Z是1或更大的整数)，并且将数目Z存储在存储器102中。

[步骤S223]纹线计算单元130将结果改变成真。

如上所述，当不存在交点时，纹线计算单元130确定移动向量dr没有与单元相交并且将返回值设置成F。相比之下，当存在至少一个交点时，纹线计算单元130确定移动向量dr与单元相交并且将返回值设置成T。在这种情况下，存储与移动向量dr具有至少一个交点的单元的单元号码和与心肌表面的交点的数目Z。

通过执行以上处理，纹线计算单元130能够适当地确定点的目的地的状态。当纹线计算单元130在该处理中确定目的地时，通过检查预测球体中的单元，纹线计算单元130能够更高效地执行处理。接下来，将详细描述用于设置预测球体的半径R的方法。

纹线计算单元130基于时间演化之前的坐标向量rk能够在时间步长Δt内移动多长来设置预测球体。在四阶龙格-库塔法的情况下，根据表达式(2)建立以下不等式。

因此，通过如表达式(24)所指示地定义半径R，半径R表示在时间步长Δt内纹线上的点P移动的最大距离。

此外，点P在时间步长Δt之后必然存在于球体的内部。此外，中间向量vi也是在具有半径R的球体内部存在的点，其将如下被指示。即，根据表达式(4)和表达式(5)来建立表达式(25)。

因此，由中间向量指示的坐标表示具有半径R的球体内部的点。通过以表达式(3)和(6)执行同样的操作，可以看出，所有中间向量vi也是具有半径R的球体内部的点。该半径R被设置为预测球体的半径。纹线计算单元130根据流体的速度场来计算最大值，预测球体半径是针对所述流体的速度场来设置的。简单地，最大值可以从在其上执行模拟的整个区域中的网格点获得。可替选地，在将模拟目标区域分割为若干适当区域之后，最大值可以从局部信息中获得。

由于从纳维叶-斯托克斯方程的角度来看，速度的最大值没有落在边界表面S(t)上，因此通过表达式(9)来对中间时间点处的速度的最大值进行插值。由于表达式(9)是相对于时间的二阶方程，因此根据通过两个时间点t＝ti和t＝ti+1处的两个点的抛物线的最大值来计算速度的最大值，其中，所述两个点的数据已经被输出。近似地，由于泰勒展开，表达式(9)中的二阶项变成二阶微量。因此，可以仅考虑一阶项，并且可以将在相应时间点t＝ti和t＝ti+1处的速度场的最大值彼此进行比较。在这种情况下，较大的最大值被用作速度场的最大值。因此，纹线计算单元130将具有计算的速度场的最大范数的点设置为最大速度|向量vmax|。接下来，当设置龙格-库塔法的时间步长Δt时，最大移动距离被设置为R＝|Δt vmax|。

然而，如果时间步长Δt太小，则如上所述设置的预测球体半径R变得小于网格宽度(相邻网格点之间的距离)。因此，在网格点的离散化中存在最小值Rmin。例如，该最小值Rmin的初始值被设置成0.001[m]作为经验最小值。最小值可以根据网格的统计分析来计算。

当预测球体的半径R变得小于计算的Rmin时，预测球体的半径R被设置成Rmin。以这种方式，避免了预测球体内部不存在单元的情况。图20示出如上所述的用于设置预测球体半径R的处理的过程。

图20是示出用于设置预测球体半径的处理的过程的示例的流程图。在下文中，将逐步描述图20所示的处理。

[步骤S231]纹线计算单元130根据单元的边长和速度场的统计分析来设置预测球体的半径的最小值Rmin。

[步骤S232]纹线计算单元130获取时间点ti处的速度向量vi的值。

[步骤S233]纹线计算单元130获取时间点ti处的速度的最大值|向量vi,max|。

[步骤S234]纹线计算单元130获取时间点ti+1处的速度的最大值|向量vi+1,max|。

[步骤S235]纹线计算单元130确定|向量vi,max|是否等于或大于|向量vi+1,max|。当|向量vi,max|等于或大于|向量vi+1,max|时，处理进行至步骤S236。当|向量vi,max|小于|向量vi+1,max|时，处理进行至步骤S237。

[步骤S236]纹线计算单元130将|向量vi,max|设置成|向量vmax|。接下来，处理进行至步骤S238。

[步骤S237]纹线计算单元130将|向量vi+1,max|设置成|向量vmax|。接下来，处理进行至步骤S238。

[步骤S238]纹线计算单元130获取龙格-库塔法中的时间步长dt。

[步骤S239]纹线计算单元130计算|向量vmax|dt并且将计算的结果设置为预测球体半径R。

[步骤S240]纹线计算单元130确定预测球体半径R是否小于预测球体半径的最小值Rmin。当预测球体半径R小于最小值Rmin时，处理进行至步骤S241。当预测球体半径R等于或大于最小值Rmin时，纹线计算单元130结束本处理。

[步骤S241]纹线计算单元130将最小值Rmin设置为预测球体半径R。接下来，纹线计算单元130结束本处理。

通过执行以上处理，纹线计算单元130能够针对预测球体设置适当的半径R。

在第二实施方式中，可以通过执行以下处理来更快地执行处理。

[通过分割方法提高计算精确度和速度]

在图9的步骤S127中，在其处给出输出文件的时间区间ti≤t≤ti+1被分成Ndiv个时间点。将详细描述如何分割时间区间。

在纹线计算中，在所有单元中只有少量单元与单个纹线上的点有关。因此，为了减小存储容量和计算时间，使用预测球体。纹线的计算成本与预测球体的半径R³成比例地增加。这将被如下进行解释。

假设空间单元号码的密度是作为坐标向量r的函数的ρ(向量r)，则作为计算目标的单元的数目Nelem由以下表达式给出。

假设密度ρ(向量r)近似均匀，则ρ(向量r)＝ρ0。因此，Nelem∝R³。

当时间区间ti≤t≤ti+1被分割成Ndiv个时间点时，基于四阶龙格-库塔法的时间演化数为Ndiv。同时，通过表达式(24)来确定预测球体半径。因此，当时间步长为1/Ndiv时，预测球体半径也变成1/Ndiv。由于单一计算量与预测球体的半径R³成比例，因此计算量由Ndiv^-3表示。由于该计算被重复Ndiv次，因此总计算量由Ndiv^-2来表示。图21示意性地示出该计算的构思。

图21示出了通过时间分割而实现计算量的减少的构思。图21假设计算作为点P的目的地的点Q的情况。当不执行时间分割并且使用长的时间步长时，分析具有以点P为其中心的大的预测球体45。然而，当进行时间分割并且使用短的时间步长时，由于点P处的粒子的可移动范围减小，因此分析比预测球体45小的预测球体46。即，纹线计算单元130并非在大的预测球体45上执行单一计算，而是通过将时间区间分割成Ndiv个时间点并且获得分割的路径来得到较小的预测球体46。因此，减少了作为计算目标的单元的数目，并且减少了计算量。

图22示出了根据时间分割数目的总计算时间的变化的示例。在图22的曲线图中，横轴表示分割数目Ndiv，而纵轴是计算时间。在图22中，虚线表示理论计算时间51的转变，而实线表示实际计算时间52的转变。实际计算时间52是时间区间ti≤t≤ti+1中的总计算时间的测量结果。当Ndiv＝0时，不使用预测球体，并且所有的单元被用作计算目标。

理论上，分割数目Ndiv越大，计算时间将越短。然而，实际上，如果分割数目太大，则每分割后单位时间执行的处理的次数增加，例如创建预测球体中的单元的列表。从而，处理时间增加。因此，存在作为分割数目Ndiv的最优值。例如，在开始纹线的计算之前，纹线计算单元130提前针对目标系统来测量最优值作为分割数目Ndiv同时改变分割数目Ndiv。

[通过统计分析实现的预测球体半径的最小值的确定方法]

流体模拟通过使用有限数目的离散点信息项来执行。因此，如果设置过小的预测球体半径，则没有离散点信息项会包括在预测球体中。这意味着存在预测球体半径的下限。同时，由于预测球体半径通过表达式(24)设置，因此存在时间步长的下限。当通过使用预测球体执行计算时，设置预测球体半径的下限Rmin。当预测球体半径等于或小于下限时，将下限Rmin设置为预测球体半径。当预测球体半径等于或小于下限时，由于使用比预测球体半径大的Rmin，因此纹线上的点没有落在预测球体的外部。因此，为了稳定地进行计算，设置下限Rmin是重要的。此外，下限Rmin的值与时间分割数目Ndiv的设置有关。由于R0＝|Δt向量vmax|是最大移动距离，当确定了下限Rmin的值时，通过使用上取整函数来根据以下表达式设置分割数目Ndiv。

由于分割数目Ndiv的值与计算速度和计算精确度有关，因此设置下限Rmin在计算速度和计算精度方面也是重要的。

然而，在设定下限Rmin时需要注意。具体地，引入概率模型，并且纹线计算单元130执行允许计算失败的推测性计算。当计算失败时，纹线计算单元130通过使用使计算必定成功的参数来执行计算。该重新计算所需的时间被认为是处罚。统计上使包括处罚的计算时间最小化的下限Rmin被设置。接下来，将详细描述如何设置下限Rmin。

首先，纹线计算单元130将半径Rw设置为用作处罚的最差值。模拟中使用的单元的最长边长被用作半径Rw。

图23示出了单元的边长的分布的示例。在图23的曲线图中，横轴表示单元的边长，并且纵轴表示指示每个边长的概率密度。在曲线图中，由虚线来表示两种心脏模拟#1和#2的概率密度函数61和62。

在图23的示例中，模拟中使用的一些单元非常粗糙，并且最长边是0.008[m]。因此，如果将这些单元直接用于计算，则分割数目Ndiv也减小，并且速度也降低。然而，大多数计算可以在不使用最长边的情况下执行。

在设置作为最差值的半径Rw之后，纹线计算单元130设置用于计算的下限Rmin。在该操作中，纹线计算单元130在假设计算失败的同时计算统计计算成本。此外，鉴于推测性计算的执行，纹线计算单元130选择使计算量最小化的半径作为下限Rmin，其中，所述计算量包括计算失败时出现的处罚。虽然这很大程度上取决于假设的概率模型，但是将通过使用简单的示例来描述推测性计算。首先，使用纹线上的点的计算假设所有单元可以是具有等概率的分析目标。在这种情况下，计算根据图23的概率密度函数61和62而成功。

以下描述假设当半径是R时计算时间是T。在这种情况下，对应的计算量与R³成比例。此外，以下描述假设通过将半径减小至βR(0<β≤1)计算成功的概率为p[％]。计算成功的情况是点的目的地落在具有半径βR的预测球体的内部的情况。当计算失败时，通过使用半径R来执行重新计算。在这种情况下，包括处罚的计算时间T'通过表达式(28)来表示。

图24示出了基于预测球体半径的计算成本的变化。在图24中，纵轴表示计算成本(T'/T)，而横轴表示预测球体半径R。如图24所示，可以看出，当预测球体半径R约为0.003[m]时的包括处罚的成本比当预测球体半径是最差值半径Rmax(0.008[m])时的包括处罚的成本小约10％。也就是说，可以看出，存在实现包括处罚的在统计上的最小计算量的半径。在实践中，由于执行有限数目的计算次数，因此所有单元并非被等概率地用作分析目标。因此，比以上值更小的值被用作半径。

在这种情况下，尽管可以指示存在最优的预测球体半径R，但是最优值本身倾向于被高估。因此，为了防止这样的高估，可以通过以下操作来获得如图25所示的移动距离的分布：获得各个离散点处的速度场的范数，并且将时间步长乘以范数。

图25示出了移动距离的概率分布。在图25中，横轴表示移动距离，并且纵轴表示出现的数目。即使当存在两个具有相同形状和尺寸的四面体时，如果一个四面体上的速度场的大小是另一四面体上的速度场的大小的两倍，则点通过目标四面体所需的时间减半。图26基于移动距离的分布通过使用表达式(28)示出作为预测球体半径R的函数的包括处罚的计算效率。

图26示出了在假设概率分布时的计算效率曲线。在图26中，横轴表示预测球体半径R，并且纵轴表示计算效率。较小的计算效率值指示较好的效率。根据图26，最佳的计算效率与近似从0.0012[m]至0.0016[m]的范围中的预测球体半径R相对应。这是支持经验最优值Rmin＝0.0015±0.0005[m]的模型。

以这种方式，获得预测球体半径R的最佳值。即，虽然当离散点落在预测球体的外部时纹线计算单元130需要执行重新计算作为处罚，但是纹线计算单元130能够鉴于该处罚来设置预测球体的最高效的半径R。因此，纹线计算的效率提高。

[具体计算示例]

在下文中，将具体描述在通过使心脏模拟器分析心肌和冠状动脉循环来实际计算纹线时测量的计算速度。

首先，将描述计算中使用的数据。输出对应于单次心跳的模拟结果，并且每0.01[秒]输出100个心脏状态数据。由于心脏搏动，因此速度场随时间变化。因此，速度场是不稳定的流。此外，由于搏动，心肌重复舒张和收缩。因此，心肌表面也移动，导致移动边界问题。为了描述心脏中的血流在该系统中的传输，将粒子生成源布置在心脏中的多个点处，并且计算纹线。

当可视化设备100从外部读取关于心肌和流体的信息时，纹线的数目M和计算的次数N被输入至可视化设备100。此外，纹线生成源的位置被输入至可视化设备100。纹线生成源被布置在心脏中的左心室和右心房中。根据纹线生成源的布置的指示，例如，可视化设备100将20个纹线生成源随机地布置在具有0.05[m]的半径的球体的内部的各个流体部分中。此外，计算次数被设置成与0.3次心跳相对应的30。通过使用四阶龙格-库塔法来执行时间演化计算。

在设置以上初始条件之后，当可视化设备100开始纹线计算时，每个步长从粒子生成源生成粒子。根据图8所示的纹线计算流程图来计算各个粒子在时间步长Δt＝0.01[秒]之后的位置。所计算的纹线被显示在监视器21上。

图27示出了纹线的显示的示例。在图27的示例中，多条纹线71被叠加在心脏70的截面图像上。

当计算纹线时，纹线计算单元130使用预测球体来减小计算量并提高计算精确度。在使用预测球体的计算中，通过分割数目Ndiv来进一步分割时间步长Δt＝0.01[秒]。将实现统计上的最小计算量的最优值自动设置为分割数目Ndiv。在该模拟中，分割数目Ndiv被设置成约3至7。

通过使用分割数目Ndiv进行优化来定量地测量计算速度，并且图22示出了测量的结果。图22的示例还示出了当没有使用预测球体时的纹线计算时间，以阐明通过使用预测球体而获得的有利效果。Ndiv＝0与当没有使用预测球体时——即当所有单元被用作计算目标时——所需的时间相对应。

当使用预测球体(Ndiv＝1)时即当基本上不执行时间分割时，需要时间来建立预测球体。因此，比不使用预测球体的情况需要更多的时间(约1.47倍)。然而，当Ndiv＝2时，与Ndiv^-2成比例的关于计算成本的降低效应比建立预测球体的固定成本具有更大的影响。因此，与不使用预测球体的情况相比，计算时间减少约22.9％。虽然计算时间随着分割数目Ndiv的增加而缩短，但是可以看出，由于建立预测球体的固定成本，所以存在性能改善的上限。在图22的示例中，最佳分割数目是4。如果分割数目Ndiv进一步增加，则由于固定成本的效应具有更大的影响，因此计算时间开始增加。当Ndiv＝10时，性能比最好的情况恶化了约53％。

如上所述，虽然分割数目Ndiv存在最佳值，但是在实际计算中最佳分割数目在计算纹线之前设置。此外，在图22的计算中，由于网格点的数量不大于约50000，因此在不使用预测球体的情况下的计算时间比使用预测球体并且分割数目Ndiv为1的情况的计算时间更短。然而，如果模拟系统增大，则不使用预测球体的情况下的计算量与网格点的数量成比例地增加。例如，如果网格数量增加10倍，则时间增加10倍。然而，当使用预测球体时，由于预测球体半径R仅基于时间步长和速度的最大值来设置并且始终恒定，因此计算量也是恒定的。因此，执行模拟的系统越大，通过使用预测球体而获得的有利效果将越大。

此外，将描述通过对预测球体进行时间分割而获得的另一有利效果——精确度的提高。虽然表达式(2)至(6)是在四阶龙格-库塔法的情况下给出，但是已知误差阶数被给出为O(Δt⁴)。因此，当时间由分割数目Ndiv分割时，单一龙格-库塔运算中的误差的结果是Ndiv^-4倍。即使当在时间区间ti≤t≤ti+1中基于龙格-库塔法的计算总计Ndiv次并且误差累积了Ndiv次时，总误差保持Ndiv^-3。当分割数目Ndiv是4时，总误差是1/64倍。

图28示出了当时间步长改变时的轨迹精确度的变化。在图28中，横轴表示模拟时间，并且纵轴表示离散点的x坐标。图28的示例示出当时间步长Δ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来设置。在实际计算中，如图25所示，鉴于单元的边长和速度场二者来假设概率模型，并且计算概率分布。此外，根据概率分布，通过使用表达式(28)来计算计算效率，并且计算实现最佳计算效率的最小值Rmin。即，如图26所示，存在作为预测球体的半径的最小值，并且最小值Rmin约为0.0015±0.0005[m]。因此，计算总是以最小的计算成本执行。

[其他实施方式]

纹线上的点彼此独立，彼此之间没有相互作用。因此，由于并行化适当地适用，因此可以通过使用消息传递接口(MPI)或开放式多处理(OpenMP)来以并行方式执行纹线的计算。以这种方式，计算速度提高。

此外，虽然在第二实施方式中已经示出通过使用心脏模拟的结果来将血流的纹线可视化的示例，但是第二实施方式也可以应用于其他流体模拟的结果。例如，当将可变翼机构布置在汽车的后部时，可以模拟和分析当可变翼机构被激活时可变翼机构周围的空气的流动。通过使用根据第二实施方式的可视化设备100，当可变翼机构被激活时可变翼机构周围的空气的流动通过纹线来可视化。此外，第二实施方式也适用于指示当飞机的翼摆动时在摆动翼周围的空气的流动的流体模拟结果。

根据一个方面，即使当结构变形时，也能够追踪纹线。