JPWO2010150758A1 - 物理量計算方法、数値解析方法、物理量計算プログラム、数値解析プログラム、物理量計算装置及び数値解析装置 - Google Patents

Abstract

物理現象を数値的に解析する数値解析方法において物理量を計算する物理量計算方法は、複数の分割領域に分割された解析領域における物理量を計算する物理量計算工程を含み、該物理量計算工程にて、前記分割領域の頂点の座標（Vertex）及び該頂点の連結情報（Connectivity）を必要としない量のみを使用すると共に重み付き残差積分法に基づいて導出された離散化された支配方程式と、各前記分割領域の体積及び隣り合う前記分割領域同士の境界面の特性を示す境界面特性量を、前記分割領域の頂点の座標（Vertex）及び該頂点の連結情報（Connectivity）を必要としない量として有する計算用データモデルと、を用いて前記物理量を計算することを特徴とする。

Description

本発明は、物理量計算方法、数値解析方法、物理量計算プログラム、数値解析プログラム、物理量計算装置及び数値解析装置に関するものである。

従来から、流速分布、応力分布及び熱分布等を数値解析によって求めるための数値解析手法として、例えば、有限要素法、有限体積法、ボクセル法及び粒子法が知られている。
このような数値解析手法は、一般的にプリ処理とソルバ処理とポスト処理とから構成されている。そして、プリ処理の中で計算用データモデルを作成し、ソルバ処理にて当該計算用データモデル及び離散化された支配方程式（以下、離散化支配方程式と呼ぶ）を用いて上記物理量の計算を行っている。

従来の有限体積法は、例えば、解析領域を複数の領域に分割し、各分割領域の体積、隣り合う分割領域の境界面の面積及び当該境界面の法線ベクトルを用いて各分割領域における物理量を計算する。
有限体積法では、プリ処理にて各分割領域の頂点の座標（Vertex）を含む計算用データモデル（通常、メッシュと呼ばれる）を作成し、ソルバ処理にて当該計算用データモデルに含まれるVertex等を用いて前述の分割領域の体積、境界面の面積及び境界面の法線ベクトルを算出し、これらの値を用いて物理量の計算を行う。Vertexは分割領域の幾何学的形状を規定するための量である。よって、有限体積法においては、ソルバ処理にて分割領域の幾何学的形状を用いて、分割領域の体積、境界面の面積及び境界面の法線ベクトルの算出を行っていると言える。
さらに、有限体積法では、隣り合う分割領域での頂点共有の条件が部分的に満足されていない部分を有してもよい。このため、有限体積法では、分割領域に対する制約が若干緩和される場合があるが、利用する解析要素タイプは例えばテトラ要素、ヘキサ要素、プリズム要素、ピラミッド要素などに限定される。
なお、特許文献１に示すように、解析要素タイプを限定しない有限体積法も提案されている。ただし、このような解析要素タイプを限定しない有限体積法であっても、前述した従来の有限体積法と同様に、プリ処理にて各分割領域の頂点の座標（Vertex）を含む計算用データモデルを作成し、ソルバ処理にて当該計算用データモデルに含まれるVertex等を用いて物理量の計算を行っている。

また、有限要素法は、周知のように補間関数を用いて各分割領域における物理量を算出する方法であるが、有限体積法と同様に、ソルバ処理にてVertex等によって規定される分割領域の幾何学的形状を用いている。

ボクセル法及び粒子法は、有限要素法や有限体積法と比較して、容易に計算用データモデルが作成可能な数値解析手法である。
ボクセル法は、基本的に同一サイズの直方体形状の複数のボクセル（直交格子）によって解析領域を定義するボクセルデータを計算用データモデルとして作成し、このボクセルデータを用いた物理量計算を行うことによって数値解析を行う方法である。ボクセル法としては、重み付き残差積分法による支配方程式を用いる重み付き残差積分法型と、例えばセル−オートマトンモデルや格子−ボルツマンモデル等を用いる非積分法型に大別される。そして、このボクセル法によれば、ボクセルデータとして、Vertex等は必要ない。
このようなボクセル法によれば、解析領域をボクセルによって分割することによって解析領域を容易に定義することができ、短時間で計算用データモデルを作成することができる。
一方、粒子法は、解析領域を複数の粒子によって定義する粒子データを計算用データモデルとして作成し、この粒子データを用いた物理量計算を行うことによって数値解析を行う方法である。粒子法は、支配方程式として、非積分法型で粒子間相互作用モデルを利用している。粒子法では、分割領域がないため、Vertex等を必要としない。このような粒子法によれば、解析領域に例えば粒子を均一に配置することによって解析領域を容易に定義することができ、計算用データモデルを短時間で作成することができる。

米国特許出願公開第２００８／００２１６８４号明細書

従来の有限要素法や有限体積法等の数値解析手法のように、ソルバ処理において分割領域の幾何学的形状を用いる場合には、当然ながら計算用データモデルに対して、分割領域の幾何学的形状を示すデータを持たせることが必須となる。
分割領域の幾何学的形状を定義するためには、Vertexに加えて、頂点の連結情報（Connectivity of Vertex、以下ではConnectivityと略記する）が必要となる。したがって、有限要素法や有限体積法では、計算用データモデルにVertexとConnectivityと持たせる必要がある。
なお、具体的には、Connectivityは、全分割領域の頂点に対して順番に付された全体ノード番号と、１つの分割領域内において頂点に対して順番に付された局所ノード番号との対応情報によって定義される。

このようなVertexとConnectivityとを有する計算用データモデルは、周知のように作成に非常に膨大な作業が必要となる。
例えば、有限要素法で用いられる計算用データモデルでは、図１に示すように、隣接する分割領域が必ずVertexを共有するという条件を満足するように計算用データモデルを作成する必要があり、全ての分割領域がこの条件を満足させるためには、非常に膨大な時間を必要とする。
一方、有限体積法で用いられる計算用データモデルは、図２に示すように、隣接する分割領域に共有されないVertexの存在を許容できる分、有限要素法と比較すればメッシュ生成の自由度が増す。しかしながら、有限体積法においても、共有されないVertexが少なくとも隣接する分割領域の辺上に存在すること、また一般的に分割領域の形状を予め設定された解析要素タイプに合わせることといった条件の中で計算用データモデルを作成する必要があり、メッシュ生成の自由度が高いとは言い難い。

また、近年においては、３次元ＣＡＤ（Computer Aided Design）データ等の３次元形状データから抽出される解析領域に対する数値解析が行われている。しかしながら、３次元形状データは、数値解析用に形成されたデータではなく、面の重なり、面の交差、面間の隙間、微小穴等を示すデータが含まれており、VertexとConnectivityとを有する計算用データモデルの作成に適さない条件を多く含んでいる。このため、これらのVertexとConnectivityとを有する計算用データモデルの作成が出来るように、３次元形状データを修正あるいは変更する必要がある。そして、VertexとConnectivityとを有する計算用データモデルが作成可能なように、３次元形状データを修正あるいは変更するためには、経験や試行錯誤を要する非常に膨大な手作業を行う必要がある。これが、有限要素法や有限体積法を実用で利用する際の大きな問題である。

また、有限体積法のようにソルバ処理において分割領域の体積、境界面の面積及び境界面の法線ベクトルの算出を行う場合には、ソルバ処理における計算量がさらに増えることとなり、ソルバ処理における計算負荷がさらに増大する。

ボクセル法では、短時間で計算用データモデルを作成することができるものの、以下の点で問題がある。ボクセル法は、基本的に解析領域が全て同一サイズのボクセル（直交格子）によって定義される。通常、有限要素法や有限体積法では、より高い解析精度を得たい領域の要素サイズ（分割領域のサイズ）を小さく設定することによって当該領域について正確な物理量計算を行い、また他の領域の要素サイズを大きく設定することによって当該領域についての計算負荷を低減させている。ところが、ボクセル法においては、全てのボクセルが基本的に同一サイズであるため、ボクセルを小さく設定した場合には計算負荷が非常に大きくなり、ボクセルを大きく設定した場合には解析精度が悪化することとなる。
また、ボクセル法では、同一サイズのボクセル（直交格子）を配列することによって解析領域を定義する必要があることから、外部領域との境界付近において解析領域を滑らかにすることができずに階段状となる場合がある。つまり、実際に解析したい領域が斜面や曲面等を有している場合であっても、ボクセルデータにおいてその領域は階段状に表されることとなる。このため、ボクセル法における解析領域形状が実際に解析したい領域形状と異なることとなり、解析精度は悪化する。

これに対して、ボクセルデータの階段状の領域を、実際に解析したい領域が有する斜面や曲面に沿って切断（境界補正）するカットセル法と呼ばれる改良方法が提案されている。しかしながら、この改良方法によると、この境界補正によって非常に小さな分割領域が生成されやすく、このような小さな分割領域が生成された場合には、解析精度の悪化を招く。
また、この改良方法では、カットセルの形成のため及びソルバ処理において、Vertexを利用することとなる。

以上のように、境界補正を行わないボクセル法においては、Vertex等が必要ないものの、ボクセルの生成、いわゆるメッシュ生成には限界がある。すなわち、十分な解析精度を得ようとすると、ボクセルの数が増え、ソルバ処理における計算負荷も増え、問題となる。さらに、境界補正を行うボクセル法の改良方法では、結果としてVertexが必要となるため、結局、分割領域の幾何学的形状の影響を受けることなり、外部領域との境界周辺での分割領域形成のための処理に、経験や試行錯誤を要する非常に膨大な手作業を伴い、形状データモデルを短時間に作成することができないことになる。

一方、粒子法では、ある特定の粒子と他の粒子との結合関係を算出する必要がある。このため、当該特定粒子の近傍に存在する粒子を探索する必要がある。そして、この粒子の近傍探索処理は、原則として全ての粒子に対して行われる。しかしながら、粒子法においては、時刻の変化と共に各粒子が移動し、これによって常に粒子同士の結合関係が変化する。このため、解析における時刻が変化するたびに近傍探索処理を行う必要があり、計算負荷の増大を招く。近傍探索の対象となる粒子を厳選する等して近傍探索処理の計算負荷を低減させる試みがなされているものの、例えば、解析精度を向上させるために粒子数を増大させた場合には、粒子数の２乗に比例して計算負荷が増大する。
このような粒子法において実用時間内における数値解析を実現するためには、大型の並列計算機で多くのＣＰＵ（Central Processing Unit）を使用する必要がある。例えば、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とを持たせることなく物理量を算出することが可能となる。したがって、本発明を用いることにより、プリ処理において、少なくとも分割領域の体積と境界面特性量（境界面の面積及び境界面の法線ベクトル）とを有する計算用データモデルを作成すれば良くなり、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を利用することとなり、計算用データモデルを作成するにあたり分割領域の幾何学的形状の影響を受けることになる。

また、従来の数値解析手法である粒子法は、分割領域という概念が存在しないため、図３に示すように、ソルバ処理において物理量を算出するにあたり、VertexとConnectivityとが必要ないものの、分割領域の代わりに計算用データモデルを定義する粒子の移動に起因し、前述のように計算負荷が増大する。また、粒子法では、保存則を満たしながら物理量の計算を行うことは難しい。

続いて、図４を参照して、本発明と従来の有限体積法とについてさらに詳しく比較を行う。
従来の有限体積法では、前述のように、プリ処理において、メッシュ分割で得られる分割領域の幾何学的な形状を定義するためのVertexとConnectivityとを有する計算用データモデルを作成する。また、一般的には、ソルバ処理において分割領域の係合情報（以下linkと呼ぶ）が必要となる。このため、プリ処理においては、VertexとConnectivityとlinkとを有する計算用データモデルが作成される。
そして、従来の有限体積法では、図４に示すように、プリ処理から、VertexとConnectivityとlinkとを有する計算用データモデルと、ソルバ処理において必要となる境界条件や初期条件等とをソルバ処理に受け渡す。ソルバ処理では、受け渡された計算用データモデルに含まれるVertexやConnectivity等を使用して離散化支配方程式を解くことによって物理量の算出を行う。

一方、本発明では、プリ処理において、任意に配置された分割領域の体積、境界面特性量（境界面の面積、境界面の法線ベクトル）及びlinkを有する計算用データモデルを作成する。また、後に詳説するが、本発明では、必要に応じて、分割領域の内部に配置されるコントロールポイントの座標を計算用データモデルに持たせる場合もある。
そして、本発明では、図４に示すように、プリ処理から、分割領域の体積、境界面特性量及びlink（必要に応じてコントロールポイントの座標）を有する計算用データモデルと、境界条件や初期条件等とをソルバ処理に受け渡す。ソルバ処理では、受け渡された計算用データモデルに含まれる分割領域の体積、境界面特性量等を使用して離散化支配方程式を解くことによって物理量の算出を行う。

