CN102804187B - 物理量计算方法、数值解析方法、物理量计算装置及数值解析装置 - Google Patents

Abstract

在对物理现象进行数值性解析的数值解析方法中计算物理量的物理量计算方法，其特征在于，包括物理量计算步骤，计算被分割为多个分割区域的解析区域中的物理量，在该物理量计算步骤中，使用离散化的控制方程式和计算用数据模型计算上述物理量，上述控制方程式以仅使用不需要上述分割区域的顶点的坐标(Vertex)及该顶点的连接信息(Connectivity)的量的方式根据加权残值法导出，在上述计算用数据模型中，具有各上述分割区域的体积及表示相邻的上述分割区域之间的边界面的特性的边界面特性量，来作为不需要上述分割区域的顶点的坐标(Vertex)及该顶点的连接信息(Connectivity)的量。

Description

物理量计算方法、数值解析方法、物理量计算装置及数值解析装置

技术领域

本发明涉及一种物理量计算方法、数值解析方法、物理量计算程序、数值解析程序、物理量计算装置及数值解析装置。

背景技术

一直以来，作为通过数值解析求出流速分布、应力分布及热分布等的数值解析方法，例如已知有限要素法、有限体积法、体素法及粒子法。

该数值解析方法一般由预处理、求解处理和后处理构成。并且，在预处理中做成计算用数据模型，在求解处理中，使用该计算用数据模型及离散化的控制方程式(以下称为离散化控制方程式)，进行上述物理量的计算。

现有的有限体积法中，例如将解析区域分割为多个区域，使用各分割区域的体积、相邻的分割区域的边界面的面积及该边界面的法线矢量，计算各分割区域中的物理量。

在有限体积法中，在预处理中做成包括各分割区域的顶点的坐标(Vertex)的计算用数据模型(通常称为网格(mesh))，在求解处理中使用该计算用数据模型中含有的Vertex等，算出上述分割区域的体积、边界面的面积及边界面的法线矢量，使用这些值进行物理量的计算。Vertex是用于规定分割区域的几何学形状的量。因此，可以说在有限体积法中，在求解处理中使用分割区域的几何学形状，进行分割区域的体积、边界面的面积及边界面的法线矢量的计算。

进一步，在有限体积法中，也可以具有部分为满足相邻的分割区域中的顶点共享的条件的部分。因此，在有限体积法中，存在对分割区域的制约略有缓和的情况，但利用的解析要素类型例如限定为四元要素、六元要素、棱柱要素、棱锥要素等。

此外，如专利文献1所示，也提出了不限定解析要素类型的有限体积法。但这种不限定解析要素类型的有限体积法中，也和上述现有的有限体积法一样，在预处理中做成包括各分割区域的顶点的坐标(Vertex)的计算用数据模型，在求解处理中使用该计算用数据模型中含有的Vertex等，进行物理量的计算。

并且，众所周知，有限要素法是使用内插函数算出各分割区域中的物理量的方法，但和有限体积法一样，在求解处理中，使用由Vertex等规定的分割区域的几何学形状。

体素法及粒子法和有限要素法、有限体积法相比，是可容易地做成计算用数据模型的数值解析方法。

体素法是如下方法：将由基本上同一尺寸的长方体形状的多个体素(正交格子)定义解析区域的体素数据，作为计算用数据模型做成，进行利用了该体素数据的物理量计算，从而进行数值解析。作为体素法从大的方面划分为：加权残值法类型，使用加权残值法下的控制方程式；非积分法类型，例如使用元胞自动机模型、格子波尔兹曼模型(latticeBoltzmannmodel)等。并且，根据该体素法，作为体素数据，不需要Vertex等。

根据该体素法，用体素分割解析区域，从而可容易地定义解析区域，可短时间内做成计算用数据模型。

另一方面，粒子法(particlemethod)是如下方法：将由多个粒子定义解析区域的粒子数据，作为计算用数据模型做成，进行使用了该粒子数据的物理量计算，从而进行数值解析。粒子法中，作为控制方程式，在非积分法类型中使用粒子间相互作用模型。在粒子法中，因没有分割区域，所以不需要Vertex等。根据该粒子法，通过对解析区域例如均匀配置粒子，可容易地定义解析区域，可短时间内做成计算用数据模型。

在先技术文献

专利文献

专利文献1：美国专利申请公开第2008/0021684号说明书

发明内容

发明要解决的问题

如现有的有限要素法、有限体积法等数值解析方法所示，在求解处理中使用分割区域的几何学形状时，对于计算用数据模型，当然需要使之具有表示分割区域的几何学形状的数据。

为了定义分割区域的几何学形状，除了Vertex外，还需要顶点的连接信息(ConnectivityofVertex、以下简称为Connectivity)。因此，在有限要素法、有限体积法中，需要使计算用数据模型具有Vertex和Connectivity。

此外，具体而言，Connectivity由对全部分割区域的顶点依次附加的全部节点代码和在一个分割区域内对顶点依次附加的局部节点代码的对应信息定义。

众所周知，这种具有Vertex和Connectivity的计算用数据模型在做成时需要非常庞大的作业。

例如，在有限要素法中使用的计算用数据模型中，如图1所示，需要满足邻接的分割区域必须共享Vertex这一条件地做成计算用数据模型，为使所有分割区域满足这一条件，需要非常庞大的时间。

另一方面，在有限体积法中使用的计算用数据模型如图2所示，允许存在不被邻接的分割区域共享的Vertex，所以和有限要素法相比，网格生成的自由度增加。但在有限体积法中，需要在不共享的Vertex至少存在于邻接的分割区域的边上、及一般情况下使分割区域的形状与预先设定的解析要素类型对应这二个条件下，做成计算用数据模型，网格生成很难说自由度较大。

并且近年来，对从三维CAD(ComputerAidedDesign：计算机辅助设计)数据等三维形状数据提取的解析区域进行数值解析。但是，三维形状数据不是在数值解析中形成的数据，含有表示面重叠、面交叉、面间间隙、小孔等的数据，包括很多不适于做成具有Vertex和Connectivity的计算用数据模型的条件。因此，为了可做成这些具有Vertex和Connectivity的计算用数据模型，需要修正或变更三维形状数据。并且，为了可做成具有Vertex和Connectivity的计算用数据模型，在修正或变更三维形状数据时，需要进行经验积累、试验错误所需的非常庞大的手工作业。这是实际利用有限要素法、有限体积法时的较大问题。

并且，如有限体积法那样在求解处理中进行分解区域的体积、边界面的面积及边界面的法线矢量的计算时，求解处理中的计算量进一步增加，求解处理中的计算负荷进一步增大。

在体素法中，虽然可短时间内做成计算用数据模型，但存在以下问题。体素法基本上由解析区域完全相同尺寸的体素(正交格子)定义。通常在有限要素法、有限体积法中，通过将想得到较高解析精度的区域的要素尺寸(分割区域的尺寸)设定得较小，从而对该区域进行正确的物理量计算，并将其他区域的要素尺寸设定得较大，从而降低对该区域的计算负荷。但在体素法中，因所有体素基本是同一尺寸，所以将体素设定得较小时，计算负荷变得非常大，将体素设定得较大时，解析精度恶化。

并且，在体素法中，需要通过排列同一尺寸的体素(正交格子)来定义解析区域，因此存在在和外部区域的边界附近无法使解析区域平滑、而变为阶梯状的情况。即，即使实际上想解析的区域具有斜面、曲面等，在体素数据中，该区域也表现为阶梯状。因此，体素法中的解析区域形状和实际想解析的区域形状不同，解析精度恶化。

对此提出了以下称为切割单元法的改良方法：将体素数据的阶梯状的区域沿着实际想解析的区域所具有的斜面、曲面切断(边界校正)。但是，根据该改良方法，通过该边界校正，易产生非常小的分割区域，当生成这种小的分割区域时，会导致解析精度的恶化。

并且，在该改良方法中，为形成切割单元及在求解处理中，使用Vertex。

如上所述，在不进行边界校正的体素法中，虽然不需要Vertex，但在体素的生成、即所谓网格生成中存在界限。即，当要获得充分的解析精度时，体素数增加、求解处理中的计算负荷也增加，成为问题。进一步，在进行边界校正的体素法的改良方法中，结果需要Vertex，因此最终受到分割区域的几何学形状的影响，在与外部区域的边界周边的分割区域形成的处理中，伴随着经验积累、试验错误所需的非常庞大的手工作业，无法短时间内做成形状数据模型。

另一方面，在粒子法中，需要算出某个特定粒子和其他粒子的结合关系。因此，需要搜索存在于该特定粒子附近的粒子。并且，该粒子的附近搜索处理原则上对全部粒子来进行。但在粒子法中，随着时间的变化，各粒子移动，因此粒子之间的结合关系总是变化。所以每当解析中的时间变化时，需要进行附近搜索处理，导致计算负荷的增加。严格挑选作为附近搜索的对象的粒子等，以尝试降低附近搜索处理的计算负荷，但例如为了提高解析精度而增大粒子数时，计算负荷与粒子数的2次方成比例地增大。

在这种粒子法中，为了实现实用时间内的数值解析，在大型的并列计算机中需要使用较多的CPU(CentralProcessingUnit：中央处理单元)。例如存在以下实例：在使用了Vertex和Connectivity的一般的有限体积法求解中，由1个CPU耗时半天的计算，在粒子法中，由使用了32个CPU的并列计算耗时1周以上。

并且，在粒子法中，紧密配置粒子时，计算负荷变得非常大，疏松配置粒子时，解析精度恶化。

进一步，在粒子法中，稍后详述，在解析基于液体、构造、热量、扩散等物理量的守恒定律的物理现象时，其守恒性未充分满足。

例如，面向解析区域和外部区域的边界面配置的粒子，没有在边界面中具有多大面积的信息。因此，即使提供了从边界面使热量进入的条件，也无法正确掌握各粒子中进入了多少热量，无法获得精度良好的定量值。

本发明鉴于作为上述现有的数值解析方法的有限要素法、有限体积法、体素法、体素法的改良方法、及粒子法的问题点而产生，其目的在于减轻计算用数据模型的做成中的作业负担，并且不伴随解析精度的恶化地降低求解处理中的计算负荷。

用于解决问题的手段

为解决上述问题，本发明提供一种物理量计算方法，在对物理现象进行数值性解析的数值解析方法中计算物理量，其采用以下的构成：包括物理量计算步骤，计算被分割为不只由正交格子形状形成的多个分割区域的解析区域中的物理量，在该物理量计算步骤中，使用离散化的控制方程式和计算用数据模型计算上述物理量，上述控制方程式以仅使用不需要上述分割区域的顶点的坐标(Vertex)及该顶点的连接信息(Connectivity)的量的方式根据加权残值法导出，在上述计算用数据模型中，具有各上述分割区域的体积及表示相邻的上述分割区域之间的边界面的特性的边界面特性量，来作为不需要上述分割区域的顶点的坐标(Vertex)及该顶点的连接信息(Connectivity)的量。

本发明中使用的离散化控制方程式不象现有技术那样以包括分割区域的几何学形状的量(Vertex和Connectivity)的形式来表现，本发明不需要规定分割区域的几何学形状的量。本发明中使用的离散化控制方程式，可以通过在根据加权残值法导出使用用于规定现有的几何学形状的量的方程式的过程中硬停留在中途而获得。这种本发明中使用的离散化控制方程式以不需要分割区域的几何学形状的量(即不需要Vertex和Connectivity和量)表现，例如可以是仅取决于分割区域的体积和边界面特性量这二个的形式。

即，在现有的有限要素法、有限体积法中，前提是将解析对象物分割成微小区域，因此以使用规定该微小区域的几何学形状的量、即Vertex和Connectivity为前提，导出离散化控制方程式，但在本发明中使用的离散化控制方程式根据和现有技术不同的新的思想而导出。

并且，本发明的特征在于，使用该基于新思想而导出的离散化控制方程式，和现有的数值解析方法不同，其不依赖几何学形状地解决现有问题，起到各种显著的效果。

在此说明分割区域的体积和边界面特性量是不需要规定分割区域的特定几何学形状的Vertex和Connectivity的量的情形。此外，不需要Vertex和Connectivity的量是指，不使用Vertex和Connectivity也可进行定义的量。

例如，考虑分割区域的体积时，用于使分割区域的体积变为某规定值的分割区域的几何学形状存在多个。即，体积取某规定值的分割区域的几何学形状可以是立方体，也可以是球。并且，例如，分割区域的体积可如下定义：在全部分割区域的总和与解析区域整体的体积一致的制约条件下，例如使分割区域的体积尽量与邻接分割区域的平均距离的3次方成比例的最优化计算来定义。因此，分割区域的体积可以是不需要分割区域的特定的几何学形状的量(不需要Vertex和Connectivity的量)。

并且，作为边界面特性量，例如包括边界面的面积、边界面的法线矢量、边界面的周长等，但这些边界面特性量取某规定值的分割区域的几何学形状(即边界面的几何学形状)存在多个。并且，例如，边界面特性量可如下定义：对包围各分割区域的所有边界面，在法线矢量的面积加权平均矢量的长度为零的制约条件下，使边界面的法线矢量的方向靠近连接邻接的二个分割区域的控制点(参照图5)的线段，且使包围分割区域的所有边界面面积的总和尽量与该分割区域的体积的3/2次方成比例的最优化计算来定义。因此，边界面特性量可以取不需要分割区域的特定的几何学形状的量(不需要Vertex和Connectivity的量)。

并且，在本发明中，“分割为不只由正交格子形状形成的多个分割区域的解析区域”是指，构成解析区域的多个分割区域的至少任意一个不取正交格子形状。即，意味着解析区域包括正交格子形状以外的形状的分割区域。

