(57)【要約】【課題】移流拡散問題においてより高精度に数値解析を行えるようにする。【解決手段】解析対象の場に設定された計算格子における計算格子点の各々について、可変な場の値を有して当該計算格子点へ移流する仮想粒子の軌跡を、時刻を遡って追跡する逆追跡を行う逆追跡ステップと、前記計算格子点の各々について、前記逆追跡の終了位置における場の値と、前記逆追跡の結果とに基づいて、拡散しながら前記逆追跡の終了位置から当該計算格子点へ到達した際に前記仮想粒子が有する場の値を、当該計算格子点における場の値として取得する場の値取得ステップと、を具備する。【選択図】図１
2014119882.tif 000002
【特許請求の範囲】
【請求項１】
解析対象の場に設定された計算格子における計算格子点の各々について、可変な場の値を有して当該計算格子点へ移流する仮想粒子の軌跡を、時刻を遡って追跡する逆追跡を行う逆追跡ステップと、
前記計算格子点の各々について、前記逆追跡の終了位置における場の値と、前記逆追跡の結果とに基づいて、拡散しながら前記逆追跡の終了位置から当該計算格子点へ到達した際に前記仮想粒子が有する場の値を、当該計算格子点における場の値として取得する場の値取得ステップと、
を具備することを特徴とする解析方法。
【請求項２】
前記逆追跡ステップでは、前記仮想粒子の軌跡を前記計算格子点の各々について求め、
前記場の値取得ステップでは、前記計算格子点の各々についての前記仮想粒子の軌跡と拡散方程式とに基づいて得られる連立方程式から、前記計算格子点の各々における場の値を取得する、
ことを特徴とする請求項１に記載の解析方法。
【請求項３】
前記逆追跡ステップでは、前記仮想粒子の軌跡を前記計算格子点の各々について求め、
前記場の値取得ステップは、
前記計算格子点の各々について、前記仮想粒子の軌跡に沿って拡散方程式を積分して、前記計算格子点の各々における場の値を算出する場の値算出ステップと、
前記場の値算出ステップにて算出された場の値が収束したか否か判定する収束性判定ステップと、
を具備し、
前記場の値算出ステップにて算出された場の値が収束したと前記収束性判定ステップにて判定するまで、前記場の値算出ステップにおける処理を繰り返す、ことを特徴とする請求項１に記載の解析方法。
【請求項４】
前記場の値取得ステップは、
前記計算格子点の各々について、当該計算格子点に到達する前記仮想粒子が前記逆追跡の終了位置において有する場の値を当該計算格子点に設定する場の値設定ステップと、
前記計算格子点の各々について、前記仮想粒子が前記逆追跡の終了位置から当該計算格子点に到達するのに要する滞留時間に応じて当該計算格子点における場の値を更新する場の値更新ステップと、
を具備することを特徴とする請求項１に記載の解析方法。
【請求項５】
解析対象の場に設定された計算格子における計算格子点の各々について、可変な場の値を有して当該計算格子点へ移流する仮想粒子の軌跡を、時刻を遡って追跡する逆追跡を行う逆追跡部と、
前記計算格子点の各々について、前記逆追跡の終了位置における場の値と、前記逆追跡の結果とに基づいて、拡散しながら前記逆追跡の終了位置から当該計算格子点へ到達した際の前記仮想粒子における場の値を、当該計算格子点における場の値として取得する場の値取得部と、
を具備することを特徴とする解析装置。
【請求項６】
コンピュータに、
解析対象の場に設定された計算格子における計算格子点の各々について、可変な場の値を有して当該計算格子点へ移流する仮想粒子の軌跡を、時刻を遡って追跡する逆追跡を行う逆追跡ステップと、
前記計算格子点の各々について、前記逆追跡の終了位置における場の値と、前記逆追跡の結果とに基づいて、拡散しながら前記逆追跡の終了位置から当該計算格子点へ到達した際に前記仮想粒子が有する場の値を、当該計算格子点における場の値として取得する場の値取得ステップと、
を実行させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、解析方法、解析装置およびプログラムに関する。
【背景技術】
【０００２】
数値解析分野の主たる解析対象の１つに移流拡散問題がある。移流拡散問題は、流体を伴う伝熱・混合問題を始めとした広い範囲に適用されており、移流拡散問題における数値解析手法の高度化や効率化が、産業利用等の観点から重要視されている。一方、一般的に移流拡散問題の解析に対する要求精度が非常に高く、ごく一部の理想的なケースを除いて、現象を正しく予測することは困難であった。
【０００３】
移流拡散問題において多用されている数値解析手法に格子法がある（例えば、特許文献１）。格子法は、オイラー的な考え方によって物理現象を捉え、空間上に設けた計算格子によって現象を離散的に表現する解析方法の総称である。
また、移流拡散問題において格子法と対比される数値解析手法に粒子法がある（例えば、特許文献２）。粒子法は、ラグランジュ的な考え方を取り入れた解析手法である。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２００４−１３４４２号公報
【特許文献２】特開平７−３３４４８４号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
格子法を用いた移流拡散問題の数値解析では、数値拡散に起因する数値誤差の発生により、解析結果の精度が低下してしまう。ここでいう数値拡散とは、移流拡散方程式の対流項の離散化に伴って付加される偽拡散のことであり、実現象の分子拡散と比較して数値拡散の方が支配的となってしまう解析結果も珍しくない。
一方、粒子法では、粒子群を用いて現象を表現することから、根本的に数値拡散が生じない。しかしながら、粒子法は、気液二相流に代表されるような自由表面を伴う現象の数値解析に優れるものの、勾配やラプラシアンなど空間における物理量の分布を定量的に評価する手段に乏しい。その解決策として用いられる粒子間相互作用モデルに対しては、移流拡散問題において分子拡散を正確に評価できないという指摘がなされている。
【０００６】
本発明は、このような事情に鑑みてなされたもので、その目的は、移流拡散問題においてより高精度に数値解析を行うことのできる解析方法、解析装置およびプログラムを提供することにある。
【課題を解決するための手段】
【０００７】
この発明は上述した課題を解決するためになされたもので、本発明の一態様による解析方法は、解析対象の場に設定された計算格子における計算格子点の各々について、可変な場の値を有して当該計算格子点へ移流する仮想粒子の軌跡を、時刻を遡って追跡する逆追跡を行う逆追跡ステップと、前記計算格子点の各々について、前記逆追跡の終了位置における場の値と、前記逆追跡の結果とに基づいて、拡散しながら前記逆追跡の終了位置から当該計算格子点へ到達した際に前記仮想粒子が有する場の値を、当該計算格子点における場の値として取得する場の値取得ステップと、を具備することを特徴とする。
【０００８】
また、本発明の一態様による解析方法は、上述の解析方法であって、前記逆追跡ステップでは、前記仮想粒子の軌跡を前記計算格子点の各々について求め、前記場の値取得ステップでは、前記計算格子点の各々についての前記仮想粒子の軌跡と拡散方程式とに基づいて得られる連立方程式から、前記計算格子点の各々における場の値を取得する、ことを特徴とする。
【０００９】
また、本発明の一態様による解析方法は、上述の解析方法であって、前記逆追跡ステップでは、前記仮想粒子の軌跡を前記計算格子点の各々について求め、前記場の値取得ステップは、前記計算格子点の各々について、前記仮想粒子の軌跡に沿って拡散方程式を積分して、前記計算格子点の各々における場の値を算出する場の値算出ステップと、前記場の値算出ステップにて算出された場の値が収束したか否か判定する収束性判定ステップと、を具備し、前記場の値算出ステップにて算出された場の値が収束したと前記収束性判定ステップにて判定するまで、前記場の値算出ステップにおける処理を繰り返す、ことを特徴とする。
【００１０】
また、本発明の一態様による解析方法は、上述の解析方法であって、前記場の値取得ステップは、前記計算格子点の各々について、当該計算格子点に到達する前記仮想粒子が前記逆追跡の終了位置において有する場の値を当該計算格子点に設定する場の値設定ステップと、前記計算格子点の各々について、前記仮想粒子が前記逆追跡の終了位置から当該計算格子点に到達するのに要する滞留時間に応じて当該計算格子点における場の値を更新する場の値更新ステップと、を具備することを特徴とする。
【００１１】
また、本発明の一態様による解析装置は、解析対象の場に設定された計算格子における計算格子点の各々について、可変な場の値を有して当該計算格子点へ移流する仮想粒子の軌跡を、時刻を遡って追跡する逆追跡を行う逆追跡部と、前記計算格子点の各々について、前記逆追跡の終了位置における場の値と、前記逆追跡の結果とに基づいて、拡散しながら前記逆追跡の終了位置から当該計算格子点へ到達した際の前記仮想粒子における場の値を、当該計算格子点における場の値として取得する場の値取得部と、を具備することを特徴とする。
【００１２】
また、本発明の一態様によるプログラムは、コンピュータに、解析対象の場に設定された計算格子における計算格子点の各々について、可変な場の値を有して当該計算格子点へ移流する仮想粒子の軌跡を、時刻を遡って追跡する逆追跡を行う逆追跡ステップと、前記計算格子点の各々について、前記逆追跡の終了位置における場の値と、前記逆追跡の結果とに基づいて、拡散しながら前記逆追跡の終了位置から当該計算格子点へ到達した際に前記仮想粒子が有する場の値を、当該計算格子点における場の値として取得する場の値取得ステップと、を実行させるためのプログラムである。
【発明の効果】
【００１３】
本発明によれば、移流拡散問題において従来方法に比べて短時間かつ少ない記憶容量でより高精度に数値解析を行うことができる。
【図面の簡単な説明】
【００１４】
【図１】本発明の第１の実施形態における解析装置の機能構成を示す概略ブロック図である。
【図２】同実施形態における逆追跡部が行う仮想粒子の軌跡の逆追跡の例を示す説明図である。
【図３】同実施形態において、解析装置が行う処理の手順を示すフローチャートである。
【図４】本発明の第２の実施形態における解析装置の機能構成を示す概略ブロック図である。
【図５】同実施形態において、解析装置が行う処理の手順を示すフローチャートである。
【図６】同実施形態における場の値算出部が行うスカラー量算出の例を示す説明図である。
【図７】本発明の第３の実施形態における解析装置の機能構成を示す概略ブロック図である。
【図８】同実施形態において解析装置が行う処理の手順を示すフローチャートである。
【図９】流体混合に用いられる典型的なＴ型ミキサーの一例を示す。
【図１０】本発明の第１の解析例での濃度場の解析における仮想粒子の逆追跡の例を示す説明図である。
【図１１】本発明の第１の解析例でのターゲット平面における混合度ＭＩの計算結果の例を、計算格子点間間隔を変えながら格子法と比較した結果を示す説明図である。
【図１２】格子法と本発明の第１の解析例における手法との、最大空間解像度における濃度分布の解析結果の例を示す説明図である。
【図１３】本発明の第２の解析例でのターゲット平面Ａにおける混合度ＭＩの計算結果の例を、計算格子点間間隔を変えながら格子法と比較した結果を示す説明図である。
【図１４】格子法と本発明の第１の解析例における手法との、濃度分布の解析結果の例を示す説明図である。
【図１５】本発明の第１の解析例における解析結果と、第２の解析例における解析結果と、格子法による解析結果との比較例を示す説明図である。
【発明を実施するための形態】
【００１５】
＜第１の実施形態＞
以下、図面を参照して、本発明の実施の形態について説明する。
なお、明細書の記載において、ベクトルないし行列を示す矢印の表記、および、２段階目以降の上付きや下付きを省略する。
図１は、本発明の第１の実施形態における解析装置の機能構成を示す概略ブロック図である。同図において、解析装置１００は、前提条件取得部１１０と、逆追跡部１２０と、場の値取得部１３０と、解析結果出力部１４０とを具備する。場の値取得部１３０は、場の値設定部１３１と、場の値更新部１３２とを具備する。
【００１６】
解析装置１００は、数値解析を行って移流拡散問題における場の値を求める。解析装置１００は、例えばコンピュータにて構築される。あるいは、解析装置１００が専用の回路にて構築されるなど、コンピュータ以外で構築されていてもよい。
【００１７】
前提条件取得部１１０は、数値解析を行う前提条件を取得する。例えば、前提条件取得部１１０は、数値解析を行う前提条件として、解析対象となる空間に設定する計算格子の規模および計算格子点間間隔や、移流速度場（場の値が移流する速度を示す速度場）のデータなどを取得する。
前提条件取得部１１０が数値解析を行う前提条件を取得する形態として様々な形態を用いることができる。あるいは、前提条件取得部１１０が、数値解析を行う前提条件を予め記憶しておくようにしてもよいし、解析装置１００が数値解析を行う際に、数値解析を行う前提条件のユーザ入力を受けるようにしてもよい。また、前提条件取得部１１０が、移流速度場を解析によって求めるようにしてもよいし、予め解析されている移流速度場のデータを取得するようにしてもよい。
【００１８】
逆追跡部１２０は、解析対象の場に設定された計算格子における計算格子点の各々について、当該計算格子点へ移流（対流または輸送ともいう）する仮想粒子の軌跡を逆追跡する。特に、逆追跡部１２０は、流入境界など場の値が既知の領域まで仮想粒子の軌跡を逆追跡する。また、逆追跡部１２０は、仮想粒子が逆追跡の終了位置から計算格子点に到達するのに要する滞留時間を求める。
逆追跡部１２０が行う処理は、逆追跡ステップにおける処理の一例に該当する。
ここでいう仮想粒子とは、可変な場の値を有して計算格子点へ移流する、質量（慣性）を持たない仮想的な粒子（流体粒子）である。また、ここでいう逆追跡とは、時刻を遡って軌跡を追跡することである。従って、仮想粒子の軌跡の逆追跡では、仮想粒子が移流する向きと逆向きに軌跡を辿っていく。
【００１９】
図２は、逆追跡部１２０が行う仮想粒子の軌跡の逆追跡の例を示す説明図である。同図において、点Ｐ１１は、計算格子点を示し、線Ｌ１１は、仮想粒子の軌跡を示している。
仮想粒子が、移流により計算格子点（点Ｐ１１）へ到達しているのに対し、逆追跡部１２０は、矢印Ｖ１１にて示すように、仮想粒子が移流する向きと逆向きに軌跡（線Ｌ１１）を辿っていく。
【００２０】
仮想粒子の軌跡を逆追跡することで、計算格子点における、対流項（移流項ともいう）による場の値を求め得る。また、仮想粒子の軌跡を逆追跡することで、仮想粒子が移流する際の拡散を算出することができる。これにより、計算格子点における場の値に対して、拡散項の影響を反映させることができる。
また、仮想粒子は仮想的な粒子なので、実在する粒子による移流拡散に限らず、熱伝導などエネルギーの移流拡散についても、仮想粒子を用いて解析することができる。
【００２１】
場の値取得部１３０は、計算格子点の各々について、逆追跡の終了位置における場の値と、逆追跡の結果とに基づいて、拡散しながら逆追跡の終了位置から計算格子点へ到達した際に仮想粒子が有する場の値を、当該計算格子点における場の値として取得する。場の値取得部１３０が行う処理は、場の値取得ステップにおける処理の一例に該当する。
【００２２】
場の値設定部１３１は、計算格子点の各々について、当該計算格子点に到達する仮想粒子が逆追跡の終了位置において有する場の値を当該計算格子点における場の値に設定する。より具体的には、逆追跡部１２０が、流入境界など場の値が既知の領域まで、仮想粒子の軌跡を逆追跡しておく。そして、場の値取得部１３０は、逆追跡の終了位置において既知となっている場の値を、計算格子点における場の値として設定する。これにより、場の値設定部１３１は、計算格子点における、対流項による場の値を設定する。
場の値設定部１３１が行う処理は、場の値設定ステップにおける処理の一例に該当する。
【００２３】
場の値更新部１３２は、計算格子点の各々について、仮想粒子が逆追跡の終了位置から当該計算格子点に到達するのに要する滞留時間に応じて、場の値設定部１３１が設定した当該計算格子点における場の値を更新する。これにより、場の値更新部１３２は、計算格子点における場の値に対して、拡散項の影響を反映させる。
場の値更新部１３２が行う処理は、場の値更新ステップにおける処理の一例に該当する。
【００２４】
解析結果出力部１４０は、場の値取得部１３０が取得した場の値を出力する。解析結果出力部１４０が場の値を出力する態様として様々な態様を用いることができる。例えば、解析結果出力部１４０が表示画面を有し、得られた場の値を、グラフ表示または数値による表示など視覚的に表示するようにしてもよい。あるいは、解析結果出力部１４０が、得られた場の値を示すデータを、解析装置に接続されている他の装置へ出力（送信）するようにしてもよい。
【００２５】
次に、図３を参照して、解析装置１００の動作について説明する。なお、以下では、解析装置１００がスカラー場の解析を行う場合を例に説明するが、解析装置１００が解析を行う場はスカラー場に限らない。また、解析装置１００は、様々な次元の空間を解析対象とすることができる。他の実施形態においても同様である。
【００２６】
以下、解析対象となる空間に速度場およびスカラー場が形成されており、スカラー場を構成するスカラー量が速度場に乗って移流しながら拡散する場合に、スカラー量の分布を求める問題を考える。より具体的には、解析装置１００は、解析対象となる空間に計算格子を設定し、計算格子点の各々におけるスカラー量を求める。なお、以下では総計算格子点数をＮとし、計算格子点を番号ｋ（０≦ｋ≦Ｎ−１）にて識別する。
【００２７】
図３は、解析装置１００が行う処理の手順を示すフローチャートである。
同図の処理において、まず前提条件取得部１１０が、計算格子の設定条件や移流速度場のデータなど、数値解析を行う前提条件を取得する（ステップＳ１０１）。
次に、逆追跡部１２０が、仮想粒子の軌跡の逆追跡を行い、逆追跡の結果に基づいて場の値設定部１３１が対流項の評価を行う（ステップＳ１０２）。以下、逆追跡部１２０が行う逆追跡および場の値設定部１３１が行う対流項の評価について、より詳細に説明する。
【００２８】
まず、スカラー量既知の領域から計算格子点ｋへ速度場に乗って移流する仮想粒子を想定する。逆追跡部１２０は、この仮想粒子の軌跡ｘｋ（ｔ）を、式（１）を時刻方向逆向きに解くことによって求める。
【００２９】
【数１】
2014119882.tif 000003
【００３０】
但し、ｔは時刻を示し、ｕは速度場における速度を示す。なお、速度場は、定常的なものであってもよいし、時刻に応じて変化してもよい。さらには、速度場の変化に伴って、スカラー場が時刻に応じて変化してもよい。すなわち、解析装置１００は、定常的なスカラー場と非定常なスカラー場とのいずれも解析することができる。非定常なスカラー場を解析する場合、解析装置１００は、ある時刻（例えば、仮想粒子が計算格子点に到達する時刻ｔ０）におけるスカラー量の分布を求める。
【００３１】
式（１）の解法として、同式を離散的に解くことが考えられる。その際、逆追跡部１２０が行う離散化スキームとして、オイラー陽解法、完全陰解法、クランクニコルソン法、ルンゲクッタ法、ルンゲクッタフェールベルグ法またはアダムスバッシュフォース法など様々なスキームを用いることができる。
【００３２】
また、速度ｕが離散的に定義されている場合、近傍の速度定義点からの速度ｕの内挿値を用いることが考えられる。その際、逆追跡部１２０が行う内挿スキームとして、一次精度内挿、バイリニア内挿、トリリニア内挿、二次精度内挿または三次精度内挿など様々なスキームを用いることができる。
【００３３】
なお、計算誤差により得られる軌跡が非現実的な軌跡となる場合、逆追跡部１２０が、それを防ぐ手段（例えば速度場の補正）を講じるようにしてもよい。例えば数値誤差の影響により、仮想粒子が、速度条件をノンスリップとする境界（壁面）に達した場合に、逆追跡部１２０が、逆追跡に用いる速度の補正を行って、仮想粒子の停止を回避するようにしてもよい。
【００３４】
逆追跡部１２０が、仮想粒子がスカラー量既知の領域（例えば流入境界）へ到達するまで逆追跡を行うことで、場の値設定部１３１は、拡散を考慮しない場合の計算格子点ｋのスカラー量を取得することができる。スカラー量既知の領域の一例として、スカラー量を予め与えられている流入境界が挙げられる。
逆追跡部１２０は、各計算格子点について上記の逆追跡を行う。そして場の値設定部１３１は、各計算格子点における逆追跡の結果に基づいて、拡散を考慮しない場合のスカラー分布を求める。
【００３５】
なお、逆追跡部１２０は、流入境界などスカラー量既知の領域を逆追跡終了位置として、当該逆追跡終了位置における仮想粒子のスカラー量と、仮想粒子が当該逆追跡終了位置から計算格子点へ到達するのに要する時間（滞留時間）を取得するものであればよい。従って、逆追跡部１２０は、必ずしも、仮想粒子の軌跡の座標を求めるなど具体的な軌跡を求める必要は無い。
【００３６】
次に、場の値更新部１３２は、モデルを用いて拡散項の評価を行う（ステップＳ１０３）。以下、場の値更新部１３２が行う拡散項の評価について、より詳細に説明する。
場の値更新部１３２は、場の値設定部１３１が取得した、拡散を考慮しない場合のスカラー分布を初期（ξ＝０）における値として、式（２）に示すモデル式をξ＝［０，１］の範囲で数値的に積分する。
【００３７】
【数２】
2014119882.tif 000004
【００３８】
但し、φはスカラー場を示す。また、τは仮想粒子の滞留時間を示す。すなわち、τは、仮想粒子が、スカラー量既知の領域から計算格子点ｋに到達するまでに要する時間を示す。また、Ｄは拡散係数を示し、▽２はラプラシアン演算子を示す。また、ξは仮に設定した変数である。
場の値更新部１３２は、積分結果のξ＝１において、拡散を考慮したスカラー分布を得る。なお、式（２）を用いたスカラー量の算出では、拡散を考慮しない場合のスカラー分布を初期状態とし、仮想粒子の軌跡における拡散を計算格子点ｋにおける拡散で近似して、仮想粒子の滞留時間τに相当する拡散を模擬している。この点において、場の値更新部１３２は、モデルを用いて拡散項を評価している。
【００３９】
なお、式（２）に対する積分を行う際、場の値更新部１３２が行う離散化スキームとして、オイラー陽解法、完全陰解法、クランクニコルソン法、ルンゲクッタ法、ルンゲクッタフェールベルグ法またはアダムスバッシュフォース法など様々な離散化スキームを用いることができる。
【００４０】
ステップＳ１０３の後、解析結果出力部１４０は、解析結果としてのスカラー場を出力する（ステップＳ１０４）。より具体的には、解析結果出力部１４０は、各計算格子点におけるスカラー量を出力する。
その後、図３の処理を終了する。
【００４１】
以上のように、逆追跡部１２０は、計算格子点へ移流する仮想粒子の軌跡を逆追跡し、場の値取得部１３０は、拡散しながら逆追跡の終了位置から計算格子点へ到達した際に仮想粒子が有する場の値を、当該計算格子点における場の値として取得する。
これにより、解析装置１００では、仮想粒子が有する場の値を計算格子点に内挿する必要が無いので数値拡散が生じず、より高精度に数値解析を行うことができる。
【００４２】
また、場の値設定部１３１は、計算格子点の各々について、当該計算格子点に到達する仮想粒子が逆追跡の終了位置において有する場の値を当該計算格子点に設定する。そして、場の値更新部１３２は、計算格子点の各々について、仮想粒子が逆追跡の終了位置から当該計算格子点に到達するのに要する滞留時間に応じて当該計算格子点における場の値を更新する。
これにより、場の値設定部１３１は、逆追跡の終了位置における場の値を計算格子点における場の値として設定するという簡単な処理によって対流項の評価を行うことができる。また、場の値更新部１３２は、例えば、仮想粒子が、当該仮想粒子の到達した計算格子点に上記の滞留時間だけ滞留した場合の場の値を算出するなど、簡単な処理によって拡散項の評価を行うことができる。
【００４３】
＜第２の実施形態＞
図４は、本発明の第２の実施形態における解析装置の機能構成を示す概略ブロック図である。同図において、解析装置２００は、前提条件取得部１１０と、逆追跡部２２０と、場の値取得部２３０と、解析結果出力部１４０とを具備する。場の値取得部２３０は、場の値算出部２３１と、収束判定部２３２とを具備する。なお、図４において、図１の各部に対応して同様の機能を有する部分には同一の符号（１１０、１４０）を付して説明を省略する。
【００４４】
解析装置２００は、解析装置１００と同様、数値解析を行って移流拡散問題における場の値を求める。但し、解析装置２００が行う具体的な処理は、解析装置１００が行う処理と異なる。解析装置２００は、例えばコンピュータにて構築される。あるいは、解析装置２００が専用の回路にて構築されるなど、コンピュータ以外で構築されていてもよい。
【００４５】
逆追跡部２２０は、解析対象の場に設定された計算格子における計算格子点の各々について、可変な場の値を有して当該計算格子点へ移流する仮想粒子の軌跡を、時刻を遡って追跡する。特に、逆追跡部２２０は、仮想粒子の軌跡の座標を求めるなど、仮想粒子の軌跡を具体的に求める。
逆追跡部２２０が行う処理は、逆追跡ステップにおける処理の一例に該当する。
【００４６】
場の値取得部２３０は、場の値取得部１３０と同様、計算格子点の各々について、逆追跡の終了位置における場の値と、逆追跡の結果とに基づいて、拡散しながら逆追跡の終了位置から計算格子点へ到達した際に仮想粒子が有する場の値を、当該計算格子点における場の値として取得する。但し、場の値取得部２３０が行う具体的な処理は、場の値取得部１３０が行う処理と異なる。
場の値取得部２３０が行う処理は、場の値取得ステップにおける処理の一例に該当する。
【００４７】
場の値算出部２３１は、計算格子点の各々について、仮想粒子の軌跡に沿って拡散方程式を積分して、計算格子点の各々における場の値を算出する。場の値算出部２３１が行う処理は、場の値算出ステップにおける処理の一例に該当する。
収束判定部２３２は、場の値算出部２３１が算出した場の値が収束したか否かを判定する。場の値算出部２３１の算出した場の値が収束したと前記収束判定部２３２が判定するまで、場の値算出部２３１は場の値の算出を繰り返す。収束判定部２３２が行う処理は、収束判定ステップにおける処理の一例に該当する。
【００４８】
次に、図５を参照して、解析装置２００の動作について説明する。
なお、第１の実施形態では、仮想粒子の移流の際には拡散を考慮していない。このため、第１の実施形態では、移流時において仮想粒子のスカラー量は変化しない。これに対して、第２の実施形態では、仮想粒子の移流の際に拡散を考慮し、仮想粒子のスカラー量を変化させる。これに伴い、第２の実施形態では、計算格子点ｋのスカラー量をφｋ（ｔ）と表記し、仮想粒子のスカラー量をφ（チルダ）ｋ（ｔ）と表記して両者を区別する。
【００４９】
図５は、解析装置２００が行う処理の手順を示すフローチャートである。図５のステップＳ２０１は、図３のステップＳ１０１と同様である。
ステップＳ２０１の後、逆追跡部２２０は、仮想粒子の軌跡を逆追跡して、当該軌跡を算出する（ステップＳ２０２）。より具体的には、逆追跡部２２０は、速度場に乗って計算格子点ｋへ移流する仮想粒子の軌跡ｘｋ（ｔ）を、式（１）を時刻方向逆向きに解くことによって求める。ここで、逆追跡部２２０が、スカラー量既知の領域へ到達するまで逆追跡を行うようにしてもよいが、これに限らず、ある程度長い時間（距離）にわたって仮想粒子を逆追跡すればよい。従って、解析装置２００は、流入境界の与えられていない系（閉じた系）における場の解析や、死水領域における場の解析を行うこともできる。
【００５０】
式（１）の解法として、同式を離散的に解くことが考えられる。その際、逆追跡部２２０が行う離散化スキームとして、オイラー陽解法、完全陰解法、クランクニコルソン法、ルンゲクッタ法、ルンゲクッタフェールベルグ法またはアダムスバッシュフォース法など様々なスキームを用いることができる。
【００５１】
また、速度ｕが離散的に定義されている場合、近傍の速度定義点からの速度ｕの内挿値を用いることが考えられる。その際、逆追跡部２２０が行う内挿スキームとして、一次精度内挿、バイリニア内挿、トリリニア内挿、二次精度内挿または三次精度内挿など様々なスキームを用いることができる。
逆追跡部２２０は、各計算格子点について上記の逆追跡を行う。
【００５２】
なお、計算誤差により得られる軌跡が非現実的な軌跡となる場合、逆追跡部２２０が、それを防ぐ手段（例えば速度場の補正）を講じるようにしてもよい。例えば数値誤差の影響により、仮想粒子が、速度条件をノンスリップとする境界（壁面）に達した場合に、逆追跡部２２０が、逆追跡に用いる速度の補正を行って、仮想粒子の停止を回避するようにしてもよい。
【００５３】
次に、場の値算出部２３１は、軌跡と拡散方程式に基づいて、スカラー場を算出する（ステップＳ２０３）。以下、場の値算出部２３１が行うスカラー場の算出について、より詳細に説明する。
まず、場の値設定部１３１（図１）と同様、場の値算出部２３１も、拡散を考慮しない場合のスカラー分布を、スカラー分布算出の初期値として求める。但し、場の値算出部２３１は、拡散を考慮しない場合のスカラー分布に限らず様々なスカラー分布を、スカラー分布算出の初期値とし得る。以下に説明するように、場の値算出部２３１は、スカラー場の算出を算出結果が収束するまで繰り返すからである。一方、場の値設定部２３１が、拡散を考慮しない場合のスカラー分布を初期値とすることで、計算が収束するまでの反復回数を削減でき、処理時間および処理負荷を低減させ得る。
次に、第１の実施形態の場合と異なり、場の値算出部２３１は、式（３）に示される拡散方程式を離散的に解く。
【００５４】
【数３】
2014119882.tif 000005
【００５５】
但し、ｔは時刻を示す。また、Ｄは拡散係数を示し、▽２はラプラシアン演算子を示し、φはスカラー場を示し、ｓは生成項を示す。
具体的には、場の値算出部２３１は、ｔ＝ｔＭからｔ＝ｔ０まで式（３）を積分して、拡散しながら計算格子点ｋに到達した仮想粒子のスカラー量φ(チルダ)ｋ（ｔ０）を求め、式（４）により、計算格子点ｋのスカラー量φｋ（ｔ０）を求める。
【００５６】
【数４】
2014119882.tif 000006
【００５７】
但し、ｔ０は、仮想粒子が、逆追跡開始位置である計算格子点ｋに到達する時刻を示す。また、ｔＭは、当該仮想粒子が逆追跡終了位置を出発する時刻を示す。ここで、軌跡に含まれる様々な位置を、逆追跡終了位置とすることができる。例えば、仮想粒子が流入境界などスカラー量既知の領域を出発する位置を逆追跡終了位置としてもよい。あるいは、スカラー量既知の領域と計算格子点ｋとの間の、軌跡上の１点を、逆追跡終了位置としてもよい。
【００５８】
式（３）を積分する際、場の値算出部２３１が行う、時刻に関する離散化スキームとして、オイラー陽解法、完全陰解法、クランクニコルソン法、ルンゲクッタ法、ルンゲクッタフェールベルグ法またはアダムスバッシュフォース法など様々なスキームを用いることができる。
このとき、逆追跡終了位置ｘｋ（ｔＭ）におけるスカラー量が境界条件等によって与えられている場合は、これを仮想粒子のスカラー量の初期条件φ（チルダ）ｋ（ｔＭ）として代入する。一方、逆追跡終了位置ｘｋ（ｔＭ）におけるスカラー量が与えられていない場合は、近傍の計算格子点から内挿する。その際、場の値算出部２３１が行う内挿スキームとして、一次精度内挿、バイリニア内挿、トリリニア内挿、二次精度内挿または三次精度内挿など様々なスキームを用いることができる。
【００５９】
また、ラプラシアンに関しては、中心位置のスカラー量として軌跡上の値を用い、近傍のスカラー分布として計算格子点のスカラー量を用いる。場の値算出部２３１が行うラプラシアンの離散化スキームとして、各種ラプラシアンモデル、二次精度中央差分法、四次精度中央差分法など様々なスキームを用いることができる。例えば、二次精度中央差分法を適用した場合、式（５）のようになる。
【００６０】
【数５】
2014119882.tif 000007
【００６１】
但し、ｅｘ、ｅｙ、ｅｚは、それぞれ、ｘ方向、ｙ方向、ｚ方向の単位ベクトル（基底ベクトル）を示す。
また、δｘ、δｙ、δｚは、それぞれ、差分法における中心位置と近傍との幅を示す。δｘ、δｙ、δｚのいずれも、計算格子点間距離と同じかそれよりも大きい値に設定する。ラプラシアンの過大評価を防止するためである。
【００６２】
また、φ（バー）（ｔ，ｘ）は、位置ｘへの計算格子点からのスカラー量の内挿を示す。その際、場の値算出部２３１が行う内挿スキームとして、一次精度内挿、バイリニア内挿、トリリニア内挿、二次精度内挿または三次精度内挿など様々なスキームを用いることができる。但し、解析精度の観点から空間精度二次以上の内挿スキームを用いることが好ましい。
【００６３】
図６は、場の値算出部２３１が行うスカラー量算出の例を示す説明図である。同図において、点Ｐ２１は、仮想粒子が到達する計算格子点を示し、線Ｌ２１は、仮想粒子の軌跡を示している。
場の値算出部２３１は、軌跡（線Ｌ２１）上の点Ｐ２２を差分法における中心位置とし、点Ｐ２３などの近傍の点のスカラー量を、図に示すように周囲の計算格子点からの内挿にて求める。そして、場の値算出部２３１は、得られた近傍の点のスカラー量を式（５）に代入する。
このように、場の値算出部２３１は、軌跡上の点と周囲とのスカラー量の差に基づいて、拡散による仮想粒子のスカラー量の変化を求める。これを軌跡に沿って積分することで、逆追跡の終了位置におけるスカラー量からの変位にて、仮想粒子が到達した計算格子点（点Ｐ２１）のスカラー量を求めることができる。
場の値算出部２３１は、各計算格子点について式（３）を解いてスカラー場を求める。
【００６４】
ステップＳ２０３の後、収束判定部２３２は、場の値算出部２３１の求めたスカラー場が収束したか否かを判定する（ステップＳ２０４）。具体的には、収束判定部２３２は、場の値算出部２３１が新たに求めたスカラー場と、場の値算出部２３１が前回求めたスカラー場との差の大きさを算出し、差の大きさが所定の閾値以下の場合に収束したと判定し、閾値より大きい場合には収束していないと判定する。
なお、場の値算出部２３１が新たに求めたスカラー場と、場の値算出部２３１が前回求めたスカラー場との差の大きさとして、収束判定部２３２は、スカラー場の変化の大きさを示す様々な値を用いることができる。例えば、収束判定部２３２が、計算格子点毎のスカラー量の変化の絶対値を全計算格子点について足し合わせた値を用いるようにしてもよい。あるいは、収束判定部２３２が、計算格子点毎のスカラー量の変化の二乗平均平方根を用いるようにしてもよい。
なお、場の値算出部２３１が最初にスカラー場を求めた場合、収束判定部２３２は、場の値算出部２３１の求めたスカラー場が収束していないと判定する。
【００６５】
ステップＳ２０４において、場の値算出部２３１の求めたスカラー場が収束していないと判定した場合（ステップＳ２０４：ＮＯ）、ステップＳ２０３へ戻る。
この場合、ステップＳ２０３において、場の値算出部２３１は、ステップＳ２０３における前回の処理で得られたスカラー場を用いて各計算格子点について式（３）を解き直す。これにより、場の値算出部２３１は、スカラー場（各計算格子点のスカラー量）を更新する。以下、場の値算出部２３１の求めたスカラー場が収束したと収束判定部２３２が判定するまで、すなわち、スカラー場の変化が充分小さくなるまで、場の値算出部２３１は、得られたスカラー場を用いて各計算格子点について式（３）を解き直す処理を繰り返す。
【００６６】
一方、ステップＳ２０４において、場の値算出部２３１の求めたスカラー場が収束したと判定した場合（ステップＳ２０４：ＹＥＳ）、ステップＳ２０５へ進む。
ステップＳ２０５は、図３のステップＳ１０４と同様である。その後、図５の処理を終了する。
【００６７】
以上のように、逆追跡部２２０は、計算格子点へ移流する仮想粒子の軌跡を逆追跡し、場の値取得部２３０は、拡散しながら逆追跡の終了位置から計算格子点へ到達した際に仮想粒子が有する場の値を、当該計算格子点における場の値として取得する。
これにより、解析装置２００では、仮想粒子が有する場の値を計算格子点に内挿する必要が無いので数値拡散が生じにくく、より高精度に数値解析を行うことができる。
【００６８】
また、場の値算出部２３１は、計算格子点の各々について、仮想粒子の軌跡に沿って拡散方程式を積分して、計算格子点の各々における場の値を算出する。
場の値算出部２３１は、仮想粒子の軌跡における拡散を、軌跡の周囲の場の値に応じて算出する点で、より正確に拡散項を評価し得る。
また、収束判定部２３２は、場の値算出部２３１の算出した場の値が収束したか否かを判定し、収束したと判定するまで、場の値算出部２３１は、場の値の算出を繰り返す。
これにより、場の値算出部２３１は、拡散の影響を充分に反映させた場の値に基づいて拡散項を評価することができ、この点において、より正確に拡散項を評価し得る。
【００６９】
なお、解析装置１００の場合と同様、解析装置２００は、定常的なスカラー場と非定常なスカラー場とのいずれも解析することができる。非定常なスカラー場を解析する場合、解析装置２００は、式（４）を参照して説明したように、時刻ｔ０におけるスカラー量の分布を求める。
【００７０】
＜第３の実施形態＞
図７は、本発明の第３の実施形態における解析装置の機能構成を示す概略ブロック図である。同図において、解析装置３００は、前提条件取得部１１０と、逆追跡部２２０と、場の値取得部３３０と、解析結果出力部１４０とを具備する。場の値取得部３３０は、方程式取得部３３１と、解析実行部３３２とを具備する。なお、図７において、図４の各部に対応して同様の機能を有する部分には同一の符号（１１０、２２０、１４０）を付して説明を省略する。
【００７１】
解析装置３００は、解析装置１００や２００と同様、数値解析を行って移流拡散問題における場の値を求める。但し、解析装置３００が行う具体的な処理は、解析装置１００や２００が行う処理と異なる。解析装置３００は、例えばコンピュータにて構築される。あるいは、解析装置３００が専用の回路にて構築されるなど、コンピュータ以外で構築されていてもよい。
【００７２】
場の値取得部３３０は、場の値取得部１３０や２３０と同様、計算格子点の各々について、逆追跡の終了位置における場の値と、逆追跡の結果とに基づいて、拡散しながら逆追跡の終了位置から計算格子点へ到達した際に仮想粒子が有する場の値を、当該計算格子点における場の値として取得する。より具体的には、場の値取得部３３０は、計算格子点の各々についての仮想粒子の軌跡と拡散方程式とに基づいて得られる連立方程式から、計算格子点の各々における場の値を取得する。
場の値取得部３３０が行う処理は、場の値取得ステップにおける処理の一例に該当する。
【００７３】
方程式取得部３３１は、計算格子点の各々についての仮想粒子の軌跡と拡散方程式とに基づいて連立方程式を取得する。
解析実行部３３２は、方程式取得部３３１が取得した連立方程式から、計算格子点の各々における場の値を取得する。
【００７４】
次に、図８を参照して、解析装置３００の動作について説明する。
以下、第２の実施形態の場合と同様、φｋ（ｔ）は、計算格子点ｋのスカラー量を示し、φ（チルダ）ｋ（ｔ）は、仮想粒子のスカラー量を示す。その他、定義を示さずに用いる変数も、第１の実施形態または第２の実施形態の場合と同様である。
【００７５】
以下では、まず、解析装置３００が非定常な場の解析を行う場合について説明し、次に、解析装置３００が定常的な場の解析を行う場合について説明する。
第３の実施形態においても、第２の実施形態の場合と同様、解析対象となる空間に速度場およびスカラー場が形成されており、スカラー場を構成するスカラー量が速度場に乗って移流しながら拡散する場合にスカラー量の分布を求める問題を考える。但し、ここでは、スカラー場が非定常である場合に、新たな時刻におけるスカラー量の分布を求める。
【００７６】
第３の実施形態においても、第２の実施形態の場合と同様、解析対象となる空間に総計算格子点数Ｎの計算格子を設定し、計算格子点の各々におけるスカラー量を求める。なお、計算格子点を番号ｋ（０≦ｋ≦Ｎ−１）にて識別する。
また、第３の実施形態においても、第２の実施形態の場合と同様、軌跡上における仮想粒子の拡散を考慮してスカラー場を求める。
一方、第３の実施形態では、各計算格子点について式（３）を解いてスカラー場を更新する処理ループを不要として処理時間を削減する。
【００７７】
図８は、解析装置３００が行う処理の手順を示すフローチャートである。図８のステップＳ３０１およびＳ３０２は、図５のステップＳ２０１およびＳ２０２と同様である。第３の実施形態においても、第２の実施形態の場合と同様、逆追跡部２２０は、速度場に乗って計算格子点ｋへ移流する仮想粒子の軌跡ｘｋ（ｔ）を、式（１）を時刻方向逆向きに解くことによって求める。第３の実施形態においても、逆追跡部２２０は、ある程度長い時間（距離）にわたって仮想粒子を逆追跡すればよい。従って、解析装置３００は、流入境界の与えられていない系（閉じた系）における場の解析や、死水領域における場の解析を行うこともできる。
【００７８】
式（１）の解法として、同式を離散的に解くことが考えられる。その際、逆追跡部２２０が行う離散化スキームとして、オイラー陽解法、完全陰解法、クランクニコルソン法、ルンゲクッタ法、ルンゲクッタフェールベルグ法またはアダムスバッシュフォース法など様々なスキームを用いることができる。
【００７９】
また、速度ｕが離散的に定義されている場合、近傍の速度定義点からの速度ｕの内挿値を用いることが考えられる。その際、逆追跡部２２０が行う内挿スキームとして、一次精度内挿、バイリニア内挿、トリリニア内挿、二次精度内挿または三次精度内挿など様々なスキームを用いることができる。
逆追跡部２２０は、各計算格子点について上記の逆追跡を行う。
【００８０】
なお、計算誤差により得られる軌跡が非現実的な軌跡となる場合、逆追跡部２２０が、それを防ぐ手段（例えば速度場の補正）を講じるようにしてもよい。例えば数値誤差の影響により、仮想粒子が、速度条件をノンスリップとする境界（壁面）に達した場合に、逆追跡部２２０が、逆追跡に用いる速度の補正を行って、仮想粒子の停止を回避するようにしてもよい。
【００８１】
ステップＳ３０２の後、方程式取得部３３１は、各計算格子点のスカラー量に関する連立方程式を示す係数行列等を算出することで、当該連立方程式を取得する（ステップＳ３０３）。そして、解析実行部３３２は、方程式取得部３３１が取得した連立方程式を数値的に解くことで、解析を実行してスカラー場を求める（ステップＳ３０４）。以下、方程式取得部３３１および解析実行部３３２が行う処理について、より詳細に説明する。
【００８２】
第３の実施形態においても、第２の実施形態の場合と同様、式（３）に示される拡散方程式を離散的に解く。但し、第２の実施形態では、仮想粒子の軌跡に沿って拡散方程式を積分したのに対し、第３の実施形態では、拡散方程式から、各計算格子点のスカラー量に関する連立方程式を示す係数行列等を取得し、当該連立方程式を解く。まず、方程式取得部３３１が取得する連立方程式について説明する。
【００８３】
第２の実施形態で説明したように、式（３）を積分する際、時刻に関する離散化スキームとして、オイラー陽解法、完全陰解法、クランクニコルソン法、ルンゲクッタ法、ルンゲクッタフェールベルグ法またはアダムスバッシュフォース法など様々なスキームを用いることができる。例えば、完全陰解法を適用した場合、式（６）のようになる。
【００８４】
【数６】
2014119882.tif 000008
【００８５】
但し、δｔは時間に関する離散化における時間幅を示す。また、φはスカラー場を示す。
また、第２の実施形態の場合と同様、ラプラシアンに関しては、中心位置のスカラー量として軌跡上の値を用い、近傍のスカラー分布として計算格子点のスカラー量を用いる。ラプラシアンの離散化スキームとして、各種ラプラシアンモデル、二次精度中央差分法、四次精度中央差分法など様々なスキームを用いることができる。例えば、二次精度中央差分法を適用した場合、式（５）のようになる。
【００８６】
ここで、任意の内挿値φ（バー）（ｔ，ｘ）が、式（７）に示されるように計算格子点上のスカラー量φの線形結合で与えられるとする。
【００８７】
【数７】
2014119882.tif 000009
【００８８】
但し、ｗｎ（ｘ）は、位置ｘに応じて計算格子点ｎ毎に設定される重み係数を示す。
すると、式（５）は、式（８）に示すように離散化される。
【００８９】
【数８】
2014119882.tif 000010
【００９０】
但し、Ｗｎｋ，ｔは、式（９）に示される。
【００９１】
【数９】
2014119882.tif 000011
【００９２】
式（８）を式（６）に代入して整理すると、式（１０）のようになる。
【００９３】
【数１０】
2014119882.tif 000012
【００９４】
但し、ａｋ，ｔは、式（１１）のように示される。
【００９５】
【数１１】
2014119882.tif 000013
【００９６】
また、Ｗ’ｎｋ，ｔは、式（１２）のように示される。
【００９７】
【数１２】
2014119882.tif 000014
【００９８】
また、δｓｋ，ｔは、式（１３）のように示される。
【００９９】
【数１３】
2014119882.tif 000015
【０１００】
仮想粒子が計算格子点ｋに到達する時刻をｔ０とすると、式（１４）が成立する。
【０１０１】
【数１４】
2014119882.tif 000016
【０１０２】
また、式（１０）に、ｔ→ｔ０、および、δｔ→ｔ０−ｔ１（ｔ０−ｔ１＞０）を代入すると、式（１４）より式（１５）のようになる。
【０１０３】
【数１５】
2014119882.tif 000017
【０１０４】
更に、時刻ｔ２（ｔ２＜ｔ１）を導入し、φ（チルダ）ｋ（ｔ１）に対して、同様の展開を施すと、式（１６）のようになる。
【０１０５】
【数１６】
2014119882.tif 000018
【０１０６】
同様の展開を、仮想粒子が逆追跡終了位置を出発する時刻ｔＭまで繰り返して、式（１７）を得る。
【０１０７】
【数１７】
2014119882.tif 000019
【０１０８】
φ（チルダ）ｋ（ｔＭ）の値が未知の場合、式（１７）を式（１８）のように近似することが考えられる。
【０１０９】
【数１８】
2014119882.tif 000020
【０１１０】
式（１９）の値が十分に小さいと仮定すれば、式（１８）の近似は妥当である。ｔ０−ｔＭの値が大きいほど、すなわち逆追跡の時間が長いほど（従って、逆追跡の距離が長いほど）、式（１９）の値が小さくなる。
【０１１１】
【数１９】
2014119882.tif 000021
【０１１２】
一方、φ（チルダ）ｋ（ｔＭ）の値が境界条件等により既知の場合は、単純に式（１７）に値を代入すればよい。
式（１７）または式（１８）を全ての計算格子点ｋ（０≦ｋ≦Ｎ−１）について考えると、計算格子点のスカラー量φｋ（ｔ）に関するＮ元連立一次方程式を得られる。
さらに、計算格子点のスカラー量φｋ（ｔ０）を成分に持つベクトルΦを式（２０）のように定義する。
【０１１３】
【数２０】
2014119882.tif 000022
【０１１４】
但し、Ｔはベクトルないし行列の転置を示す。
式（１７）または式（１８）に基づいて得られたＮ元連立一次方程式は、Φを用いて式（２１）のように表すことができる。
【０１１５】
【数２１】
2014119882.tif 000023
【０１１６】
但し、Ａは、式（２２）に示される成分を有するＮ次正方行列である。
【０１１７】
【数２２】
2014119882.tif 000024
【０１１８】
また、ｂはＮ次元ベクトルを示す。
φ（チルダ）ｋ（ｔＭ）の値が既知の場合、従って式（１７）に基づく場合、ｂの成分は式（２３）のようになる。
【０１１９】
【数２３】
2014119882.tif 000025
【０１２０】
一方、φ（チルダ）ｋ（ｔＭ）の値が未知の場合、従って式（１８）に基づく場合、ｂの成分は式（２４）のようになる。
【０１２１】
【数２４】
2014119882.tif 000026
【０１２２】
方程式取得部３３１は、式（２１）における係数行列Ａおよびベクトルｂを求めることで、式（２１）に示される連立方程式を取得する。具体的には、方程式取得部３３１は、式（２２）〜式（２４）に基づいて、係数行列Ａやベクトルｂの各要素を求める。
【０１２３】
Φ（ｔ）：［ｔＭ，ｔ１］が既知であれば、ＳＯＲ法（Successive Over-relaxation Method）またはＢｉＣＧＳｔａｂ法（Biconjugate Gradient Stabilized Method）など公知の解法を用いて、新しい（未だ算出していない）時刻ｔ０における全ての計算格子点のスカラー量Φが求まる。但し、Φ（ｔ）は、時刻ｔにおけるスカラー場を示す。
そこで、解析実行部３３２は、方程式取得部３３１が求めた連立方程式と、取得済みのスカラー量Φ（ｔ）とに基づいて、新しい時刻ｔ０におけるスカラー場（全ての計算格子点のスカラー量）を求める。
【０１２４】
ステップＳ３０４の後、ステップＳ３０５へ進む。ステップＳ３０５は、図５のステップＳ２０５と同様である。
その後、図８の処理を終了する。
【０１２５】
次に、スカラー場が定常の場合について説明する。
第２の実施形態の場合と同様、解析対象となる空間に速度場およびスカラー場が形成されており、スカラー場を構成するスカラー量が速度場に乗って移流しながら拡散する場合に、定常状態となったスカラー量の分布を求める問題を考える。
なお、以下では、計算格子点ｋのスカラー量をφｋと表記する。
【０１２６】
スカラー場が非定常の場合と同様、解析対象となる空間に総計算格子点数Ｎの計算格子を設定し、計算格子点の各々におけるスカラー量を求める。なお、計算格子点を番号ｋ（０≦ｋ≦Ｎ−１）にて識別する。
また、軌跡上における仮想粒子の拡散を考慮してスカラー場を求める点、および、各計算格子点について式（３）を解いてスカラー場を更新する処理ループを不要として処理時間を削減する点も、スカラー場が非定常の場合と同様である。
【０１２７】
また、スカラー場が非定常の場合と同様、逆追跡部２２０は、速度場に乗って計算格子点ｋへ移流する仮想粒子の軌跡ｘｋ（ｔ）を、式（１）を時刻方向逆向きに解くことによって求める。スカラー場が非定常の場合と同様、逆追跡部２２０は、ある程度長い時間（距離）にわたって仮想粒子を逆追跡すればよい。従って、解析装置３００は、流入境界の与えられていない系（閉じた系）における場の解析や、死水領域における場の解析を行うこともできる。
また、式（１）の解法や、様々なスキームを用いることができる点や、計算誤差により得られる軌跡が非現実的な軌跡となることを防ぐ手段等も、スカラー場が非定常の場合と同様である。
【０１２８】
また、式（３）に示される拡散方程式を離散的に解く点および解法や、様々なスキームを用いることができる点も、スカラー場が非定常の場合と同様である。式（３）を積分する際、時刻に関する離散化スキームとして完全陰解法を適用した場合、式（６）のようになる点も、スカラー場が非定常の場合と同様である。
【０１２９】
また、非定常の場合と同様、ラプラシアンに関しては、中心位置のスカラー量として軌跡上の値を用い、近傍のスカラー分布として計算格子点のスカラー量を用いる。ラプラシアンの離散化スキームとして、各種ラプラシアンモデル、二次精度中央差分法、四次精度中央差分法など様々なスキームを用いることができる。例えば、二次精度中央差分法を適用した場合、式（２５）のようになる。
【０１３０】
【数２５】
2014119882.tif 000027
【０１３１】
式（２５）では、定常状態にあるスカラー場を扱うため、φ（バー）の直接の引数に時刻ｔが含まれていない。それ以外は、式（５）と同様である。
ここで、任意の内挿値φ（バー）（ｘ）が、式（２６）に示されるように計算格子点上のスカラー量φの線形結合で与えられるとする。
【０１３２】
【数２６】
2014119882.tif 000028
【０１３３】
すると、式（２５）は、式（２７）に示すように離散化される。
【０１３４】
【数２７】
2014119882.tif 000029
【０１３５】
但し、非定常の場合と同様、Ｗｎｋ，ｔは、式（９）に示される。
式（２７）を式（６）に代入して整理すると、式（２８）のようになる。
【０１３６】
【数２８】
2014119882.tif 000030
【０１３７】
但し、非定常の場合と同様、ａｋ，ｔ、Ｗ’ｎｋ，ｔ、δｓｋ，ｔは、それぞれ式（１１）、（１２）、（１３）のように示される。
仮想粒子が計算格子点ｋに到達する時刻をｔ０とすると、式（２９）が成立する。
【０１３８】
【数２９】
2014119882.tif 000031
【０１３９】
式（２８）に、ｔ→ｔ０、および、δｔ→ｔ０−ｔ１（ｔ０−ｔ１＞０）を代入すると、式（３０）のようになる。
【０１４０】
【数３０】
2014119882.tif 000032
【０１４１】
更に、時刻ｔ２（ｔ２＜ｔ１）を導入し、φ（チルダ）ｋ（ｔ１）に対して、同様の展開を施すと、式（３１）のようになる。
【０１４２】
【数３１】
2014119882.tif 000033
【０１４３】
同様の展開を時刻ｔＭまで繰り返して、式（３２）を得る。
【０１４４】
【数３２】
2014119882.tif 000034
【０１４５】
φ（チルダ）ｋ（ｔＭ）の値が未知の場合、式（３２）を式（３３）のように近似することが考えられる。
【０１４６】
【数３３】
2014119882.tif 000035
【０１４７】
非定常の場合と同様、式（１９）の値が十分に小さいと仮定すれば、式（３３）の近似は妥当である。
一方、φ（チルダ）ｋ（ｔＭ）の値が境界条件等により既知の場合は、単純に式（３２）に値を代入すればよい。
式（３２）または式（３３）を全ての計算格子点ｋ（ｋ∈［０，Ｎ−１］）について考えると、計算格子点のスカラー量φｋに関するＮ元連立一次方程式を得られる。
さらに、計算格子点のスカラー量φｋを成分に持つベクトルΦを式（３４）のように定義する。
【０１４８】
【数３４】
2014119882.tif 000036
【０１４９】
但し、Ｔはベクトルないし行列の転置を示す。
式（３２）または式（３３）に基づいて得られたＮ元連立一次方程式は、Φを用いて式（３５）のように表すことができる。
【０１５０】
【数３５】
2014119882.tif 000037
【０１５１】
但し、ＡはＮ次正方行列の係数行列を示し、ｂはＮ次元ベクトルを示す。
φ（チルダ）ｋ（ｔＭ）の値が既知の場合、従って式（３２）に基づく場合、Ａの成分は式（３６）のようになる。
【０１５２】
【数３６】
2014119882.tif 000038
【０１５３】
また、式（３２）に基づく場合、ｂの成分は式（３７）のようになる。
【０１５４】
【数３７】
2014119882.tif 000039
【０１５５】
一方、φ（チルダ）ｋ（ｔＭ）の値が未知の場合、従って式（３３）に基づく場合、Ａの成分は式（３８）のようになる。
【０１５６】
【数３８】
2014119882.tif 000040
【０１５７】
また、式（３３）に基づく場合、ｂの成分は式（３９）のようになる。
【０１５８】
【数３９】
2014119882.tif 000041
【０１５９】
方程式取得部３３１は、式（３５）における係数行列Ａおよびベクトルｂを求めることで、式（３５）に示される連立方程式を取得する。具体的には、方程式取得部３３１は、式（３６）、式（２３）、式（３８）および式（２４）に基づいて、係数行列Ａやベクトルｂの各要素を求める。
【０１６０】
ＳＯＲ法またはＢｉＣＧＳｔａｂ法など公知の解法を式（３５）に用いて、全ての計算格子点のスカラー量Φが求まる。
そこで、解析実行部３３２は、方程式取得部３３１が求めた連立方程式を解いて、スカラー場（全ての計算格子点のスカラー量）を求める。
【０１６１】
以上のように、逆追跡部２２０は、計算格子点へ移流する仮想粒子の軌跡を逆追跡し、場の値取得部３３０は、拡散しながら逆追跡の終了位置から計算格子点へ到達した際に仮想粒子が有する場の値を、当該計算格子点における場の値として取得する。
これにより、解析装置３００では、仮想粒子が有する場の値を計算格子点に内挿する必要が無いので数値拡散が生じにくく、より高精度に数値解析を行うことができる。
【０１６２】
また、場の値取得部３３０は、計算格子点の各々についての仮想粒子の軌跡と拡散方程式とに基づいて得られる連立方程式から、計算格子点の各々における場の値を取得する。
これにより、場の値取得部３３０は、場の値の更新を繰り返す必要が無い。この点において、処理時間や、処理に必要な記憶容量を削減することができる。
【０１６３】
なお、本発明者は、上記の実施形態と同様の数値解析にて、実際に場の値を求めた。その結果本発明の有効性を確認することができた。以下、この点について説明する。
【０１６４】
＜第１の解析例＞
第１の解析例では、図９に示すＴ型ミキサーにおいて濃度場が定常状態にある場合の二液の混合を、第１の実施形態と同様の数値解析にて解析した。
図９は、流体混合に用いられる典型的なＴ型ミキサーの一例を示す。同図において、各寸法は流路高さｈに対する比によって無次元化されている。
２つの流入口Ａ３１またはＡ３２から入った流体は、流入流路を通り混合流路で合流する。合流した２つの流体は混合流路内で混合されて排出口Ａ３３から排出される。
以下では、液体を混合する場合を例に説明する。流入する二液について、ニュートン流体、非圧縮性、および、二液の物資同一性を仮定する。
この場合、無次元形式の支配方程式は、連続の式、ナビエ・ストークス方程式、濃度に関する移流拡散方程式の三式であり、連続の式は、式（４０）のように示される。
【０１６５】
【数４０】
2014119882.tif 000042
【０１６６】
但し、ｕは速度を示す。また、▽は、▽２にてラプラシアン演算子を示す。
また、ナビエ・ストークス方程式は、式（４１）のように示される。
【０１６７】
【数４１】
2014119882.tif 000043
【０１６８】
但し、ｔは時刻を示し、ｐは圧力を示す。また、Ｒｅはレイノルズ数を示す。レイノルズ数Ｒｅは、式（４２）によって定義される無次元数である。
【０１６９】
【数４２】
2014119882.tif 000044
【０１７０】
但し、Ｕは入口平均流速を示し、νは動粘性係数を示し、Ｄは拡散係数を示す。
また、濃度に関する移流拡散方程式は、式（４３）のように示される。
【０１７１】
【数４３】
2014119882.tif 000045
【０１７２】
但し、ｃは濃度を示し、Ｓｃはシュミット数を示す。シュミット数Ｓｃは、式（４４）によって定義される無次元数である。
【０１７３】
【数４４】
2014119882.tif 000046
【０１７４】
本解析では、二液を物性値の等しい流体と仮定しており、速度場の支配方程式（連続の式およびナビエ・ストークス方程式）は濃度を含んでいない。このため、濃度場と分離して速度場を解くことができる。
そこで、濃度場の解析を行う前に、速度場の解析を行った。速度場の解析には、公知の解析手法である格子法を用い、総計算格子点数を約２４０万とした。また、レイノルズ数Ｒｅ＝１５０とした。
速度場の解析において、速度場に対する大きな格子依存性は確認されなかった。また、流れは定常状態となっていた。
【０１７５】
濃度場の解析では、速度場の解析に用いた計算格子とは異なる計算格子を設定した。濃度場に対する計算格子は、図９に示すｙ座標にてｙ＝７．６７の位置における、混合流路の断面に設定した。
各計算格子点に到達する仮想粒子の軌跡ｘｋ（ｔ）は、仮想粒子の移流に関する関係式である式（１）を、時刻を遡る方向にオイラー陽解法的に離散化した式（４５）を数値的に解くことによって求める。
【０１７６】
【数４５】
2014119882.tif 000047
【０１７７】
式（４５）を解く処理は、式（１）を解く処理の一例に該当する。
図１０は、濃度場の解析における仮想粒子の逆追跡の例を示す説明図である。同図において、２つの流入流路に、それぞれ初期平面Ａ４１、Ａ４２が設定されている。初期平面Ａ４１、Ａ４２は、スカラー値既知の領域の一例に該当し、濃度ｃの値は、それぞれ０、１となっている。また、混合流路には、計算格子の設定されているターゲット平面Ａ４３が示されている。
軌跡の算出では、ターゲット平面Ａ４３上の計算格子点の１つである逆追跡開始位置Ｐ４２から、矢印Ｖ４１で示すように時刻を遡る方向に、初期平面Ａ４１またはＡ４２に到達するまで軌跡Ｌ４１を求める。
【０１７８】
なお、軌跡を求めるにあたって、速度ｕの内挿スキームとして速度の定義点（速度場の解析にて速度を求めた位置）からのトリリニア内挿を適用した。また、クーラン数を０．１に設定し、この条件を満たすように時間幅δｔを決定した。
得られる軌跡は数値誤差の影響を含むため、速度条件をノンスリップとする境界（壁面）に達してしまう場合がある。この場合の対処として、壁面への到達を検知した前の時間ステップにおいて、壁面に対して垂直方向の速度がゼロとなるように、逆追跡に用いる速度の補正を行った。
【０１７９】
逆追跡において、仮想粒子はやがて２つの流入口のいずれかに向かう。その際、仮想粒子の位置の座標値のうち、図１０に示すｘ座標の値に基づいて、流入時における初期濃度値を判定するようにした。具体的には、ｘ≦−２となった場合、濃度ｃ＝０（初期平面Ａ４１側）と判定する。一方、ｘ≧２となった場合、濃度ｃ＝１（初期平面Ａ４２側）と判定する。
【０１８０】
全ての計算格子点について逆追跡を行うことで、各計算格子点のスカラー値が０または１に設定される。このようにして得られたスカラー分布を初期（ξ＝０）における値として、拡散方程式をモデル化した式（４６）に基づいて、拡散効果を反映したスカラー分布を得た。
【０１８１】
【数４６】
2014119882.tif 000048
【０１８２】
但し、τは仮想粒子の滞留時間を示す。また、ξは仮に設定した変数である。
具体的には、式（４６）を完全陰解法で離散化した式（４７）に対してＳＯＲ法を適用し、ξ＝１におけるスカラー分布を求めた。
【０１８３】
【数４７】
2014119882.tif 000049
【０１８４】
なお、シュミット数Ｓｃ＝３６００とした。なお式（４６）は、式（２）の一例に該当する。
図１１は、ターゲット平面Ａ４３（図１０）における混合度ＭＩの計算結果の例を、計算格子点間間隔を変えながら格子法と比較した結果を示す説明図である。図１１において、丸で示す各点は、本実施形態の解析方法（以下、本実施形態において「本手法」と称する）における計算結果を示す。一方、三角で示す各点は、格子法における計算結果を示す。
混合度ＭＩは、式（４８）によって定義される。
【０１８５】
【数４８】
2014119882.tif 000050
【０１８６】
但し、ｗは面Ａ（式（４８）における積分の対象となっている面）に対する垂直方向速度を示す。また、ｃ（バー）は、完全に混合したときの濃度（すなわち、ｃ（バー）＝０．５）を示す。
格子法による解析には市販のＣＦＤ（Computational Fluid Dynamics）コードを使用し、本手法による解析には本発明者がＣ＋＋にて生成したコードを使用した。また、格子法における最大計算格子点数は、約６６００万であり、その際の空間解像度は、一方向の単位長さあたりの点数で１２０であった。一方、本手法における最大計算格子点数は約５２万であり、その際の空間解像度は、一方向の単位長さあたりの点数で５１２であった。
【０１８７】
図１２は、格子法と本手法との、最大空間解像度における濃度分布の解析結果の例を示す説明図である。上側の図Ｆ１１は、格子法による濃度分布を示し、下側の図Ｆ１２は、本手法による濃度分布を示す。
図Ｆ１１とＦ１２とを比較すると、図Ｆ１１のほうが混合の度合いが大きい。これは、格子法において数値拡散により拡散の過大評価を引き起こしたためと考えられる。
【０１８８】
格子法に使用した計算格子では、一般的な計算機の限界ともいえるほど大規模な最大計算格子点数となっている。それにもかかわらず格子法では、数値拡散により拡散の過大評価を引き起こしており未だ計算格子依存性が強く残っている。その結果、格子法では正しい解析結果を得られていない。
これに対し、本手法では、計算格子点に到達する仮想粒子のスカラー場を求めるので、仮想粒子のスカラー量を計算格子点のスカラー量に内挿する必要がない。また、逆追跡終了位置においてスカラー量既知なので、計算格子点のスカラー量を仮想粒子のスカラー量に内挿する必要がない。これらの点において、本手法は、根本的に数値拡散が生じにくいという特徴を有する。これにより、本手法は、圧倒的に小さい計算格子依存性を示している。
【０１８９】
なお、最大空間解像度における所要解析時間は、格子法で約１０時間であったのに対し、本手法では３０分程度であり、約２０分の１の時間で解析を行うことができた。また、使用した記憶容量は、格子法で約４５ギガバイト（ＧＢ）であったのに対し、本手法では１ギガバイト未満であり、４５分の１程度の記憶容量で解析を行うことができた。
【０１９０】
＜第２の解析例＞
第２の解析例では、図９に示すＴ型ミキサーにおいて濃度場が定常状態にある場合の二液の混合を、第３の実施形態と同様の数値解析にて解析した。
第１の解析例と同様、液体を混合する場合を例に説明する。流入する二液について、ニュートン流体、非圧縮性、および、二液の物資同一性を仮定する。
また、第１の解析例と同様、支配方程式は、式（４０）に示す連続式と、式（４１）に示すナビエ・ストークス方程式と、式（４３）に示す濃度に関する移流拡散方程式との三式である。
【０１９１】
また、第１の解析例と同様、濃度場の解析を行う前に、速度場の解析を行った。第１の解析例と同様、速度場の解析には、公知の解析手法である格子法を用い、総計算格子点数を約２４０万とし、レイノルズ数Ｒｅ＝１５０とした。
【０１９２】
一方、第１の解析例と異なり、本解析例では、濃度場に対する計算格子を、速度場に対する計算格子と同じ領域に設定した。但し、速度場に対する計算格子と比較すると空間解像度を変更し、また、着目位置から見て遠方に存在する領域の計算省略を適宜行った。
【０１９３】
濃度場に対する各計算格子点に到達する仮想粒子の軌跡ｘｋ（ｔ）は、第１の解析例と同様、仮想粒子の移流に関する関係式である式（１）を、時刻を遡る方向にオイラー陽解法的に離散化した式（４５）を数値的に解くことによって求める。
なお、軌跡を求めるにあたって、速度ｕの内挿スキームとして速度の定義点からのトリリニア内挿を適用した。また、クーラン数を０．１、拡散数を０．０７５に設定し、この条件を満たすように時間幅δｔを決定した。
【０１９４】
第１の解析例と同様、速度条件をノンスリップとする境界（壁面）に仮想粒子が達した場合の対処として、壁面への到達を検知した前の時間ステップにおいて、壁面に対して垂直方向の速度がゼロとなるように、逆追跡に用いる速度の補正を行った。
また、逆追跡の終了条件は、仮想粒子が濃度既知の位置へ到達するか、あるいは、式（４９）を満たすこととした。
【０１９５】
【数４９】
2014119882.tif 000051
【０１９６】
ここで、許容値ε（ε≧０）を小さく設定するほど解析精度が高くなるが、解析にかかるコスト（解析時間や計算負荷など）は増大する。本解析では、経験的にε＝０．１とした。
上記の設定のもとで解析を行い、仮想粒子の軌跡ｘｋ（ｔ）：［ｔＭ，ｔ０］を求めた。但し、Ｎを総計算格子点数として、ｋは、０≦ｋ≦Ｎ−１の正整数である。また、ｔＭ≦ｔ０である。
【０１９７】
濃度場の解析では、式（３５）に示すＮ元連立一次方程式を数値的に解いて濃度場を求めた。その際、式（５０）に示すように、式（３４）におけるスカラー量φ０，φ１，・・・，φＮ−１の具体例として濃度ｃ０，ｃ１，・・・，ｃＮ−１を用いる。
【０１９８】
【数５０】
2014119882.tif 000052
【０１９９】
また、位置ｘｋ（ｔＭ）における仮想粒子の濃度ｃ（チルダ）ｋ（ｔＭ）が未知の場合、係数行列Ａを、式（３８）に示すようにする。
但し、ａｋ，ｔｍは、式（５１）のように示される。
【０２００】
【数５１】
2014119882.tif 000053
【０２０１】
また、Ｗ’ｎｋ，ｔｍは、式（５２）のように示される。
【０２０２】
【数５２】
2014119882.tif 000054
【０２０３】
また、Ｗｎｋ，ｔｍは、式（５３）のように示される。
【０２０４】
【数５３】
2014119882.tif 000055
【０２０５】
また、ｗｎ（ｘ）は、計算格子点上のスカラー分布から位置ｘへの内挿値ｃ（バー）（ｘ）を与える重み係数である。内挿値ｃ（バー）（ｘ）は、式（５４）のように示される。
【０２０６】
【数５４】
2014119882.tif 000056
【０２０７】
また、位置ｘｋ（ｔＭ）における仮想粒子の濃度ｃ（チルダ）ｋ（ｔＭ）が未知の場合、ベクトルｂを式（５５）に示すようにする。
【０２０８】
【数５５】
2014119882.tif 000057
【０２０９】
一方、位置ｘｋ（ｔＭ）における仮想粒子の濃度ｃ（チルダ）ｋ（ｔＭ）が既知の場合、係数行列Ａを、式（３６）に示すようにする。
但し、位置ｘｋ（ｔＭ）における仮想粒子の濃度ｃ（チルダ）ｋ（ｔＭ）が未知の場合と同様、ａｋ，ｔｍ、Ｗ’ｎｋ，ｔｍ、Ｗｎｋ，ｔｍは、それぞれ、式（５０）、式（５１）、式（５２）に示すようにする。また、ｗｎ（ｘ）についても、位置ｘｋ（ｔＭ）における仮想粒子の濃度ｃ（チルダ）ｋ（ｔＭ）が未知の場合と同様である。
また、位置ｘｋ（ｔＭ）における仮想粒子の濃度ｃ（チルダ）ｋ（ｔＭ）が既知の場合、ベクトルｂを式（５６）に示すようにする。
【０２１０】
【数５６】
2014119882.tif 000058
【０２１１】
なお、式（５１）、式（５２）、式（５３）、式（５４）、式（５５）、式（５６）は、それぞれ、式（１１）、式（１２）、式（９）、式（２６）、式（２４）、式（２３）の一例に該当する。
式（３５）に示すＮ元連立一次方程式の解法として、ＢｉＣＧＳｔａｂ法を用いた。また、シュミット数Ｓｃ＝１０とした。なお、第２の解析例では、濃度場の移流拡散方程式の時間項を完全陰解法的に離散化し、ラプラシアンを二次精度中央差分法的に離散化した場合の例を示している。
【０２１２】
図１３は、ターゲット平面Ａ４３（図１０）における混合度ＭＩの計算結果の例を、計算格子点間間隔を変えながら格子法と比較した結果を示す説明図である。図１３において、黒丸で示す各点は、第３の実施形態の解析方法（以下、第２の解析例において「本手法」と称する）における計算結果を示す。一方、白丸で示す各点は、格子法における計算結果を示す。なお、混合度ＭＩの定義は式（４８）に示している。
【０２１３】
格子法による解析には市販のＣＦＤコードを使用し、本手法による解析には本発明者がＣ＋＋にて生成したコードを使用した。
また、格子法における最大計算格子点数は、約７６００万であり、その際の空間解像度は、一方向の単位長さあたりの点数で１２６であった。格子法に使用した計算格子では、一般的な計算機の限界ともいえるほど大規模な最大計算格子点数となっている。それにもかかわらず格子法では、数値拡散により拡散の過大評価を引き起こしており未だ計算格子依存性が強く残っている。その結果、格子法では正しい解析結果を得られていない。
【０２１４】
これに対し、本手法では、計算格子点に到達する仮想粒子のスカラー場を求めるので、仮想粒子のスカラー量を計算格子点のスカラー量に内挿する必要がない。また、逆追跡終了位置においてスカラー量既知の場合は、計算格子点のスカラー量を仮想粒子のスカラー量に内挿する必要がない。これらの点において、本手法は、根本的に数値拡散が生じにくいという特徴を有する。これにより、本手法は、圧倒的に小さい計算格子依存性を示し、解析は収束に至っている。すなわち、計算格子点間隔をより小さく（短く）しても、解があまり変化しない状態となっている。
なお、第２の実施形態における解析方法でも同様の解析を行い、本手法と比較した。その結果、本手法では、第２の実施形態における解析方法と比較して、１０分の１以下の解析時間、および、１０分の１程度の記憶容量で解析を行うことができた。
【０２１５】
図１４は、格子法と本手法との、濃度分布の解析結果の例を示す説明図である。上側の図Ｆ２１は、格子法の最大空間解像度において得られた濃度分布を示し、下側の図Ｆ２２は、本手法にて得られた濃度分布を示す。
図Ｆ２１とＦ２２とを比較すると、図Ｆ２２のほうが荒い計算格子にて図Ｆ２１と同程度の濃度分布を得られている。
【０２１６】
格子法の最大空間解像度における所要解析時間が約１０時間であったのに対し、本手法では、混合度でそれと同程度の精度を与える空間解像度（計算格子点数約２６万、単位長さあたり３２点）において５分間程度であった。従って、本手法では、格子法の約１２０分の１の時間で同程度の解析結果を得ることができた。また、使用した記憶容量は、格子法で約５０ギガバイトであったのに対し、本手法では２ギガバイト未満であった。従って、本手法では、格子法の２５分の１程度の記憶容量で同程度の解析結果を得ることができた。
【０２１７】
図１５は、第１の解析例における解析結果と、第２の解析例における解析結果と、格子法による解析結果との比較例を示す説明図である。同図に示す各図のうち、図Ｆ３１１〜Ｆ３１４は格子法による解析結果の例を示す。また、図Ｆ３２１〜Ｆ３２４は第２の解析例における解析結果の例を示す。また、図Ｆ３３１〜Ｆ３３４は第１の解析例における解析結果の例を示す。
【０２１８】
また、左側の図ほど拡散の度合いが大きくなっている。具体的には、図Ｆ３１１、Ｆ３２１，Ｆ３３１は、シュミット数Ｓｃ＝１に設定されている。また、図Ｆ３１２、Ｆ３２２，Ｆ３３２は、Ｓｃ＝１０に設定されている。また、図Ｆ３１３、Ｆ３２３，Ｆ３３３は、Ｓｃ＝１００に設定されている。また、図Ｆ３１４、Ｆ３２４，Ｆ３３４は、Ｓｃ＝１０００に設定されている。
【０２１９】
図Ｆ３１３と図Ｆ３１４とでは、拡散の度合の設定が異なるにもかかわらず、同様の解析結果となっている。これは、格子法における数値拡散の影響が拡散の度合の変化よりも大きくなったためと考えられる。すなわち、図Ｆ３１４においては、格子法では数値拡散の影響を受けて高精度な解析を行えていない。一方、図Ｆ３２３、Ｆ３２４、Ｆ３３３、Ｆ３３４では、拡散の度合いに応じた解析結果を得られている。このように、第１の実施形態における解析方法（第１の解析例における解析結果）や、第３の実施形態における解析方法（第２の解析例における解析結果）では、数値拡散の影響をほとんど受けずに高精度な解析を行うことができている。
【０２２０】
また、図Ｆ３２１〜Ｆ３２４では、いずれも良好な解析結果を得られている。このように、第３の実施形態における解析方法では、拡散の度合いが大きい場合と小さい場合とのいずれにおいても、良好な解析結果を得ることができる。
一方、図Ｆ３３１〜Ｆ３３４に関しては、拡散の度合いが大きい図Ｆ３３１では、やや解析精度が低下しているが、拡散の度合いが小さい図Ｆ３３３やＦ３３４では良好な解析結果を得られている。第１の実施形態における解析方法では、拡散の度合いが比較的小さい場合に、短い処理時間および少ない記憶容量にて高精度な解析結果を得ることができる。
【０２２１】
なお、解析装置１００や２００や３００の全部または一部の機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することで各部の処理を行ってもよい。なお、ここでいう「コンピュータシステム」とは、ＯＳや周辺機器等のハードウェアを含むものとする。
また、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。
また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ−ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムを送信する場合の通信線のように、短時間の間、動的にプログラムを保持するもの、その場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリのように、一定時間プログラムを保持しているものも含むものとする。また上記プログラムは、前述した機能の一部を実現するためのものであっても良く、さらに前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるものであっても良い。
【０２２２】
以上、本発明の実施形態を図面を参照して詳述してきたが、具体的な構成はこの実施形態に限られるものではなく、この発明の要旨を逸脱しない範囲の設計変更等も含まれる。
【符号の説明】
【０２２３】
１００、２００、３００  解析装置
１１０  前提条件取得部
１２０、２２０  逆追跡部
１３０、２３０、３３０  場の値取得部
１３１  場の値設定部
１３２  場の値更新部
１４０  解析結果出力部
２３１  場の値算出部
２３２  収束判定部
３３１  方程式取得部
３３２  解析実行部
【図１】
2014119882.tif 000059
【図２】
2014119882.tif 000060
【図３】
2014119882.tif 000061
【図４】
2014119882.tif 000062
【図５】
2014119882.tif 000063
【図６】
2014119882.tif 000064
【図７】
2014119882.tif 000065
【図８】
2014119882.tif 000066
【図９】
2014119882.tif 000067
【図１０】
2014119882.tif 000068
【図１１】
2014119882.tif 000069
【図１３】
2014119882.tif 000070
【図１２】
2014119882.tif 000071
【図１４】
2014119882.tif 000072
【図１５】
2014119882.tif 000073