そして、図４から分かるように、本発明では、ソルバ処理において、VertexとConnectivityとを使用しないで物理量を計算している点が従来の有限体積法と大きく異なり、この点が本発明の大きな特徴である。このような特徴は、ソルバ処理にて、VertexとConnectivityとを必要としない量のみを使用する離散化支配方程式を用いることによって得られるものである。
この結果、図４に示すように、本発明では、ソルバ処理にVertexとConnectivityとを受け渡す必要がなくなり、プリ処理において、VertexとConnectivityとを持たない計算用データモデルを作成すれば良いこととなる。したがって、従来の有限体積法と比較して本発明では、遥かに容易に計算用データモデルを作成することが可能となり、計算用データモデルの作成における作業負担を軽減することが可能となる。

また、数値解析を行う解析領域の形状が時系列的に変化する場合、すなわち解析領域が移動境界を含む場合がある。このような場合、移動境界に合わせて分割領域を移動及び変形させる必要がある。
従来の有限体積法では、移動境界の移動ごとにおけるVertexを予めストアする方法か、分割領域のいびつな変形により計算不可となった場合に領域分割を再実行する方法によって、移動境界を含む場合の物理量の計算を行っている。これに対して、本発明では、Vertexに替えて分割領域の体積や境界面特性値等を予め求めておいてストアしておく方法か、領域分割の再実行によって移動境界を含む場合の物理量の計算を行うことができる。
従来の有限体積法または本発明において、前述のいずれの方法を採用する場合であっても、複数の計算用データモデルを作成する必要がある。しかしながら、従来の有限体積法では、１つでも膨大な作業量が必要となる計算用データモデルを複数作成するとなると、その作業量は、多くの場合において現実的に負担可能な範囲で行えるものではなくなる。
一方、本発明では、計算用データモデルがVertexとConnectivityとを持つ必要がなく、領域分割過程でVertexやConnectivityの整合性を考慮する必要がないため、とても高速に計算用データモデルを作成することができ、移動境界を含む場合の物理量の計算を容易に行うことができる。

ここで、前述のlinkについて補足をする。linkは、物理量のやり取りを行う分割領域同士を関連付ける情報である。そして、このlinkによって関連付けられる分割領域は、必ずしも空間的に隣接されている必要はなく、空間的に離間していても構わない。このようなlinkは、VertexやConnectivityと関連するものではなく、Vertex及びConnectivityと比較すると、極めて短時間で作成することができるものである。

次に、本発明を用いた数値解析手法（以下、本数値解析手法と称する）の原理、すなわち重み付き残差積分法に基づいて導出された離散化支配方程式と、分割領域の体積と境界面特性量とによって物理量を算出可能となる原理について詳細に説明する。なお、以下の説明において、［］にて挟まれた文字は、図面において太字で記されたベクトルを示す。
まず、本数値解析手法における計算用データモデルは、解析領域を分割して得られる各分割領域の体積と、隣り合う分割領域同士の境界面の特性を示す境界面特性量とを用いて定義されている。
図５は、このような本数値解析手法の計算用データモデルの一例を示す概念図である。
この図において、セルＲ_１，Ｒ_２，Ｒ_３・・・は、解析領域を分割して得られる分割領域であり、各々が体積Ｖ_ａ，Ｖ_ｂ，Ｖ_ｃ・・・を有している。また、境界面Ｅは、セルＲ_１とセルＲ_２との間において物理量の交換が行われる面であり、本発明における境界面に相当するものである。また、面積Ｓ_ａｂは、境界面Ｅの面積を示し、本発明における境界面特性量の１つである。また、［ｎ］_ａｂは、境界面Ｅの法線ベクトルを示し、本発明における境界面特性量の１つである。また、コントロールポイントａ，ｂ，ｃ・・・は、各セルＲ_１，Ｒ_２，Ｒ_３の内部に配置されおり、図５においては各セルＲ_１，Ｒ_２，Ｒ_３・・・の重心位置に配置されている。ただし、コントロールポイントａ，ｂ，ｃ・・・は、必ずしも各セルＲ_１，Ｒ_２，Ｒ_３・・・の重心位置に配置される必要はない。また、αは、コントロールポイントａからコントロールポイントｂまでの距離を１とした場合におけるコントロールポイントａから境界面Ｅまでの距離を示し、境界面Ｅがコントロールポイントａとコントロールポイントｂとを結ぶ線分のどの内分点に存在するかを示す比率である。
なお、境界面は、セルＲ_１とセルＲ_２との間のみに限らず、隣り合う全てのセル間に存在する。そして、境界面の法線ベクトル及び境界面の面積も、境界面ごとに与えられるものである。

そして、実際の計算用データモデルは、各コントロールポイントａ，ｂ，ｃ・・・の配置データと、各コントロールポイントａ，ｂ，ｃ・・・が存在するセルＲ_１，Ｒ_２，Ｒ_３・・・の体積Ｖ_ａ，Ｖ_ｂ，Ｖ_ｃ・・・を示す体積データと、各境界面の面積を示す面積データと、各境界面の法線ベクトルを示す法線ベクトルデータとを有するデータ群として構築されている。
すなわち、本数値解析手法の計算用データモデルは、セルＲ_１，Ｒ_２，Ｒ_３・・・の体積Ｖ_ａ，Ｖ_ｂ，Ｖ_ｃ・・・と、隣り合うセルＲ_１，Ｒ_２，Ｒ_３・・・同士の境界面の特性を示す境界面特性量である境界面の面積と、隣り合うセルＲ_１，Ｒ_２，Ｒ_３・・・同士の境界面の特性を示す境界面特性量である境界面の法線ベクトルとを有して定義されている。
なお、各セルＲ_１，Ｒ_２，Ｒ_３・・・は、コントロールポイントａ，ｂ，ｃ・・・を有している。このため、セルＲ_１，Ｒ_２，Ｒ_３・・・の体積Ｖ_ａ，Ｖ_ｂ，Ｖ_ｃ・・・は、コントロールポイントａ，ｂ，ｃ・・・が仮想的に占める空間（コントロールボリューム）の体積として捉えることができる。
また、本数値解析手法の計算用データモデルは、必要に応じて、境界面が挟まれたコントロールポイント同士を結ぶ線分のどの内分点に存在するかの比率αを示す比率データを有する。

以下では、前述の計算用データモデルを用いて解析領域の各セル（分割領域）における流速を求める物理量計算例について説明する。なお、ここでは、各コントロールポイントにおける流速を各セルにおける流速として求める。
まず、本物理量計算において本数値解析手法は、流体解析の場合、下式（１）で示すナビエ・ストークスの式と、下式（２）で示す連続の式とを用いる。

なお、式（１），（２）において、ｔが時間を示し、ｘ_ｉ（ｉ＝１，２，３）がカーテシアン系における座標を示し、ρが流体密度を示し、ｕ_ｉ（ｉ＝１，２，３）が流体の流速成分を示し、Ｐが圧力を示し、μが流体の粘性係数を示し、添字ｉ（ｉ＝１，２，３），ｊ（ｊ＝１，２，３）がカーテシアン座標系における各方向成分を示している。また、添字ｊに関しては総和規約に従うものとする。

そして、式（１），（２）を、重み付き残差積分法に基づいて、コントロールボリュームの体積に対して積分して示すと、式（１）が下式（３）のように示され、式（２）が下式（４）のように示される。

なお、式（３），（４）において、Ｖがコントロールボリュームの体積を示し、∫_ＶｄＶが体積Ｖに関する積分を示し、Ｓがコントロールボリュームの面積を示し、∫_ＳｄＳが面積Ｓに関する積分を示し、［ｎ］がＳの法線ベクトルを示し、ｎ_ｉ（ｉ＝１，２，３）が法線ベクトル［ｎ］の成分を示し、∂／∂ｎが法線方向微分を示している。

ここで、説明を簡単化するために、流体の密度ρと粘性係数μを定数とする。ただし、以下の定数化は、流体の物性値が、時間、空間、温度等によって変化する場合に対して拡張可能である。
そして、図５のコントロールポイントａについて、境界面Ｅの面積Ｓ_ａｂについて離散化し、代数方程式による近似式に変換すると、式（３）が下式（５）、式（４）が下式（６）のように示される。

ここで、添字ａｂが付く、［ｎ］_ａｂ，［ｕ］_ａｂ，ｕ_ｉａｂ，ｎ_ｉａｂ，Ｐ_ａｂ，（∂ｕ_ｉ／∂ｎ）_ａｂは、コントロールポイントａとコントロールポイントｂとの間の境界面Ｅ上における物理量であることを示す。また、ｎ_ｉａｂは［ｎ］_ａｂの成分である。また、ｍは、コントロールポイントａと結合関係（境界面を挟む関係）にある全てのコントロールポイントの数である。

そして、式（５），（６）をＶ_ａ（コントロールポイントａのコントロールボリュームの体積）で割ると、式（５）が下式（７）のように示され、式（６）が下式（８）のように示される。

ここで、下式（９）とする。

すると、式（７）が下式（１０）のように示され、式（８）が下式（１１）のように示される。

式（１０），（１１）において、［ｕ］_ａｂ、ｕ_ｉａｂ、Ｐ_ａｂ、（∂ｕ_ｉ／∂ｎ）_ａｂは、コントロールポイントａとコントロールポイントｂ上の物理量の重み付け平均（移流項については、風上性を考慮した重み付け平均）により近似的に求められ、コントロールポイントａ，ｂ間の距離及び向きと、その間に存在する境界面Ｅとの位置関係（上記比率α）と、境界面Ｅの法線ベクトルの向きに依存して決定される。ただし、［ｕ］_ａｂ、ｕ_ｉａｂ、Ｐ_ａｂ、（∂ｕ_ｉ／∂ｎ）_ａｂは、境界面Ｅの幾何学的な形状には無関係な量（すなわちセル形状を規定するVertexとConnectivityを必要としない量）である。
また、式（９）で定義されるφ_ａｂも（面積／体積）という量であり、コントロールボリュームの幾何学的形状には無関係な量（すなわちセル形状を規定するVertexとConnectivityを必要としない量）である。
つまり、このような式（１０），（１１）は、セル形状を規定するVertexとConnectivityを必要としない量のみを使用して物理量が算出可能な、重み付き残差積分法に基づく演算式である。

このため、物理量計算（ソルバ処理）に先立って前述の計算用データモデルを作成し、物理量計算において当該計算用データモデルと、式（１０），（１１）の離散化支配方程式とを用いることによって、物理量計算においてコントロールボリュームの幾何学的形状（すなわちセルの形状を規定するVertexとConnectivity）を全く使用せずに、流速の計算を行うことが可能である。

このように、物理量計算においてVertexとConnectivityとを全く使用せずに流速の計算が行えることから、計算用データモデルにVertexとConnectivityとを持たせる必要がなくなる。よって、計算用データモデルの作成にあたり、セルの幾何学的形状に縛られる必要がなくなるため、セルの形状を任意に設定することができる。このため、本数値解析手法によれば、前述のように３次元形状データの修正作業に対する規制を大幅に緩和することができる。

なお、実際に式（１０），（１１）を解くにあたり、［ｕ］_ａｂやＰ_ａｂ等の境界面Ｅ上の物理量は、通常、線形補間によって補間される。例えば、コントロールポイントａの物理量をψ_ａ、コントロールポイントｂの物理用をψ_ｂとすると、境界面Ｅ上の物理量ψ_ａｂは、下式（１２）によって求めることができる。

また、物理量ψ_ａｂは、境界面が挟まれたコントロールポイント同士を結ぶ線分のどの内分点に存在するかの比率αを用いることによって、下式（１３）によって求めることもできる。

したがって、計算用データモデルが比率αを示す比率データを有している場合には、式（１３）を用いて境界面Ｅ上の物理量をコントロールポイントａとコントロールポイントｂとからの離間する距離に応じた重み付け平均を用いて算出することができる。

また、連続体モデルの方程式（ナビエ・ストークスの式等）には、式（１）に示すように、１階の偏導関数（偏微分）が含まれる。
ここで、連続体モデルの方程式の微係数を部分積分、ガウスの発散定理、あるいは一般化されたグリーンの定理を利用して、体積分を面積分に変換し、微分の次数を下げる。これによって１次微分は０次微分（スカラー量またはベクトル量）とすることができる。
例えば一般化されたグリーンの定理では、物理量をψとすると、下式（１４）という関係が成り立つ。

なお、式（１４）において、ｎ_ｉ（ｉ＝１，２，３）は、表面Ｓ上の単位法線ベクトル［ｎ］のｉ方向の成分である。
連続体モデルの方程式の１次微分項は、体積分から面積分の変換により、境界面上ではスカラー量またはベクトル量として取り扱われる。そして、これらの値は、前述の線形補間等によって、各コントロールポイント上の物理量から補間することができる。

また、後述するが、連続体モデルの方程式によっては、２階の偏導関数が含まれる場合もある。
式（１４）の被積分関数をさらに１階微分した式は下式（１５）となり、連続体モデルの方程式の２次微分項は、体積分から面積分の変換により境界面Ｅ上では下式（１６）となる。