并且，在本发明中，“仅使用不需要Vertex和Connectivity的量”是指，代入到离散化控制方程式的值仅是不需要Vertex和Connectivity的量。

接着参照图3的概念图，对比使用了本发明的数值解析方法和现有的数值解析方法中的预处理及求解处理的同时，对本发明的显著效果进行更详细的说明。

使用本发明的数值解析方法中，如图3所示，在求解处理(本发明的物理量计算步骤)中，使用仅利用了不需要Vertex和Connectivity的量的离散化控制方程式，进行分割区域中的物理量的计算。因此，在解离散化控制方程式时，在预处理中做成的计算用数据模型中，无需含有Vertex和Connectivity。

并且，在使用本发明时，作为不需要Vertex和Connectivity的量，使用分割区域的体积和边界面特性量。因此，在预处理中做成的计算用数据模型不具有Vertex和Connectivity，具有分割区域的体积、边界面特性量、其他辅助数据(例如下述分割区域的配合信息、控制点坐标等)。

因此，在使用本发明时，如上所述，根据分割区域的体积和上述边界面特性量、即不需要分割区域的几何学形状的量，可计算各分割区域中的物理量。因此，计算用数据模型中可不具有分割区域的几何学形状、即不具有Vertex和Connectivity地算出物理量。因此，通过使用本发明，在预处理中，做成至少具有分割区域的体积和边界面特性量(边界面的面积及边界面的法线矢量)的计算用数据模型即可，可不做成具有Vertex和Connectivity的计算用数据模型地进行物理量的计算。

不具有Vertex和Connectivity的计算用数据模型无需分割区域的几何学形状，因此可不受限于分割区域的几何学形状地做成。

因此，对三维形状数据的修正作业的限制也大幅降低。因此，不具有Vertex和Connectivity的计算用数据模型和具有Vertex和Connectivity的计算用数据模型相比，可非常容易地做成。因此，根据本发明，可减轻计算用数据模型的做成中的作业负担。

并且，在使用本发明时，在预处理中也可使用Vertex和Connectivity。即，在预处理中，可使用Vertex和Connectivity算出分割区域的体积、边界面特性值等。这种情况下，在求解处理中只要有分割区域的体积、边界面特性值，就可计算物理量，因此在预处理中即使利用Vertex和Connectivity，也不会产生对分割区域的几何学形状的限制、例如源自分割区域的应变、扭曲等造成的制约，可减轻计算用数据模型的做成中的作业负担。

并且，通过使用本发明，在预处理中不再有对分割区域的几何学形状的制约，因此可将分割区域变更为任意的形状。因此，可不增加分割区域的个数，使解析区域容易地与实际想解析的区域对应，不增加计算负荷地提高解析精度。

进一步，通过使用本发明，分割区域的分布密度也可任意地变更，因此在必要范围内可允许计算负荷增大的同时，进一步提高解析精度。

并且，通过使用本发明，和现有的数值解析方法不同，在求解处理中，无需使用Vertex和Connectivity来算出分割区域的体积及边界面特性量。因此，可降低求解处理中的计算负荷。

并且，在本发明中，当解析区域的形状不变化时，无需分割区域的移动，因此无需如粒子法那样，每当时间变化时就进行附近搜索处理，计算负荷较小。并且稍后详述，通过使用本发明，和粒子法不同，可在满足物理量的守恒定律的同时进行物理量的计算。

另一方面，作为现有的数值解析方法的有限体积法中，在预处理中做成具有表示分割区域的几何学形状的Vertex和Connectivity的计算用数据模型，在求解处理中使用计算用数据模型中含有的Vertex和Connectivity，算出分割区域的体积和边界面特性量(边界面的面积及边界面的法线矢量)后，计算各分割区域中的物理量。此时，要求对几何学形状的制约、即Vertex和Connectivity的关系上不存在问题。因此，需要在分割区域的应变、扭曲等的制约中做成计算用数据模型(即网格)，如上所述，存在计算用数据模型做成中需要庞大的手工作业的问题。

并且，有限要素法中，在求解处理中也使用计算用数据模型中含有的Vertex和Connectivity来计算物理量，因此需要在预处理中做成具有表示分割区域的几何学形状的Vertex和Connectivity的计算用数据模型，在计算用数据模型的做成中产生庞大的手工作业。

并且，作为现有的数值解析方法的体素法如图3所示，在求解处理中算出物理量时，虽然不需要Vertex和Connectivity，但分割区域的形状限定为体素，因此如上所述，存在与外部区域的边界变为阶梯状的问题。因此，如上所述，为了获得充分的解析精度，体素数也增加，求解处理中的计算负荷增大，成为问题。并且，在进行边界校正的体素法中，最终算出分割区域的体积等时，利用Vertex，在做成计算用数据模型时，受到分割区域的几何学形状的影响。

并且，在作为现有的数值解析法的粒子法中，因不存在分割区域的概念，所以如图3所示，在求解处理中算出物理量时，虽然不需要Vertex和Connectivity，但替代分割区域，因定义计算用数据模型的粒子的移动，如上所述，计算负荷增大。并且在粒子法中，难以在满足守恒定律的同时进行物理量计算。

接着参照图4，对本发明和现有的有限体积法进一步详细进行比较。

在现有的有限体积法中，如上所述，在预处理中，做成具有用于定义通过网格分割获得的分割区域的几何学形状的Vertex和Connectivity的计算用数据模型。并且一般情况下，在求解处理中需要分割区域的配合信息(以下称为link)。因此在预处理中，做成具有Vertex和Connectivity的计算用数据模型。

并且，在现有的有限体积法中，如图4所示，将具有Vertex、Connectivity和link的计算用数据模型、及在求解处理中必须的边界条件、初始条件等，从预处理交接到求解处理。在求解处理中，使用交接的计算用数据模型中含有的Vertex、Connectivity等，解离散化控制方程式，从而进行物理量的计算。

另一方面，在本发明中，在预处理中，做成具有任意配置的分割区域的体积、边界面特性量(边界面的面积、边界面的法线矢量)及link的计算用数据模型。并且稍后详述，在本发明中，也存在以下情况：根据需要，使计算用数据模型中具有配置在分割区域内部的控制点的坐标。

并且在本发明中，如图4所示，将具有分割区域的体积、边界面特性量及link(根据需要也可是控制点的坐标)的计算用数据模型、和边界条件、初始条件等，从预处理交接到求解处理。在求解处理中，使用交接的计算用数据模型中含有的分割区域的体积、边界面特性量等，解离散化控制方程式，从而进行物理量的计算。

并且从图4可知，在本发明中，在求解处理中，不使用Vertex和Connectivity地计算物理量，这一点和现有的有限体积法大为不同，这一点是本发明的最大特征。这一特征在求解处理中，使用仅利用不需要Vertex和Connectivity的量的离散化控制方程式来获得。

其结果如图4所示，在本发明中，无需将Vertex和Connectivity交接到求解处理，在预处理中，做成不具有Vertex和Connectivity的计算用数据模型即可。因此，和现有的有限体积法相比，本发明可极容易地做成计算用数据模型，可减轻计算用数据模型的做成中的作业负担。

并且，存在进行数值解析的解析区域的形状以时间序列变化的情况，即解析区域包括移动边界的情况。这种情况下，需要对应移动边界使分割区域移动及变形。

在现有的有限体积法中，通过以下方法进行包括移动边界时的物理量的计算：每当移动边界移动时预先存储Vertex的方法；因分割区域的扭曲变形而不可计算时，再执行区域分割的方法。与之相对，在本发明中，可通过替代Vertex而预先求出分割区域的体积、边界面特性值等并存储的方法，或通过区域分割的再执行，来进行包括移动边界时的物理量的计算。

在现有的有限体积法或本发明中，采用上述任意一种方法时，均需要做成多个计算用数据模型。但是在现有的有限体积法中，做成多个至少一个需要庞大的作业量的计算用数据模型时，其作业量大多情况下并非在实际上可负担的范围下进行。

另一方面，在本发明中，计算用数据模型无需具有Vertex和Connectivity，在区域分割过程中无需考虑Vertex、Connectivity的匹配性，因此可非常快速地做成计算用数据模型，容易地进行含有移动边界时的物理量计算。

在此补充说明上述link。Link是使进行物理量的交换的分割区域之间建立关联的信息。并且，通过该link建立关联的分割区域不必在空间上邻接，可在空间上分离。这种link与Vertex、Connectivity不关联，和Vertex及Connectivity相比，可以极短的时间做成。

接着详细说明使用了本发明的数值解析方法(以下称为本数值解析方法)的原理、即通过根据加权残值法导出的离散化控制方程式、分割区域的体积和边界面特性量可算出物理量的原理。并且在以下说明中，[]中的文字表示附图中粗体字所示的矢量。

首先，本数值解析方法中的计算用数据模型使用分割解析区域所获得的各分割区域的体积、表示相邻的分割区域之间的边界面的特性的边界面特性量来定义。

图5是表示该本数值解析方法的计算用数据模型的一例的概念图。

在该图中，单元R₁、R₂、R₃……是分割解析区域而获得的分割区域，分别具有体积V_a、V_b、V_c……。并且，边界面E是在单元R₁和单元R₂之间进行物理量交换的面，相当于本发明中的边界面。并且，面积S_ab表示边界面E的面积，是本发明中的边界面特性量之一。并且，[n]_ab表示边界面E的法线矢量，是本发明中的边界面特性量之一。并且，控制点a、b、c……配置在各单元R₁、R₂、R₃的内部，在图5中，配置在各单元R₁、R₂、R₃……的重心位置上。但控制点a、b、c……并非必须配置在各单元R₁、R₂、R₃……的重心位置上。并且，α表示从控制点a到控制点b的距离为1时的、从控制点a到边界面E为止的距离，是表示边界面E存在于连接控制点a和控制点b的线段的哪个内分点的比率。

此外，边界面不仅存在于单元R₁和单元R₂之间，而且存在于相邻的所有单元之间。并且，边界面的法线矢量及边界面的面积也按照各边界面提供。

并且，实际的计算用数据模型作为具有以下数据的数据组构建：各控制点a、b、c……的配置数据；表示存在各控制点a、b、c……的单元R₁、R₂、R₃……的体积V_a、V_b、V_c……的体积数据；表示各边界面的面积的面积数据；表示各边界面的法线矢量的法线矢量数据。

即，本数值解析方法的计算用数据模型的定义中具有：单元R₁、R₂、R₃……的体积V_a、V_b、V_c……；作为表示相邻的单元R₁、R₂、R₃……之间的边界面的特性的边界面特性量的、边界面的面积；作为表示相邻的单元R₁、R₂、R₃……之间的边界面的特性的边界面特性量的、边界面的法线矢量。

此外，各单元R₁、R₂、R₃……具有控制点a、b、c……。因此，单元R₁、R₂、R₃……的体积V_a、V_b、V_c……可作为控制点a、b、c……虚拟占据的空间(控制量)的体积。

并且，本数值解析方法的计算用数据模型根据需要可具有比率数据，其表示边界面存在于连接位于该边界面两侧的控制点之间的线段的哪个内分点上的比率α。

以下说明使用上述计算用数据模型求出解析区域的各单元(分割区域)中的流速的物理量计算示例。并且，在此将各控制点中的流速作为各单元中的流速求出。

首先，在本物理量计算中，本数值解析方法在液体解析时，使用下式(1)所示的纳维-斯托克斯(Navier-Stokes)式、及下式(2)所示的连续式。

(数式1)

$\frac{\∂}{\∂t}\left({\mathrm{\ρu}}_i\right)+\frac{\∂}{\∂x_j}\left({\mathrm{\ρu}}_j{u}_i\right)=-\frac{\∂P}{\∂x_i}+\frac{\∂}{\∂x_j}\[\mathrm{\μ}\frac{\∂u_i}{\∂x_j}\]\mathrm{...}\left(1\right)$

(数式2)

$\frac{\∂\left({\mathrm{\ρu}}_j\right)}{\∂x_j}=\mathrm{0...}\left(2\right)$

并且在公式(1)、(2)中，t表示时间，x_i(i＝1、2、3)表示笛卡尔(Cartesian)系中的坐标，ρ表示流体密度，u_i(i＝1、2、3)表示流体的流速分量，P表示压力，μ表示液体的粘性系数，下标i(i＝1、2、3)、j(j＝1、2、3)表示笛卡尔坐标系中的各方向分量。并且，下标j基于求和约定。

并且，将公式(1)、(2)根据加权残值法对控制量的体积积分表示时，公式(1)如下式(3)所示，公式(2)如下式(4)所示。

(数式3)

${\∫}_V\frac{\∂}{\∂t}\left({\mathrm{\ρu}}_i\right)dV+{\∫}_S(n\·u)u_{i}dS=-{\∫}_S{n}_{i}PdS+{\∫}_S\mathrm{\μ}\frac{\∂u_i}{\∂n}dS\mathrm{...}\left(3\right)$

(数式4)

∫_Sρn·udS＝0…(4)

并且在公式(3)、(4)中，V表示控制量的体积，∫_vdV表示和体积V相关的积分，S表示控制量的面积，∫_sdS表示和面积S相关的积分，[n]表示S的法线矢量，n_i(i＝1、2、3)表示法线矢量[n]的分量，表示法线方向微分。

在此为简化说明，设流体密度ρ和粘性系数μ为常数。但以下常数化可扩展到流体的物性值根据时间、空间、温度等变化的情况。

并且，对于图5的控制点a，对边界面E的面积S_ab离散化，变换为代数方程式的近似式时，公式(3)如下式(5)所示，公式(4)如下式(6)所示。

(数式5)

$V_a\·\mathrm{\ρ}\frac{\∂u_i}{\∂t}+\underset{b=1}{\overset{m}{\Σ}}\[S_{ab}\·(n_{ab}\·u_{ab})u_{iab}\]=-\underset{b=1}{\overset{m}{\Σ}}\[S_{ab}\·n_{iab}\·P_{ab}\]+\underset{b=1}{\overset{m}{\Σ}}\[S_{ab}\·\mathrm{\μ}\·{\left(\frac{\∂u_i}{\∂n}\right)}_{ab}\]\mathrm{...}\left(5\right)$

(数式6)

$\underset{b=1}{\overset{m}{\Σ}}\[S_{ab}\·(n_{ab}\·u_{ab})\]=\mathrm{0...}\left(6\right)$

其中，加了下标ab的[n]_ab、[u]_ab、u_iab、n_iab、P_ab、表示控制点a和控制点b之间的边界面E上的物理量。并且，n_iab是[n]_ab的分量。并且，m是与控制点a为结合关系(夹着边界面的关系)的所有控制点的个数。

并且，公式(5)、(6)除以V_a(控制点a的控制量的体积)时，公式(5)如下式(7)所示，公式(6)如下式(8)所示。

(数式7)

$\mathrm{\ρ}\frac{\∂u_i}{\∂t}+\underset{b=1}{\overset{m}{\Σ}}\[\frac{S_{ab}}{V_a}\·(n_{ab}\·u_{ab})u_{iab}\]=-\underset{b=1}{\overset{m}{\Σ}}\[\frac{S_{ab}}{V_a}\·n_{iab}\·P_{ab}\]+\underset{b=1}{\overset{m}{\Σ}}\[\frac{S_{ab}}{V_a}\·\mathrm{\μ}\·{\left(\frac{\∂u_i}{\∂n}\right)}_{ab}\]\mathrm{...}\left(7\right)$

(数式8)

$\underset{b=1}{\overset{m}{\Σ}}\[\frac{S_{ab}}{V_a}\·(n_{ab}\·u_{ab})\]=\mathrm{0...}\left(8\right)$

在此设下式(9)。

(数式9)

${\mathrm{\φ}}_{ab}\≡\frac{S_{ab}}{V_a}\mathrm{...}\left(9\right)$

这样一来，公式(7)如下式(10)所示，公式(8)如下式(11)所示。

(数式10)

$\mathrm{\ρ}\frac{\∂u_i}{\∂t}+\underset{b=1}{\overset{m}{\Σ}}\[{\mathrm{\φ}}_{ab}\·(n_{ab}\·u_{ab})u_{iab}\]=-\underset{b=1}{\overset{m}{\Σ}}\[{\mathrm{\φ}}_{ab}\·n_{iab}\·P_{ab}\]+\underset{b=1}{\overset{m}{\Σ}}\[{\mathrm{\φ}}_{ab}\·\mathrm{\μ}\·{\left(\frac{\∂u_i}{\∂n}\right)}_{ab}\]\mathrm{...}\left(10\right)$

(数式11)

$\underset{b=1}{\overset{m}{\Σ}}\[{\mathrm{\φ}}_{ab}\·(n_{ab}\·u_{ab})\]=\mathrm{0...}\left(11\right)$

在公式(10)、(11)中，[u]_ab、u_iab、P_ab、通过控制点a和控制点b上的物理量的加权平均(对于对流项，是考虑到迎风性的加权平均)近似求出，依赖控制点a、b间的距离及方向、与存在于其间的边界面E的位置关系(上述比率α)、边界面E的法线矢量的方向来决定。但[u]_ab、u_iab、P_ab、是与边界面E的几何学形状无关的量(即是不需要规定单元形状的Vertex和Connectivity的量)。

并且，公式(9)定义的φ_ab也是(面积/体积)的量，是和控制量的几何学形状无关的量(即是不需要规定单元形状的Vertex和Connectivity的量)。

也就是说，该公式(10)、(11)是仅使用无需规定单元形状的Vertex和Connectivity的量即可算出物理量的、基于加权残值法的运算公式。

因此，在物理量计算(求解处理)前，做成上述计算用数据模型，在物理量计算中通过使用该计算用数据模型、及公式(10)、(11)的离散化控制方程式，在物理量计算中可完全不使用控制量的几何学形状(即规定单元形状的Vertex和Connectivity)地进行流速的计算。

这样一来，在物理量计算中，可完全不使用Vertex和Connectivity地进行流速的计算，因此计算用数据模型中不再需要具有Vertex和Connectivity。所以在做成计算用数据模型时，无需受限于单元的几何学形状，可任意设定单元形状。因此，根据本数值解析方法，如上所述，可大幅降低对三维形状数据的修正作业。

此外，实际解公式(10)、(11)时，[u]_ab、P_ab等边界面E上的物理量通常通过线性内插来内插。例如设控制点a的物理量为ψ_a、控制点b的物理量为ψ_b时，边界面E上的物理量ψ_ab可通过下式(12)求出。

(数式12)

${\mathrm{\ψ}}_{ab}=\frac{1}{2}({\mathrm{\ψ}}_a+{\mathrm{\ψ}}_b)\mathrm{...}\left(12\right)$

并且，物理量ψ_ab也可以通过使用比率α而由下式(13)求出，上述比率α表示边界面存在于连接其两侧的控制点之间的线段的哪个内分点上。

(数式13)

Ψ_ab＝(1-α)·Ψ_a+α·Ψ_b…(13)

因此，计算用数据模型具有表示比率α的比率数据时，可使用公式(13)，对边界面E上的物理量使用和距控制点a和控制点b分离的距离对应的加权平均来算出。

并且，连续体模型的方程式(纳维-斯托克斯式等)中，如公式(1)所示，包括1阶的偏导函数(偏微分)。

其中，对连续体模型的方程式的微系数使用部分积分、高斯散度定理、或者广义格林定理，将体积分(volumeintegration)变换为面积分(surfaceintegration)，降低微分的阶次。这样一来，1次微分可成为0次微分(纯量或矢量)。

例如，在广义格林定理中，物理量为ψ时，下式(14)的关系成立。

(数式14)

${\∫}_V\frac{\∂\mathrm{\ψ}}{\∂x_i}dV={\∫}_S{\mathrm{\ψn}}_{i}dS\mathrm{...}\left(14\right)$

此外，在公式(14)中，n_i(i＝1、2、3)是表面S上的单位法线矢量[n]的i方向的分量。

连续体模型的方程式的1次微分项通过体积分变换为面积分，在边界面上作为纯量或矢量处理。并且，这些值通过上述线性内插等，可由各控制点上的物理量内插。

并且稍后论述，根据连续体模型的方程式的不同，也包括含有2阶的偏导函数的情况。

进一步将公式(14)的被积分函数进行1阶微分后，公式变为下式(15)，连续体模型的方程式的2次微分项通过体积分变换为面积分，在边界面E上变为下式(16)。

(数式15)

${\∫}_V\frac{{\∂}^2\mathrm{\ψ}}{\∂x_i\∂x_j}dV={\∫}_S\frac{\∂\mathrm{\ψ}}{\∂x_j}{n}_{i}dS={\∫}_S\frac{\∂\mathrm{\ψ}}{\∂n}{n}_i{n}_{j}dS\mathrm{...}\left(15\right)$

(数式16)

${\∫}_{S_{ab}}\frac{\∂\mathrm{\ψ}}{\∂n_{ab}}\·n_{iab}\·n_{jab}dS\mathrm{...}\left(16\right)$

此外，在公式(15)中，表示法线方向微分，在公式(16)中，表示[n]_ab方向微分。

即，连续体模型的方程式的2次微分项通过体积分到面积分的变换，变为向物理量ψ的法线方向微分(S_ab的法线[n]_ab方向上的微分)乘以[n]的分量n_iab、n_jab的形式。

其中，公式(16)中的与下式(17)近似。

(数式17)

$\frac{\∂\mathrm{\ψ}}{\∂n_{ab}}=\frac{{\mathrm{\ψ}}_b-{\mathrm{\ψ}}_a}{r_{ab}\·n_{ab}}\mathrm{...}\left(17\right)$

此外，控制点a和控制点b的控制点间矢量[r]_ab根据控制点a的位置矢量[r]_a和控制点b的位置矢量[r]_b如下式(18)定义。

(数式18)

r_ab≡r_b-r_a…(18)

因此，边界面E的面积为Sab，所以公式(16)变为下式(19)，利用它可算出公式(16)。

(数式19)

${\∫}_{S_{ab}}\frac{\∂\mathrm{\ψ}}{\∂n_{ab}}\·n_{iab}\·n_{jab}dS=S_{ab}\·\frac{{\mathrm{\ψ}}_b-{\mathrm{\ψ}}_a}{r_{ab}\·n_{ab}}\·n_{iab}\·n_{jab}\mathrm{...}\left(19\right)$

此外，在导出公式(16)时可获知以下内容。

所有的线性偏微分方程式以对常数和1次、2次、其他偏导函数乘以系数的项的线性和来表示。在公式(15)到公式(18)中，将物理量ψ置换为ψ的1次偏导函数时，对较高次的偏导函数的体积分，可如公式(14)所示，通过低次的偏导函数的面积分求出。从低次的偏微分依次重复该步骤时，构成线性偏微分方程式的所有项的偏导函数可全部根据以下求出：控制点的物理量ψ；作为通过公式(12)或公式(13)计算的边界面上的ψ的ψ_ab；根据由公式(18)定义的控制点间矢量求出的控制点间距离；公式(5)所示的边界面E的面积S_ab、公式(16)所示的法线矢量的分量n_iab和n_jab。接着在非线性的偏微分方程式中，例如作为公式(20)所示的非线性项的ψ和ψ的一次偏导函数相乘的项、一次偏导函数的平方，可通过反复计算进行数值计算。即，将各自的项的ψ、一次偏导函数作为反复计算中的一次反复前的计算值近似，并反复计算即可。通过该方法，偏微分方程式中的非线性项均可数值计算。

以上特别说明了连续体模型的方程式，但可知对其以外的任意偏微分方程式均可进行无需Vertex和Connectivity的离散化。但对于守恒定律需要其他成立条件。对此稍后论述。

(数式20)

$\mathrm{\φ}\·\frac{\∂\mathrm{\φ}}{\∂x_i},{\left(\frac{\∂\mathrm{\φ}}{\∂x_i}\right)}^2\mathrm{...}\left(20\right)$

以上论述了在本数值解析方法中，在物理量计算时无需Vertex和Connectivity的情况。因此，在计算用数据模型做成(预处理)时，如不使用Vertex和Connectivity地求出控制量的体积、边界面的面积及法线矢量，则可使用公式(10)、(11)的离散化控制方程式，完全不使用控制量的几何学形状(即单元的几何学形状)地进行流速的计算。

但在本数值解析方法中，并非必须不使用控制量的具体几何学形状来求出控制量的体积、边界面的面积及法线矢量。即，在求解处理中不使用Vertex和Connectivity，因此即使使用控制量的具体几何学形状、具体而言使用Vertex和Connectivity，也不存在现有的有限要素法、有限体积法这样的分割区域相关的制约，即对分割区域的应变、弯曲的制约，因此如上所述，可容易地做成计算用数据模型。

在本数值解析方法中，根据条件不同，可将上述法线矢量置换为连接控制量之间的距离矢量。以下说明其理由。

图5中所示的边界面E的法线矢量[n]_ab朝向和连接控制点a和控制点b的距离矢量[r]_ab相同的方向时，法线矢量[n]可如下式(21)所示。

(数式21)

$n_{ab}=\frac{r_{ab}}{\left|r_{ab}\right|}\mathrm{...}\left(21\right)$

因此，边界面E的法线矢量[n]_ab朝向和距离矢量[r]_ab相同的方向时，通过向公式(10)、(11)所示的离散化控制方程式代入公式(21)，可获得法线矢量[n]_ab朝向控制点a、b间的方向时的离散化控制方程式。即，边界面和连接位于该边界面两侧的控制点的矢量正交时，可将离散化控制方程式中的法线矢量置换为距离矢量。

根据该离散化控制方程式，可仅根据控制点的位置坐标决定边界面E的法线矢量[n]_ab。

并且，通过使边界面E和距离矢量尽量接近正交，计算物理量时的精度提高。因此，通过将法线矢量置换为距离矢量，可提高计算精度。

进一步，法线矢量的任意性可在连接控制点的方向上固定化。

但是，因法线矢量的任意性消失，存在无法使边界面的姿态具有自由度的情况。此时，与法线矢量存在任意性时相比，做成计算用数据模型时，控制量的体积和边界面的设定中产生制约。并且，将法线矢量置换为连接控制点的矢量时，距离矢量变得必需，为计算该距离矢量，需要控制点的坐标，因此计算用数据模型中需要具有控制点的坐标数据或距离矢量数据。但这种情况下，在本数值解析方法中，在物理量计算中无需Vertex和Connectivity。

接着为了明确本发明和粒子法的不同点，说明物理量计算时满足物理量的守恒定律的条件。

例如如图6所示，考虑将L字形的流路分割为四个的单元R_a、R_b、R_c、R_d的情况。此外，配置在各单元R_a、R_b、R_c、R_d内部的控制点a、b、c、d设置在单元R_a、R_b、R_c、R_d的中心。并且，各边界面中的流速矢量相对边界面垂直。