なお、式（１５）において∂／∂ｎは法線方向微分を示し、式（１６）において∂／∂ｎ_ａｂは、［ｎ］_ａｂ方向微分を示す。
つまり、連続体モデルの方程式の２次微分項は、体積分から面積分の変換により、物理量ψの法線方向微分（Ｓ_ａｂの法線［ｎ］_ａｂ方向への微分）に、［ｎ］の成分ｎ_ｉａｂ、ｎ_ｊａｂを乗じた形となる。

ここで式（１６）中の∂ψ／∂ｎ_ａｂは、下式（１７）と近似される。

なお、コントロールポイントａとコントロールポイントｂとのコントロールポイント間ベクトル［ｒ］_ａｂは、コントロールポイントａの位置ベクトル［ｒ］_ａとコントロールポイントｂの位置ベクトル［ｒ］_ｂから下式（１８）のように定義される。

したがって、境界面Ｅの面積がＳ_ａｂであるため、式（１６）は下式（１９）となり、これを利用して式（１６）を計算することができる。

なお、式（１６）の導出にあたり、次のことがわかる。
すべての線形偏微分方程式は、定数と、１次、２次、その他の偏導関数に係数を乗じた項の線形和で表わされる。式（１５）から式（１８）において、物理量ψをψの１次偏導関数に置き換えると、より高次の偏導関数の体積分を、式（１４）のように低次の偏導関数の面積分により求めることができる。この手順を、低次の偏微分から順次繰り返すと、線形偏微分方程式を構成するすべての項の偏導関数は、コントロールポイントの物理量ψと、式（１２）又は式（１３）で計算される境界面上のψであるψ_ａｂと、式（１８）で定義されるコントロールポイント間ベクトルから求められるコントロールポイント間距離と、式（５）に示される境界面Ｅの面積Ｓ_ab、式（１６）に示される法線ベクトルの成分ｎ_ｉａｂとｎ_ｊａｂから、すべて求めることができる。
次に非線形の偏微分方程式では、例えば、式（２０）に示す非線形項であるψとψの一次偏導関数とを乗算した項や、一次偏導関数の二乗は、反復計算により数値計算することができる。すなわち、それぞれの項のψ、一次偏導関数の方を反復計算における一反復前の計算値として、近似して、反復計算をすればよい。このような方法により、偏微分方程式における非線形項はすべて数値計算が可能になる。
以上から、前述では特に連続体モデルの方程式について説明したが、それ以外のいかなる偏微分方程式に対しても、VertexとConnectivityとを必要としない離散化が可能であることがわかった。ただし、保存則については別の成立条件が必要である。これについては後述する。

さて、本数値解析手法において物理量計算にあたりVertexとConnectivityとを必要としないことは前述した。このため、計算用データモデルの作成（プリ処理）にあたり、コントロールボリュームの体積と、境界面の面積及び法線ベクトルとを、VertexとConnectivityと使用しないで求めれば、式（１０），（１１）の離散化支配方程式を用いて、コントロールボリュームの幾何学的形状（すなわちセルの幾何学的形状）を全く使用せずに、流速の計算を行うことが可能である。
ただし、本数値解析手法においては、必ずしも、コントロールボリュームの体積と、境界面の面積及び法線ベクトルとを、コントロールボリュームの具体的な幾何学的形状を使用しないで求める必要はない。すなわち、ソルバ処理においてVertexとConnectivityを利用しないので、コントロールボリュームの具体的な幾何学的形状、具体的にはVertexとConnectivityを利用するとしても、従来の有限要素法、有限体積法のような分割領域に関わる制約、すなわち分割領域の歪みや捩じれに対する制約がないため、前述のように容易に計算用データモデルの作成が可能となる。

さて、本数値解析手法においては、条件によっては、前述の法線ベクトルを、コントロールボリューム同士を結ぶ距離ベクトルに置き換えることも可能である。以下にその理由について説明する。
図５において示す境界面Ｅの法線ベクトル［ｎ］_ａｂが、コントロールポイントａとコントロールポイントｂとを結ぶ距離ベクトル［ｒ］_ａｂと同じ方向を向く場合には、法線ベクトル［ｎ］は、下式（２１）のように示すことができる。

したがって、境界面Ｅの法線ベクトル［ｎ］_ａｂが距離ベクトル［ｒ］_ａｂと同じ方向を向く場合には、式（１０），（１１）で示す離散化支配方程式に式（２１）を代入することによって、法線ベクトル［ｎ］_ａｂがコントロールポイントａ，ｂ間方向を向いている場合の離散化支配方程式が得られる。すなわち、境界面と該境界面を挟むコントロールポイントを結ぶベクトルが直交する場合には、離散化支配方程式における法線ベクトルを距離ベクトルに置き換えることができる。
このような離散化支配方程式によれば、コントロールポイントの位置座標だけから、境界面Ｅの法線ベクトル［ｎ］_ａｂを決定することができる。
また、境界面Ｅと距離ベクトルとを出来るだけ直交に近づけることによって物理量を計算する際の精度が向上する。したがって、法線ベクトルを、距離ベクトルに置き換えることによって、計算精度の向上を図ることができる。
さらに、法線ベクトルの任意性を、コントロールポイントを結ぶ方向に固定化することができる。

ただ、法線ベクトルの任意性がなくなることで境界面の姿勢に自由度を持たせることができなくなる場合がある。このような場合には、法線ベクトルの任意性が存在する場合と比較して、計算用データモデルを作成する際にコントロールボリュームの体積と境界面の設定に制約が生じることとなる。また、法線ベクトルを、コントロールポイントを結ぶベクトルに置き換える場合には、距離ベクトルが必要となり、該距離ベクトルを算出するためにコントロールポイントの座標が必要となるため、計算用データモデルにコントロールポイントの座標データあるいは距離ベクトルデータを持たせる必要がある。ただし、この場合でも、本数値解析手法においては、物理量計算にてVertexとConnectivityとは必要ない。

次に、本発明と粒子法との違いを明確にするために、物理量計算にあたり、物理量の保存則が満足される条件について説明する。

例えば、図６に示すように、Ｌ字形の流路を４つに分割するセルＲ_ａ，Ｒ_ｂ，Ｒ_ｃ，Ｒ_ｄについて考える。なお、各セルＲ_ａ，Ｒ_ｂ，Ｒ_ｃ，Ｒ_ｄの内部に配置されるコントロールポイントａ，ｂ，ｃ，ｄは、セルＲ_ａ，Ｒ_ｂ，Ｒ_ｃ，Ｒ_ｄの中心に設置されている。また、各境界面における流速ベクトルは、境界面に対して垂直であるものとする。
なお、図６において、Ｖ_ａがセルＲ_ａの体積（コントロールポイントａのコントロールボリュームの体積）、Ｖ_ｂがセルＲ_ｂの体積（コントロールポイントｂのコントロールボリュームの体積）、Ｖ_ｃがセルＲ_ｃの体積（コントロールポイントｃのコントロールボリュームの体積）、Ｖ_ｄがセルＲ_ｄの体積（コントロールポイントｄのコントロールボリュームの体積）、ρ_ａがセルＲ_ａにおける密度、ρ_ｂがセルＲ_ｂにおける密度、ρ_ｃがセルＲ_ｃにおける密度、ρ_ｄがセルＲ_ｄにおける密度、Ｓ_ａがセルＲ_ａと外部領域との境界面の面積、Ｓ_ｃがセルＲ_ｃと外部領域との境界面の面積、Ｓ_ｄがセルＲ_ｄと外部領域との境界面の面積、Ｓ_ａｂがセルＲ_ａとセルＲ_ｂとの境界面の面積、Ｓ_ａｃがセルＲ_ａとセルＲ_ｃとの境界面の面積、Ｓ_ｂｄがセルＲ_ｂとセルＲ_dとの境界面の面積、Ｓ_ｃｄがセルＲ_ｃとセルＲ_ｄとの境界面の面積、ｕ_ａがセルＲ_ａと外部領域との境界面における流速、ｕ_ｃがセルＲ_ｃと外部領域との境界面における流速、ｕ_ｄがセルＲ_ｄと外部領域との境界面における流速、ｕ_ａｂがセルＲ_ａとセルＲ_ｂとの境界面における流速、ｕ_ａｃがセルＲ_ａとセルＲ_ｃとの境界面における流速、ｕ_ｂｄがセルＲ_ｂとセルＲ_dとの境界面における流速、ｕ_ｃｄがセルＲ_ｃとセルＲ_ｄとの境界面における流速、ρ_ａがセルＲ_ａにおける密度、ρ_abがセルＲ_ａとセルＲ_ｂとの境界面における密度、ρ_acがセルＲ_ａとセルＲ_cとの境界面における密度、ρ_bdがセルＲ_bとセルＲ_dとの境界面における密度、を示している。

図６に示すような４つのコントロールポイントａ，ｂ，ｃ，ｄ（４つのセルＲ_ａ，Ｒ_ｂ，Ｒ_ｃ，Ｒ_ｄ）上で質量保存の式を離散化することを考える。なお、質量保存の式を離散化した離散化支配方程式は、後述の式（４３）に示す。
コントロールポイントａにおける離散化支配方程式は、下式（２２）に示される。

コントロールポイントｂにおける離散化支配方程式は、下式（２３）に示される。

コントロールポイントｃにおける離散化支配方程式は、下式（２４）に示される。

コントロールポイントｄにおける離散化支配方程式は、下式（２５）に示される。

そして、式（２２）〜（２５）を全て足し加えると、下式（２６）を得る。

各コントロールポイントａ，ｂ，ｃ，ｄのコントロールボリュームの体積Ｖ_ａ，Ｖ_ｂ，Ｖ_ｃ，Ｖｄが時間について一定であることを利用して時間微分項の中にくり入れると、式（２６）は、下式（２７）となる。

今、各コントロールポイントａ，ｂ，ｃ，ｄのコントロールボリュームが占める全体積をＶ_ａｂｃｄ、平均密度ρ（オーバーライン）_ａｂｃｄとすると、全体積をＶ_ａｂｃｄが下式（２８）で示され、平均密度ρ（オーバーライン）_ａｂｃｄが下式（２９）で示される。

したがって式（２７）は、下式（３０）のように示される。

式（３０）は、コントロールポイントａ，ｂ，ｃ，ｄのコントロールボリュームが占める全領域に流入する質量流束と流出する質量流束の差が単位時間にコントロールポイントａ，ｂ，ｃ，ｄのコントロールボリュームが占める全領域の平均密度の時間変化（質量の時間変化）に等しいことを意味する。つまり、各コントロールポイントａ，ｂ，ｃ，ｄごとに離散化した質量保存の式は、全コントロールポイントのコントロールボリュームが占める領域に対しても成立する。
すなわち、コントロールポイントが示すコントロールボリューム領域についての離散化支配方程式は、全てのコントロールポイントについて足し加えると、計算対象である解析領域の全領域に関する保存則を満足する方程式にならなくてはならない。

続いて、コントロールポイントの全数をＮとし、式（４３）に示される質量保存の式を足し加えると、下式（３１）が得られる。

式（３１）において、各コントロールポイント間の境界面の面積は、コントロールポイントａ側から見てもコントロールポイントｂ側から見ても等しいとすると、各コントロールポイント間の質量流束（ρ［ｎ］・［ｕ］）・Ｓは、コントロールポイントａ側とコントロールポイントｂ側で正負が逆で絶対値が等しくなるので差し引きゼロとなりキャンセルされる。つまり、式（３１）は、計算する全領域に対して、流入する質量と流出する質量との差が全領域での質量の単位時間変化に等しいことを示す。したがって、式（３１）は、解析領域全体における質量保存の式となる。
よって、式（３１）が、計算する全領域についての質量保存則を満足するためには、２つのコントロールポイント間の境界面の面積と一致するという条件及び法線ベクトルが一方のコントロールポイント側から見た場合と他方のコントロールポイント側から見た場合とで絶対値が一致するという条件が必要である。

また、質量保存則を満足するためには、下式（３２）に示される全コントロールポイントのコントロールボリュームの占める体積が解析領域の全体積と一致するという条件が必要である。

これは、連続体の密度ρが、ρ_１＝ρ_２＝・・・＝ρと１つの変数で表されることを考えれば容易に理解される。
このように質量保存則を満足するためには、全コントロールポイントのコントロールボリュームの体積の総和が解析領域の体積と一致するという条件が必要である。

なお、ここでは、質量保存の式に対して説明を行ったが、保存則は、連続体の運動量やエネルギに対しても成立しなければならない。これらの物理量に対しても、後述の式（５０）、（５５）において全コントロールポイントに対して足し加えることによって、保存側が満足されるためには、全コントロールポイントのコントロールボリュームの占める体積が解析領域の全体積と一致するという条件と、２つのコントロールポイント間の境界面の面積が一致する条件及び法線ベクトルが一方のコントロールポイント側から見た場合と他方のコントロールポイント側から見た場合とで絶対値が一致する（正負逆符号）という条件とが必要であることが分かる。

また、保存則を満たすためには、図７に示すように、コントロールポイントａの占めるコントロールボリュームを考えた場合に、コントロールポイントａを通り、任意の向きの単位法線ベクトル［ｎ］_ｐを持つ無限に広い投影面Ｐを考えたときに下式（３３）が成り立つという条件が必要である。