此外在图6中，V_a表示单元R_a的体积(控制点a的控制量的体积)，V_b表示单元R_b的体积(控制点b的控制量的体积)，V_c表示单元R_c的体积(控制点c的控制量的体积)，V_d表示单元R_d的体积(控制点d的控制量的体积)，ρ_a表示单元R_a中的密度，ρ_b表示单元R_b中的密度，ρ_c表示单元R_c中的密度，ρ_d表示单元R_d中的密度，S_a表示单元R_a和外部区域的边界面的面积，S_c表示单元R_c和外部区域的边界面的面积，S_d表示单元R_d和外部区域的边界面的面积，S_ab表示单元R_a和单元R_b的边界面的面积，S_ac表示单元R_a和单元R_c的边界面的面积，S_bd表示单元R_b和单元R_d的边界面的面积，S_cd表示单元R_c和单元R_d的边界面的面积，u_a表示单元R_a和外部区域的边界面中的流速，u_c表示单元R_c和外部区域的边界面中的流速，u_d表示单元R_d和外部区域的边界面中的流速，u_ab表示单元R_a和单元R_b的边界面中的流速，u_ac表示单元R_a和单元R_c的边界面中的流速，u_bd表示单元R_b和单元R_d的边界面中的流速，u_cd表示单元R_c和单元R_d的边界面中的流速，ρ_a表示单元R_a中的密度，ρ_ab表示单元R_a和单元R_b的边界面中的密度，ρ_ac表示单元R_a和单元R_c的边界面中的密度，ρ_bd表示单元R_b和单元R_d的边界面中的密度。

考虑在图6所示的四个控制点a、b、c、d(四个单元R_a、R_b、R_c、R_d)上，使质量守恒的公式离散化。此外，使质量守恒的公式离散化的离散化控制方程式如下式(43)所示。

控制点a中的离散化控制方程式如下式(22)所示。

(数式22)

$V_a\·\frac{\∂{\mathrm{\ρ}}_a}{\∂t}+\[-S_a\·{\mathrm{\ρ}}_a\·u_a+S_{ab}\·{\mathrm{\ρ}}_{ab}\·u_{ab}+S_{ac}\·{\mathrm{\ρ}}_{ac}\·u_{ac}\]=\mathrm{0...}\left(22\right)$

控制点b中的离散化控制方程式如下式(23)所示。

(数式23)

$V_b\·\frac{\∂{\mathrm{\ρ}}_b}{\∂t}+\[S_{ab}\·{\mathrm{\ρ}}_{ab}\·u_{ab}+S_{bd}\·{\mathrm{\ρ}}_{bd}\·u_{bd}\]=\mathrm{0...}\left(23\right)$

控制点c中的离散化控制方程式如下式(24)所示。

(数式24)

$V_c\·\frac{\∂{\mathrm{\ρ}}_c}{\∂t}+\[-S_c\·{\mathrm{\ρ}}_c\·u_c-S_{ac}\·{\mathrm{\ρ}}_{ac}\·u_{ac}+S_{cd}\·{\mathrm{\ρ}}_{cd}\·u_{cd}\]=\mathrm{0...}\left(24\right)$

控制点d中的离散化控制方程式如下式(25)所示。

(数式25)

$V_d\·\frac{\∂{\mathrm{\ρ}}_d}{\∂t}+\[-S_{cd}\·{\mathrm{\ρ}}_{cd}\·u_{cd}-S_{bd}\·{\mathrm{\ρ}}_{bd}\·u_{bd}+S_d\·{\mathrm{\ρ}}_d\·u_d\]=\mathrm{0...}\left(25\right)$

并且，全部加上公式(22)～(25)，获得下式(26)。

(数式26)

$\[V_a\frac{\∂{\mathrm{\ρ}}_a}{\∂t}+V_b\frac{\∂{\mathrm{\ρ}}_b}{\∂t}+V_c\frac{\∂{\mathrm{\ρ}}_c}{\∂t}+V_d\frac{\∂{\mathrm{\ρ}}_d}{\∂t}\]+\[-S_a\·{\mathrm{\ρ}}_a\·u_a-S_c\·{\mathrm{\ρ}}_c\·u_c+S_d\·{\mathrm{\ρ}}_d\·u_d\]=\mathrm{0...}\left(26\right)$

利用各控制点a、b、c、d的控制量的体积V_a、V_b、V_c、V_d相对时间一定，导入到时间微分项中，公式(26)变为下式(27)。

(数式27)

$\frac{\∂}{\∂t}({\mathrm{\ρ}}_a{V}_a+{\mathrm{\ρ}}_b{V}_b+{\mathrm{\ρ}}_c{V}_c+{\mathrm{\ρ}}_d{V}_d)+\[-S_a\·{\mathrm{\ρ}}_a\·u_a-S_c\·{\mathrm{\ρ}}_c\·u_c+S_d\·{\mathrm{\ρ}}_d\·u_d\]=\mathrm{0...}\left(27\right)$

现在设各控制点a、b、c、d的控制量所占的全部体积为V_abcd、平均密度为ρ(上划线)_abcd时，全部体积V_abcd由下式(28)所示，平均密度ρ(上划线)_abcd由下式(29)所示。

(数式28)

V_abca≡V_a+V_b+V_c+V_d…(28)

(数式29)

${\stackrel{\‾}{\mathrm{\ρ}}}_{abcd}\≡\frac{{\mathrm{\ρ}}_a{V}_a+{\mathrm{\ρ}}_b{V}_b+{\mathrm{\ρ}}_c{V}_c+{\mathrm{\ρ}}_d{V}_d}{V_{abcd}}\mathrm{...}\left(29\right)$

因此，公式(27)如下式(30)所示。

(数式30)

$V_{abcd}\·\frac{\∂{\stackrel{\‾}{\mathrm{\ρ}}}_{abcd}}{\∂t}+\[-S_a\·{\mathrm{\ρ}}_a\·u_a-S_c\·{\mathrm{\ρ}}_c\·u_c+S_d\·{\mathrm{\ρ}}_d\·u_d\]=\mathrm{0...}\left(30\right)$

公式(30)意味着，控制点a、b、c、d的控制量所占的全部区域中流入的质量流束和流出的质量流束的差在单位时间下，与控制点a、b、c、d的控制量所占的全部区域的平均密度的时间变化(质量的时间变化)相等。即，按照各控制点a、b、c、d离散化的质量守恒的公式，对于全部控制点的控制量所占的区域也成立。

即，控制点所示的控制量区域的离散化控制方程式对所有控制点相加时，必须变为满足和作为计算对象的解析区域的全部区域相关的守恒定律的方程式。

接着，设控制点的全数为N，加上公式(43)所示的质量守恒的公式，获得下式(31)。

(数式31)

$\underset{a=1}{\overset{N}{\Σ}}\[V_a\frac{\∂{\mathrm{\ρ}}_a}{\∂t}+\underset{b=1}{\overset{m}{\Σ}}\{S_{ab}\·{\mathrm{\ρ}}_{ab}\·(n_{ab}\·u_{ab})\}\]=\mathrm{0...}\left(31\right)$

在公式(31)中，各控制点间的边界面的面积从控制点a侧观察时、从控制点b侧观察时均相等时，各控制点间的质量流束(ρ[n]·[u])·S在控制点a侧和控制点b侧正负相反，绝对值相等，因此相减为零，被消去。即，公式(31)表示，对计算的全部区域，流入的质量和流出的质量的差与全部区域中的质量的单位时间变化相等。因此，公式(31)变解析区域整体的质量守恒的公式。

因此，公式(31)为满足计算的全部区域的质量守恒定律，需要以下条件：二个控制点间的边界面的面积一致的条件；对法线矢量从一个控制点侧观察时和从另一个控制点侧观察时，绝对值一致的条件。

并且，为满足质量守恒定律，需要以下条件：下式(32)所示的全部控制点的控制量所占的体积与解析区域的全部体积一致。

(数式32)

$V_{total}=\underset{a=1}{\overset{N}{\Σ}}{V}_a\mathrm{...}\left(32\right)$

如果考虑到连续体的密度ρ由ρ₁＝ρ₂＝……＝ρ和一个变量表示则容易理解。

因此，为了满足质量守恒定律，需要以下条件：全部控制点的控制量的体积的总和与解析区域的体积一致。

此外，在此对质量守恒的公式进行了说明，守恒定律对于连续体的运动量、能量也必须成立。对这些物理量，在下述公式(50)、(55)中对全部控制点相加，从而满足守恒定律，为此需要以下条件：全部控制点的控制量所占的体积和解析区域的全部体积一致的条件；二个控制点间的边界面的面积一致的条件、及对法线矢量从一个控制点侧观察时和从另一个控制点侧观察时，绝对值一致(正负相反符号)的条件。

并且，为满足守恒定律，如图7所示，考虑控制点_a所占的控制量时，在考虑通过控制点a、具有任意方向的单位法线矢量[n]_p的无限大的投影面P时，需要下式(33)成立的条件。

(数式33)

$\underset{i=1}{\overset{m}{\Σ}}\[(n_i\·n_P)\·S_i\]=\mathrm{0...}\left(33\right)$

此外，在图7及公式(33)中，S_i表示边界面E_i的面积，[n]_i表示边界面E_i的单位法线矢量，m表示控制量的面的总数。

公式(33)表示构成控制量的多面体构成闭包空间。该公式(33)在构成控制量的多面体的一部分凹陷时也成立。

此外如图8所示，对于二维的三角形，公式(33)也成立。并且，当取设多面体的一个面为微小面dS、m为∞的极限时，变为下式(34)，对图9所示的闭包曲面体也成立。

(数式34)

∫_Sn·n_PdS＝0…(34)

公式(33)成立这一条件，是高斯的散度定理、公式(14)所示的广义格林定理成立所需的条件。

并且，广义格林定理是作为连续体离散化的基本的定理。因此，根据格林定理使体积分变化为面积分并离散化时，为满足守恒定律，需要公式(33)成立这一条件。

因此，使用上述计算用数据模型及物理量计算进行数值解析时，为满足物理量的守恒定律，需要以下三个条件。

(a)全部控制点的控制量的体积(全部分割区域的体积)的总和与解析区域的体积一致。

(b)二个控制点间的边界面的面积一致，及对法线矢量从一个控制点侧(位于边界面两侧的一个分割区域)观察时和从另一个控制点侧(位于边界面两侧的另一个分割区域)观察时，绝对值一致。

(c)考虑通过控制点(通过分割区域)、具有任意方向的单位法线矢量[n]_p的无限大的投影面P时，公式(33)成立。

即，满足守恒定律时，需要满足这些条件地做成计算用数据模型。但如上所述，在本数值解析方法中，在做成计算用数据模型时，可将单元形状任意变形，所以可容易地满足上述三个条件地做成计算用数据模型。

接着详细说明作为现有的粒子法的MPS(MovingParticleSemi-Implicit：移动粒子半隐)法无法满足守恒定律的理由、计算负荷增大的理由、及本数值解析方法相对粒子法的优越性。

MPS法是，检测出适当设定的半径r_e的球内存在的粒子，与这些粒子建立结合关系并进行计算的方法，例如如图10所示，粒子i周围存在多个粒子j时，粒子i中的拉普拉斯算子(Laplacian)(▽²ψ)_i如下式(35)所示近似。

(数式35)

${\left({\▿}^2\mathrm{\φ}\right)}_i=\frac{2d}{m}\·\underset{j\≠i}{\overset{m}{\Σ}}\[\frac{{\mathrm{\ψ}}_j-{\mathrm{\ψ}}_i}{{\left|r_{ij}\right|}^2}\·\mathrm{\ω}\left(r\right)\]\mathrm{...}\left(35\right)$

此外，在图10中，ψ_i表示粒子i中的物理量，ψ_j表示粒子j中的物理量，[r]_ij表示粒子i到粒子j的距离矢量。

并且，公式(35)中的d是表示维度数的常数，当是三维时变为3。并且，公式(35)中的ω(r)是权重函数，如下式(36)所示。并且，公式(35)中的m表示处于结合关系下的粒子个数。

(数式36)

$\mathrm{\ω}\left(r\right)=\left\{\begin{array}{cc}\frac{r_e}{e}-1& (0\≤r\≤r_e)\\ 0& (r_e\≤r)\end{array}\right\}\mathrm{...}\left(36\right)$

另一方面，如图11所示，粒子i、j作为控制点处理，粒子i的控制量为V_i、粒子i和粒子j之间的边界面的面积为S_ij、粒子i和粒子j之间的边界面的法线矢量为[n]_ij、粒子i到粒子j的距离矢量为[r]_ij时，粒子i中的拉普拉斯算子(▽²ψ)₁如下式(37)所示近似。

(数式37)

${\left({\▿}^2\mathrm{\ψ}\right)}_i=\frac{1}{V_i}\·\underset{j\≠i}{\overset{m}{\Σ}}\[\frac{{\mathrm{\ψ}}_j-{\mathrm{\ψ}}_i}{\left|r_{ij}\·n_{ij}\right|}\·S_{ij}\]\mathrm{...}\left(37\right)$

并且，比较公式(36)和公式(37)，如果[r]_ij和[n]_ij为相同方向时，可获得下式(38)、(39)。

(数式38)

$\frac{S_{ij}}{V_i}\·\frac{1}{\left|r_{ij}\·n_{ij}\right|}=\frac{2d}{m}\·\frac{\mathrm{\ω}\left(r\right)}{{\left|r_{ij}\right|}^2}\mathrm{...}\left(38\right)$

(数式39)

$\frac{S_{ij}}{V_i}=\frac{2d}{m}\·\frac{\mathrm{\ω}\left(r\right)}{\left|r_{ij}\right|}\mathrm{...}\left(39\right)$

公式(39)中，左边和右边的维度以(1/距离)一致。因此，右边的MPS法公式化可解释为，仅根据2个粒子i、j间的距离，算出下式(40)所示的(面积/体积)的量(即公式(9)定义的比率)的方法。

(数式40)

${\mathrm{\φ}}_{ij}=\frac{S_{ij}}{V_i}\mathrm{...}\left(40\right)$