なお、図７及び式（３３）において、Ｓ_ｉが境界面Ｅ_ｉの面積、［ｎ］_ｉが境界面Ｅ_ｉの単位法線ベクトル、ｍがコントロールボリュームの面の総数を示す。
式（３３）は、コントロールボリュームを構成する多面体が、閉包空間を構成することを示す。この式（３３）は、コントロールボリュームを構成する多面体の一部が凹んでいる場合であっても成立する。
なお、図８に示すように、２次元における三角形についても式（３３）が成り立つ。また、多面体の１つの面を微小面ｄＳとし、ｍを∞とする極限を取ると、下式（３４）となり、図９に示すような閉包曲面体についても成り立つことが分かる。

式（３３）が成り立つという条件は、ガウスの発散定理や、式（１４）に示す一般化されたグリーンの定理が成り立つために必要な条件である。
そして、一般化されたグリーンの定理は、連続体の離散化のための基本となる定理である。したがって、グリーンの定理にしたがって体積分を面積分に変化し離散化させる場合において、保存則を満足させるためには式（３３）が成り立つという条件は必須となる。

このように、前述の計算用データモデル及び物理量計算を用いて数値解析を行う際に、物理量の保存則が満足されるためには、以下の３つの条件が必要となる。
（ａ）全コントロールポイントのコントロールボリュームの体積（全分割領域の体積）の総和が解析領域の体積と一致する。
（ｂ）２つのコントロールポイント間の境界面の面積が一致する及び法線ベクトルが一方のコントロールポイント側（境界面を挟む一方の分割領域）から見た場合と他方のコントロールポイント側（境界面を挟む他方の分割領域）から見た場合とで絶対値が一致する。
（ｃ）コントロールポイントを通り（分割領域を通り）、任意の向きの単位法線ベクトル［ｎ］_ｐを持つ無限に広い投影面Ｐを考えたときに式（３３）が成り立つ。

つまり、保存則を満足させる場合には、これらの条件を満足するように計算用データモデルを作成する必要がある。ただし、前述のように本数値解析手法においては、計算用データモデルの作成にあたり、セル形状を任意に変形することができることから、容易に上記３つの条件を満足するように計算用データモデルを作成することができる。

次に、従来の粒子法であるＭＰＳ（Moving Particle Semi-Implicit）法が保存則を満足できない理由、計算負荷が増大する理由、及び本数値解析手法の粒子法に対する優位性について詳説する。

ＭＰＳ法は、適切に設定された半径ｒ_ｅの球内に存在する粒子を検出し、それらと結合関係を結んで計算する手法であり、例えば図１０に示すように、粒子ｉ周りに複数の粒子ｊが存在する場合に、粒子ｉにおけるラプラシアン（∇^２ψ）_ｉは、下式（３５）のように近似される。

なお、図１０において、ψ_ｉが粒子ｉにおける物理量を示し、ψ_ｊが粒子ｊにおける物理量を示し、［ｒ］_ｉｊが粒子ｉから粒子ｊへの距離ベクトルを示している。
また、式（３５）におけるｄは次元数を示す定数であり、３次元である場合には３となる。また式（３５）におけるω（ｒ）は重み関数であり、下式（３６）のように示される。また、式（３５）におけるｍは、結合関係にある粒子の数を示す。

一方、図１１に示すように、粒子ｉ，ｊをコントロールポイントとして捉え、粒子ｉのコントロールボリュームがＶ_ｉ、粒子ｉと粒子ｊとの間の境界面の面積がＳ_ｉｊ、粒子ｉと粒子ｊとの間の境界面の法線ベクトルが［ｎ］_ｉｊ、粒子ｉから粒子ｊへの距離ベクトルが［ｒ］_ｉｊである場合には、粒子ｉにおけるラプラシアン（∇^２ψ）_ｉは、下式（３７）のように近似される。

そして、式（３６）と式（３７）とを比較すると、仮に［ｒ］_ｉｊと［ｎ］_ｉｊが同じ向きであるとした場合には、下式（３８），（３９）が得られる。

式（３９）は、左辺と右辺の次元が（１／距離）で一致している。したがって、右辺のＭＰＳ法定式化は、２つの粒子ｉ，ｊ間の距離のみから、下式（４０）に示す（面積／体積）という量（すなわち式（９）で定義された比率）を算出する手法と解釈出来る。

しかし、ｍ個の粒子の結合関係のみから、境界面の面積Ｓ_ｉｊと体積Ｖ_ｉを求めるには、関係式が不足しており、具体的な数値は決定できず、あくまでも、式（４０）の比率のみが決定される。
このため、仮に、ＭＰＳ法の離散化支配方程式から、境界面の面積Ｓ_ｉｊと体積Ｖ_ｉを求めたとしても、前述した保存則が満足される条件（ａ）〜（ｃ）を満足する保証は全くない。このことは、ＭＰＳ法が、保存則の満足性という点で大きな問題を有していることを示している。
数値解析を工学上の問題、特に機械設計問題やプラント設計問題に適用した場合において、定量値（圧力、温度、熱量等）の評価は極めて重要であるが、数値解析が保存則を満足しない場合には、定量性が保証されなくなる。
つまり、ＭＰＳ法では、質量保存、運動量保存、エネルギ保存という保存則が満足される保証がなく、定量性が保証されない。
これに対して、本数値解析手法によれば、保存則を満足させることができ、定量性を保証することができる。

なお、ＭＰＳ法は、前述のように、時刻の変化に伴って粒子が移動するため、例えば、その都度、前述の半径ｒ_ｅの球内に存在する粒子を検出する近傍探索処理を行う必要があり、物理量計算における計算負荷が大きくなる。
これに対して、本数値解析手法においてコントロールボリューム及びコントロールポイントは、時刻が変化した場合であっても移動しない。このため、予めコントロールボリュームやコントロールポイントの配置関係が分かっている場合は近傍探索処理を行うことなく物理量計算を行うことができる。このため、物理量計算における計算負荷をＭＰＳ法と比較して小さくすることができる。なお、コントロールボリューム及びコントロールポイントの配置関係が予め分かっていない場合であっても、最初に一度のみ、コントロールボリュームやコントロールポイントの配置関係を定める処理を行えば良い。

なお、前述の説明においては、ナビエ・ストークスの式及び連続の式から重み付き残差積分法に基づいて導出した離散化支配方程式を用いる物理量の計算例について説明したが、本数値解析手法において用いられる離散化支配方程式はこれに限られるものではない。
つまり、種々の方程式（質量保存の方程式、運動量保存の方程式、エネルギ保存の方程式、移流拡散方程式、及び波動方程式等）から重み付き残差積分法に基づいて導出されると共に、VertexとConnectivityとを必要としない量のみを使用して物理量を算出可能な離散化支配方程式であれば本数値解析手法に用いることができる。
そして、このような離散化支配方程式の特性によって、従来の有限要素法や有限体積法のようにいわゆるメッシュを必要としない、メッシュレスでの計算が可能となる。また、たとえ、プリ処理において、セルの幾何学的形状を規定するVertexとConnectivityとを利用するとしても、従来の有限要素法、有限体積法、ボクセル法のようなメッシュに対する制約がないため、計算用データモデルの作成に伴う作業負荷を低減できる。
以下に、質量保存の方程式、運動量保存の方程式、エネルギ保存の方程式、移流拡散方程式、及び波動方程式から、重み付き残差積分法に基づいて、VertexとConnectivityとを必要としない量のみを使用する離散化支配方程式が導出可能であること、すなわち本数値解析手法において他の支配方程式を用いることができることについて説明する。

（１）質量保存の方程式
オイラー座標系における質量保存の方程式は、下式（４１）のように微分形式で示される。

なお、式（４１）においてｔが時間、ｘ_ｉ（ｉ＝１，２，３）がカーテシアン系における座標、ρが密度、ｕ_ｉ（ｉ＝１，２，３）が変形速度の成分、添字ｉ（ｉ＝１，２，３）がカーテシアン座標系における各方向成分を示している。また、添字ｉに関しては総和規約に従うものとする。
そして、式（４１）を、重み付き残差積分法に基づいて、コントロールポイントのコントロールボリュームの体積Ｖに対して積分すると下式（４２）と示される。

さらに、図５に示すコントロールポイントａについて離散化し、代数方程式に変換すると、下式（４３）のように示される。

ここで、添字ａｂが付く、ρ_ａｂ、［ｎ］_ａｂは、コントロールポイントａとコントロールポイントｂとの間の境界面Ｅ上における物理量であることを示す。また、ｍは、コントロールポイントａと結合関係（境界面を挟む関係）にある全てのコントロールポイントの数である。

そして、式（４３）をコントロールポイントａのコントロールボリュームの体積であるＶ_ａで割ると、下式（４４）が得られ、さらに下式（４５）とすると、質量保存の方程式が離散化された下式（４６）を得る。

このような式（４６）は、重み付き残差積分法に基づいて導出されると共に、VertexとConnectivityとを必要としない量のみを使用する方程式であるため、本数値解析手法の離散化支配方程式として用いることができる。

なお、本数値解析手法で用いられる離散化支配方程式は、従来の幾何学的形状を示す量を使用する方程式を重み付き残差積分法に基づいて導出する過程で敢えて途中にて留めることによって得られることは前述した。
つまり式（４６）は、質量保存の方程式を、重み付き残差積分法に基づいて、Vertex等を使用する方程式を導出する過程で得られるものである。
ここで、図１２は、２次元三角形状のセルを示す模式図である。図１２における三角形aの面積、辺の長さ、法線ベクトルを下表に示す。なお、下表において記号×は、外積を示している。

図１２に示すように、セルが２次元三角形状である場合において、質量保存の方程式を、重み付き残差積分法に基づいて、Vertex等を使用する離散化支配方程式とすると、下式（４７）となる。

ただし、式（４７）において、［ｒ］_ｉが頂点（Vertex）_ｉの位置ベクトル（座標）であり、記号×がベクトルの外積を示している。また、ρ_ａｂ及びρを一定とし、［ｒ］_ｉｊを［ｒ］_ｊ−［ｒ］_ｉとし、［ｒ］_４を［ｒ］_１とする。

（２）運動量保存の方程式
オイラー座標系における運動量保存の方程式は、下式（４８）のように微分形式で示される。

なお、式（４８）においてσ_ｉｊ（ｉ，ｊ＝１，２，３）が連続体内部応力、ｆ_ｉ（ｉ＝１，２，３）が連続体に作用する外力（例えば重力）を示し、他の量は式（４１）と同様である。また、添字ｊに関しては総和規約に従うものとする。
この式（４７）は、構造体や材料、流体等の応力場の基礎方程式である。

そして、式（４７）を、重み付き残差積分法に基づいて、コントロールポイントのコントロールボリュームの体積Ｖに対して積分すると下式（４９）と示される。

さらに、図５に示すコントロールポイントａについて離散化し、代数方程式に変換すると、下式（５０）のように示される。

式（５０）をコントロールポイントａのコントロールボリュームの体積であるＶ_ａで割り、さらに式（４５）を導入すると、運動量保存の方程式が離散化された下式（５１）が得られる。

なお、運動量保存の方程式において応力テンソルの対称性を考慮すると、角運動量保存の方程式も、運動量保存の方程式と同様に離散化できる。

（３）移流拡散方程式
ある物質Ｃの連続体内への移流拡散現象は、下式（５２）の移流拡散方程式で示される。

なお、式（５２）において、Ｃが物質Ｃの濃度、μ_Ｃが物質Ｃの拡散係数、ｑ_Ｃが物質Ｃのソース（シンク）項、ρが連続体の密度、ｕ_ｉが連続体の変形速度を示す。
そして、式（５２）を、重み付き残差積分法に基づいて積分し、さらに離散化して離散化支配方程式に変換すると、下式（５３）を得る。

なお、Ｃａ等の添字ａが付く量がコントロールポイントａの物理量であり、Ｃ_ａｂ等の添字ａｂが付く量がコントロールポイントａとコントロールポイントｂとの間の境界面における物理量である。

（４）エネルギ保存の方程式
エネルギ保存則は、熱エネルギ保存と運動エネルギ保存に分けられるが、運動エネルギの保存は前述の運度量保存に含まれるため、ここでは熱エネルギ保存の方程式の一般形を下式（５４）に示す。

なお、式（５４）において、Ｕが連続体の内部エネルギ、ｑ_ｉが熱流束ベクトル、ｒが熱源，熱エネルギのソース（シンク）項、σ_ｉｊが応力テンソル、Ｄ_ｉｊが変形速度テンソルを示す。なお、σ_ｉｊ・Ｄ_ｉｊのテンソル２重積の項は、応力仕事率と呼ばれる。また、添字ｉ，ｊについては、総和規約を適用する。
そして、式（５４）を、重み付き残差積分法に基づいて積分し、さらに離散化して離散化支配方程式に変換すると、下式（５５）を得る。