但是，仅根据m个粒子的结合关系求出边界面的面积S_ij和体积V_i，关系式是不够的，无法决定具体的数值，最终仅决定公式(40)的比率。

因此，即使根据MPS法的离散化控制方程式求出了边界面的面积S_ij和体积V_i，也完全不能保证满足上述可满足守恒定律的条件(a)～(c)。这意味着，MPS法在守恒定律满足性上有较大的问题。

将数值解析适用于工学问题、尤其是机械设计问题、设备设计问题时，定量值(压力、温度、热量等)的评估是极为重要的，但当数值解析不满足守恒定律时，不能保证定量性。

即，在MPS法中，质量守恒、动量守恒、能量守恒这些守恒定律不能保证满足，不能保证定量性。

与之相对，根据本数值解析方法，可满足守恒定律，可保证定量性。

此外，MPS法如上所述，随着时间的变化粒子移动，因此例如届时需要进行检测出存在于上述半径r_e的球内的粒子的附近搜索处理，物理量计算中的计算负荷变大。

与之相对，在本数值解析方法中，控制量及控制点在时间变化时也不会移动。因此，预先知晓控制量、控制点的配置关系时，可不进行附近搜索处理地进行物理量计算。因此，可使物理量计算中的计算负荷和MPS法相比较小。此外，控制量及控制点的配置关系预先不知道的情况下，在最初进行一次确定控制量、控制点的配置关系的处理即可。

此外，在上述说明中，说明了由纳维-斯托克斯式及连续式根据加权残值法导出的离散化控制方程式的物理量的计算例，但在本数值解析方法中使用的离散化控制方程式不限于此。

即，只要是由各种方程式(质量守恒的方程式、动量守恒的方程式、能量守恒的方程式、对流扩散方程式、及波动方程式等)根据加权残值法导出，并仅使用不需要Vertex和Connectivity的量就可算出物理量的离散化控制方程式，均可用于本数值解析方法。

并且，根据该离散化控制方程式的特性，无需象现有的有限要素法、有限体积法中的网格，可进行无网格的计算。并且，即使在预处理中利用规定单元的几何学形状的Vertex和Connectivity，也不存在现有的有限要素法、有限体积法、体素法中的对网格的制约，因此可降低伴随着计算用数据模型的做成的作业负荷。

以下说明由质量保持的方程式、动量守恒的方程式、能量守恒的方程式、对流扩散方程式、及波动方程式，根据加权残值法，可导出仅使用不需要Vertex和Connectivity的量的离散化控制方程式的情况，即在本数值解析方法中可使用其他控制方程式的情况。

(1)质量守恒的方程式

欧拉(Euler)坐标系中的质量守恒的方程式如下式(41)所示，以微分形式表示。

(数式41)

$\frac{\∂\mathrm{\ρ}}{\∂t}+\frac{\∂{\mathrm{\ρu}}_i}{\∂x_i}=\mathrm{0...}\left(41\right)$

此外在公式(41)中，t表示时间，x_i(i＝1、2、3)表示笛卡尔系中的坐标，ρ表示密度，u_i(i＝1、2、3)表示变形速度的分量，下标i(i＝1、2、3)表示笛卡尔系坐标中的各方向分量。并且，下标i基于求和约定。

并且，根据加权残值法使公式(41)对控制点的控制点的体积V积分后，如下式(42)所示。

(数式42)

${\∫}_V\frac{\∂\mathrm{\ρ}}{\∂t}dV+{\∫}_S\mathrm{\ρ}n\·udS=\mathrm{0...}\left(42\right)$

进一步，对图5所示的控制点a离散化，变换为代数方程式后，如下式(43)所示。

(数式43)

$V_a\·\frac{\∂{\mathrm{\ρ}}_a}{\∂t}+\underset{b=1}{\overset{m}{\Σ}}\[S_{ab}\·{\mathrm{\ρ}}_{ab}\·(n_{ab}\·u_{ab})\]=\mathrm{0...}\left(43\right)$

其中，附加下标ab的ρ_ab、[n]_ab表示控制点a和控制点b之间的边界面E上的物理量。并且，m是和控制点a处于结合关系(位于边界面两侧的关系)的所有控制点的个数。

并且，公式(43)除以作为控制点a的控制量的体积的V_a后，获得下式(44)，进一步形成下式(45)，从而获得质量守恒的方程式离散化的下式(46)。

(数式44)

$\frac{\∂{\mathrm{\ρ}}_a}{\∂t}+\underset{b=1}{\overset{m}{\Σ}}\[\frac{S_{ab}}{V_a}\·{\mathrm{\ρ}}_{ab}\·(n_{ab}\·u_{ab})\]=\mathrm{0...}\left(44\right)$

(数式45)

${\mathrm{\φ}}_{ab}=\frac{S_{ab}}{V_a}\mathrm{...}\left(45\right)$

(数式46)

$\frac{\∂{\mathrm{\ρ}}_a}{\∂t}+\underset{b=1}{\overset{m}{\Σ}}\[{\mathrm{\φ}}_{ab}\·{\mathrm{\ρ}}_{ab}\·(n_{ab}\·u_{ab})\]=\mathrm{0...}\left(46\right)$

该公式(46)是根据加权残值法导出、并仅使用不需要Vertex和Connectivity的量的方程式，因此可作为本数值解析方法的离散化控制方程式使用。

此外，本数值解析方法中使用的离散化控制方程式，在根据加权残值法导出使用表示现有的几何学形状的量的方程式的过程中，暂时停留在中途而获得，对此以上进行了说明。

即，公式(46)是在对质量守恒方程式根据加权残值法导出使用Vertex等的方程式的过程中获得的。

其中，图12是表示二维三角形的单元的示意图。图12中的三角形a的面积、边长、法线矢量如下表所示。此外，下表中记号×表示叉积。

(表1)

如图12所示，当单元是二维三角形时，使质量守恒的方程式根据加权残值法成为使用Vertex等的离散化控制方程式时，变为下式(47)。

(数式47)

$\frac{\∂{\mathrm{\ρ}}_a}{\∂t}+\underset{b=1}{\overset{3}{\Σ}}\[\frac{\left|r_{bb+1}\right|}{\left|r_{12}\×r_{13}\right|}\·\frac{r_{bb+1}\×(r_b\×r_{b+1})}{|r_{bb+1}\×(r_b\×r_{b+1})|}\·\mathrm{\ρ}(u_a+u_b)\]=\mathrm{0...}\left(47\right)$

其中在公式(47)中，[r]_i是顶点(Vertex)_i的位置矢量(坐标)，记号×表示矢量的叉积。并且，使ρ_ab及ρ一定，使[r]_ij为[r]_i-[r]_j，[r]₄为[r]₁。

(2)动量守恒的方程式

在欧拉坐标系中的动量守恒的方程式如下式(48)所示，以微分形式表示。

(数式48)

$\frac{\∂}{\∂t}\left({\mathrm{\ρu}}_i\right)+\frac{\∂}{\∂x_j}\left({\mathrm{\ρu}}_j{u}_i\right)=\frac{\∂{\mathrm{\σ}}_{ij}}{\∂x_j}+{\mathrm{\ρf}}_i\mathrm{...}\left(48\right)$

此外在公式(48)中，σ_ij(i、j＝1、2、3)表示连续体内部应力，f_i(i＝1、2、3)表示作用于连续体的外力(例如重力)，其他量和公式(41)相同。并且，下标j基于求和约定。

该公式(47)是构造体、材料、流体等的应力场的基础方程式。

并且，将公式(47)根据加权残值法对控制点的控制量的体积V积分后，如下式(49)所示。

(数式49)

${\∫}_V\frac{\∂\mathrm{\ρ}}{\∂t}\left({\mathrm{\ρu}}_i\right)dV+{\∫}_S\mathrm{\ρ}(n\·u)u_{i}dS={\∫}_S{\mathrm{\σ}}_{ij}\·n_{j}dS+{\∫}_V{\mathrm{\ρf}}_i\·dV\mathrm{...}\left(49\right)$

进一步，对图5所示的控制点a离散化，变换为代数方程式后，如下式(50)所示。

(数式50)

$\begin{array}{c}{V}_a\·\frac{\∂{\mathrm{\ρ}}_a{u}_{ia}}{\∂t}+\underset{b=1}{\overset{m}{\Σ}}\[S_{ab}\·\mathrm{\ρ}\·(n_{ab}\·u_{ab})u_{iab}\]\\ =\underset{b=1}{\overset{m}{\Σ}}\[S_{ab}\·{\mathrm{\σ}}_{ijab}{n}_{jab}\]+V_a\·f_{ia}\·{\mathrm{\ρ}}_a\end{array}\mathrm{...}\left(50\right)$

用公式(50)除以作为控制点a的控制量的体积的V_a，进一步导入公式(45)中后，获得动量守恒的方程式离散化的下式(51)。

(数式51)

$\begin{array}{c}\frac{\∂{\mathrm{\ρ}}_a{u}_{ia}}{\∂t}+\underset{b=1}{\overset{m}{\Σ}}\[{\mathrm{\φ}}_{ab}\·\mathrm{\ρ}\·(n_{ab}\·u_{ab})u_{iab}\]\\ =\underset{b=1}{\overset{m}{\Σ}}\[{\mathrm{\φ}}_{ab}\·{\mathrm{\σ}}_{ijab}{n}_{jab}\]+f_{ia}\·\mathrm{\ρ}\end{array}\mathrm{...}\left(51\right)$

此外，在动量守恒的方程式中考虑到应力张量的对称性时，角动量守恒方程式也可和动量守恒方程式一样离散化。

(3)对流扩散方程式

某个物质C到连续体内的对流扩散现象如下式(52)的对流扩散方程式所示。

(数式52)

$\frac{\∂}{\∂t}\left(\mathrm{\ρ}C\right)+\frac{\∂}{\∂x_j}\left({\mathrm{\ρu}}_{j}C\right)=\frac{\∂}{\∂x_j}\[{\mathrm{\μ}}_C\frac{\∂C}{\∂x_j}\]+\mathrm{\ρ}\·q_C\mathrm{...}\left(52\right)$

并且在公式(52)中，C表示物质C的浓度，μ_c表示物质C的扩散系数，q_c表示物质C的源(汇)项，ρ表示连续体的密度，u_i表示连续体的变形速度。

并且，对公式(52)根据加权残值法积分，进一步离散化并变换为离散化控制方程式后，获得下式(53)。

(数式53)

$\begin{array}{c}\frac{\∂{\mathrm{\ρ}}_a{C}_a}{\∂t}+\underset{b=1}{\overset{m}{\Σ}}\[{\mathrm{\φ}}_{ab}\·\mathrm{\ρ}\·(n_{ab}\·u_{ab})\·C_{ab}\]\\ =\underset{b=1}{\overset{m}{\Σ}}\[{\mathrm{\φ}}_{ab}\·{\mathrm{\μ}}_C\·{\left(\frac{\∂C}{\∂n}\right)}_{ab}\]+V_a\·{\mathrm{\ρ}}_a\·q_{C_a}\end{array}\mathrm{...}\left(53\right)$

此外，C_a等的添加了下标a的量是控制点a的物理量，C_ab等添加了下标ab的量是控制点a和控制点b之间的边界面中的物理量。

(4)能量守恒的方程式

能量守恒定律分为热能守恒和动能守恒，但动能的守恒包含在上述动量守恒中，因此在此将热能守恒的方程式的通式用下式(54)表示。

(数式54)

$\frac{\∂\mathrm{\ρ}U}{\∂t}+\frac{\∂}{\∂x_j}({\mathrm{\ρu}}_i\·U)=-\frac{\∂q_j}{\∂x_j}+\mathrm{\ρ}\·r+{\mathrm{\σ}}_{ij}\·D_{ij}\mathrm{...}\left(54\right)$

此外在公式(54)中，U表示连续体的内部能量，q_i表示热流束矢量，r表示热源、热能的源(汇)项，σ_ij表示应力张量，D_ij表示变形速度张量。此外，σ_ij·D_ij的张量2重积的项称为应力功率。并且，下标i、j适用求和约定。

并且，将公式(54)根据加权残值法积分，进一步离散化并变换为离散化控制方程式后，获得下式(55)。

(数式55)

$\begin{array}{c}\frac{\∂{\mathrm{\ρ}}_a{U}_a}{\∂t}+\underset{b=1}{\overset{n}{\Σ}}\[{\mathrm{\φ}}_{ab}\·\mathrm{\ρ}\·(n_{ab}\·u_{ab})\·U_{ab}\]\\ =-\underset{b=1}{\overset{m}{\Σ}}\[{\mathrm{\φ}}_{ab}\·(n_{ab}\·q_{ab})\]+U_a\·{\mathrm{\ρ}}_a\·r_a+U_a\·{(\mathrm{\σ}:D)}_a\end{array}\mathrm{...}\left(55\right)$

并且，对热流束矢量[q]除了热流束外加上电气或化学能量等所有非力学的能量的流束时，能量守恒的方程式变为表示非常大范围的能量守恒的方程式。

(5)波动方程式

上述质量守恒的方程式、动量守恒扩方程式、能量守恒的方程式、对流扩散方程式等以守恒形式表示的物理定律的方程式，是同时具有称为“抛物型”和“椭圆型”的偏微分方程式的性质的方程式。与之相对，表示波的传送、振动的传送的波动方程式称为“双曲型”，通式如下式(56)所示。

(数式56)

$\frac{{\∂}^{2}u}{\∂t^2}={\mathrm{\α}}^2\·\[\frac{{\∂}^{2}u}{\∂x^2}+\frac{{\∂}^{2}u}{\∂y^2}+\frac{{\∂}^{2}u}{\∂z^2}\]\mathrm{...}\left(56\right)$

并且在公式(56)中，u表示振幅、位移、α表示波的传播速度。

并且，根据加权残值法对公式(56)积分，进一步离散化并变换为离散化控制方程式后，获得下式(57)。

(数式57)

$\frac{{\∂}^2{u}_a}{\∂t^2}=\underset{b=1}{\overset{m}{\Σ}}\[{\mathrm{\φ}}_{ab}\·{\mathrm{\α}}^2\·{\left(\frac{\∂u}{\∂n}\right)}_{ab}\]\mathrm{...}\left(57\right)$

由公式(56)可知，在空间方向上可直接适用不需要控制量的几何学形状的离散化方法。但在时间方向上是2阶的导函数，因此使用高精度的时间积分法。

上述公式(46)、(51)、(53)、(55)、(57)，是根据加权残值法导出并且仅使用无需Vertex和Connectivity的量而可算出物理量的方程式，因此可作为本数值解析方法的离散化控制方程式使用。

并且，本数值解析方法通过使用上述离散化控制方程式，可适用于常规及非常规的流体力学、热传导、对流扩散、构造力学、波动、及这些物理现象组成的现象中的数值解析。

其次，在使用本数值解析方法时，可在将通过不同坐标系设计的配件容易地组装的状态下，进行数值解析。这是因为，在本数值解析方法中，和有限要素法中使用的要素不同，不需要控制量的具体的几何学形状，并且二个邻接的控制点间的“距离”无需是绝对坐标系中的距离，只要是“计算上的距离”即可。

因此，根据本数值解析方法，如图13所示，在组装了以不同坐标系设计的配件A、B、C的状态下想进行数值解析时，可不使各配件的坐标系一致地进行数值解析。

此外，称为一般的有限要素法的数值解析方法，是图14、图15所示的完全不允许要素交叉的方法。因此，在有限要素法中使用的应用程序软件中检测出要素交叉时，形成错误输出。并且，因这一情况，在有限要素中做成计算用数据模型需要非常大的作业负担，对此之前已经论述。

另一方面，根据本数值解析方法，如图14、图15所示，控制点之间的结合中允许交叉。各控制点占据的控制量的体积(分割区域的体积)、边界面的面积、及边界面的法线这些信息量，不需要具体的几何学形状。因此，在控制点间的结合交叉时，也可执行物理量计算。

因此，做成计算用数据模型时的制约减少，计算用数据模型做成的自由度大幅上升。但考虑到流体等中的物理现象时，根据来自与一个控制点接近的控制点的信息不同，具有该控制点的物理性质被更新的性质，因此与远方的控制点结合时，计算精度可能恶化。因此在本数值解析方法中优选，在适当的范围内形成控制点间的结合。

附图说明

图1是表示现有的有限要素法下的计算用数据模型的一例的概念图。

图2是表示现有的有限体积法下的计算用数据模型的一例的概念图。

图3是用于比较本发明和现有的数值解析方法的表。

图4是详细比较本发明和现有的有限体积法的图。

图5是表示本数值解析方法的计算用数据模型的一例的概念图。

图6是表示用于说明在本数值解析方法中满足物理量的守恒定律的条件的多个分割区域的示意图。

图7是表示通过控制点、具有任意方向的单位法线矢量的无限大的投影面的示意图。

图8是说明考虑到二维的三角形的控制量时满足物理量的守恒定律的条件的示意图。

图9是说明考虑到球的控制量时满足物理量的守恒定律的条件的示意图。

图10是表示粒子法下的附近搜索的示意图。

图11是表示本数值解析方法中的控制量的示意图。

图12是表示二维三角形的单元的示意图。

图13是表示组成了以不同坐标系设计的配件的状态的示意图。

图14是表示要素交叉的形态的示意图。

图15是表示控制点之间的结合交叉的形态的示意图。

图16是概要表示本发明的一个实施方式中的数值解析装置的硬件构成的框图。

图17是表示本发明的一个实施方式中的数值解析方法的流程图。

图18是表示本发明的一个实施方式中的数值解析方法中进行的预处理的流程图。

图19是表示本发明的一个实施方式中的数值解析方法中进行的求解处理的流程图。

图20是表示对三维形状数据中含有的车辆的室内空间用微小的闭合曲面搭叠(lapping)的形态的示意图。

图21是表示通过同一正交格子形状的分割区域形成包括三维形状数据中含有的车辆的室内空间的解析区域的形态的示意图。

图22是表示删除了从室内空间溢出的分割区域的形态的示意图。

图23是表示将新的分割区域填充到在解析区域的边界面和室内空间的边界面之间形成的间隙中的形态的示意图。

图24是表示解析区域包括移动边界时的本发明的一个实施方式中的数值解析方法的流程图。

标号说明

A……数值解析装置(物理量计算装置)

P……数值解析程序

P1……预处理程序

P2……求解处理程序(物理量计算程序)

P3……后处理程序

M……计算用数据模型

V_i……粒子i虚拟所占体积

S_ij……粒子i和粒子j的边界面积

r_ij……i到j的距离矢量

n_ij……S_ij的法线矢量

1……CPU

2……存储装置

2a……程序存储部

2b……数据存储部

3……DVD驱动器

4……输入装置

4a……键盘

4b……鼠标

5……输出装置

5a……显示器

5b……打印机

具体实施方式

以下参照附图说明本发明涉及的物理量计算方法、数值解析方法、物理量计算程序、数值解析程序、物理量计算装置及数值解析装置。

并且在以下说明中，说明包括本发明涉及的物理量计算方法的数值解析方法、包括本发明涉及的物理量计算程序的数值解析程序、及包括本发明涉及的物理量计算装置的数值解析装置的实施方式。

并且，在以下实施方式中，说明通过数值解析求出车辆的室内空间中的空气流速的情况。

图16是概要表示本实施方式的数值解析装置A的硬件构成的框图。

如该图所示，本实施方式的数值解析装置A由个人计算机、工作站等计算机构成，具有CPU1、存储装置2、DVD(DigitalVersatileDisc)驱动器3、输入装置4、输出装置5、及通信装置6。

CPU1与存储装置2、DVD驱动器3、输入装置4、输出装置5、及通信装置6电连接，处理从这些各种装置输入的信号，并输出处理结果。

存储装置2由存储器等内部存储装置及硬盘驱动器等外部存储装置构成，存储从CPU1输入的信息，并且根据从CPU1输入的指令输出存储的信息。

并且，在本实施方式中，存储装置2具有程序存储部2a和数据存储部2b。

程序存储部2a存储数值解析程序P。该数值解析程序P是在规定的OS中执行的应用程序，使由计算机构成的本实施方式的数值解析装置A进行数值解析地作用。并且，数值解析程序P使本实施方式的数值解析装置A例如作为计算用数据模型做成单元和物理量计算单元作用。

并且如图16所示，数值解析程序P具有：预处理程序P1、求解处理程序P2、后处理程序P3。

预处理程序P1使用于执行求解处理的前处理(预处理)，由本实施方式的数值解析装置A执行，使本实施方式的数值解析装置A作为计算用数据模型做成单元作用，从而做成计算用数据模型。并且，预处理程序P1使本实施方式的数值解析装置A在执行求解处理时，执行必要的条件设定，进一步执行上述计算用数据模型、汇总了设定的条件的求解输入数据文件F的做成。

并且，预处理程序P1使本实施方式的数值解析装置A作为计算用数据模型做成单元作用时，首先对本实施方式的数值解析装置A，取得包括车辆的室内空间的三维形状数据，执行表示该取得的三维形状数据中含有的车辆的室内空间的解析区域的做成。

此外，稍后详述，在本实施方式中，在求解处理中，使用在利用了上述本发明的数值解析方法中说明的离散化控制方程式(仅使用不需要Vertex和Connectivity的量，并且根据加权残值法导出的离散化控制方程式)。因此，在做成计算用数据模型时，在满足守恒定律的条件下，可将分割区域的形状及解析区域的形状任意变更。因此，三维形状数据中含有的车辆的室内空间的修正或变更作业中，通过搭叠处理等简单的处理就变得充分。因此，在本实施方式中，预处理程序P1使本实施方式的数值解析装置A执行下述搭叠处理：对取得的三维形状数据中含有的车辆的室内空间以微小的闭合曲面搭叠，从而修缮车辆的室内空间中存在的孔、间隙。

之后，预处理程序P1使本实施方式的数值解析装置A形成同一尺寸的分割区域，并执行包括进行了搭叠处理的室内空间的全部区域的解析区域的做成。接着，预处理程序P1使本实施方式的数值解析装置A，在做成为同一尺寸的分割区域中切除从室内空间溢出的区域，从而执行表示室内空间的解析区域的做成。其中，在求解处理中使用上述离散控制方程式，因此可容易地切割解析区域中从室内空间溢出的区域。这样一来，不会象体素法那样与外部空间的边界变为阶梯状，并且对于体素法的切割单元法这样的外部空间的边界附近的解析区域的形成，无需伴随着经验积累、试验错误所需的非常庞大的手工作业的修正或处理。即，在本实施方式中，在体素中成为问题的与外部空间的边界处理不会产生问题。

此外，在本实施方式中，如下所述，通过向室内空间和切割的区域之间的间隙填充新的任意形状的分割区域，以不仅由正交格子形状形成的分割区域构成解析区域，进一步在解析区域中不重叠分割区域地填充。

并且，预处理程序P1使本实施方式的数值解析装置A作为计算用数据模型做成单元作用时，使本实施方式的数值解析装置A，执行将一个控制点虚拟配置到表示做成的室内空间的解析区域中含有的分割区域各自的内部的处理，存储控制点的配置信息、及各控制点占据的控制量的体积数据。

并且，预处理程序P1使本实施方式的数值解析装置A作为计算用数据模型做成单元作用时，使本实施方式的数值解析装置，执行作为上述分割区域之间的边界面的边界面面积及法线矢量的计算，存储这些边界面面积及法线矢量。

并且，预处理程序P1使本实施方式的数值解析装置A作为计算用数据模型做成单元作用时，做成控制量(控制点)的结合信息(Link)，并存储该Link。

并且，预处理程序P1使本实施方式的数值解析装置A，汇总上述各控制点占据的控制量的体积、边界面面积及法线矢量、控制点(即分割区域)的配置信息(坐标)、Link，做成计算用数据模型。

并且，预处理程序P1使本实施方式的数值解析装置A在执行上述求解处理时进行必要的条件设定的情况下，进行物性值的设定、边界条件的设定、初始条件的设定、计算条件的设定。

其中，物性值是指室内空间中的空气密度、粘性系数等。

边界条件用来规定控制点间的物理量的交换定律，在本实施方式中，是基于上述公式(10)所示的纳维-斯托克斯式的离散化控制方程式、及基于公式(11)所示的连续式的离散化控制方程式。

并且，边界条件中包括表示面向室内空间和外部空间的边界面的分割区域的信息。

初始条件表示执行求解处理时的最初的物理量，是各分割区域的流速的初始值。

计算条件是求解处理中的计算条件，例如是反复次数、收敛基准。

并且，预处理程序P1使本实施方式的数值解析装置A形成GUI(GraphicalUserInterface：图像化用户界面)。具体而言，预处理程序P1使输出装置5具有的显示器5a显示图像，并且成为能够通过输入装置4具有的键盘4a、鼠标4b进行操作的状态。

求解处理程序P2(物理量计算程序)用于使本实施方式的数值解析装置A执行求解处理，使本实施方式的数值解析装置A作为物理量计算装置作用。

并且，求解处理程序P2使本实施方式的数值解析装置A作为物理量计算单元作用时，使用包括计算用数据模型具有的控制量的体积和边界面的面积及法线矢量的求解输入数据文件F，来物理量计算解析区域中的物理量。

并且，求解处理程序P2使本实施方式的数值解析装置A作为物理量计算单元作用时，使本实施方式的数值解析装置A，执行求解输入数据文件F中含有的纳维-斯托克斯式及连续式的离散化系数矩阵的做成，并且执行矩阵形成用的数据表格的做成。

并且，求解处理程序P2使本实施方式的数值解析装置A作为物理量计算单元作用时，使本实施方式的数值解析装置A，根据基于上述公式(10)所示的纳维-斯托克斯式的离散化控制方程式、及基于上述公式(11)所示的连续式的离散化控制方程式，执行下式(58)所示的矩阵计算用的大型稀疏矩阵方程式的形成。

(数式58)

A·X＝B…(58)

并且在公式(58)中，[A]表示大型稀疏矩阵，[B]表示边界条件矢量，[X]表示流速的解。

并且，求解处理程序P2在上述离散化控制方程式中存在非压缩性等附加条件时，使本实施方式的数值解析装置A执行该附加条件向矩阵方程式的组合。

并且，求解处理程序P2使本实施方式的数值解析装置A，执行基于CG法(共轭梯度法)等的矩阵方程式的解的计算、该解的使用了下式(59)的解的更新、收敛条件的判断，取得最终的计算结果。

(数式59)

A(Xⁿ)·Xⁿ⁺¹＝B(Xⁿ)…(59)

后处理程序P3使本实施方式的数值解析装置A执行后处理，使本实施方式的数值解析装置A执行基于在求解处理中取得的计算结果的处理。

具体而言，后处理程序P3使本实施方式的数值解析装置A，执行计算结果的可视化处理、提取处理。

其中，可视化处理例如是，将截面轮廓显示、矢量显示、等值面显示、动画显示输出到输出装置5的处理。并且，提取处理是如下处理：提取作业人员指定的区域的定量值，作为数值、图表输出到输出装置5，或者提取作业人员指定的区域的定量值，输出文件化的内容。

并且，后处理程序P3使本实施方式的数值解析装置A，执行自动报告做成、计算残值的显示及分析。