そして、熱流束ベクトル［ｑ］に、熱流束に加えて電気的または化学的エネルギ等の全ての非力学的エネルギの流束を加えると、エネルギ保存の方程式は、非常に広い範囲のエネルギの保存を示す方程式となる。

（５）波動方程式
前述の質量保存の方程式、運動量保存の方程式、エネルギ保存の方程式、移流拡散方程式等の保存形で表される物理法則の方程式は、「放物型」と「楕円型」と呼ばれる偏微分方程式の性質を合わせ持った方程式である。これに対して、波の伝わりや振動の伝わりを表す波動方程式は「双曲型」と呼ばれ、一般形として下式（５６）と示される。

なお、式（５６）において、ｕが振幅，変位、αが波の伝播速度を示している。
そして、式（５６）を、重み付き残差積分法に基づいて積分し、さらに離散化して離散化支配方程式に変換すると、下式（５７）が得られる。

式（５６）から、空間方向にはコントロールボリュームの幾何学的形状を必要としない離散化手法がそのまま適用できることが分かる。ただし、時間方向には２階の導関数であるので、高精度な時間積分法が用いられる。

以上のような式（４６），（５１），（５３），（５５），（５７）は、重み付き残差積分法に基づいて導出されると共に、VertexとConnectivityとを必要としない量のみを使用して物理量を算出可能な方程式であるため、本数値解析手法の離散化支配方程式として用いることができる。
そして、本数値解析手法は、上記離散化支配方程式を用いることによって、定常及び非定常の、流体力学、熱伝導、移流拡散、構造力学、波動、及びこれらの物理現象が連成した現象における数値解析に適用することができる。

次に、本数値解析手法を用いる場合には、異なる座標系で設計された部品を容易に組み上げた状態で数値解析を行うことができる。これは、本数値解析手法においては、有限要素法で用いられる要素と異なり、コントロールボリュームの具体的な幾何学的形状を必要とせず、また２つの隣接するコントロールポイント間の「距離」は、絶対座標系における距離である必要はなく、「計算上の距離」であれば良いためである。
したがって、本数値解析手法によれば、図１３に示すように、異なる座標系で設計された部品Ａ，Ｂ，Ｃを組み上げた状態で数値解析したい場合において、各部品の座標系を一致させることなく数値解析を行うことができる。

なお、一般的な有限要素法と呼ばれる数値解析手法は、図１４や図１５に示すような要素の交差を全く許さない手法である。このため、有限要素法において用いられるアプリケーションソフトにて要素の交差が検出されるとエラー出力がなされる。そして、このようなことから、有限要素において計算用データモデルの作成に非常に大きな作業負担がかかることは前述の通りである。
一方、本数値解析手法によれば、図１４や図１５に示すように、コントロールポイント間の結合に交差を許すことが可能である。これは、各コントロールポイントが占めるコントロールボリュームの体積（分割領域の体積）、境界面の面積、及び境界面の法線ベクトルという情報量は、具体的な幾何学的形状を有している必要がない。したがって、コントロールポイント間の結合が交差する場合であっても物理量計算を実行することが可能である。
したがって、計算用データモデルを作成する際の制約が減り、計算用データモデルの作成の自由度が大幅に拡大する。ただし、流体などにおける物理現象を考えた場合には、１つのコントロールポイントに近接するコントロールポイントからの情報によって、そのコントロールポイントの物理用がアップデートされるという性質を有しているため、遠方のコントロールポイントと結合させると計算精度が悪化する虞がある。このため、本数値解析手法においても適切な範囲においてコントロールポイント間結合を形成することが好ましい。

従来の有限要素法による計算用データモデルの一例を示す概念図である。 従来の有限体積法による計算用データモデルの一例を示す概念図である。 本発明と従来の数値解析手法とを比較するための表である。 本発明と従来の有限体積法とを詳細に比較するための図である。 本数値解析手法の計算用データモデルの一例を示す概念図である。 本数値解析手法において物理量の保存則が満足される条件について説明するための複数の分割領域を示す模式図である。 コントロールポイントを通り、任意の向きの単位法線ベクトルを持つ無限に広い投影面を示す模式図である。 ２次元における三角形のコントロールボリュームを考えた場合において物理量の保存則が満足される条件について説明する模式図である。 球のコントロールボリュームを考えた場合において物理量の保存則が満足される条件について説明する模式図である。 粒子法による近傍探索について示す模式図である。 本数値解析手法におけるコントロールボリュームを示す模式図である。 ２次元三角形状のセルを示す模式図である。 異なる座標系で設計された部品を組み上げた状態を示す模式図である。 要素が交差する様子を示す模式図である。 コントロールポイント間の結合が交差する様子を示す模式図である。 本発明の一実施形態における数値解析装置のハードウェア構成を概略的に示すブロック図である。 本発明の一実施形態における数値解析方法を示すフローチャートである。 本発明の一実施形態における数値解析方法にて行うプリ処理を示すフローチャートである。 本発明の一実施形態における数値解析方法にて行うソルバ処理を示すフローチャートである。 ３次元形状データに含まれる車両の室内空間を微小な閉曲面によってラッピングした様子を示す模式図である。 同一直交格子形状の分割領域によって３次元形状データに含まれる車両の室内空間を含む解析領域を形成した様子を示す模式図である。 室内空間から食み出した分割領域を削除した様子を示す模式図である。 新たな分割領域を解析領域の境界面と室内空間の境界面との間に形成される隙間に充填した様子を示す模式図である。 解析領域が移動境界を含む場合における本発明の一実施形態における数値解析方法を示すフローチャートである。

以下、図面を参照して、本発明に係る物理量計算方法、数値解析方法、物理量計算プログラム、数値解析プログラム、物理量計算装置及び数値解析装置について説明する。
なお、以下の説明においては、本発明に係る物理量計算方法を含む数値解析方法と、本発明に係る物理量計算プログラムを含む数値解析プログラムと、本発明に係る物理量計算装置を含む数値解析装置との実施形態について説明する。
また、以下の実施形態においては、車両の室内空間における空気の流速を数値解析によって求める場合について説明する。

図１６は、本実施形態の数値解析装置Ａのハードウェア構成を概略的に示すブロック図である。
この図に示すように、本実施形態の数値解析装置Ａは、パーソナルコンピュータやワークステーション等のコンピュータによって構成されるものであり、ＣＰＵ１、記憶装置２、ＤＶＤ（Digital Versatile Disc）ドライブ３、入力装置４、出力装置５、及び通信装置６を備えている。

ＣＰＵ１は、記憶装置２、ＤＶＤドライブ３、入力装置４、出力装置５、及び通信装置６と電気的に接続されており、これらの各種装置から入力される信号を処理すると共に、処理結果を出力するものである。

記憶装置２は、メモリ等の内部記憶装置及びハードディスクドライブ等の外部記憶装置によって構成されており、ＣＰＵ１から入力される情報を記憶すると共にＣＰＵ１から入力される指令に基づいて記憶した情報を出力するものである。
そして、本実施形態において記憶装置２は、プログラム記憶部２ａとデータ記憶部２ｂとを備えている。

プログラム記憶部２ａは、数値解析プログラムＰを記憶している。この数値解析プログラムＰは、所定のＯＳにおいて実行されるアプリケーションプログラムであり、コンピュータから構成される本実施形態の数値解析装置Ａを、数値解析を行うように機能させるものである。そして、数値解析プログラムＰは、本実施形態の数値解析装置Ａを、例えば計算用データモデル作成手段と物理量計算手段として機能させるものである。
そして、図１６に示すように、数値解析プログラムＰは、プリ処理プログラムＰ１と、ソルバ処理プログラムＰ２と、ポスト処理プログラムＰ３とを有している。

プリ処理プログラムＰ１は、ソルバ処理を実行するための前処理（プリ処理）を本実施形態の数値解析装置Ａに実行させるものであり、本実施形態の数値解析装置Ａを計算用データモデル作成手段として機能させることによって計算用データモデルを作成させる。また、プリ処理プログラムＰ１は、本実施形態の数値解析装置Ａに、ソルバ処理を実行するにあたり必要となる条件の設定を実行させ、さらには上記計算用データモデルや設定された条件を纏めたソルバ入力データファイルＦの作成を実行させる。

そして、プリ処理プログラムＰ１は、本実施形態の数値解析装置Ａを計算用データモデル作成手段として機能させる場合に、まず本実施形態の数値解析装置Ａに対して、車両の室内空間を含む３次元形状データを取得させ、この取得させた３次元形状データに含まれる車両の室内空間を示す解析領域の作成を実行させる。
なお、後に詳説するが、本実施形態においては、ソルバ処理において、前述の本発明を用いた数値解析手法にて説明した離散化支配方程式（VertexとConnectivityとを必要としない量のみを使用すると共に重み付き残差積分法に基づいて導出された離散化支配方程式）を用いる。このため、計算用データモデルの作成にあたり、保存則を満たす条件の下、分割領域の形状及び解析領域の形状を任意に変更することができる。よって、３次元形状データに含まれる車両の室内空間の修正あるいは変更作業は、ラッピング処理等の簡易的なもので充分となる。そこで、本実施形態においてプリ処理プログラムＰ１は、本実施形態の数値解析装置Ａに対して、取得させた３次元形状データに含まれる車両の室内空間を微小な閉曲面によってラッピングすることによって車両の室内空間に存在する穴や隙間を修繕するラッピング処理を実行させる。
その後、プリ処理プログラムＰ１は、本実施形態の数値解析装置Ａに対して、同一サイズの分割領域を形成し、ラッピング処理された室内空間の全領域を含む解析領域の作成を実行させる。続いて、プリ処理プログラムＰ１は、本実施形態の数値解析装置Ａに対して、同一サイズに作成した分割領域のうち室内空間から食み出した領域をカットすることによって、室内空間を示す解析領域の作成を実行させる。ここでも、ソルバ処理において前述の離散化支配方程式を用いることから、解析領域のうち室内空間から食み出した領域を容易にカットすることができる。これにより、ボクセル法のように、外部空間との境界が階段状になることがなく、また、ボクセル法のカットセル法のような外部空間の境界付近の解析領域の形成に対して、経験や試行錯誤を要する非常に膨大な手作業を伴う特別な修正または処理を必要としない。すなわち、本実施形態では、ボクセル法で問題となる外部空間との境界の処理に関わる問題がない。
なお、本実施形態においては、後述のように室内空間とカットした領域との隙間に新たな任意形状の分割領域を充填することによって、直交格子形状のみによらない分割領域で解析領域が構成されるようにし、さらには解析領域に分割領域を重なることなく充填させている。

また、プリ処理プログラムＰ１は、本実施形態の数値解析装置Ａを計算用データモデル作成手段として機能させる場合に、本実施形態の数値解析装置Ａに対して、作成させた室内空間を示す解析領域に含まれる分割領域の各々の内部に対して１つのコントロールポイントを仮想的に配置する処理を実行させ、コントロールポイントの配置情報、及び各コントロールポイントが占めるコントロールボリュームの体積データを記憶させる。
また、プリ処理プログラムＰ１は、本実施形態の数値解析装置Ａを計算用データモデル作成手段として機能させる場合に、本実施形態の数値解析装置Ａに対して、上記分割領域同士の境界面である境界面の面積及び法線ベクトルの算出を実行させ、これらの境界面の面積及び法線ベクトルを記憶させる。
また、プリ処理プログラムＰ１は、本実施形態の数値解析装置Ａを計算用データモデル作成手段として機能させる場合に、コントロールボリューム（コントロールポイント）の結合情報（Link）を作成させ、このLinkを記憶させる。
そして、プリ処理プログラムＰ１は、本実施形態の数値解析装置Ａに対して、上記各コントロールポイントが占めるコントロールボリュームの体積と、境界面の面積及び法線ベクトルと、コントロールポイント（すなわち分割領域）の配置情報（座標）と、Linkとを纏めさせて計算用データモデルを作成させる。

また、プリ処理プログラムＰ１は、本実施形態の数値解析装置Ａに、前述のソルバ処理を実行するにあたり必要となる条件の設定を行わせる場合には、物性値の設定、境界条件の設定、初期条件の設定、計算条件の設定を行わせる。
ここで、物性値とは、室内空間における空気の密度、粘性係数等である。
境界条件とは、コントロールポイント間の物理量の交換の法則を規定するものであり、本実施形態においては前述した式（１０）で示されるナビエ・ストークスの式に基づく離散化支配方程式、及び式（１１）で示される連続の式に基づく離散化支配方程式である。
また、境界条件には、室内空間と外部空間との境界面に臨む分割領域を示す情報が含まれる。
初期条件とは、ソルバ処理を実行する際の最初の物理量を示すものであり、各分割領域の流速の初期値である。
計算条件とは、ソルバ処理における計算の条件であり、例えば反復回数や収束基準である。

また、プリ処理プログラムＰ１は、本実施形態の数値解析装置Ａに、ＧＵＩ（Graphical User Interface）を形成させる。より詳細には、プリ処理プログラムＰ１は、出力装置５が備えるディスプレイ５ａに対してグラフィックを表示させると共に、入力装置４が備えるキーボード４ａやマウス４ｂによって操作が可能な状態とさせる。

ソルバ処理プログラムＰ２（物理量計算プログラム）は、本実施形態の数値解析装置Ａにソルバ処理を実行させるものであり、本実施形態の数値解析装置Ａを物理量計算装置として機能させる。
そして、ソルバ処理プログラムＰ２は、本実施形態の数値解析装置Ａを物理量計算手段として機能させる場合に、計算用データモデルが有するコントロールボリュームの体積と境界面の面積及び法線ベクトルとを含むソルバ入力データファイルＦを用いて、解析領域における物理量を物理量計算させる。

そして、ソルバ処理プログラムＰ２は、本実施形態の数値解析装置Ａを物理量計算手段として機能させる場合に、本実施形態の数値解析装置Ａに対して、ソルバ入力データファイルＦに含まれるナビエ・ストークスの式及び連続の式の離散化係数行列の作成を実行させると共に、マトリックス形成用のデータテーブルの作成を実行させる。
また、ソルバ処理プログラムＰ２は、本実施形態の数値解析装置Ａを物理量計算手段として機能させる場合に、本実施形態の数値解析装置Ａに対して、前述した式（１０）で示されるナビエ・ストークスの式に基づく離散化支配方程式、及び前述した式（１１）で示される連続の式に基づく離散化支配方程式から、下式（５８）で示すマトリックス計算用の大規模粗行列方程式の組み上げを実行させる。

なお、式（５８）において［Ａ］が大規模粗行列を示し、［Ｂ］が境界条件ベクトルを示し、［Ｘ］が流速の解を示す。

また、ソルバ処理プログラムＰ２は、上記離散化支配方程式に非圧縮性等の付帯条件が存在する場合には、本実施形態の数値解析装置Ａに対して、この付帯条件の行列方程式への組み上げを実行させる。
そして、ソルバ処理プログラムＰ２は、本実施形態の数値解析装置Ａに対して、ＣＧ法（共役勾配法）等による行列方程式の解の計算、当該解の下式（５９）を用いた解のアップデート、収束条件の判定を実行させ、最終的な計算結果を取得させる。

ポスト処理プログラムＰ３は、本実施形態の数値解析装置Ａにポスト処理を実行させるものであり、本実施形態の数値解析装置Ａに対して、ソルバ処理において取得された計算結果に基づく処理を実行させる。
より詳細には、ポスト処理プログラムＰ３は、本実施形態の数値解析装置Ａに対して、計算結果の可視化処理、抽出処理を実行させる。
ここで、可視化処理とは、例えば、断面コンタ表示、ベクトル表示、等値面表示、アニメーション表示を出力装置５に出力させる処理である。また、抽出処理とは、作業者が指定する領域の定量値を抽出して数値やグラフとして出力装置５に出力させる、あるいは作業者が指定する領域の定量値を抽出してファイル化したものの出力を実行させる処理である。
また、ポスト処理プログラムＰ３は、本実施形態の数値解析装置Ａに対して、自動レポート作成、計算残差の表示及び分析を実行させる。

データ記憶部２ｂは、図１６に示すように、計算用データモデルＭ、境界条件を示す境界条件データＤ１、計算条件を示す計算条件データＤ２、物性値を示す物性値データＤ３、及び初期条件を示す初期条件データＤ４を有するソルバ入力データファイルＦと、３次元形状データＤ５と、計算結果データＤ６等を記憶するものである。また、データ記憶部２ｂは、ＣＰＵ１の処理過程において生成される中間データを一時的に記憶する。

ＤＶＤドライブ３は、ＤＶＤメディアＸを取り込み可能に構成されており、ＣＰＵ１から入力される指令に基づいて、ＤＶＤメディアＸに記憶されるデータを出力するものである。そして、本実施形態においては、ＤＶＤメディアＸに数値解析プログラムＰが記憶されており、ＤＶＤドライブ３は、ＣＰＵ１から入力される指令に基づいて、ＤＶＤメディアＸに記憶される数値解析プログラムＰを出力する。

入力装置４は、本実施形態の数値解析装置Ａと作業者とのマンマシンインターフェイスであり、ポインティングデバイスであるキーボード４ａやマウス４ｂを備えている。
出力装置５は、ＣＰＵ１から入力される信号を可視化して出力するものであり、ディスプレイ５ａ及びプリンタ５ｂを備えている。
通信装置６は、本実施形態の数値解析装置ＡとＣＡＤ装置Ｃ等の外部装置との間においてデータの受け渡しを行うものであり、社内ＬＡＮ（Local Area Network）等のネットワークＢに対して電気的に接続されている。

次に、このように構成された本実施形態の数値解析装置Ａを用いた数値解析方法（本実施形態の数値解析方法）について、図１７〜図１９のフローチャートを参照して説明する。

図１７のフローチャートに示すように、本実施形態の数値解析方法は、プリ処理（ステップＳ１）と、ソルバ処理（ステップＳ２）と、ポスト処理（ステップＳ３）とから構成されている。
なお、本実施形態の数値解析方法を行うより前に、ＣＰＵ１は、ＤＶＤドライブ３に取り込まれたＤＶＤメディアＸに記憶された数値解析プログラムＰをＤＶＤメディアＸから取り出し、記憶装置２のプログラム記憶部２ａに記憶させる。
そして、ＣＰＵ１は、入力装置４から数値解析の開始を指示する信号が入力されると、記憶装置２に記憶された数値解析プログラムＰに基づいて数値解析を実行する。より詳細には、ＣＰＵ１は、プログラム記憶部２ａに記憶されたプリ処理プログラムＰ１に基づいてプリ処理（ステップＳ１）を実行し、プログラム記憶部２ａに記憶されたソルバ処理プログラムＰ２に基づいてソルバ処理（ステップＳ２）を実行し、プログラム記憶部２ａに記憶されたポスト処理プログラムＰ３に基づいてポスト処理（ステップＳ３）を実行する。なお、このようにＣＰＵ１がプリ処理プログラムＰ１に基づくプリ処理（ステップＳ１）を実行することによって、本実施形態の数値解析装置Ａが計算用データモデル作成手段として機能される。また、ＣＰＵ１がソルバ処理プログラムＰ２に基づくソルバ処理（ステップＳ２）を実行することによって、本実施形態の数値解析装置Ａが物理量計算手段として機能される。

図１８は、プリ処理（ステップＳ１）を示すフローチャートである。この図に示すように、プリ処理（ステップＳ１）が開始されると、ＣＰＵ１は、通信装置６に、ネットワークＢを介してＣＡＤ装置Ｃから車両の室内空間を含む３次元形状データＤ５を取得させる（ステップＳ１ａ）。ＣＰＵ１は、取得した３次元形状データＤ５を記憶装置２のデータ記憶部２ｂに記憶させる。

続いて、ＣＰＵ１は、ステップＳ１ａで取得した３次元形状データＤ５の分析を行い（ステップＳ１ｂ）、データ記憶部２ｂに記憶された３次元形状データＤ５に含まれる、曲面の重なり、交差した曲面、曲面間の隙間、微小穴等を検出する。

続いて、ＣＰＵ１は、ステップＳ１ａで取得した３次元形状データＤ５の修正あるいは変更処理を実行する（ステップＳ１ｃ）。
より詳細には、ＣＰＵ１は、図２０に示すように、３次元形状データＤ５に含まれる室内空間Ｋを微小な閉曲面によってラッピング処理することによって、曲面の重なり、交差した曲面、曲面間の隙間、微小穴等の存在が排除された室内空間の３次元形状データＤ５とする。
なお、ＣＰＵ１は、当該修正あるいは変更処理（ステップＳ１ｃ）において、ＧＵＩを形成し、ＧＵＩから指令（例えばラッピング領域を示す指令）が入力された場合には、当該指令を反映させた修正あるいは変更処理を実行する。

続いて、ＣＰＵ１は、計算用データモデルの作成を実行する（ステップＳ１ｄ）。
より詳細には、まずＣＰＵ１は、図２１に示すように、ステップＳ１ｃにおいて修正あるいは変更された３次元形状データＤ５から、室内空間の全領域を含むと共に同一形状（直交格子）の分割領域にて均等に分割された解析領域Ｋ１の作成を実行する。なお、ここでは、時間の短縮のために、解析領域Ｋ１が同一形状の分割領域で分割するようにしているが、解析領域Ｋ１を構成する分割領域は、必ずしも同一形状である必要はなく、任意の形状とすることができる。
次に、ＣＰＵ１は、図２２に示すように、室内空間Ｋから食み出した分割領域を削除することで、解析領域Ｋを室内空間Ｋに食み出すことなく収容する。この結果、図２２に示す拡大図のように、解析領域Ｋ１の境界面Ｋ１Ｆと室内空間Ｋの境界面ＫＦとの間に隙間Ｙが形成される。

次に、ＣＰＵ１は、図２３に示すように、解析領域の一部となる新たな分割領域ＲＡを隙間Ｙに充填する。
本実施形態の数値解析方法では、前述の原理説明で詳説したように、VertexとConnectivityとを持たない計算用データモデルを作成するため、分割領域の幾何学的形状に制約を課すことなく計算用データモデル作成することができる。つまり、計算用データモデルを作成するにあたり、解析領域Ｋ１を構成する分割領域は、任意の形状を取ることができる。したがって、ＣＰＵ１は、新たな分割領域ＲＡを隙間Ｙに充填するにあたり、図２３の拡大図のように、分割領域ＲＡの形状を任意に設定することができる。このため、極めて容易に隙間Ｙを分割領域ＲＡで充填することができ、例えば、ＧＵＩにより作業者が作業しなくとも、自動で分割領域ＲＡを形成することも充分に可能である。仮に、従来の有限体積法において隙間Ｙを分割領域で充填する場合には、前述のように、分割領域の幾何学的形状に対する制約の下、許容外の歪みや捩れが生じないように分割領域を配置する必要がある。この作業は、作業者の手作業となり、結果、作業者に膨大な負担を強いることとなると共に、解析時間の長期化を招くこととなる。

次にＣＰＵ１は、室内空間を示す解析領域に含まれる各分割領域内に１つのコントロールポイントを仮想的に配置する。ここでは、ＣＰＵ１は、分割領域の重心を算出し、各々の重心に対して１つのコントロールポイントを仮想的に配置する。そして、ＣＰＵ１は、コントロールポイントの配置情報、各コントロールポイントが占めるコントロールボリュームの体積（コントロールポイントが配置される分割領域の体積）を算出し、記憶装置２のデータ記憶部２ｂに一時的に記憶させる。
また、ＣＰＵ１は、分割領域同士の境界面である境界面の面積及び法線ベクトルを算出し、これらの境界面の面積及び法線ベクトルを記憶装置２のデータ記憶部２ｂに一時的に記憶させる。
また、ＣＰＵ１は、Linkを作成し、このLinkを記憶装置２のデータ記憶部２ｂに一時的に記憶させる。
そして、ＣＰＵ１は、データ記憶部２ｂに記憶された、コントロールポイントの配置情報と、各コントロールポイントが占めるコントロールボリュームの体積と、境界面の面積及び法線ベクトルと、Linkとをデータベース化することによって計算用データモデルＭを作成し、作成した計算用データモデルを記憶装置２のデータ記憶部２ｂ内に記憶させる。

なお、このようなステップＳ１ｄにおいては、室内空間を含む解析領域を同一形状の分割領域にて均等に分割し、さらに室内空間から食み出した分割領域を削除し、さらにその結果生じた解析領域と室内空間との隙間に新たな分割領域を充填することによって最終的な解析領域Ｋ２が作成される。このため、室内空間の全領域が重ならない分割領域によって充填された状態とされる。
したがって、計算用データモデルは、前述した保存則を満足するための３つの条件（ａ）〜（ｃ）を満たすものとされている。

また、本実施形態では、ステップＳ１ｄにおいて、先に分割領域を形成し、その後コントロールポイントを配置し、各コントロールポイントに対して、自らが配置された分割領域の体積を割り当てる構成を採用している。
しかしながら、本発明においては、先にコントロールポイントを解析領域に配置し、各コントロールポイントに対して後から体積を割り当てることも可能である。
具体的には、例えば、異なるコントロールポイントにぶつかるまでの半径や、結合関係にある（Linkで関連付けられた）コントロールポイントまでの距離に基づいて、各コントロールポイントに対して重み付けを行う。
ここでコントロールポイントｉの重みをｗ_ｉ、基準体積をＶ^＋とし、コントロールポイントｉに割り当てられる体積Ｖ_ｉを下式（６０）とする。

さらに、各コントロールポイントの体積Ｖ_ｉの総和は、解析領域の体積Ｖ_{ｔｏｔａｌ}と等しいため、下式（６１）が成り立つ。

この結果、基準体積Ｖ^＋は下式（６２）で求めることができる。

したがって、各コントロールポイントに割り当てる体積は、式（６１），（６２）から求めることができる。
このような方法を用いれば、プリ処理において、VertexとConnectivityとを用いることなく、計算用データモデルに持たせる分割領域の体積を求めることができる。