数据存储部2b如图16所示存储：求解输入数据文件F，其具有计算用数据模型M、表示边界条件的边界条件数据D1、表示计算条件的计算条件数据D2、表示物性值的物性值数据D3、及表示初始条件的初始条件数据D4；三维形状数据D5；计算结果数据D6等。并且，数据存储部2b暂时存储在CPU1的处理过程中生成的中间数据。

DVD驱动器3构成为可取入DVD媒体X，根据从CPU1输入的指令，输出存储在DVD媒体X中的数据。并且在本实施方式中，DVD媒体X中存储有数值解析程序P，DVD驱动器3根据从CPU1输入的指令，输出DVD媒体X中存储的数值解析程序P。

输入装置4是本实施方式的数值解析装置A和作业人员的人机界面，具有作为指示设备的键盘4a、鼠标4b。

输出装置5使CPU1输入的信号可视化并输出，具有显示器5a及打印机5b。

通信装置6在本实施方式的数值解析装置A和CAD装置C等外部装置之间进行数据的收发，相对内部LAN(局域网)等网络B电连接。

接着参照图17～图19的流程图说明使用了上述构成的本实施方式的数值解析装置A的数值解析方法(本实施方式的数值解析方法)。

如图17的流程图所示，本实施方式的数值解析方法由预处理(步骤S1)、求解处理(步骤S2)、后处理(步骤S3)构成。

此外，在进行本实施方式的数值解析方法前，CPU1将取入到DVD驱动器3中的DVD媒体X中存储的数值解析程序P，从DVD媒体X取出，存储到存储装置2的程序存储部2a中。

并且，CPU1在从输入装置4输入了指示开始数值解析的信号时，根据存储装置2中存储的数值解析程序P执行数值解析。具体而言，CPU1根据程序存储部2a中存储的预处理程序P1执行预处理(步骤S1)，根据程序存储部2a中存储的求解处理程序P2执行求解处理(步骤S2)，根据程序存储部2a中存储的后处理程序P3执行后处理(步骤S3)。此外，通过CPU1执行基于预处理程序P1的处理(步骤S1)，本实施方式的数值解析装置A作为计算用数据模型做成单元作用。并且，通过CPU1执行基于求解处理程序P2的求解处理(步骤S2)，本实施方式的数值解析装置A作为物理量计算单元作用。

图18是表示预处理(步骤S1)的流程图。如该图所示，开始预处理(步骤S1)时，CPU1使通信装置6经由网络B从CAD装置C取得包括车辆的室内空间的三维形状数据D5(步骤S1a)。CPU1将取得的三维形状数据D5存储到存储装置2的数据存储部2b。

接着，CPU1进行在步骤S1a中取得的三维形状数据D5的分析(步骤S1b)，检测出数据存储部2b中存储的三维形状数据D5中含有的曲面的重叠、交叉的曲面、曲面间的间隙、小孔等。

接着，CPU1进行在步骤S1a中取得的三维形状数据D5的修正或变更处理(步骤S1c)。

具体而言，CPU1如图20所示，对三维形状数据D5中含有的室内空间K通过微小的闭合曲面进行搭叠处理，从而形成排除了曲面的重叠、交叉的曲面、曲面间的间隙、小孔等存在的室内空间的三维形状数据D5。

并且，CPU1在该修正或变更处理(步骤S1c)中，形成GUI，从GUI输入指令(例如表示搭叠区域的指令)时，执行反映了该指令的修正或变更处理。

接着，CPU1执行计算用数据模型的做成(步骤S1d)。

具体而言，首先，CPU1如图21所示，根据在步骤S1c中修正或变更的三维形状数据D5，执行含有室内空间的全部区域并且以同一形状(正交格子)的分割区域均等分割的解析区域K1的做成。其中，为缩短时间，解析区域K1以同一形状的分割区域分割，但构成解析区域K1的分割区域不一定是同一形状，可以是任意形状。

接着，CPU1如图22所示，通过删除从室内空间K溢出的分割区域，使解析区域K不溢出地收容到室内空间K。其结果如图22的放大图所示，在解析区域K1的边界面K1F和室内空间K的边界面KF之间，形成间隙Y。

接着，CPU1如图23所示，将作为解析区域的一部分的新的分割区域RA填充到间隙Y。

在本实施方式的数值解析方法中，如在上述原理说明中所详述，做成不具有Vertex和Connectivity的计算用数据模型，因此可不对分割区域的几何学形状形成制约地做成计算用数据模型。即，在做成计算用数据模型时，构成解析区域K1的分割区域可以是任意的形状。因此，CPU1将新的分割区域RA填充到间隙Y时，如图23的放大图所示，可任意设定分割区域RA的形状。因此，可极容易地将间隙Y用分割区域RA填充，例如即使作业人员不通过GUI进行作业，也可充分自动地形成分割区域RA。假如在现有的有限体积法中用分割区域填充间隙Y时，如上所述，在对分割区域的几何学形状的制约下，需要不产生允许范围外的应变、弯曲地配置分割区域。该作业变为作业人员的手工作业，结果强化了作业人员的庞大负担，并且导致解析时间的长期化。

接着，CPU1在表示室内空间的解析区域中含有的各分割区域内，虚拟配置一个控制点。其中，CPU1算出分割区域的重心，对各自的重心虚拟配置一个控制点。并且，CPU1算出控制点的配置信息、各控制点占据的控制量的体积(配置控制点的分割区域的体积)，暂时存储到存储装置2的数据存储部2b中。

并且，CPU1算出作为分割区域之间的边界面的边界面面积及法线矢量，将这些边界面面积及法线矢量暂时存储到存储装置2的数据存储部2b中。

并且，CPU1做成Link，将该Link暂时存储到存储装置2的数据存储部2b中。

并且，CPU1使数据存储部2b中存储的控制点的配置信息、各控制点占据的控制量的体积、边界面面积及法线矢量、及Link数据库化，从而做成计算用数据模型M，将做成的计算用数据模型存储到存储装置2的数据存储部2b内。

此外，在该步骤S1d中，将包括室内空间的解析区域以同一形状的分割区域均等地分割，进一步删除从室内空间溢出的分割区域，向产生的解析区域和室内空间的间隙填充新的分割区域，从而做成最终的解析区域K2。因此，室内空间的全部区域变为由不重叠的分割区域填充的状态。

因此，计算用数据模型满足上述用于满足守恒定律的三个条件(a)～(c)。

并且在本实施方式中，在步骤S1d中采用如下构成：先形成分割区域，之后配置控制点，并对各控制点分配配置了自身的分割区域的体积。

但是在本发明中也可是，先将控制点配置到解析区域，对各控制点稍后分配体积。

具体而言，例如根据和不同控制点接触为止的半径、处于结合关系的(以Link建立关联的)控制点为止的距离，对各控制点进行加权。

在此设控制点i的权重为w_i、基准体积为V⁺，分配到控制点i的体积V_i如下式(60)所示。

(数式60)

V_i＝w_i·V⁺…(60)

进一步，各控制点的体积V_i的总和与解析区域的体积V_total相等，因此下式(61)成立。

(数式61)

$\underset{i}{\Σ}{V}_i=V^{+}\·\underset{i}{\Σ}{w}_i=V_{total}\mathrm{...}\left(61\right)$

其结果是，基准体积V⁺可通过下式(62)求出。

(数式62)

$V^{+}=V_{total}/\underset{i}{\Σ}{w}_i\mathrm{...}\left(62\right)$

因此，分配到各控制点的体积可由公式(61)、(62)求出。

如使用该方法，则在预处理中，可不使用Vertex和Connectivity地求出计算用数据模型中具有的分割区域的体积。

并且，在该计算用数据模型的做成(步骤S1d)中，CPU1形成GUI，从GUI输入了指令(例如表示分割区域的密度的指令、表示分割区域的形状的指令)时，执行反映了该指令的处理。因此，作业人员通过操作GUI，可任意调节控制点的配置、分割区域的形状。

CPU1对照满足数值解析程序中存储的守恒定律的三个条件，在从GUI输入的指令脱离该条件时，将这一消息显示到显示器5a中。

接着，CPU1进行物性值数据的设定(步骤S1e)。具体而言，CPU1使用GUI，在显示器5a上显示物性值的输入画面，将从键盘4a或鼠标4b输入的表示物性值的信号作为物性值数据D3暂时存储到数据存储部2b，从而进行物性值的设定。此外，这里所述的物性值是室内空间中的流体(即空气)的特性值，是空气密度、粘性系数等。

接着，CPU1进行边界条件数据的设定(步骤S1f)。具体而言，CPU1使用GUI在显示器5a上显示边界条件的输入画面，将从键盘4a或鼠标4b输入的表示边界条件的信号，作为边界条件数据D1，暂时存储到数据存储部2b，从而进行边界条件数据的设定。此外，这里所述的边界条件表示：控制室内空间的物理现象的离散化控制方程式、面向室内空间和外部空间的边界面的控制点的特定信息、及室内空间和外部空间之间的热的导热条件等。

此外，在本实施方式的数值解析方法中，其目的在于通过数值解析求出室内空间中的流速，因此作为上述离散化控制方程式，使用基于上述纳维-斯托克斯式的离散化控制方程式(10)及基于连续式的离散化控制方程式(11)。

此外，这些离散化控制方程式例如通过将数值解析程序P中预先存储的多个离散化控制方程式显示到显示器5a上、并由作业人员通过使用键盘4a、鼠标4b而从该多个离散化控制方程式中选择。

接着，CPU1进行初始条件数据的设定(步骤S1g)。具体而言，CPU1使用GUI在显示器5a上显示初始条件的输入画面，将从键盘4a或鼠标4b输入的表示初始条件的信号作为初始条件数据D4，暂时存储到数据存储部2b，从而进行初始条件数据的设定。此外，这里所述的初始条件是指各控制点(各分割区域)中的初始流速。

接着，CPU1使用GUI进行计算条件数据的设定(步骤S1h)。具体而言，CPU1使用GUI在显示器5a上显示计算条件的输入画面，将从键盘4a或鼠标4b输入的表示计算条件的信号作为计算条件数据D2，暂时存储到数据存储部2b，从而进行计算条件数据的设定。此外，这里所述的计算条件是求解处理(步骤S2)中的计算条件，例如表示反复次数、收敛基准。

接着，CPU1进行求解输入数据文件F的做成(步骤S1i)。

具体而言，CPU将在步骤S1d中做成的计算用数据模型M、步骤S1e中设定的物性值数据D3、步骤S1f中设定的边界条件数据D1、步骤S1g中设定的初始条件数据D4、步骤S1h中设定的计算条件数据D2收容到求解输入数据文件F，从而做成求解输入数据文件F。此外，该求解输入数据文件F存储到数据存储部2b中。

以上预处理(步骤S1)完成后，CPU1根据求解处理程序P2，执行图19的流程图所示的求解处理(步骤S2)。

如图19所示，当开始求解处理(步骤S2)时，CPU1取得通过预处理(步骤S1)做成的求解输入数据文件F(步骤S2a)。此外，如本实施方式所示的数值解析方法所示，通过单一装置(本实施方式的数值解析装置A)执行预处理及求解处理时，已经在数据存储部2b中存储了求解输入数据文件F，因此可省略步骤S2a。但当预处理(步骤S1)和求解处理(步骤S2)在不同的装置中执行时，需要取得通过网络、可移动磁盘传送的求解输入数据文件F，因此需要进行该步骤S2a。

接着，CPU1判断求解输入数据的匹配性(步骤S2b)。此外，求解输入数据表示求解输入数据文件F中收容的数据，是计算用数据模型M、边界条件数据D1、计算条件数据D2、物性值数据D3及初始条件数据D4。

具体而言，CPU1在求解处理中分析可执行物理量计算的求解输入数据是否全部收容在求解输入数据文件F中，从而进行求解输入数据的匹配性的判断。

并且，CPU1判断求解输入数据不匹配时，在显示器5a中显示错误(步骤S2b)，进一步显示用于输入不匹配的部分的数据的画面。之后，CPU11根据从GUI输入的信号，进行求解输入数据的调整(步骤S2c)，再次执行步骤S2a。

另一方面，CPU1在步骤S2b中判断求解输入数据有匹配性时，执行初始计算处理(步骤S2e)。

具体而言，CPU1根据边界条件数据D1中存储的离散化控制方程式(即基于纳维-斯托克斯式的离散化控制方程式(10)、及基于连续式的离散化控制方程式(11))，做成离散化系数矩阵，进一步通过进行矩阵计算用的数据表格的做成，进行初始计算处理。

接着，CPU1进行大型稀疏矩阵方程式的形成(步骤S2f)。具体而言，CPU1根据基于纳维-斯托克斯式的离散化控制方程式(10)、及基于连续式的离散化控制方程式(11)，进行由上述公式(58)所示的矩阵计算用的大型稀疏矩阵方程式的形成。

接着，CPU11判断离散化控制方程式中是否存在非压缩性、接触等附加条件。该附加条件例如作为边界条件数据收容在求解输入数据文件F中。

并且，CPU1判断离散化控制方程式中存在附加条件时，在执行了该附加条件的大型矩阵方程式的形成(步骤S2h)后，执行大型矩阵方程式的计算(步骤S2i)。另一方面，CPU1判断离散化控制方程式中不存在附加条件时，不执行附加条件的大型矩阵方程式的形成(步骤S2h)，执行大型矩阵方程式的计算(步骤S2i)。

并且，CPU1例如通过CG法(共轭坡度法)解出大型矩阵方程式，使用上述公式(59)进行解的更新(步骤S2j)。

接着，CPU1进行公式(59)的残值是否达到收敛条件的判断(步骤S2g)。具体而言，CPU1计算公式(59)的残值，与计算条件数据D2中含有的收敛条件比较，从而进行公式(59)的残值是否达到收敛条件的判断。