また、当該計算用データモデルの作成（ステップＳ１ｄ）において、ＣＰＵ１は、ＧＵＩを形成し、ＧＵＩから指令（例えば分割領域の密度を示す指令や分割領域の形状を示す指令）が入力された場合には、当該指令を反映させた処理を実行する。したがって、作業者は、ＧＵＩを操作することによって、コントロールポイントの配置や分割領域の形状を任意に調節することが可能とされている。
ただし、ＣＰＵ１は、数値解析プログラムに記憶された保存則を満足するための３つの条件に照らし合わせ、ＧＵＩから入力される指令が、当該条件から外れる場合には、その旨をディスプレイ５ａに表示させる。

続いて、ＣＰＵ１は、物性値データの設定（ステップＳ１ｅ）を行う。具体的には、ＣＰＵ１は、ＧＵＩを用いて、ディスプレイ５ａ上に物性値の入力画面を表示し、キーボード４ａあるいはマウス４ｂから入力される物性値を示す信号を物性値データＤ３としてデータ記憶部２ｂに一時的に記憶させることで物性値の設定を行う。なお、ここで言う物性値とは、室内空間における流体（すなわち空気）の特性値であり、空気の密度、粘性係数等である。

続いて、ＣＰＵ１は、境界条件データの設定（ステップＳ１ｆ）を行う。具体的には、ＣＰＵ１は、ＧＵＩを用いて、ディスプレイ５ａ上に境界条件の入力画面を表示し、キーボード４ａあるいはマウス４ｂから入力される境界条件を示す信号を境界条件データＤ１としてデータ記憶部２ｂに一時的に記憶させることで境界条件データの設定を行う。なお、ここで言う境界条件とは、室内空間の物理現象を支配する離散化支配方程式や、室内空間と外部空間との境界面に臨むコントロールポイントの特定情報、及び室内空間と外部空間との間における熱の伝熱条件等を示す。
なお、本実施形態の数値解析方法においては、室内空間における流速を数値解析により求めることを目的とするものであるため、上記離散化支配方程式として、前述のナビエ・ストークスの式に基づく離散化支配方程式（１０）及び連続の式に基づく離散化支配方程式（１１）が用いられる。
なお、これらの離散化支配方程式は、例えば、数値解析プログラムＰに予め記憶された複数の離散化支配方程式をディスプレイ５ａ上に表示された複数の離散化支配方程式から作業者がキーボード４ａやマウス４ｂを用いることによって選択される。

続いて、ＣＰＵ１は、初期条件データの設定（ステップＳ１ｇ）を行う。具体的には、ＣＰＵ１は、ＧＵＩを用いて、ディスプレイ５ａ上に初期条件の入力画面を表示し、キーボード４ａあるいはマウス４ｂから入力される初期条件を示す信号を初期条件データＤ４としてデータ記憶部２ｂに一時的に記憶させることで初期条件データの設定を行う。なお、ここ言う初期条件とは、各コントロールポイント（各分割領域）における初期流速である。

続いて、ＣＰＵ１は、計算条件データの設定（ステップＳ１ｈ）を行う。具体的には、ＣＰＵ１は、ＧＵＩを用いて、ディスプレイ５ａ上に計算条件の入力画面を表示し、キーボード４ａあるいはマウス４ｂから入力される計算条件を示す信号を計算条件データＤ２としてデータ記憶部２ｂに一時的に記憶させることで計算条件データの設定を行う。なお、ここで言う計算条件とは、ソルバ処理（ステップＳ２）における計算の条件であり、例えば、反復回数や収束基準を示す。

続いて、ＣＰＵ１は、ソルバ入力データファイルＦの作成（ステップＳ１ｉ）を行う。
具体的には、ＣＰＵ１は、ステップＳ１ｄにて作成された計算用データモデルＭと、ステップＳ１ｅで設定された物性値データＤ３と、ステップＳ１ｆで設定された境界条件データＤ１と、ステップＳ１ｇで設定された初期条件データＤ４と、ステップＳ１ｈで設定された計算条件データＤ２とをソルバ入力データファイルＦに格納することによってソルバ入力データファイルＦを作成する。なお、このソルバ入力データファイルＦは、データ記憶部２ｂに記憶される。

以上のようなプリ処理（ステップＳ１）が完了すると、ＣＰＵ１は、ソルバ処理プログラムＰ２に基づいて、図１９のフローチャートに示すソルバ処理（ステップＳ２）を実行する。
図１９に示すように、ソルバ処理（ステップＳ２）が開始されると、ＣＰＵ１は、プリ処理（ステップＳ１）で作成されたソルバ入力データファイルＦを取得する（ステップＳ２ａ）。なお、本実施形態に示す数値解析方法のように、単一の装置（本実施形態の数値解析装置Ａ）によってプリ処理及びソルバ処理を実行する場合には、既にデータ記憶部２ｂにソルバ入力データファイルＦが記憶されているため、ステップＳ２ａを省略することができる。ただし、プリ処理（ステップＳ１）とソルバ処理（ステップＳ２）とが異なる装置において実行される場合には、ネットワークやリムーバルディスクによって搬送されるソルバ入力データファイルＦを取得する必要があるため、本ステップＳ２ａを行う必要がある。

続いて、ＣＰＵ１は、ソルバ入力データの整合性を判定する（ステップＳ２ｂ）。なお、ソルバ入力データとは、ソルバ入力データファイルＦに格納されたデータを示し、計算用データモデルＭ、境界条件データＤ１、計算条件データＤ２、物性値データＤ３及び初期条件データＤ４である。
具体的には、ＣＰＵ１は、ソルバ処理において物理量計算を実行可能なソルバ入力データがソルバ入力データファイルＦに全て格納されているかを分析することによってソルバ入力データの整合性の判定を行う。
そして、ＣＰＵ１は、ソルバ入力データが不整合であると判定した場合には、ディスプレイ５ａにエラーを表示させ（ステップＳ２ｂ）、さらには不整合である部分のデータを入力するための画面を表示させる。その後、ＣＰＵ１は、ＧＵＩから入力される信号に基づいてソルバ入力データの調整を行い（ステップＳ２ｃ）、再度ステップＳ２ａを実行する。

一方、ＣＰＵ１は、ステップＳ２ｂにおいてソルバ入力データの整合性があると判定した場合には、初期計算処理（ステップＳ２ｅ）を実行する。
具体的には、ＣＰＵ１は、境界条件データＤ１に記憶された離散化支配方程式（すなわちナビエ・ストークスの式に基づく離散化支配方程式（１０）及び連続の式に基づく離散化支配方程式（１１））から離散化係数行列を作成し、さらにマトリクス計算用のデータテーブルの作成を行うことによって初期計算処理を行う。

続いて、ＣＰＵ１は、大規模粗行列方程式の組み上げ（ステップＳ２ｆ）を行う。具体的には、ＣＰＵ１は、ナビエ・ストークスの式に基づく離散化支配方程式（１０）及び連続の式に基づく離散化支配方程式（１１）から、前述の式（５８）で示すマトリックス計算用の大規模粗行列方程式の組み上げを行う。

続いて、ＣＰＵ１は、離散化支配方程式に、非圧縮性や接触等の付帯条件が存在するかの判定を行う。この付帯条件は、例えば、境界条件データとしてソルバ入力データファイルＦに格納されている。
そして、ＣＰＵ１は、離散化支配方程式に付帯条件が存在すると判定した場合には当該付帯条件の大規模行列方程式の組み込み（ステップＳ２ｈ）を実行した後に大規模行列方程式の計算（ステップＳ２ｉ）を実行する。一方、ＣＰＵ１は、離散化支配方程式に付帯条件が存在しないと判定した場合には付帯条件の大規模行列方程式の組み込み（ステップＳ２ｈ）を実行することなく大規模行列方程式の計算（ステップＳ２ｉ）を実行する。
そしてＣＰＵ１は、大規模行列方程式を例えば、ＣＧ法（共役勾配法）によって解き、前述の式（５９）を用いて解のアップデート（ステップＳ２ｊ）を行う。

続いて、ＣＰＵ１は、式（５９）の残差が収束条件に達したか否かの判定（ステップＳ２ｇ）を行う。具体的には、ＣＰＵ１は、式（５９）の残差を計算し、計算条件データＤ２に含まれる収束条件と比較し、これによって式（５９）の残差が収束条件に達したか否かの判定を行う。
そして、残差が収束条件に達していないと判定した場合には、ＣＰＵ１は、物性値のアップデートを行った後、再度ステップＳ２ｇを実行する。つまり、ＣＰＵ１は、式（５９）の残差が収束条件に達するまで、物性値のアップデートを行いながらステップＳ２ｆ〜Ｓ２ｇを繰り返し行う。

一方、残差が収束条件に達したと判定した場合には、ＣＰＵ１は、計算結果の取得を行う（ステップＳ２ｌ）。具体的には、ＣＰＵ１は、直前のステップＳ２ｉにおいて計算された物理量の解を計算結果データとしてデータ記憶部２ｂに記憶させることによって計算結果の取得を行う。

このようなソルバ処理（ステップＳ２）によって、室内空間における空気の流速が求められる。なお、このようなソルバ処理（ステップＳ２）は、本発明の物理量計算方法に相当するものである。

以上のようなソルバ処理（ステップＳ２）が完了すると、ＣＰＵ１は、ポスト処理プログラムＰ３に基づいてポスト処理（ステップＳ３）を実行する。
具体的には、例えばＣＰＵ１は、ＧＵＩから入力される指令に基づいて、計算結果データから、例えば断面コンタデータ、ベクトルデータ、等値面データ、アニメーションデータを生成し、当該データを、出力装置５に可視化させる。
また、ＣＰＵ１は、ＧＵＩから入力される指令に基づいて、室内空間の一部における定量値（計算結果）を抽出して数値やグラフとし、この数値やグラフを出力装置５に可視化させ、さらには数値やグラフをファイルとして纏めて出力する。また、ＣＰＵ１は、ＧＵＩから入力される指令に基づいて、例えば計算結果データから自動レポート作成、計算残差の表示及び分析を行ってその結果を出力する。

以上のような本実施形態の数値解析装置Ａ、数値解析方法及び数値解析プログラムによれば、プリ処理にてコントロールボリュームの体積と境界面の面積及び法線ベクトルとを有する計算用データモデルＭが作成され、ソルバ処理にて計算用データモデルＭに含まれるコントロールボリュームの体積と境界面の面積及び法線ベクトルとを用いて各コントロールボリュームにおける物理量が計算される。
つまり、前述の数値解析手法にて説明したように、本実施形態の数値解析装置Ａ、数値解析方法及び数値解析プログラムによれば、VertexとConnectivityを有する計算用データモデルを作成することなく数値解析を行うことが可能となる。このため、３次元形状データの修正あるいは変更作業に対する規制も大幅に緩和され、VertexとConnectivityを有する計算用データモデルと比較して計算用データモデルＭを遥かに容易に作成することが可能となる。したがって、計算用データモデルＭの作成における作業負担を軽減することが可能となる。
また、本実施形態の数値解析装置Ａ、数値解析方法及び数値解析プログラムにおいては、従来の数値解析手法と異なり、ソルバ処理においてVertexとConnectivityを用いてコントロールボリュームの体積及び境界面の面積及び法線ベクトルを算出する必要がない。したがって、ソルバ処理における計算負荷を低減させることが可能となる。
したがって、本実施形態の数値解析装置Ａ、数値解析方法及び数値解析プログラムによれば、計算用データモデルの作成における作業負担を軽減すると共に、ソルバ処理における計算負荷の低減を図ることが可能となる。

また、本実施形態の数値解析装置Ａ、数値解析方法及び数値解析プログラムにおいては、解析領域に分割領域を重なることなく充填させている。このため、前述した保存側を満足するための３つの条件（ａ）〜（ｃ）が満たされることとなり、保存則を満足して流速を計算することができる。

なお、本実施形態における計算用データモデルＭは、従来の有限体積法や有限要素法に用いられるメッシュデータや、従来の粒子法に用いられる粒子データ、または単なる点群データから容易に計算用データモデルを作成することができる。この場合でも、ボクセル法のような外部空間との境界領域における分割領域の修正は必要ない。
例えば、メッシュデータから計算用データモデルを作成する場合には、各要素を本実施形態における分割領域（コントロールポイントが占めるコントロールボリューム）として捉えることによって容易に行うことができる。また、粒子データから計算用データモデルを作成する場合には、解析領域に配置されたコントロールポイントを囲う閉空間を解析領域に充填して得られる閉空間を本実施形態における分割領域（コントロールポイントが占めるコントロールボリューム）として捉えることによって容易に行うことができる。

また、本実施形態の数値解析装置Ａ、数値解析方法及び数値解析プログラムによれば、前述のように、プリ処理における計算用データモデルの作業負担が大幅に減少し、またソルバ処理における計算負荷を低減させることが可能となる。
したがって、解析領域の形状が時系列的に変化する場合、すなわち解析領域が移動境界を含む場合であっても、本発明によれば、図２４のフローチャートに示すように、解析領域が形状変化するたびにプリ処理とソルバ処理とを繰り返し行うことによって、現実的な時間内で物理量の計算が可能となる。