并且，当判断残值未达到收敛条件时，CPU1在进行了物性值的更新后，再次执行步骤S2g。即，CPU1在公式(59)的残值达到收敛条件为止，进行物性值的更新的同时，重复进行步骤S2f～S2g。

另一方面，当判断残值达到收敛条件时，CPU1进行计算结果的取得(步骤S21)。具体而言，CPU1将在之前的步骤S2i中计算的物理量的解，作为计算结果数据存储到数据存储部2b，从而进行计算结果的取得。

通过该求解处理(步骤S2)，求出室内空间的空气流速。此外，该求解处理(步骤S2)相当于本发明的物理量计算方法。

上述求解处理(步骤S2)完成后，CPU1根据后处理程序P3执行后处理(步骤S3)。

具体而言，例如，CPU1根据从GUI输入的指令，由计算结果数据生成截面轮廓数据、矢量数据、等值面数据、动画数据，将该数据可视化于输出装置5中。

并且，CPU1根据从GUI输入的指令，提取室内空间的一部分中的定量值(计算结果)作为数值、图表，使该数据、图表在输出装置5中可视化，进一步将数值、图表作为文件汇总输出。并且，CPU1根据从GUI输入的指令，例如根据计算结果数据进行自动报告做成、计算残值的显示及分析，并输出其结果。

根据上述本实施方式的数值解析装置A、数值解析方法及数值解析程序，在预处理中做成具有控制量的体积、边界面的面积及法线矢量的计算用数据模型M，在求解处理中，使用计算用数据模型M中含有的控制量的体积、边界面的面积及法线矢量，计算各控制量中的物理量。

即，如在上述数值解析方法中所说明的，根据本实施方式的数值解析装置A、数值解析方法及数值解析程序，可不做成具有Vertex和Connectivity的计算用数据模型地进行数值解析。因此，对三维形状数据的修正或变更作业的限制也大幅降低，与具有Vertex和Connectivity的计算用数据模型相比，可极为容易地做成计算用数据模型M。因此，可减轻计算用数据模型M的做成中的作业负担。

并且，在本实施方式的数值解析装置A、数值解析方法及数值解析程序中，和现有的数值解析方法不同，在求解处理中，无需使用Vertex和Connectivity算出控制量的体积及边界面的面积及法线矢量。因此，可降低求解处理中的计算负荷。

因此，根据本实施方式的数值解析装置A、数值解析方法及数值解析程序，可减轻计算用数据模型的做成中的作业负担，并降低求解处理中的计算负荷。

并且，在本实施方式的数值解析装置A、数值解析方法及数值解析程序中，将分割区域不重叠地填充到解析区域。因此，满足了上述用于满足守恒定律的三个条件(a)～(c)，可满足守恒定律地计算流速。

此外，本实施方式中的计算用数据模型M可根据现有的有限体积法、有限要素法中使用的网格数据、现有的粒子法中使用的粒子数据、或单纯的点群数据，容易地做成计算用数据模型。这种情况下，也无需体素法中的与外部空间的边界区域中的分割区域的修正。

例如，由网格数据做成计算用数据模型时，将各要素作为本实施方式中的分割区域(控制点占据的控制量)处理，从而可容易地进行。并且，由粒子数据做成计算用数据模型时，将包围配置在解析区域中的控制点的闭合空间填充到解析区域而获得的闭合空间，作为本实施方式中的分割区域(控制点占据的控制量)处理，从而可容易地进行。

并且，根据本实施方式的数值解析装置A、数值解析方法及数值解析程序，如上所述，预处理中的计算用数据模型的作业负担大幅减少，并且可降低求解处理中的计算负荷。

因此，当解析区域的形状以时间序列变化时，即解析区域包括移动边界时，根据本发明，如图24的流程图所示，每当解析区域形状变化时，通过反复进行预处理和求解处理，可在现实的时间内进行物理量的计算。

并且，在本实施方式中，在物理量计算中算出边界面上的流速时，可将对位于该边界面两侧的控制点中的流速根据从该控制点到边界面为止的距离的比率进行加权平均而得到的值，作为边界面上的流速。

此时，与边界面上的流速单纯是位于该边界面两侧的控制点的流速的平均时相比，可更正确地求出室内空间的流速。但这种情况下，计算用数据模型作为数据需要具有表示边界面存在于连接控制点和控制点的线段的哪个内分点的比率α。因此，在预处理中做成计算用数据模型M，其进一步含有配置在相邻的控制量的内部的各控制点、到位于该控制点之间的边界面为止的距离的比率，并根据该比率，在求解处理中计算边界面中的流速。

此外，在物理量计算中，当连接控制点之间的矢量和位于该控制点之间的边界面正交时，可将法线矢量置换为连接控制点之间的距离矢量的单位矢量，并进行物理量计算。

但这种情况下，计算用数据模型中需要具有表示上述距离矢量或用于算出该距离矢量的控制点的位置坐标的数据。因此，在预处理中做成计算用数据模型，其进一步含有连接配置在相邻的控制量内部的各连接点的距离矢量或用于算出该距离矢量的控制点坐标，在距离矢量与位于该控制点之间的边界面正交时，将距离矢量置换为边界面的法线矢量，并计算流速。

以上参照附图的同时说明了本发明的优选实施方式，但本发明当然不限于上述实施方式。在上述实施方式中所示的各构成部件的形状、组合等仅是一例，在不脱离本发明主旨的范围内，可根据设计要求等进行各种变更。

在上述实施方式中，说明了使用由作为动量守恒方程式的变形例的纳维-斯托克斯式及连续式导出的离散化控制方程式，并通过数值解析求出空气流速的构成。

但本发明不限于此，可使用由质量守恒的方程式、动量守恒的方程式、角动量守恒的方程式、能量守恒的方程式、对流扩散方程式及波动方程式的至少任意一个导出的离散化控制方程式，通过数值解析求出物理量。

并且，在上述实施方式中，作为本发明的边界面特性量，说明了使用边界面的面积和边界面的法线矢量的构成。

但是本发明不限于此，作为边界面特性量也想中使用其他量(例如边界面的周长)。

并且在上述实施方式中，说明了为满足守恒定律，满足上述三个条件地做成计算用数据模型的构成。

但本发明不限于此，无需满足守恒定律时，无需满足上述三个条件地做成计算用数据模型。

并且在上述实施方式中，说明了将分割区域的体积作为配置在该分割区域内部的控制点所占的控制量的体积来处理的构成。

但是本发明不限于此，对分割区域的内部不一定配置控制点。此时，可通过将控制点占据的控制量的体积置换为分割区域的体积，来进行数值解析。

并且，在上述实施方式中，说明了数值解析程序P存储在DVD媒体X中并可传送的构成。

但是本发明不限于此，也可采用将数值解析程序P存储到其他可移动介质并传送的构成。

并且，也可将预处理程序P1和求解处理程序P2存储到不同的可移动介质中并传送。并且，数值解析程序P可通过网络传送。

并且在上述实施方式中，说明了边界面是平面的构成。但本发明不限于此，边界面也可是曲面。

对本申请详细且参照特定的实施方式进行了说明，对本领域技术人员而言，在不脱离本发明的精神和范围的情况下可进行各种变更、修正，这是不言而喻的。

本申请基于2009年6月25日申请的日本专利申请(特愿2009-150945)，其内容作为参照援引入其中。

Claims (18)

1.一种物理量计算方法，在对物理现象进行数值性解析的数值解析方法中计算物理量，其特征在于，

包括物理量计算步骤，计算被分割为多个分割区域的解析区域中的物理量，

在该物理量计算步骤中，使用离散化的控制方程式和计算用数据模型计算上述物理量，上述控制方程式以仅使用不需要上述分割区域的顶点的坐标Vertex及该顶点的连接信息Connectivity的量的方式根据加权残值法导出，在上述计算用数据模型中，具有各上述分割区域的体积及表示相邻的上述分割区域之间的边界面的特性的边界面特性量，来作为不需要上述分割区域的顶点的坐标Vertex及该顶点的连接信息Connectivity的量。

2.根据权利要求1所述的物理量计算方法，其中，上述控制方程式根据质量守恒的方程式、动量守恒的方程式、能量守恒的方程式、对流扩散方程式及波动方程式导出。

3.根据权利要求1或2所述的物理量计算方法，其中，上述边界面特性量是上述边界面的面积和上述边界面的法线矢量。

4.一种数值解析方法，具有以下处理：做成计算用数据模型的预处理；和求解处理，通过基于该计算用数据模型及控制方程式的物理量计算方法，计算解析区域中的物理量，

所述数值解析方法的特征在于，

在上述求解处理中，使用权利要求1～3的任意一项所述的物理量计算方法。

5.根据权利要求4所述的数值解析方法，其中，

在上述预处理中，以满足下述条件的方式形成上述分割区域：

全部分割区域的体积的总和与解析区域的体积一致的条件；

上述边界面的面积一致的条件；

在从位于该边界面两侧的一个分割区域观察时和从另一个分割区域观察时法线矢量的绝对值一致的条件；和

通过上述分割区域的无限大的投影面P的任意方向的单位法线矢量为[n]_p、边界面的面积为S_i、边界面的单位法线矢量为[n]_i、分割区域的面的总数为m、上述[n]_i和[n]_p分别对应于下式(1)中的n_i和n_p时，下式(1)成立的条件：

$\underset{i=1}{\overset{m}{\Σ}}\[(n_i\·n_P)\·S_i\]=0...\left(1\right).$

6.根据权利要求5所述的数值解析方法，其中，通过对上述解析区域不重叠地填充上述分割区域，而满足上述条件。

7.根据权利要求4～6的任意一项所述的数值解析方法，其中，通过上述预处理做成上述计算用数据模型，该计算用数据模型还包括从配置在相邻的上述分割区域内部的控制点到位于该分割区域之间的上述边界面为止的距离的比率，使用上述比率，通过上述求解处理计算上述边界面中的物理量。

8.根据权利要求4～6的任意一项所述的数值解析方法，其中，通过上述预处理做成上述计算用数据模型，该计算用数据模型还包括连接配置在相邻的上述分割区域内部的各控制点的距离矢量或用于算出该距离矢量的上述控制点的坐标，当上述距离矢量与位于该控制点之间的上述边界面正交时，将上述距离矢量作为上述边界面特性量之一来计算上述物理量。

9.根据权利要求4～6的任意一项所述的数值解析方法，其中，当上述解析区域包括移动边界时，每当上述移动边界移动时，进行上述预处理及上述求解处理。

10.一种物理量计算装置，在对物理现象进行数值性解析的数值解析方法中计算物理量，其特征在于，

包括物理量计算模块，计算被分割为多个分割区域的解析区域中的物理量，

该物理量计算模块使用离散化的控制方程式和计算用数据模型计算上述物理量，上述控制方程式以仅使用不需要上述分割区域的顶点的坐标Vertex及该顶点的连接信息Connectivity的量的方式根据加权残值法导出，在上述计算用数据模型中，具有各上述分割区域的体积及表示相邻的上述分割区域之间的边界面的特性的边界面特性量，来作为不需要上述分割区域的顶点的坐标Vertex及该顶点的连接信息Connectivity的量。

11.根据权利要求10所述的物理量计算装置，其中，上述控制方程式根据质量守恒的方程式、动量守恒的方程式、能量守恒的方程式、对流扩散方程式及波动方程式导出。

12.根据权利要求10或11所述的物理量计算装置，其中，上述边界面特性量是上述边界面的面积和上述边界面的法线矢量。

13.一种数值解析装置，包括：做成计算用数据模型的预处理模块；和求解处理模块，通过基于该计算用数据模型及控制方程式的物理量计算方法，计算解析区域中的物理量，

所述数值解析装置的特征在于，

上述求解处理模块使用权利要求1～3的任意一项所述的物理量计算方法。

14.根据权利要求13所述的数值解析装置，其中，

上述预处理模块以满足下述条件的方式形成上述分割区域：

全部分割区域的体积的总和与解析区域的体积一致的条件；

上述边界面的面积一致的条件；

在从位于该边界面两侧的一个分割区域观察时和从另一个分割区域观察时法线矢量的绝对值一致的条件；和

通过上述分割区域的无限大的投影面P的任意方向的单位法线矢量为[n]_p、边界面的面积为Si、边界面的单位法线矢量为[n]_i、分割区域的面的总数为m、上述[n]_i和[n]_p分别对应于下式(1)中的n_i和n_p时，下式(1)成立的条件：

$\underset{i=1}{\overset{m}{\Σ}}\[(n_i\·n_P)\·S_i\]=0...\left(1\right).$

15.根据权利要求14所述的数值解析装置，其中，通过对上述解析区域不重叠地填充上述分割区域，而满足上述条件。

16.根据权利要求13～15的任意一项所述的数值解析装置，其中，通过上述预处理模块做成上述计算用数据模型，该计算用数据模型还包括从配置在相邻的上述分割区域内部的控制点到位于该分割区域之间的上述边界面为止的距离的比率，使用上述比率，通过上述求解处理模块计算上述边界面中的物理量。

17.根据权利要求13～15的任意一项所述的数值解析装置，其中，通过上述预处理模块做成上述计算用数据模型，该计算用数据模型还包括连接配置在相邻的上述分割区域内部的各控制点的距离矢量或用于算出该距离矢量的上述控制点的坐标，当上述距离矢量与位于该控制点之间的上述边界面正交时，将上述距离矢量作为上述边界面特性量之一来计算上述物理量。

18.根据权利要求13～15的任意一项所述的数值解析装置，其中，当上述解析区域包括移动边界时，每当上述移动边界移动时，进行预处理及求解处理。