また、本実施形態においては、物理量計算において、境界面上における流速を算出する場合に、当該境界面を挟むコントロールポイントにおける流速を当該コントロールポイントから境界面までの距離の比率に応じて重み付け平均した値を境界面上における流速としても良い。
このような場合には、境界面における流速を、単純に当該境界面を挟むコントロールポイントの流速の平均とする場合と比較して、より正確に室内空間の流速を求めることが可能となる。ただし、このような場合には、計算用データモデルが、境界面がコントロールポイントとコントロールポイントとを結ぶ線分のどの内分点に存在するかを示す比率αをデータとして有している必要がある。このため、隣り合うコントロールボリュームの内部に配置された各コントロールポイントから当該コントロールボリュームに挟まれる境界面までの距離の比率をさらに含む計算用データモデルＭをプリ処理にて作成し、この比率に基づいて境界面における流速をソルバ処理にて計算する。

なお、物理量計算において、コントロールポイント間を結ぶベクトルと当該コントロールポイントに挟まれた境界面が直交する場合には、法線ベクトルをコントロールポイント間を結ぶ距離ベクトルの単位ベクトルに置き換えて物理量計算することができる。
ただし、このような場合には、計算用データモデルに上記距離ベクトルあるいは当該距離ベクトルを算出するためのコントロールポイントの位置座標を示すデータを備えさせる必要がある。このため、隣り合うコントロールボリュームの内部に配置された各コントロールポイントを結ぶ距離ベクトルあるいは該距離ベクトルを算出するためのコントロールポイントの座標をさらに含む計算用データモデルＭをプリ処理にて作成し、距離ベクトルが当該コントロールポイントに挟まれた境界面と直交する場合に、距離ベクトルを境界面の法線ベクトルに置き換えて流速を計算する。

以上、添付図面を参照しながら本発明の好適な実施形態について説明したが、本発明は、上記実施形態に限定されないことは言うまでもない。前述した実施形態において示した各構成部材の諸形状や組み合わせ等は一例であって、本発明の主旨から逸脱しない範囲において設計要求等に基づき種々変更可能である。

上記実施形態においては、運動量保存の方程式の変形例であるナビエ・ストークスの式及び連続の式から導出した離散化支配方程式を用いて空気の流速を数値解析によって求める構成について説明した。
しかしながら、本発明はこれに限定されるものではなく、質量保存の方程式、運動量保存の方程式、角運動量保存の方程式、エネルギ保存の方程式、移流拡散方程式及び波動方程式の少なくともいずれかから導出した離散化支配方程式を用いて物理量を数値解析によって求めることが可能である。

また、上記実施形態においては、本発明の境界面特性量として、境界面の面積と境界面の法線ベクトルとを用いる構成について説明した。
しかしながら、本発明はこれに限定されるものではなく、境界面特性量として他の量（例えば境界面の周長）を用いることもできる。

また、上記実施形態においては、保存則を満足するために前述の３つの条件を満たすように計算用データモデルを作成する構成について説明した。
しかしながら、本発明はこれに限定されるものではなく、保存則を満足させる必要がない場合には、計算用データモデルを必ずしも前述の３つの条件を満たすように作成する必要はない。

また、上記実施形態においては、分割領域の体積を、当該分割領域の内部に配置されるコントロールポイントが占めるコントロールボリュームの体積として捉えた構成について説明した。
しかしながら、本発明はこれに限定されるものではなく、分割領域の内部に対してコントロールポイントを配置する必要は必ずしもない。このような場合には、コントロールポイントが占めるコントロールボリュームの体積を分割領域の体積に置き換えることによって数値解析を行うことができる。

また、上記実施形態においては、数値解析プログラムＰがＤＶＤメディアＸに記憶されて搬送可能な構成について説明した。
しかしながら、本発明はこれに限定されるものではなく、数値解析プログラムＰを他のリムーバブルメディアに記憶させて搬送可能とする構成を採用することもできる。
また、プリ処理プログラムＰ１とソルバ処理プログラムＰ２とを別々のリムーバブルメディアに記憶させて搬送可能とすることもできる。また、数値解析プログラムＰは、ネットワークを介して伝達することも可能である。

また、上記実施形態においては、境界面が平面である構成について説明した。しかしながら、本発明はこれに限定されるものではなく、境界面が曲面であっても良い。

本出願を詳細にまた特定の実施態様を参照して説明したが、本発明の精神と範囲を逸脱することなく様々な変更や修正を加えることができることは当業者にとって明らかである。
本出願は、２００９年６月２５日出願の日本特許出願（特願２００９-１５０９４５）に基づくものであり、その内容はここに参照として取り込まれる。

Ａ……数値解析装置（物理量計算装置）
Ｐ……数値解析プログラム
Ｐ１……プリ処理プログラム
Ｐ２……ソルバ処理プログラム（物理量計算プログラム）
Ｐ３……ポスト処理プログラム
Ｍ……計算用データモデル
Ｖ_i……粒子iが仮想的に占める体積
Ｓ_ｉｊ……粒子iと粒子jの境界面積
ｒ_ｉｊ……ｉからjの距離ベクトル
ｎ_ｉｊ……Ｓ_ｉｊの法線ベクトル
１……ＣＰＵ
２……記憶装置
２ａ……プログラム記憶部
２ｂ……データ記憶部
３……ＤＶＤドライブ
４……入力装置
４ａ……キーボード
４ｂ……マウス
５……出力装置
５ａ……ディスプレイ
５ｂ……プリンタ

上記課題を解決するために、本発明の一態様は、物理現象を数値的に解析する数値解析方法においてコンピュータが物理量を計算する物理量計算方法であって、直交格子形状のみによらない複数の分割領域に分割された解析領域における物理量をコンピュータが計算する物理量計算工程を含み、該物理量計算工程にて、コンピュータが、前記分割領域の頂点の座標（Vertex）及び該頂点の連結情報（Connectivity）を必要としない量のみを使用すると共に重み付き残差積分法に基づいて導出された離散化された支配方程式を記憶している記憶装置から前記支配方程式を読み出し、コンピュータが、各前記分割領域の体積及び隣り合う前記分割領域同士の境界面の特性を示す境界面特性量を、前記分割領域の頂点の座標（Vertex）及び該頂点の連結情報（Connectivity）を必要としない量として有する計算用データモデルと、前記支配方程式と、を用いて前記物理量を計算するという構成を採用する。
また、本発明の一態様は、コンピュータが計算用データモデルを作成するプリ処理と、コンピュータが該計算用データモデル及び支配方程式に基づく物理量計算方法によって解析領域における物理量を計算するソルバ処理とを有する数値解析方法であって、コンピュータが、前記ソルバ処理において上記物理量計算方法を用いることを特徴とする数値解析方法である。
また、本発明の一態様は、物理現象を数値的に解析する数値解析方法において物理量を計算する物理量計算方法をコンピュータに実行させる物理量計算プログラムであって、複数の分割領域に分割された解析領域における物理量を計算する物理量計算工程を含み、該物理量計算工程にて、前記分割領域の頂点の座標（Vertex）及び該頂点の連結情報（Connectivity）を必要としない量のみを使用すると共に重み付き残差積分法に基づいて導出された離散化された支配方程式と、各前記分割領域の体積及び隣り合う前記分割領域同士の境界面の特性を示す境界面特性量を、前記分割領域の頂点の座標（Vertex）及び該頂点の連結情報（Connectivity）を必要としない量として有する計算用データモデルと、を用いて前記物理量を計算することを特徴とする物理量計算方法を前記コンピュータに実行させる物理量計算プログラムである。
また、本発明の一態様は、計算用データモデルを作成するプリ処理と、該計算用データモデル及び支配方程式に基づく物理量計算方法によって解析領域における物理量を計算するソルバ処理とを有する数値解析方法をコンピュータに実行させる数値解析プログラムであって、前記ソルバ処理において上記の物理量計算プログラムによる処理を用いることを特徴とする数値解析方法を前記コンピュータに実行させる数値解析プログラムである。
また、本発明の一態様は、物理現象を数値的に解析する数値解析方法において物理量を計算する物理量計算装置であって、複数の分割領域に分割された解析領域における物理量を計算する演算手段と、前記分割領域の頂点の座標（Vertex）及び該頂点の連結情報（Connectivity）を必要としない量のみを使用すると共に重み付き残差積分法に基づいて導出された離散化された支配方程式を記憶する記憶装置と、を備え、前記演算手段は、各前記分割領域の体積及び隣り合う前記分割領域同士の境界面の特性を示す境界面特性量を、前記分割領域の頂点の座標（Vertex）及び該頂点の連結情報（Connectivity）を必要としない量として有する計算用データモデルと、前記記憶装置に記憶されている前記支配方程式と、を用いて前記物理量を計算することを特徴とする物理量計算装置である。
また、本発明の一態様は、計算用データモデルを作成するプリ処理を行う計算用データモデル作成手段と、上記の物理量計算装置を含み、該計算用データモデル及び支配方程式に基づく前記物理量計算装置による物理量計算によって解析領域における物理量を計算するソルバ処理を行う物理量計算手段と、を備える数値解析装置である。

Claims (13)

1. 物理現象を数値的に解析する数値解析方法において物理量を計算する物理量計算方法であって、
   複数の分割領域に分割された解析領域における物理量を計算する物理量計算工程を含み、
   該物理量計算工程にて、前記分割領域の頂点の座標（Vertex）及び該頂点の連結情報（Connectivity）を必要としない量のみを使用すると共に重み付き残差積分法に基づいて導出された離散化された支配方程式と、各前記分割領域の体積及び隣り合う前記分割領域同士の境界面の特性を示す境界面特性量を、前記分割領域の頂点の座標（Vertex）及び該頂点の連結情報（Connectivity）を必要としない量として有する計算用データモデルと、を用いて前記物理量を計算することを特徴とする物理量計算方法。
2. 前記支配方程式は、質量保存の方程式、運動量保存の方程式、エネルギ保存の方程式、移流拡散方程式、及び波動方程式から導出されている請求項１記載の物理量計算方法。
3. 前記境界面特性量は、前記境界面の面積と前記境界面の法線ベクトルである請求項１または２記載の物理量計算方法。
4. 計算用データモデルを作成するプリ処理と、該計算用データモデル及び支配方程式に基づく物理量計算方法によって解析領域における物理量を計算するソルバ処理とを有する数値解析方法であって、
   前記ソルバ処理において請求項１〜３のいずれか１項に記載の物理量計算方法を用いることを特徴とする数値解析方法。
5. 前記プリ処理にて、
   全分割領域の体積の総和が解析領域の体積と一致するという条件と、
   前記境界面の面積が一致するという条件及び法線ベクトルが当該境界面を挟む一方の分割領域から見た場合と他方の分割領域から見た場合とで絶対値が一致するという条件と、 前記分割領域を通る無限に広い投影面Ｐの任意の向きの単位法線ベクトルが［ｎ］_ｐ、境界面の面積がＳ_ｉ、境界面の単位法線ベクトルが［ｎ］_ｉ、分割領域の面の総数がｍ、上記［］で囲まれた文字がベクトルを示す太字であるときに下式（１）が成り立つという条件と、
   が満足されるように前記分割領域が形成される請求項４記載の数値解析方法。
   

   
6. 前記解析領域に前記分割領域を重なることなく充填させることによって前記条件を満足させる請求項５記載の数値解析方法。
7. 隣り合う前記分割領域の内部に配置されるコントロールポイントから当該分割領域に挟まれる前記境界面までの距離の比率をさらに含む前記計算用データモデルを前記プリ処理にて作成し、前記比率を用いて前記境界面における物理量を前記ソルバ処理にて計算する請求項４〜６のいずれか１項に記載の数値解析方法。
8. 隣り合う前記分割領域の内部に配置される各コントロールポイントを結ぶ距離ベクトルあるいは該距離ベクトルを算出するための前記コントロールポイントの座標をさらに含む前記計算用データモデルを前記プリ処理にて作成し、前記距離ベクトルが当該コントロールポイントに挟まれた前記境界面と直交する場合に、前記距離ベクトルを前記境界面特性量の１つとして前記物理量を計算する請求項４〜７のいずれか１項に記載の数値解析方法。
9. 前記解析領域が移動境界を含む場合に、前記移動境界が移動するたびに前記プリ処理及び前記ソルバ処理を行う請求項４〜８のいずれか１項に記載の数値解析方法。
10. 請求項１〜３のいずれか１項に記載の物理量計算方法をコンピュータに実行させることを特徴とする物理量計算プログラム。
11. 請求項４〜９のいずれか１項に記載の数値解析方法をコンピュータに実行させることを特徴とする数値解析プログラム。
12. 請求項１〜３のいずれか１項に記載の物理量計算方法を行うことを特徴とする物理量計算装置。
13. 請求項４〜９のいずれか１項に記載の数値解析方法を行うことを特徴とする数値解析装置。

Images