(57)【要約】  優れたシミュレーション結果を得ることができるシミュレーション装置を提供すること。  移流項を含む方程式の解を算出するシミュレーション装置であって、速度場を取得する速度場取得部と、方程式に対して、少なくとも２つの操作を加えて得られる離散化方程式に、取得した速度場を適用して、離散化方程式における解を算出する解算出部とを具備し、上述の２つの操作は、移流項を保存型に変形することと、移流項の保存型の誤差を修正する項を含めることであることを特徴とする。
2014069421.tif 000002
【特許請求の範囲】
【請求項１】
移流項を含む方程式の解を算出するシミュレーション装置であって、
速度場を取得する速度場取得部と、
前記方程式に対して、少なくとも２つの操作を加えて得られる離散化方程式に、前記速度場取得部で取得した速度場を適用して、前記離散化方程式における解を算出する解算出部と
を具備し、
前記２つの操作は、前記移流項を保存型に変形することと、前記移流項の保存型の誤差を修正する項を含めることであること
を特徴とするシミュレーション装置。
【請求項２】
前記移流項の保存型の誤差を修正する項が、前記速度場が連続の式を満たさないことにより生じる誤差を修正する項である、請求項１に記載のシミュレーション装置。
【請求項３】
前記誤差を修正する項が、前記連続の式に対して、前記方程式に対する積分と同じ積分を少なくとも行うことで得られる項である、請求項２に記載のシミュレーション装置。
【請求項４】
前記解算出部は、少なくとも連続の式を用いて、次の時間ステップにおける速度場を算出し、
前記速度場取得部が取得する速度場は、前記解算出部が算出した速度場である、
請求項２に記載のシミュレーション装置。
【請求項５】
前記速度場取得部は、予め記憶された速度場を取得する、請求項２に記載のシミュレーション装置。
【請求項６】
移流項を含む方程式の解を算出するシミュレーション方法であって、
速度場を取得する第１の過程と、
移流項を含む方程式に対して、少なくとも２つの操作を加えて得られる離散化方程式に、前記第１の過程で取得した速度場を適用して、前記離散化方程式における解を算出する第２の過程と
を有し、
前記２つの操作は、前記移流項を保存型に変形することと、前記移流項の保存型の誤差を修正する項を含めることであること
を特徴とするシミュレーション方法。
【請求項７】
コンピュータを、
速度場を取得する速度場取得部、
移流項を含む方程式に対して、少なくとも２つの操作を加えて得られる離散化方程式に、前記速度場取得部で取得した速度場を適用して、前記離散化方程式における解を算出する解算出部
として機能させるためのプログラムであって、
前記２つの操作は、前記移流項を保存型に変形することと、前記移流項の保存型の誤差を修正する項を含めることであること
を特徴とするプログラム。
【請求項８】
複数個の観測器からの情報に基づき、流体の発生源情報を推定する拡散物質の発生源の推定をするシミュレーション装置であって、
前記観測器の位置情報および計測した濃度情報を入手する観測情報入手手段と、
拡散物質の濃度の推定を予定する領域に複数の仮想放出地点を設定する仮想放出地点設定手段と、
前記観測器の位置情報、計測した濃度情報、および前記仮想放出地点の情報に基づいて前記仮想放出地点での拡散物質の濃度の推定をする請求項１に記載の装置と、
を有するシミュレーション装置。
【請求項９】
複数個の観測器からの情報に基づき、流体の発生源情報を推定する拡散物質の発生源の推定をするシミュレーション方法であって、
前記観測器の位置情報および計測した濃度情報を入手する観測情報入手ステップと、
拡散物質の濃度の推定を予定する領域に複数の仮想放出地点を設定する仮想放出地点設定ステップと、
前記観測器の位置情報、計測した濃度情報および前記仮想放出地点の情報に基づいて請求項６に記載の方法によって前記仮想放出地点での拡散物質の濃度の推定をするステップと、
を有するシミュレーション方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーション装置、シミュレーション方法およびプログラムに関する。
【背景技術】
【０００２】
従来、流体などの連続体の数値シミュレーションに用いられる方法として、有限体積法がある。この有限体積法では、シミュレーションする領域をセルあるいはコントロールボリュームと呼ばれる領域に分割しておく。そして、各セルについて、支配方程式を体積積分し、さらに離散化した式を用いて、解を算出する。通常、流体をシミュレーションするときは、支配方程式であるナビエ・ストークス方程式を体積積分した後、離散化した式を用いる。この体積積分をする際に、移流項と呼ばれる項については、保存型と呼ばれる形式に変換しておく。これにより、移流項の体積積分の結果は表面積分となる（例えば、非特許文献１のＰ．２７参照）。なお、この非特許文献１においては、移流項のことを対流項と呼んでいる。また、非圧縮性流れに対する完全保存型については、離散化の差分スキームについて研究がされている（非特許文献２、３）。【０００３】
上述の移流項について説明する。図１９に示すように連続体中に検査体積Ｖ（ｔ）を取る。連続体が移動、変形するときは、検査体積Ｖ（ｔ）も連続体の移動、変形に伴って移動、変形する。図１９は、時刻ｔ→ｔ＋Δｔと、Δｔだけ時間進行するとき、原点Ｏ’（ｔ）の位置にあった検査体積Ｖ（ｔ）が、Ｏ’（ｔ＋Δｔ）の位置へ移動し、検査体積Ｖ（ｔ＋Δｔ）、表面Ｓ（ｔ＋Δｔ）の状態へ、移動、変形する様子を示す。この移動、変形は、Ｖ（ｔ）を囲む表面Ｓ（ｔ）が速度ベクトルｕｓ（ｔ）の速度で移動することにより、生じている。【０００４】
Ｖ→０、Δｔ→０の極限操作により、連続体を表す偏微分方程式が得られる。この偏微分方程式の記述方法には、ラグランジュ座標系による記述方法と、オイラー座標系による記述方法がある。ラグランジュ座標系は、図１９の原点Ｏ’（ｔ）のように、連続体とともに移動する原点を持つ座標系である。オイラー座標系は、図１９の原点Ｏのように、連続体とは別に設定された座標系である。オイラー座標系では、原点Ｏは、連続体とともに移動しない。【０００５】
連続体として表される物理量をφとし、φの偏微分方程式が、ラグランジュ座標系で以下の式（１−１）のように表されるものとする。式（１−１）をオイラー座標系で表すと、式（１−２）となる。【０００６】
【数１】
2014069421.tif 000003
【０００７】
ここで、ｕ、ｖ、ｗは、連続体内の座標（ｘ、ｙ、ｚ）の位置の速度ベクトルｕのｘ、ｙ、ｚ成分であり、速度ベクトルｕは、式（１−４）のように表される。【０００８】
【数２】
2014069421.tif 000004
【０００９】
式（１−３）の左辺にあるＤ／Ｄｔは、全微分あるいは物質微分と呼ばれる。式（１−３）の右辺第２項から第４項（式（１−５））が、移流項と呼ばれる。【００１０】
【数３】
2014069421.tif 000005
【００１１】
添え字とアインシュタインの総和規則を用いると、式（１−２）、（１−３）は、式（１−６）、（１−７）となる。ここで添え字ｉ、ｊは、ｉ＝１、２、３、ｊ＝１、２、３のインデックスを示し、１，２，３→ｘ、ｙ、ｚに対応する。同じ添え字が２度表れたときは、総和規則に従う。【００１２】
【数４】
2014069421.tif 000006
【００１３】
式（１−６）に示されている物質微分Ｄ／Ｄｔは、次に示すように、連続体として表現される全ての物理量の、時間進行による連続体の移動、変形を表す項として、オイラー座標系で表された偏微分方程式の中に現れる。式（１−３）、（１−７）で示したように、物質微分は、移流項を含んでいるので、物質微分が現れる偏微分方程式は、移流項を含んでいる。【００１４】
物質微分を含む偏微分方程式の例をいくつか以下に示す。
１）流体力学（ナビエ・ストークス方程式）
流体力学の支配方程式であるナビエ・ストークス方程式は、式（１−８）の様に表される。式（１−８）において、ｔ：時間、ｘｉ（ｉ＝１、２、３）：カーテシアン系における座標、ρ：流体の密度、μ：流体の粘性係数、ｕｉ（ｉ＝１、２、３）：流体の速度成分、ｐ：圧力である。また、添え字ｉ（ｉ=１、２、３）、ｊ（ｊ＝１、２、３）はカーテシアン座標系における各方向成分を示す。式（１−８）では、左辺に物質微分が現れている。すなわち、移流項を含んでいる。【００１５】
【数５】
2014069421.tif 000007
【００１６】
２）応力場の方程式
オイラー座標系における応力場の方程式は、式（１−９）である。式（１−９）において、ρ：連続体の密度、ｕｉ（ｉ＝１、２、３）：連続体の変形速度、σｉｊ（ｉ、ｊ＝１、２、３）：連続体の内部応力、ｆｉ＝連続体に作用する外力（例えば、重力）である。式（１−９）でも、左辺に物質微分が現れている。すなわち、移流項を含んでいる。【００１７】
【数６】
2014069421.tif 000008
【００１８】
３）エネルギー保存の方程式
エネルギー保存則は、熱エネルギー保存と運動エネルギー保存とに分けられる。運動エネルギーの保存は、式（１−９）と同じであるので、ここでは熱エネルギー保存の方程式を式（１−１０）に示す。式（１−１０）において、Ｕ：連続体の内部エネルギー、ｑｊ（ｊ＝１、２、３）：熱流束ベクトル成分、ρ：連続体の密度、γ：熱源（熱エネルギーのソース／シンク）、σｉｊ（ｉ、ｊ＝１、２、３）：連続体の内部応力テンソル、Ｄｉｊ（ｉ、ｊ＝１、２、３）：連続体の変形速度テンソルである。式（１−１０）でも、左辺に物質微分が現れている。すなわち、移流項を含んでいる。【００１９】
【数７】
2014069421.tif 000009
【００２０】
４）移流拡散方程式
ある物質Ｃの連続体内での移流拡散現象は、以下の式（１−１１）である。式（１−１１）において、Ｃ：物質Ｃの濃度、μｃ：物質Ｃの拡散係数、ｑｃ：物質Ｃのソース・シンク、ρ：連続体の密度、ｕｊ（ｊ＝１、２、３）：連続体の変形速度成分である。式（１−１１）でも、左辺に物質微分が現れている。すなわち、移流項を含んでいる。【００２１】
【数８】
2014069421.tif 000010
【先行技術文献】
【非特許文献】
【００２２】
【非特許文献１】Ｈ．Ｋ．Ｖｅｒｓｔｅｅｇ、Ｗ．Ｍａｌａｌａｓｅｋｅｒａ（原著）、松下洋介、齋藤泰洋、青木秀之、三浦隆利（共訳）、“数値流体力学［第２版］”、（Ａｎ  Ｉｎｔｒｏｄｕｃｔｉｏｎ  ｔｏ  Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ｆｌｕｉｄ  Ｄｙｎａｍｉｃｓ  Ｓｅｃｏｎｄ  Ｅｄｉｔｉｏｎ）、第２版、森北出版株式会社、２０１１年５月３０日．
【非特許文献２】森西洋平、〔竜門賞受賞記念解説〕非圧縮性流体解析のための完全保存形高次精度差分スキーム、ながれ １９（２０００）、Ｐ．１６４−１７０．
【非特許文献３】Ｍｏｒｉｎｉｓｈｉ，Ｙ．，Ｌｕｎｄ，Ｔ．Ｓ．，Ｖａｓｉｌｙｅｖ，Ｏ．Ｖ．，Ｍｏｉｎ，Ｐ．，“Ｆｕｌｌｙ  ｃｏｎｓｅｒｖａｔｉｖｅ  ｈｉｇｈｅｒ  ｏｒｄｅｒ  ｆｉｎｉｔｅ  ｄｉｆｆｅｒｅｎｃｅ  ｓｃｈｅｍｅｓ  ｆｏｒ  ｉｎｃｏｍｐｒｅｓｓｉｂｌｅ  ｆｌｏｗ，Ｊ．Ｃｏｍｐｕｔ．Ｐｈｙｓ．，Ｖｏｌ．１４３（１９９８），９０−１２４．
【発明の概要】
【発明が解決しようとする課題】
【００２３】
上述したように移流項は、ナビエ・ストークス方程式だけでなく、物質微分を含む方程式（以下、輸送方程式という）に現れる項である。本発明が解決しようとする課題を説明する前に、移流項を保存型にして解く方法の一例として、非特許文献１に記載の有限体積法の場合に用いられている式を説明する。変数φに対する輸送方程式は、式（２−１）あるいは式（２−２）のように表される。上述したように、この式（２−２）の左辺第２項が、移流項である。なお、ナビエ・ストークス方程式は、変数φを速度ｕあるいは運動量ρｕとした輸送方程式と言うことが出来る。【００２４】
【数９】
2014069421.tif 000011
【００２５】
密度ρが一定であると仮定すると、一般に連続体は、連続の式と呼ばれる式（２−３）を満たしている。【００２６】
【数１０】
2014069421.tif 000012
【００２７】
そこで、この式（２−３）の添え字ｉをｊに変えて、両辺にφを乗じた式を、式（２−２）に加えると、式（２−４）が得られる。さらに、式（２−４）の左辺第２項と第３項との和を変形すると式（２−５）が得られる。式（２−５）の左辺第２項が、密度が一定のときの保存型の移流項である。【００２８】
【数１１】
2014069421.tif 000013
【００２９】
一方、連続体の密度ρが一定でなく、ρ＝ρ（ｔ，ｘ，ｙ，ｚ）である場合は、まず、式（２−１）、式（２−２）の両辺にρを乗じて、式（２−６）、式（２−７）を得る。また、密度ρが一定でない場合の連続の式は、式（２−８）である。【００３０】
【数１２】
2014069421.tif 000014
【００３１】
式（２−４）のときと同様に、式（２−８）の添え字ｉをｊに変えて、両辺にφを乗じた式を、式（２−７）に加えると、式（２−９）が得られる。さらに、式（２−９）の左辺を変形すると式（２−１０）が得られる。式（２−１０）の左辺第２項が、密度が一定でない（密度が変化する）ときの保存型の移流項である。【００３２】
【数１３】
2014069421.tif 000015
【００３３】
式（２−１０）において、ρφ→φ’、ρＦ（φ）→Ｆ’（φ’）と置き換えると、式（２−１０）は、式（２−５）と同じ形で表すことができるので、以降では、式（２−５）の保存型の輸送方程式を基に説明する。【００３４】
式（２−５）を、コントロールボリュームの体積で積分すると、式（２−１１）が得られる。また、連続の式（２−３）についても、同様にコントロールボリュームの体積で積分すると、式（２−１２）が得られる。【００３５】
【数１４】
2014069421.tif 000016
【００３６】
積分方程式である式（２−１１）、式（２−１２）の各々を、ｍ個の境界面Ｅａｂ（ｂは１〜ｍ）を持つコントロームボリュームａについて離散化すると、以下の式（２−１３）、式（２−１４）が得られる。なお、式（２−１３）の左辺第３項である[Ｆ（φ）の離散化項］は、Ｆ（φ）を体積積分して離散化した項であるが、詳細な説明は省略する。【００３７】
【数１５】
2014069421.tif 000017
【００３８】
本発明が解決しようとする課題を分かり易くするために、式（２−１３）において、Ｆ（φ）＝０とした式（２−１５）を用いて説明する。【００３９】
【数１６】
2014069421.tif 000018
【００４０】
上述したように、式（２−１５）は、変数φに対する保存型の輸送方程式を、体積積分し、離散化した式である。すなわち、物理量φが、連続体の移動、変形に伴って、移流する様子を表した式である。ある時点で、連続体の全ての領域でφが同じ値であれば、連続体が移動、変形していても、式（２−１５）でシミュレーションした結果は、連続体の全ての領域でφが変化せず、同じ値となる必要がある（条件Ａ）。例えば、物理量φが水の塩分濃度であり、全ての領域で塩分濃度φが同じであるときには、時間が経過して、どのように水が流れても、塩分濃度φは変化しない。この条件Ａについて、図２０に示す１次元流れ問題を例に考える。【００４１】
図２０に示す１次元流れでは、流体はｘ軸方向の正の向きに流れており、この流れにより物理量φが移流している。この１次元流れの流路を、図２０に示すように、ｘ軸に沿ってコントロールボリュームａ、ｂ、ｃに分割する。なお、コントロールボリュームｃは、コントロールボリュームａの上流側に隣接しており、コントロールボリュームｂは、コントロールボリュームａの下流側に隣接している。流路の断面積は一定であり、１とする。すなわち、式（２−１６）が成り立つ。
Ｓａｂ＝Ｓａｃ＝１  ・・・（２−１６）
また、コントロールボリューム間の境界面は、ｘ軸に垂直であるとすると、式（２−１７）が成り立つ。各コントロールボリュームの幅は同じで、Δｘであるとすると、式（２−１８）が成り立つ。また、コントロールボリュームａの境界面における速度ベクトルｕａｂ、ｕａｃは、式（２−１９）で表されるとする。【００４２】
【数１７】
2014069421.tif 000019
【００４３】
これらを用いて、式（２−１５）を表すと、式（２−２０）となる。同様に、これらを用いて連続の式（２−１４）を表すと、式（２−２１）となる。【００４４】
【数１８】
2014069421.tif 000020
【００４５】
したがって、ｕａｂ＝ｕａｃ＝Ｕ０として、式（２−２０）の両辺をΔｘで割ると、式（２−２２）が得られる。【００４６】
【数１９】
2014069421.tif 000021
【００４７】
ここで、流路全体で、φが一定である、すなわち、φａｂ＝φａｃとすると、式（２−２２）から、式（２−２３）が得られる。この式（２−２３）は、φａは、時間が進行しても増減せずに一定の値となることを示す。これは、式（２−１５）の説明で述べた条件Ａを満たしている。【００４８】
【数２０】
2014069421.tif 000022
【００４９】
ところが、実際のシミュレーション計算の結果は、数値計算上の誤差を含んでいる。また、測定値は測定誤差を含んでいる。そのため、数値計算により得られた速度であっても、測定により得られた速度であっても、離散化された連続の式（２−１４）を完全に満たしていない。すなわち、式（２−２１）は満たされず、ｕａｂ−ｕａｃ≠０である。
そこで、その誤差をδＵ０とし、ｕａｂ−ｕａｃ＝δＵ０とする。これを、式（２−２０）に代入すると、φａｂ＝φａｃ＝φ＝一定としても、式（２−２３）とは異なる式（２−２４）が得られる。【００５０】
【数２１】
2014069421.tif 000023
【００５１】
この式（２−２４）は、δＵ0が正の値であるときは、コントロールボリュームａには、｜δＵ0・φ／Δｘ｜の大きさの擬似的なシンク（sink）があり、時間が進行するとφａは減少していくことを示す。また、δＵ0が負の値であるときは、コントロールボリュームａには、｜δＵ0・φ／Δｘ｜の大きさの擬似的なソース（source）があり、時間が進行するとφａは増加していくことを示す。【００５２】
すなわち、速度が連続の式を満たしていないときは、式（２−１５）を用いても条件Ａは満たされない。
ここでは、１次元流れ問題を例に考えたが、２次元、３次元の輸送方程式でも同様である。つまり、連続の式（２−１４）が完全には満足されず、式（２−２５）であるときは、計算領域全てでφが一定であっても、式（２−２６）となる。φａは、時間が進行すると増減してしまい、条件Ａは満たされない。【００５３】
【数２２】
2014069421.tif 000024
【００５４】
以上を纏めると、以下の１）、２）が言える。
１)数値シミュレーションによって算出された連続体の速度は、完全に連続の式を満たすことは無い。有効桁数を大きくし、十分に高精度な計算を行っても、誤差は発生する。また、計測により得られた連続体の速度も、計測誤差を必ず含んでいるので、完全に連続の式を満たすことは無い。
２）移流項を保存型にした輸送方程式では、速度が連続の式を満足していないと、その誤差に比例した大きさの擬似的なソース／シンク項が発生してしまう。【００５５】
さらに、１）で述べたように、擬似的なソース／シンク項の大きさは、誤差に比例していることから以下の３）、４）が言える。
３）誤差が小さいときは、擬似的なソース／シンク項の大きさも小さくなる。この場合、対象の方程式に含まれている拡散項や、メッシュ分割による数値拡散があるため、擬似的なソース／シンク項により発生した増減は拡散され、目立たなくなる。したがって、良好な解が得られているように見える。
４）解析領域が複雑な形状の場合は、メッシュ分割をすると、歪んだ形状のコントロールボリュームが生成されることが多い。コントロールボリュームの形状が歪んでいると、数値シミュレーションによって速度を算出する際に、誤差が大きくなりやすい。このため、擬似的なソース／シンク項により発生した増減も大きくなり、数値拡散では抑えきれずに、解の精度を悪化させて異常な解を生じさせる。さらに、著しい場合には、計算を発散させてしまうこともある。【００５６】
このように、移流項を保存型にした輸送方程式を用いたシミュレーションでは、速度が連続の式を満足していないと、異常な解が生じてしまうことがあるという問題がある。
非特許文献２、３では、積分方程式が解析的には保存則を満足していても、これらを離散化した時点で保存則が満足されない問題点について提起している。しかし、非特許文献２、３では、離散化式において保存則を満たすように変形することによって、この保存則が満足されない問題についての解決方法を提示している。すなわち、非特許文献２、３は、数値計算上の誤差の問題については何ら解決策を提示していない。このため、非特許文献２、３の方法によって離散化段階で保存則を満足することによって精度を向上しても、離散化方程式の入力データとなる速度場が誤差を持っている場合には、依然として異常な解が生じてしまう問題がある。【００５７】
本発明は、このような事情に鑑みてなされたもので、その目的は、良好なシミュレーション結果を得ることができるシミュレーション装置、シミュレーション方法およびプログラムを提供することにある。
【課題を解決するための手段】
【００５８】
上述したように、従来、輸送方程式を保存型に変形する際には、連続の式（２−３）、（２−８）が完全に満たされていることを前提にして、連続の式（２−３）、（２−８）にφを乗じたものを、輸送方程式に加算している。しかし、上述したように、数値シミュレーションによって算出したものでも、測定により得たものでも、速度は誤差を含んでおり、連続の式を完全には満足していない。すなわち、この前提は成り立っていないために、保存型にした移流項に誤差が生じてしまう。そこで、本発明では、この保存型にした移流項で生じる誤差を修正する項を、方程式に含めておき、該方程式の解を算出する。
従来の研究は、離散化上の問題に着目して行われてきた。しかし、工業的に数値シミュレーションを利用する場合には、通常、莫大な容量の計算で、計算対象の形状が複雑で、複雑な現象の再現が求められる。このため、本発明者らは、数値シミュレーションを利用する場合に、単に離散化上の精度向上の対策では不十分であり数値計算上の誤差の存在を前提にした対策が重要であることに着目して、この発明をした。【００５９】
（１）この発明は上述した課題を解決するためになされたもので、本発明の一態様は、移流項を含む方程式の解を算出するシミュレーション装置であって、速度場を取得する速度場取得部と、前記方程式に対して、少なくとも２つの操作を加えて得られる離散化方程式に、前記速度場取得部で取得した速度場を適用して、前記離散化方程式における解を算出する解算出部とを具備し、前記２つの操作は、前記移流項を保存型に変形することと、前記移流項の保存型の誤差を修正する項を含めることであることを特徴とする。【００６０】
（２）また、本発明の他の態様は、上述のシミュレーション装置であって、前記移流項の保存型の誤差を修正する項が、前記速度場が連続の式を満たさないことにより生じる誤差を修正する項である。【００６１】
（３）また、本発明の他の態様は、上述のシミュレーション装置であって、前記誤差を修正する項が、前記連続の式に対して、前記方程式に対する積分と同じ積分を少なくとも行うことで得られる項である。【００６２】
（４）また、本発明の他の態様は、上述のシミュレーション装置であって、前記解算出部は、少なくとも連続の式を用いて、次の時間ステップにおける速度場を算出し、前記速度場取得部が取得する速度場は、前記解算出部が算出した速度場である。【００６３】
（５）また、本発明の他の態様は、上述のシミュレーション装置であって、前記速度場取得部は、予め記憶された速度場を取得する。【００６４】
（６）また、本発明の他の態様は、移流項を含む方程式の解を算出するシミュレーション方法であって、速度場を取得する第１の過程と、移流項を含む方程式に対して、少なくとも２つの操作を加えて得られる離散化方程式に、前記第１の過程で取得した速度場を適用して、前記離散化方程式における解を算出する第２の過程とを有し、前記２つの操作は、前記移流項を保存型に変形することと、前記移流項の保存型の誤差を修正する項を含めることであることを特徴とする。【００６５】
（７）また、本発明の他の態様は、コンピュータを、速度場を取得する速度場取得部、移流項を含む方程式に対して、少なくとも２つの操作を加えて得られる離散化方程式に、前記速度場取得部で取得した速度場を適用して、前記離散化方程式における解を算出する解算出部として機能させるためのプログラムであって、前記２つの操作は、前記移流項を保存型に変形することと、前記移流項の保存型の誤差を修正する項を含めることであることを特徴とするプログラムである。
（８）また、本発明の他の態様は、複数個の観測器からの情報に基づき、流体の発生源情報を推定する拡散物質の発生源の推定をするシミュレーション装置であって、前記観測器の位置情報および計測した濃度情報を入手する観測情報入手手段と、拡散物質の濃度の推定を予定する領域に複数の仮想放出地点を設定する仮想放出地点設定手段と、前記観測器の位置情報、計測した濃度情報、および前記仮想放出地点の情報に基づいて前記仮想放出地点での拡散物質の濃度の推定をする上述のシミュレーション装置と、を有することを特徴とする。
（９）また、本発明の他の態様は、複数個の観測器からの情報に基づき、流体の発生源情報を推定する拡散物質の発生源の推定をするシミュレーション方法であって、前記観測器の位置情報および計測した濃度情報を入手する観測情報入手ステップと、拡散物質の濃度の推定を予定する領域に複数の仮想放出地点を設定する仮想放出地点設定ステップと、前記観測器の位置情報、計測した濃度情報および前記仮想放出地点の情報に基づいて上述のシミュレーション方法によって前記仮想放出地点での拡散物質の濃度の推定をするステップと、を有することを特徴とする。
【発明の効果】
【００６６】
この発明によれば、良好なシミュレーション結果を得ることができる。
【図面の簡単な説明】
【００６７】
【図１】簡単なシミュレーション例におけるコントロールボリュームを示す図である。
【図２】簡単なシミュレーション例における連続の式の値を示す図である。
【図３】従来の方法による簡単なシミュレーション例を示す図である。
【図４】本発明の簡単なシミュレーション例を示す図である。
【図５】この発明の第１の実施形態による流体シミュレーション装置１０の構成を示す概略ブロック図である。
【図６】同実施形態におけるセルを説明する概念図である。
【図７】同実施形態におけるセル構成記憶部１０２が記憶するセル情報の例を示す図である。
【図８】同実施形態におけるセル構成記憶部１０２が記憶するセル間情報の例を示す図である。
【図９】同実施形態における速度場記憶部１０４が記憶する速度場情報の例を示す図である。
【図１０】同実施形態における圧力場記憶部１０７が記憶する圧力場情報の例を示す図である。
【図１１】同実施形態における流体シミュレーション装置１０の動作を説明するフローチャートである。
【図１２】この発明の第２の実施形態による流体シミュレーション装置１０ａの構成を示す概略ブロック図である。
【図１３】この発明の第３の実施形態による連続体シミュレーション装置１０ｂの構成を示す概略ブロック図である。
【図１４】この発明の第４の実施形態による連続体シミュレーション装置１０ｃの構成を示す概略ブロック図である。
【図１５】シミュレーション対象の速度場を示すコンター図である。
【図１６】シミュレーション対象の速度場における速度の発散の大きさを示すコンター図である。
【図１７】温度分布を、従来の方法を用いて算出した場合の例を示すコンター図である。
【図１８】温度分布を、本発明を用いて算出した場合の例を示すコンター図である。
【図１９】オイラー座標系とラグランジュ座標系とを説明する図である。
【図２０】一次元流れにおける各パラメータを説明する図である。
【発明を実施するための形態】
【００６８】
各実施形態を説明する前に、本発明の概略を説明する。まず、連続体の密度ρが一定の場合について説明する。上述したように、従来、輸送方程式（２−１）、（２−２）を、保存型の式（２−５）に変形する際には、式（２−３）でも示した連続の式（３−１）を用いている。
【００６９】
【数２３】
2014069421.tif 000025
【００７０】
しかし、上述したように、シミュレーションや測定により得られた速度は、誤差を含んでいるために、式（３−１）を完全には満たさず、右辺はゼロにならない。そこで、ここでは、式（３−１）の右辺をゼロからＤに置き換えた式（３−２）を用いて、保存型の輸送方程式を導出する。すなわち、連続の式が完全には満たされず、式（３−１）の右辺がゼロにはならないことを前提として、保存型の輸送方程式を導出する。
【００７１】
【数２４】
2014069421.tif 000026
【００７２】
まず、輸送方程式（２−２）の左辺第２項と第３項の間に、+φ・Ｄ−φ・Ｄを挿入すると、式（３−３）が得られる。次に、式（３−３）の左辺第３項のＤに、式（３−２）を代入すると、式（３−４）が得られる。この式（３−４）の左辺第２項と第３項とは、式（２−４）と同じであるので、同様に変形すると、式（３−５）が得られる。
【００７３】
【数２５】
2014069421.tif 000027
【００７４】
本発明では、密度ρが一定のときは、この式（３−５）を保存型の輸送方程式として用いる。すなわち、本発明では、式（３−５）の解を算出する。なお、この（３−５）は、左辺第３項（−φ・Ｄ）を有することを除けば、従来の保存型の輸送方程式（２−５）と同一である。この左辺第３項を有することにより、式（３−２）に示したように、連続の式の右辺がゼロとならず、誤差Ｄを含むことを前提とした、保存型の輸送方程式となっている。この式（３−５）を用いたシミュレーションでは、式（２−２４）に現れたような擬似的なソース／シンクを、−φ・Ｄにより修正することができる。すなわち、擬似的なソース／シンクが現れないようにすることができる。
【００７５】
式（３−５）を用いたシミュレーションの例として、式（３−５）を有限体積法で離散化して、シミュレーションする場合を説明する。なお、以下では、有限体積法で離散化する場合を説明しているが、差分法など、保存型の移流項を持つ方程式を離散化する方法であれば、その他の離散化の方法であってもよい。
【００７６】
まず、式（２−１１）を導出したときと同様に、式（３−５）を、コントロールボリュームの体積Ｖに対して積分すると、式（３−６）が得られる。式（３−６）における各記号は、式（２−１１）と同様である。また、式（３−２）も、同様にコントロールボリュームの体積Ｖに対して積分すると、式（３−７）が得られる。
【００７７】
【数２６】
2014069421.tif 000028
【００７８】
これらの式（３−６）、（３−７）を、式（２−１３）、（２−１４）を導出した際と同様に、ｍ個の境界面Ｅａｂ（ｂは１〜ｍ）を持つコントロームボリュームａについて離散化すると、以下の式（３−８）、式（３−９）が得られる。
【００７９】
【数２７】
2014069421.tif 000029
【００８０】
さらに、式（３−８）に、式（３−９）を代入すると、式（３−１０）が得られる。この離散化した輸送方程式（３−１０）を解くことで、良好なシミュレーション結果を得ることが出来る。
【００８１】
【数２８】
2014069421.tif 000030
【００８２】
次に、連続体の密度ρが一定でない場合を説明する。密度ρが一定でない場合は、密度ρが一定のときと同様に、密度ρが一定でないときの連続の式の右辺をＤとした式（３−１１）を用いて、保存型の輸送方程式を導出する。
【００８３】
【数２９】
2014069421.tif 000031
【００８４】
密度ρが一定でないときの輸送方程式（２−７）の左辺第１項と第２項の間に、+φ・Ｄ−φ・Ｄを挿入すると、式（３−１２）が得られる。次に、式（３−１２）の左辺第３項のＤに、式（３−１１）を代入すると、式（３−１３）が得られる。この式（３−１３）の左辺第１項、第２項は、式（２−９）と同様であるので、同様に変形すると、式（３−１４）が得られる。
【００８５】
【数３０】
2014069421.tif 000032
【００８６】
本発明では、密度ρが一定でないときは、この式（３−１４）を保存型の輸送方程式として用いる。すなわち、式（３−１４）の解を算出する。なお、この（３−１４）は、左辺第３項（−φ・Ｄ）を有することを除けば、従来の保存型の輸送方程式（２−１０）と同一である。この左辺第３項を有することにより、式（３−１１）に示したように、連続の式の右辺がゼロとならず、誤差Ｄを含むことを前提とした、保存型の輸送方程式となっている。この式（３−１４）を用いたシミュレーションでは、式（２−２４）に現れたような擬似的なソース／シンクを、−φ・Ｄにより修正することができる。すなわち、密度ρが一定でないときでも、擬似的なソース／シンクが現れないようにすることができる。
【００８７】
次に、式（３−１４）を用いたシミュレーションの例として、密度ρが一定のときと同様に、式（３−１４）を有限体積法で離散化して、シミュレーションする場合を説明する。式（２−１１）を導出したときと同様に、式（３−１４）を、コントロールボリュームの体積Ｖに対して積分すると、式（３−１５）が得られる。式（３−１５）における各記号は、式（２−１１）と同様である。また、式（３−１１）も、同様にコントロールボリュームの体積Ｖに対して積分すると、式（３−１６）が得られる。
【００８８】
【数３１】
2014069421.tif 000033
【００８９】
これらの式（３−１５）、（３−１６）を、式（３−８）、（３−９）を導出した際と同様に、ｍ個の境界面Ｅａｂ（ｂは１〜ｍ）を持つコントロームボリュームａについて離散化すると、以下の式（３−１７）、式（３−１８）が得られる。
【００９０】
【数３２】
2014069421.tif 000034
【００９１】
さらに、式（３−１７）に、式（３−１８）を代入すると、式（３−１９）が得られる。
【００９２】
【数３３】
2014069421.tif 000035
【００９３】
ここで、式（３−２０）から、式（３−１９）は、式（３−２１）となる。この離散化した輸送方程式（３−２１）を解くことで、密度ρが一定でない場合でも、良好なシミュレーション結果を得ることが出来る。
【００９４】
【数３４】
2014069421.tif 000036
【００９５】
ところで、式（３−１０）、（３−２１）をそれぞれ変形すると、以下の式（３−２２）、（３−２３）が得られる。
【００９６】
【数３５】
2014069421.tif 000037
【００９７】
式（３−２２）、（３−２３）の各々の左辺第２項のΣの中にある、（φａｂ−φａ）に注目する。このφａｂは、コントロールボリュームａ、ｂの境界面Ｅａｂにおけるφの値である。φａｂは、φａとφｂとを用いた補間により算出する。すなわち、φａｂは、式（３−２４）で表すことができる。ここで、ｘａ、ｘｂは、コントロールボリュームａ、ｂ各々の代表点であるコントロールポイントＣa、Ｃｂの座標ベクトルである。
【００９８】
【数３６】
2014069421.tif 000038
【００９９】
このφａｂを算出するための補間関数は、様々なものが従来から知られている。いずれの補間関数（スキームともいう）を用いるかは、シミュレーションの目的、対象とする連続体の種類、数値計算方法などにより決められるが、いずれを用いるようにしてもよい。ここでは、代表的なスキームをいくつか示す。
【０１００】
（１）一次風上スキーム
流体を対象とするときに良く使用されるスキームである。φａｂを導出するための補間関数として、以下の式（３−２５）、（３−２６）を用いる。単位法線ベクトルｎａｂは、コントロールボリュームａに対して外向きのベクトルであるので、式（３−２５）は、コントロールボリュームａに対して、流入の場合の補間関数であり、式（３−２６）は、流出の場合の補間関数である。
【０１０１】
【数３７】
2014069421.tif 000039
【０１０２】
これらを、式（３−２２）、（３−２３）の（φａｂ−φａ）に代入すると、以下の式（３−２７）、（３−２８）となる。
【０１０３】
【数３８】
2014069421.tif 000040
【０１０４】
従って、式（３−２２）、（３−２３）のシグマのうち、流出する境界面の項目は０になる。そこで、流入する境界面についての総和を式（３−２９）のように記載すると、式（３−２２）、（３−２３）は、式（３−３０）、（３−３１）となる。
【０１０５】
【数３９】
2014069421.tif 000041
【０１０６】
これら、式（３−３０）、（３−３１）におけるチルダ付きのシグマは、コントロールボリュームａと隣接するコントロールボリュームｂからφｂが流入する際の連続体の流量Ｓａｂ・（ｎａｂ・ｕａｂ）と同じ流量でφａが流出することを示す式となっている。すなわち、コントロールボリュームａに対する流入/流出の流量の収支が合う式となっている。本発明の有効性を示す一例となっている。
【０１０７】
一次風上スキームを用いた場合の簡単なシミュレーション例を示す。１次元流れで、密度は一定であるとする。図１に示すように、流れ方向に４つのコントロールボリュームａ１、ａ２、ａ３、ａ４に分割する。各境界面における速度は、コントロールボリュームａ１の上流側の境界面から順に、１．０、１．１、０．９、０．８、１．０であったとする。各コントロールボリュームの体積を１、境界面を１とし、この速度分布を各コントロールボリュームにおいて連続の式に代入すると、図２のようにコントロールボリュームａ１から順に、−０．１、０．２、０．１、−０．２となり、いずれも０とならない。
【０１０８】
式（３−３０）は、簡単のためにＦ（φ）＝０、時間微分を１次差分で近似し、Δｔ＝１とすると、式（３−３０’）となる。
φｉ＋ｕｉ（φｉ−１−φｉ）＝
（１−ｕｉ）φａ＋ｕｉφｉ−１＝φｉ’  ・・・（３−３０’）
ここで、ｕｉは、コントロールボリュームａｉの上流側のコントロールボリュームａｉ−１との境界面での速度、φｉは、当該コントロールボリュームａｉの物理量（温度）、φｉ’は、前の時間ステップにおける当該コントロールボリュームの物理量、φｉ―１は、上流側のコントロールボリュームａｉ−１の物理量である。
【０１０９】
一方、従来は、式（２−２０）に同様の条件を当てはめて、式（２−２０’）となる。
φｉ＋ｕｉ・φｉ−１−ｕｉ＋１・φｉ＝
（１−ｕｉ＋１）φｉ＋ｕｉ・φｉ−１＝φｉ’  ・・・（２−２０’）
【０１１０】
ここで、初期条件として全てのコントロールボリュームで温度が１００℃、すなわちφｉ’＝１００℃（ｉ＝１〜４）として、式（２−２０’）、（３−３０’）を用いて、次の時間ステップにおける各コントロールボリュームの温度を算出する。式（２−２０’）を用いたとき、すなわち従来の算出結果例を図３に、式（３−３０’）を用いたとき、すなわち本発明の算出結果例を図４に示す。このように、一様に１００℃であったので、本来は算出結果も一様に１００℃となるべきであるが、従来は、一様な分布にならないが、本発明の算出結果例では、一様に１００℃となる。
【０１１１】
（２）二次風上スキーム、三次風上スキーム
二次風上スキーム、三次風上スキームとも、流体計算で良く使用されるスキームである。これらのスキームの詳細は、ここでは省略するが、φａｂは、式（３−３２）となる。すなわち、φａ、φｂのｇｒａｄ（勾配ベクトル）を使用した補間関数を用いる。
【０１１２】
【数４０】
2014069421.tif 000042
【０１１３】
（１）の風上スキームと異なり、式（３−２８）のように、流出する場合に移流項が単純にゼロにはならない。したがって、式（３−２２）、（３−２３）におけるシグマの項は、コントロールボリュームａに隣接する全てのコントロールボリュームｂに関する移流項の総和となる。
【０１１４】
（３）その他のスキーム
一次、二次、三次風上スキーム以外にも、指数関数ｅｘなど、非線形の関数を利用したスキームがある。また、固体や流れの遅い流体などでは、中心差分と同じ考え方から、式（３−３３）が用いられることがある。
【０１１５】
【数４１】
2014069421.tif 000043
【０１１６】
以上、いくつかのφａｂを導出するためのスキームの例を示したが、どのようなスキームを用いても、連続体の移動、変形速度が連続の式を満たさず誤差を含んだ状態であっても、良好なシミュレーション結果を得ることができる。
【０１１７】
［第１の実施形態］
以下、図面を参照して、本発明の第１の実施形態について説明する。図５は、この発明の第１の実施形態による流体シミュレーション装置１０の構成を示す概略ブロック図である。本実施形態における流体シミュレーション装置１０は、有限体積法により、非圧縮性流体をシミュレーションして、シミュレーションする領域における速度と圧力とを算出するシミュレーション装置である。すなわち、本実施形態では、運動量に関する輸送方程式であるナビエ・ストークス方程式の解を算出するシミュレーション装置を説明する。
【０１１８】
図５に示すように流体シミュレーション装置１０は、セル構成設定部１０１、セル構成記憶部１０２、初期条件設定部１０３、速度場記憶部１０４、速度場取得部１０５、速度・圧力算出部１０６、圧力場記憶部１０７を含んで構成される。セル構成設定部１０１は、シミュレーションする領域（空間）を分割したセルに関する情報を設定し、セル構成記憶部１０２に記憶させる。セル構成記憶部１０２は、セル構成設定部１０１が設定したセルに関する情報を記憶する。なお、セルに関する情報は、各セルの情報であるセル情報と、隣接するセル間の情報であるセル間情報とを含む。
【０１１９】
初期条件設定部１０３は、初期条件として、各セルの速度ベクトルを設定し、速度場記憶部１０４に記憶させる。速度場記憶部１０４は、初期条件設定部１０３が設定した各セルの速度ベクトルを記憶する。また、速度場記憶部１０４は、各時間ステップについて、速度・圧力算出部１０６が算出した各セルの速度ベクトルも記憶する。速度場取得部１０５は、速度・圧力算出部１０６が時間ステップｎに関する解を算出するときは、時間ステップｎ−１における各セルの速度ベクトルを、速度場記憶部１０４から読み出す。
【０１２０】
速度・圧力算出部１０６は、速度場取得部１０５が速度場記憶部１０４から読み出した時間ステップｎ−１における各セルの速度ベクトルと、セル構成記憶部１０２が記憶するセル情報とを用いて、対象時間ステップｎにおける各セルの速度ベクトルと圧力とを算出する。速度・圧力算出部１０６は、算出した各セルの速度ベクトルを、速度場記憶部１０４に記憶させ、圧力を圧力場記憶部１０７に記憶させる。なお、速度・圧力算出部１０６は、速度ベクトルと圧力とを算出する際に、非圧縮性流体のナビエ・ストークス方程式と、非圧縮性流体の連続の式（質量保存の式ともいう）を離散化した離散化方程式を用いる。本実施形態における非圧縮性流体のナビエ・ストークス方程式の離散化方程式は、移流項を後述する保存型に変形することと、移流項の保存型の誤差を修正する項を含める操作を非圧縮性流体のナビエ・ストークス方程式に対して行っている。これらの詳細については後述する。圧力場記憶部１０７は、各時間ステップについて、速度・圧力算出部１０６が算出した各セルの圧力を記憶する。
【０１２１】
図６は、本実施形態におけるセルを説明する概念図である。図６では、符号Ｒ１〜Ｒ１５を付した領域の各々が１〜１５番目のセルである。図６に示すように、各セルは、多面体である。また、シミュレーションする領域を隙間なく埋めるように、セルは配置されている。同図において、代表点Ｃｎ（ｎ＝１〜４）は、それぞれｎ番目のセルＲｎの代表点である。また、境界面Ｅ１２は、セルＲ１とＲ２との境界面である。単位法線ベクトルｎ１２は、境界面Ｅ１２の単位法線ベクトルである。距離α１２は、代表点Ｃ１からＣ２までの距離を１としたときの、点Ｃ１とＣ２を結ぶ線の境界面Ｅ１２との交点から代表点Ｃ１までの距離である。
【０１２２】
図７は、セル構成記憶部１０２が記憶するセル情報の例を示す図である。図７に示すようにセル情報は、全てのセルのセルＩＤ、代表点座標、体積を含む。図７に示す例では、セルＩＤが「０００１」のセルのセル情報は、セルＩＤ「０００１」、代表点座標「１５，２２５，２１４」、体積「２４．３」である。同様に、セルＩＤが「０００２」のセルのセル情報は、セルＩＤ「０００２」、代表点座標「１５，２２５，２１４」、体積「２４．３」であり、セルＩＤが「２４５３」のセルのセル情報は、セルＩＤ「２４５３」、代表点座標「２１５、２５、４４」、体積「３１．２」である。なお、本実施形態において、代表点座標の３つの数値は、ｘ座標と、ｙ座標と、ｚ座標である。以降、Ｘｎは、セルＩＤがｎのセルの代表点座標を示し、Ｖｎは、セルＩＤがｎのセルの体積を示す。
【０１２３】
図８は、セル構成記憶部１０２が記憶するセル間情報の例を示す図である。図８に示すようにセル間情報は、基準とするセルのセルＩＤ、該セルに隣接するセルのセルＩＤ、これら２つのセル間の境界面の面積、該境界面の単位法線ベクトル、これら２つのセルの代表点を通る直線と境界面との交点から基準とするセルの代表点までの距離αを含む。なお、単位法線ベクトルは、長さが「１」の単位ベクトルである。また、該距離は、２つのセルの代表点間の距離を「１」としたときの距離であり、図６のα１２に対応する。
【０１２４】
図８に示す例では、セルＩＤ「０００１」のセルとセルＩＤ「０００５」のセルとのセル間に関するセル間情報は、セルＩＤ「０００１」、セルＩＤ「０００５」、面積「４．２１４５８」、単位法線ベクトル「０．７００２４３，０．５５２４６６，０．８９１９４１」、距離「０．４２２５６１」である。同様に、セルＩＤ「０００１」のセルとセルＩＤ「０１０２」のセルとのセル間に関するセル間情報は、セルＩＤ「０００１」、セルＩＤ「０１０２」、面積「３．１３３５２」、単位法線ベクトル「０．５８４２７６，−０．６１７１２３，０．５２７０４９」、距離「０．５５４８７８」であり、セルＩＤ「２４５３」のセルとセルＩＤ「２４５２」のセルとのセル間に関するセル間情報は、セルＩＤ「２４５３」、セルＩＤ「２４５２」、面積「３．３１２１５」、単位法線ベクトル「０．３３５７８６，０．８１４２６６，０．４７３５１７」、距離「０．５１５４８６」である。
【０１２５】
図９は、速度場記憶部１０４が記憶する速度場情報の例を示す図である。速度場記憶部１０４は、初期条件として設定された初期時間ステップにおける速度場情報と、速度・圧力算出部１０６が算出した各時間ステップにおける速度場情報とを記憶する。図９に示す例は、それらのうちの１つの時間ステップにおける速度場情報である。１つの時間ステップにおける速度場情報は、各セルのセルＩＤと該セルの速度ベクトルとを含む。図９に示す速度場情報の例は、セルＩＤ「０００１」と該セルの速度ベクトル「１．２２１３５４，３．２５５４６，０．２５１３５６」と、セルＩＤ「０００２」と該セルの速度ベクトル「１．２４０１３４，３．１５２２４，０．６６２１５８」と、セルＩＤ「２４５３」と該セルの速度ベクトル「３．２２１５１，３．００１２５，１．６５８５５」とを含む。
【０１２６】
図１０は、圧力場記憶部１０７が記憶する圧力場情報の例を示す図である。圧力場記憶部１０７は、速度・圧力算出部１０６が算出した各時間ステップにおける圧力場情報を記憶する。図１０に示す例は、それらのうちの１つの時間ステップにおける圧力場情報である。１つの時間ステップにおける圧力場情報は、各セルのセルＩＤと該セルの圧力とを含む。図１０に示す圧力場情報の例は、セルＩＤ「０００１」と該セルの圧力「１．２８５６８」と、セルＩＤ「０００２」と該セルの圧力「１．２２３２５」と、セルＩＤ「２４５３」と該セルの圧力「１．５６８４６」とを含む。
【０１２７】
図１１は、本実施形態における流体シミュレーション装置１０の動作を説明するフローチャートである。まず、セル構成設定部１０１が、外部から取得したセルに関する情報（セル情報およびセル間情報）をセル構成記憶部１０２に記憶させる。また、初期条件設定部１０３が、外部から取得した初期条件である初期時間ステップの速度場情報を、速度場記憶部１０４に記憶させる（Ｓ１）。次に、速度場取得部１０５は、対象時間ステップを初期値「１」とする（Ｓ２）。次に、速度場取得部１０５が、速度場記憶部１０４から、対象時間ステップの１ステップ前の速度場情報を取得する（Ｓ３）。なお、速度場取得部１０５は、対象時間ステップが「１」のときは、初期時間ステップを１ステップ前とし、初期時間ステップの速度場情報を、速度場記憶部１０４から取得する。
【０１２８】
次に、速度・圧力算出部１０６は、非圧縮性流体のナビエ・ストークス方程式を離散化した式と、非圧縮性流体の連続の式を離散化した式とからなる連立一次方程式を、Ａ・φ＝ｂとしたときの行列Ａとベクトルｂとを算出する（Ｓ４）。速度・圧力算出部１０６は、ステップＳ３にて速度場取得部１０５が取得した速度場情報を用いて、行列Ａおよびベクトルｂの各要素を算出する。
【０１２９】
なお、ベクトルφは、対象時間ステップにおける全てのセルの速度ベクトルの各要素と圧力とを要素とするベクトルである。すなわち、ベクトルφは、セルの数×４個の要素を持つ。一方、前記連立一次方程式は、各セルに関する、各軸方向のナビエ・ストークス方程式を離散化した式と、各セルに関する非圧縮性流体の連続の式を離散化した式とからなる。すなわち、前記連立一次方程式は、ベクトルφの要素数と同じ数の式からなるため、解の算出が可能である。なお、ステップＳ４の行列Ａと、ベクトルｂとの算出方法についての詳細は、後述する。
【０１３０】
次に、速度・圧力算出部１０６は、前記Ａ・φ＝ｂの解であるベクトルφを算出する。前記Ａ・φ＝ｂは、ベクトルφの各要素を未知数とする連立一次方程式であるので、速度・圧力算出部１０６は、ベクトルφの算出には、例えば、共役勾配法などの公知の連立一次方程式の解法を用いる。速度・圧力算出部１０６は、算出したベクトルφの要素のうちの速度ベクトルを、対象時間ステップの速度ベクトルとして速度場記憶部１０４に記憶させる。同様に、速度・圧力算出部１０６は、算出したベクトルφの要素のうちの圧力を、対象時間ステップの圧力として圧力場記憶部１０７に記憶させる（Ｓ５）。
【０１３１】
次に、速度・圧力算出部１０６は、終了条件が満たされているか否かを判定する（Ｓ６）。この終了条件としては、例えば、対象時間ステップが、予め設定された値と一致することなどが挙げられる。終了条件が満たされていないと判定したときは（Ｓ６−Ｎｏ）、ステップＳ７に遷移して、速度・圧力算出部１０６は、対象時間ステップを「１」増加させ（Ｓ７）、ステップＳ３に戻る。これにより、次の時間ステップについての処理を行う。また、ステップＳ６にて、終了条件が満たされていると判定したときは（Ｓ６−Ｙｅｓ）、処理を終了する。
【０１３２】
次に、速度・圧力算出部１０６による行列Ａ、ベクトルｂの算出処理、すなわち図１１のステップＳ４の詳細を説明する。前述のように、Ａ・φ＝ｂは、非圧縮性流体のナビエ・ストークス方程式を離散化した式と、非圧縮性流体の連続の式を離散化した式との連立一次方程式である。一般に、第ｉ軸方向の非圧縮性流体のナビエ・ストークス方程式は、式（４−１）のように表される。また、一般に、非圧縮性流体の連続の式は、式（４−２）のように表される。なお、本実施形態では、座標系として直交座標系（カーテシアン座標系）を用いる。また、式（４−１）（４−２）を含む以降の各式では、アインシュタインの総和則を用いる。
【０１３３】
【数４２】
2014069421.tif 000044
【０１３４】
式（４−１）において、ｔは時間である。ｕｉ（ｉ＝１，２，３）は、速度の各軸方向の成分である。ｘｉ（ｉ＝１，２，３）は、座標の各軸成分である。pは圧力である。ρは密度である。νは、動粘性係数である。Ｆｉ（ｉ＝１，２，３）は、外力の各軸方向の成分である。式（４−１）では、左辺の第２項が移流項である。なお、式（４−１）は第ｉ軸方向の運動量保存の式であることと、非圧縮性流体であるため流体の密度ρは一定であることとから、式（４−１）における移流項では、移流されるのは第ｉ軸方向の速度ｕｉとなっている。以降では、説明の簡易のために、式（４−１）の右辺の第２項（粘性項）と第３項（外力）との和を−Ｆ（ｕ）と表記し、右辺の全ての項を左辺に移動させた式（４−３）を用いて説明する。
【０１３５】
【数４３】
2014069421.tif 000045
【０１３６】
従来の有限体積法では、ナビエ・ストークス方程式をコントロールボリュームで体積積分し、さらに離散化した式を解いている。有限体積法では、この体積積分の際に、連続の式である式（４−２）の両辺に速度ｕｉを乗じた式（４−４）をナビエ・ストークス方程式に加えて、式（４−５）のように移流項を保存型と呼ばれる形に変形したものを用いる。この移流項の保存型への変形は、速度ベクトルが連続の式（４−２）を満たしていること前提とした変形である。
【０１３７】
【数４４】
2014069421.tif 000046
【０１３８】
一方、本実施形態では、速度ベクトルｕを、連続の式の左辺に代入しても「０」にはならず、式（４−６）に示すように誤差Ｄが存在すると仮定する。例えば、計算容量などの問題から、反復計算により速度を求める際に充分に収束するまで計算できないときなどに、誤差Ｄの大きさが無視できないものとなることがある。このような誤差は、セルの形状が、特定の軸方向に長いときなどに発生しやすい。
【０１３９】
次に、式（４−６）の両辺に速度ｕｉを乗じた式（４−７）をナビエ・ストークス方程式（４−１）に加えると、式（４−８）が得られる。この式（４−８）の移流項を保存型に変形し、さらに式（４−８）の右辺を左辺に移動させると、式（４−９）が得られる。この式（４−９）は、左辺の第３項があることを除けば、従来の有限体積法で用いている式（４−５）と同じであり、第３項に、保存型の移流項における誤差を修正する項（すなわち、−ｕｉＤ）がある。以降、該誤差を修正する項を含んだままで、体積積分および離散化を行う。
【０１４０】
【数４５】
2014069421.tif 000047
【０１４１】
この式（４−９）を、コントロールボリュームで体積積分すると、式（４−１０）が得られる。また、式（４−６）についても、コントロールボリュームで体積積分すると、式（４−１１）が得られる。なお、式（４−１０）、（４−１１）において、Ｖはコントロールボリュームの体積であり、Ｓはコントロールボリュームの表面である。また、式（４−１０）の左辺第２項および式（４−１１）の左辺第１項における太字体のｎとｕは、それぞれコントロールボリュームの表面Ｓの単位法線ベクトルと、速度ベクトルである。単位法線ベクトルｎは、コントロールボリュームの表面Ｓに対して垂直で、外向き、かつ、長さが「１」のベクトルである。
【０１４２】
【数４６】
2014069421.tif 000048
【０１４３】
これらの式（４−１０）、（４−１１）を、例えば、図６のセルＲ１のコントロールボリュームＶ１、境界面Ｅ１ｋ（ｋ∈Ｎ１であり、集合Ｎ１は{２、３、４、１１、１２、１３、１４}である）の面積Ｓ１ｋ、時間Δｔについて離散化し、代数方程式による近似式に変換すると、式（４−１２）、式（４−１３）が得られる。なお、境界面Ｅ１ｋは、セルＲ１とセルＲｋとの間の境界面である。集合Ｎ１は、セルＲ１に隣接するセルのセルＩＤの集合である。また、Δｔは、時間ステップの幅である。
【０１４４】
【数４７】
2014069421.tif 000049
【０１４５】
面積Ｓ１ｋは、セルＲ１とセルＲｋとの間の境界面の面積である。また、Ｄ１は、セルＲ１における連続の式の誤差である。ｕ１ｉは、セルＲ１における速度のｉ番目の軸方向の成分である。ｕ１ｉ（ｔ−Δｔ）は、対象時間ステップの前の時間ステップでのセルＲ１における速度のｉ番目の軸方向の成分である。ｕ１ｋｉは、境界面Ｅ１ｋにおける速度のｉ番目の軸方向の成分である。ｎ１ｋは、境界面Ｅ１ｋにおけるセルＲｋに向いた単位法線ベクトルである。u１ｋは、境界面Ｅ１ｋにおける速度ベクトルである。ｐｋは、セルＲｋにおける圧力である。
次に、式（４−１３）を式（４−１２）の第３項に代入し、Σの項（第２項と第３項）を纏めると、式（４−１４）が得られる。
【０１４６】
【数４８】
2014069421.tif 000050
【０１４７】
本実施形態では、この式（４−１４）を、前記連立一次方程式を構成する一次方程式の一つとしている。ただし、式（４−１４）は、第２項に未知数である速度同士の乗算を含んでおり、一次方程式となっていない。そこで、第２項における速度ベクトルｕ１ｋとして、前の時間ステップの速度ベクトルを用いるという近似を行う。また、Ｆ（ｕ）の離散化項については、説明を簡単にするため、前の時間ステップの速度ベクトルを用いるという近似を行う。
【０１４８】
また、境界面Ｅ１ｋにおける速度ｕ１ｋｉを、セルＲ１における速度ｕ１ｉと、セルＲｋにおける速度ｕｋｉとから補間する。補間方法としては、様々な方法があり、いずれの方法を用いてもよいが、本実施形態では一次風上スキームを用いる。一次風上スキームでは、境界面の速度として、風上側のセルにおける速度を用いる。単位法線ベクトルｎ１ｋは、セルＲ１から見て外向きのベクトルである。このため、単位法線ベクトルｎ１ｋと境界面の速度ベクトルｖ１ｋとの内積が負のときは、境界面Ｅ１ｋにおいてはセルＲｋが風上であり、正のときはセルＲ１が風上である。この一次風上スキームを式に表すと、式（４−１５）となる。なお、前記内積が負のときをセルＲ１へ流入、正のときをセルＲ１から流出ともいう。
【０１４９】
【数４９】
2014069421.tif 000051
【０１５０】
この式（４−１５）を式（４−１４）に代入すると、セルＲ１から流出のときはΣの中は「０」になる。さらに、一次方程式とするために、境界面Ｅ１ｋの速度ベクトルｕ１ｋを、前の時間ステップにおける各セルの速度ベクトルから補間して求めた速度ベクトルｕ１ｋ’とする。本実施形態では、補間方法としては、式（４−１６）に示すように、セル構成記憶部１０２が記憶する距離α１ｋを用いて、前の時間ステップにおける速度ベクトルｕ１’とｕ１’とを比例配分する方法を用いる。
【０１５１】
【数５０】
2014069421.tif 000052
【０１５２】
これらにより、式（４−１４）は、式（４−１７）のようになる。ただし、集合Ｎ１’は、セルＲ１に隣接するセルのうち、セルＲ１へ流入しているセルのセル番号を要素とする集合である。この式（４−１７）を変形すると式（４−１８）となり、セルＲ１の速度のｉ軸成分ｕ１ｉと、隣接する各セルＲｋの速度のｉ軸成分ｕｋｉと、セルＲ１の圧力ｐとを未知数とする一次方程式となる。なお、３次元であれば、ｉ＝１，２，３の各々について、式（４−１８）が立てられる。
【０１５３】
【数５１】
2014069421.tif 000053
【０１５４】
また、連続の式である式（４−２）も、セルＲ１のコントロールボリュームＶ１、境界面Ｅ１ｋ（ｋ∈Ｎ１であり、集合Ｎ１は{２、３、４、１１、１２、１３、１４}である）の面積Ｓ１ｋについて離散化し、代数方程式による近似式に変換すると、式（４−１９）が得られる。前記式（４−１６）と同様の補間を用いると、式（４−１９）は式（４−２０）となり、セルＲ１の速度の各成分と、隣接セルＲｋの速度の各成分を未知数とする一次方程式となる。
【０１５５】
【数５２】
2014069421.tif 000054
【０１５６】
これらの式（４−１８）、（４−２０）は、セルＲ１に関する式であるが、その他のセル各々についても同様の式を立てることができる。速度・圧力算出部１０６は、全てのセルについての式（４−１８）、（４−２０）における各未知数の係数を、行列Ａの要素とし、右辺をベクトルｂの要素とする。
【０１５７】
上述のように、非圧縮性流体のナビエ・ストークス方程式を離散化する際に、連続の式を加えて移流項を保存型に変形している。このときに、本実施形態では、連続の式が誤差Ｄを有すると仮定し、該誤差Ｄを修正する項を含めて離散化している。よって、各時間ステップの速度が、連続の式を完全には満たしていなくても、その誤差による影響を抑えることができるので、優れたシミュレーション結果を得ることができる。
【０１５８】
なお、本実施形態では、ナビエ・ストークス方程式と、連続の式とを連立させて解く場合を説明したが、これらに加えて他の方程式を連立させて解くようにしてもよい。その場合、他の方程式も移流項を含むときは、ナビエ・ストークス方程式と同様に移流項の誤差Ｄを修正する項を含めて離散化するようにしてもよい。これにより、該方程式についても誤差による影響を抑えることができるので、優れたシミュレーション結果を得ることができる。
【０１５９】
また、速度・圧力算出部１０６が各時間ステップの解を算出する際に、共役勾配法などの反復法を用いる場合、解の初期値が必要になるが、初期値として前の時間ステップの解を用いるようにしてもよい。その場合、速度場取得部１０５は、圧力場記憶部１０７から前の時間ステップの圧力を取得し、速度・圧力算出部１０６に入力する。
【０１６０】
［第２の実施形態］
以下、図面を参照して、本発明の第２の実施形態について説明する。第２の実施形態では、圧縮性流体をシミュレーションする場合を説明する。図１２は、この発明の第２の実施形態による流体シミュレーション装置１０ａの構成を示す概略ブロック図である。本実施形態における流体シミュレーション装置１０ａは、有限体積法により、圧縮性流体をシミュレーションして、シミュレーションする領域における速度と圧力と密度とを算出するシミュレーション装置である。
【０１６１】
図１２に示すように流体シミュレーション装置１０ａは、セル構成設定部１０１、セル構成記憶部１０２、初期条件設定部１０３ａ、密度・速度場記憶部１０４ａ、密度・速度場取得部１０５ａ、解算出部１０６ａ、圧力場記憶部１０７を含んで構成される。同図において、図５の各部に対応する部分には、同一の符号（１０１、１０２、１０７）を付し、説明を省略する。
【０１６２】
初期条件設定部１０３ａ、初期条件として、各セルの速度ベクトルと密度とを設定し、密度・速度場記憶部１０４ａに記憶させる。密度・速度場記憶部１０４ａは、初期条件設定部１０３ａが設定した各セルの速度ベクトルと密度とを記憶する。また、密度・速度場記憶部１０４ａは、各時間ステップについて、解算出部１０６ａが算出した各セルの速度ベクトルと密度も記憶する。密度・速度場取得部１０５ａは、解算出部１０６ａが時間ステップｎに関する解を算出するときは、時間ステップｎ−１における各セルの速度ベクトルと密度を、密度・速度場記憶部１０４ａから読み出す。
【０１６３】
解算出部１０６ａは、密度・速度場取得部１０５ａが密度・速度場記憶部１０４ａから読み出した時間ステップｎ−１における各セルの速度ベクトルおよび密度と、セル構成記憶部１０２が記憶するセル情報とを用いて、対象時間ステップｎにおける各セルの速度ベクトルと圧力と密度とを算出する。解算出部１０６ａは、算出した各セルの速度ベクトルと密度とを、密度・速度場記憶部１０４に記憶させ、圧力を圧力場記憶部１０７に記憶させる。
【０１６４】
解算出部１０６ａは、圧縮性流体のナビエ・ストークス方程式を離散化した離散化方程式と、圧縮性流体の連続の式を離散化した離散化方程式と、密度と圧力の関係を表す状態方程式を離散化した離散化方程式とを含む連立一次方程式を解き、各セルの速度ベクトル、圧力、密度を算出する。解算出部１０６ａは、連立一次方程式を解く際に、例えば、共役勾配法などの公知の連立一次方程式の解法を用いる。なお、状態方程式としては、例えば、ｐ・１／ρ＝ＲＴを用いる。ここで、Ｒは気体定数であり、Ｔは温度（熱力学温度）である。本実施形態では、温度Ｔは一定である。なお、温度Ｔは、一定でなく、既知であってもよい。また、温度Ｔが、既知でないときは、温度Ｔを変数とする離散化方定式を、上述の連立一次方程式に含めてもよい。
【０１６５】
上述の圧縮性流体のナビエ・ストークス方程式を離散化した離散化方程式を説明する。一般に、第ｉ軸方向の圧縮性流体のナビエ・ストークス方程式は、式（４−２１）のように表される。また、一般に、非圧縮性流体の連続の式は、式（４−２２）のように表される。なお、本実施形態でも、座標系として直交座標系（カーテシアン座標系）を用いる。式（４−２１）においてμは、粘性係数である。
【０１６６】
【数５３】
2014069421.tif 000055
【０１６７】
本実施形態でも、連続の式である式（４−２２）に速度ベクトルｕを代入しても、誤差Ｄが存在すると仮定する。第１の実施形態と同様に、式（４−２２）の右辺を誤差Ｄとし、両辺にｕｉを乗じると、式（４−２３）が得られる。そして、この式（４−２３）を、ナビエ・ストークス方程式（４−２１）に加えると、式（４−２４）が得られる。
【０１６８】
【数５４】
2014069421.tif 000056
【０１６９】
この式（４−２４）を整理すると、保存型である式（４−２５）となる。この式（４−２５）は、左辺の第３項があることを除けば、従来の有限体積法で用いている式と同様であり、第３項に、保存型の移流項における誤差を修正する項（すなわち、−ｕｉＤ）がある。
【０１７０】
【数５５】
2014069421.tif 000057
【０１７１】
この式（４−２５）を、コントロールボリュームで体積積分すると、式（４−２６）が得られる。また、式（４−２２）の右辺を誤差Ｄとした式についても、コントロールボリュームで体積積分すると、式（４−２７）が得られる。第１の実施形態と同様に、Ｖはコントロールボリュームの体積であり、Ｓはコントロールボリュームの表面である。また、式（４−２６）の左辺第２項における太字体のｎとｕは、それぞれコントロールボリュームの表面Ｓの単位法線ベクトルと、速度ベクトルである。単位法線ベクトルｎは、コントロールボリュームの表面Ｓに対して垂直で、外向き、かつ、長さが「１」のベクトルである。
【０１７２】
【数５６】
2014069421.tif 000058
【０１７３】
これらの式（４−２６）、（４−２７）を、例えば、図６のセルＲ１のコントロールボリュームＶ１、境界面Ｅ１ｋ（ｋ∈Ｎ１であり、集合Ｎ１は{２、３、４、１１、１２、１３、１４}である）の面積Ｓ１ｋ、時間Δｔについて離散化し、代数方程式による近似式に変換すると、式（４−２８）、式（４−２９）が得られる。
【０１７４】
【数５７】
2014069421.tif 000059
【０１７５】
式（４−２８）の第３項に式（４−２９）を代入すると、式（４−３０）が得られる。この式（４−３０）に第１の実施形態と同様に、一次風上スキームなどを適用して、速度ベクトルと、圧力とを未知数とする一次方程式にする。
【０１７６】
【数５８】
2014069421.tif 000060
【０１７７】
この式（４−３０）は、セルＲ１に関する式であるが、その他のセル各々についても同様の式を立てることができる。解算出部１０６ａは、全てのセルについての式（４−３０）における各未知数の係数を、行列Ａの要素とし、右辺をベクトルｂの要素とする。
【０１７８】
上述のように、圧縮性流体のナビエ・ストークス方程式を離散化する際に、連続の式を加えて移流項を保存型に変形している。このときに、本実施形態では、連続の式が誤差Ｄを有すると仮定し、該誤差Ｄを修正する項を含めて離散化している。よって、各時間ステップの速度が、連続の式を完全には満たしていなくても、その誤差による影響を抑えることができるので、優れたシミュレーション結果を得ることができる。
【０１７９】
なお、本実施形態では、ナビエ・ストークス方程式と、連続の式と、状態方程式とを連立させて解く場合を説明したが、これらに加えて他の方程式を連立させて解くようにしてもよい。その場合、他の方程式も移流項を含むときは、ナビエ・ストークス方程式と同様に移流項の誤差Ｄを修正する項を含めて離散化するようにしてもよい。これにより、該方程式についても誤差による影響を抑えることができるので、優れたシミュレーション結果を得ることができる。
【０１８０】
［第３の実施形態］
以下、図面を参照して、本発明の第３の実施形態について説明する。第１および第２の実施形態では、移流項を含む方程式として、ナビエ・ストークス方程式を挙げ、これを解く流体シミュレーション装置１０、１０ａを説明した。本実施形態では、移流項を含む方程式として連続体の熱エネルギー保存の式を挙げ、これを有限体積法により解く連続体シミュレーション装置１０ｂを説明する。なお、本実施形態の連続体シミュレーション装置１０ｂは、シミュレーションする領域における各時間ステップの速度ベクトルは、外部から入力された値を用いる。この外部から入力された速度ベクトルはシミュレーションにより、予め算出されたものでもよいし、観測データ、実験データであってもよい。
【０１８１】
図１３は、この発明の第３の実施形態による連続体シミュレーション装置１の構成を示す概略ブロック図である。本実施形態における連続体シミュレーション装置１０ｂは、セル構成設定部１０１、セル構成記憶部１０２、速度場設定部１０３ｂ、速度場記憶部１０４ｂ、速度場取得部１０５ｂ、解算出部１０６ｂ、解記憶部１０７ｂを含んで構成される。同図において、図５の各部に対応する部分には、同一の符号（１０１、１０２）を付し、説明を省略する。
【０１８２】
速度場設定部１０３ｂは、外部から入力された各セルの速度ベクトルを、速度場記憶部１０４ｂに記憶させる。なお、本実施形態では、速度ベクトルとして、シミュレーション対象となる全ての時間ステップにおける速度ベクトルが入力される。速度場記憶部１０４ｂは、各セルの速度ベクトルを記憶する。速度場取得部１０５ｂは、解算出部１０６ｂが時間ステップｎに関する解を算出するときは、時間ステップｎ−１における各セルの速度ベクトルを、速度場記憶部１０４ｂから読み出す。
【０１８３】
解算出部１０６ｂは、速度場取得部１０５ｂが速度場記憶部１０４ｂから読み出した時間ステップｎ−１における各セルの速度ベクトルと、セル構成記憶部１０２が記憶するセル情報と、解記憶部１０７ｂが記憶する時間ステップｎ−１における解とを用いて、対象時間ステップｎにおける熱エネルギー保存の式の解を算出する。解算出部１０６ｂは、算出した解を、対象時間ステップｎの解として解記憶部１０７ｂに記憶させる。
【０１８４】
なお、解算出部１０６ｂは、解を算出する際に、式（４−３１）に示すエネルギー保存の式を離散化した離散化方程式を用いる。式（４−３１）に示すように、エネルギー保存の式には、左辺の第２項に移流項がある。密度ρは一定ではないため、解算出部１０６ｂは、第２の実施形態におけるナビエ・ストークス方程式の離散化と同様にして、式（４−３１）を離散化した離散化方程式を用いる。
【０１８５】
【数５９】
2014069421.tif 000061
【０１８６】
上述のように、エネルギー保存の式を離散化する際に、連続の式を加えて移流項を保存型に変形している。このときに、本実施形態では、連続の式が誤差Ｄを有すると仮定し、該誤差Ｄを修正する項を含めて離散化している。よって、外部から入力された各時間ステップの速度が、離散化された連続の式を完全には満たしていなくても、その誤差による影響を抑えることができるので、優れたシミュレーション結果を得ることができる。
【０１８７】
なお、連続体中の物質の移流拡散現象を記述した移流拡散方程式も、式（４−３２）に示すように物質の濃度が移流される移流項を有する。このため、この移流拡散方程式を、エネルギー保存の式と同様に離散化した離散化方程式を解くようにしてもよい。
【０１８８】
【数６０】
2014069421.tif 000062
【０１８９】
また、連続体における応力場を記述した応力場の方程式も、式（４−３３）に示すように連続体の変形速度が移流される移流項を有する。このため、この応力場の方程式を、ナビエ・ストークス方程式と同様に離散化した離散化方程式を解くようにしてもよい。
【０１９０】
【数６１】
2014069421.tif 000063
【０１９１】
なお、連続体の移動・変形を伴う現象を表す方程式には、移流項を含む物質微分が用いられる。該方程式を、上述の各実施形態と同様に離散化した離散化方程式を解くようにしてもよい。
【０１９２】
［第４の実施形態］
以下、図面を参照して、本発明の第４の実施形態について説明する。第４の実施形態では、複数の移流項を含む方程式を有限体積法により解く流体シミュレーション装置１０ｃを説明する。なお、本実施形態では、移流項を持つ方程式を複数同時に解くシミュレーション装置の一例として、ナビエ・ストークス方程式と、熱エネルギー保存の式とを同時に解く流体シミュレーション装置１０ｃを説明する。
【０１９３】
図１４は、この発明の第４の実施形態による流体シミュレーション装置１０ｃの構成を示す概略ブロック図である。本実施形態における流体シミュレーション装置１０ｃは、セル構成設定部１０１、セル構成記憶部１０２、初期条件設定部１０３、速度場記憶部１０４、速度場取得部１０５、解算出部１０６ｃ、解記憶部１０７ｃを含んで構成される。同図において、図５の各部に対応する部分には、同一の符号（１０１〜１０５）を付し、説明を省略する。
【０１９４】
解算出部１０６ｃは、速度場取得部１０５が速度場記憶部１０４から読み出した時間ステップｎ−１における各セルの速度ベクトルと、セル構成記憶部１０２が記憶するセル情報と、解記憶部１０７ｃが記憶する時間ステップｎ−１における解とを用いて、対象時間ステップｎにおけるナビエ・ストークス方程式、連続の式、熱エネルギー保存の式の解を算出する。解算出部１０６ｃは、算出した解を、対象時間ステップｎの解として解記憶部１０７ｃに記憶させる。なお、解算出部１０６ｃは、算出した解のうち速度については、速度場記憶部１０４に記憶させる。
【０１９５】
なお、解算出部１０６ｃは、第１の実施形態におけるナビエ・ストークス方程式および連続の式を離散化した離散化方程式と、第３の実施形態における熱エネルギー保存の式を離散化した離散化方程式を用いる。
なお、本実施形態では、移流項を持つ方程式を複数同時に解くシミュレーション装置の一例として、ナビエ・ストークス方程式と、熱エネルギー保存の式とを同時に解く場合を示したがその他の組み合わせであってもよい。
【０１９６】
本実施形態においても、移流項を持つ方程式各々を離散化する際に、連続の式を加えて移流項を保存型に変形している。このときに、第１から第３の実施形態と同様に連続の式が誤差Ｄを有すると仮定し、該誤差Ｄを修正する項を含めて離散化している。よって、算出した各時間ステップの速度が、離散化された連続の式を完全には満たしていなくても、その誤差による影響を抑えることができるので、優れたシミュレーション結果を得ることができる。
【０１９７】
次に、本発明を用いたシミュレーション結果の例として、第３の実施形態の連続体シミュレーション装置１０ｂにより温度分布を算出した例を説明する。図１５は、シミュレーション対象の速度場を示すコンター図（等高線図）である。すなわち、図１５は、速度場設定部１０３ｂが、速度場記憶部１０４ｂに記憶させる各セルの速度ベクトルの大きさのコンター図である。なお、速度場は、時間変化がないものとして計算を行った。図１５において、横軸はｘ軸、縦軸はｙ軸であり、該速度場の大きさはｘ軸方向が４ｍ、ｙ軸方向が１ｍである。図１５のコンター図は、流体が一様に１ｍ／ｓで左から右へ（ｘ軸の正方向）と流れているが、速度場中に「ｅｒｒ」という形状において１．２ｍ／ｓで左から右へと流れている領域があることを示している。なお、ｙ軸方向の速度は０である。
【０１９８】
図１６は、図１５に示す速度場における速度の発散の大きさを示すコンター図である。図１６においても、横軸はｘ軸、縦軸はｙ軸である。速度の発散は、式（２−３）に示した連続の式の左辺と一致する。すなわち、図１５の速度場が連続の式を完全に満たしていれば、図１６のコンター図は一様に「０」となるはずである。しかし、図１５に示す速度場では、速度が１．２ｍ／ｓの「ｅｒｒ」という形状の領域があるため、該領域の左右端において「０」とはならない。
【０１９９】
図１７は、図１５に示す速度場の温度分布を、従来の方法を用いて算出した場合の例を示すコンター図である。なお、温度分布の初期値は、一様に１００℃としており、図１７に示す温度分布は、定常状態における温度分布である。図１７において、横軸はｘ軸、縦軸はｙ軸であり、等高線は８３℃から１０７℃を１３段階に分割した線である。図１７に示すように、従来の方法にて算出した場合、温度分布の初期値が一様に１００℃であるにも関わらず、「ｅｒｒ」という形状の左右端から下流に、１００℃を超える、あるいは、下回る領域が出来てしまう。これは、「ｅｒｒ」という形状の左右端において、速度の発散が「０」になっていないため、これらの箇所に擬似的な熱のシンク／ソースが発生しているからである。
【０２００】
図１８は、図１５に示す速度場の温度分布を、第３の実施形態の連続体シミュレーション装置１０ｂにて算出した場合の例を示すコンター図である。図１８の例も、温度分布の初期値は、一様に１００℃としており、定常状態における温度分布である。図１８においても、横軸はｘ軸、縦軸はｙ軸であり、等高線は８３℃から１０７℃を１３段階に分割した線である。図１８に示すように、連続体シミュレーション装置１０ｂにより算出された温度分布は、擬似的な熱のシンク／ソースが発生していないため、一様に１００℃となる。一様に１００℃の流体が流れていれば、時間が経過しても一様に１００℃であるので、連続体シミュレーション装置１０ｂにより良好なシミュレーション結果が得られていることがわかる。
【０２０１】
また、図５、図１２、図１３、図１４における各部の機能もしくはその一部の機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することにより流体シミュレーション装置１０、１０ａ、連続体シミュレーション装置１０ｂを実現するようにしてもよい。なお、ここでいう「コンピュータシステム」とは、ＯＳや周辺機器等のハードウェアを含むものとする。
【０２０２】
また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ−ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムを送信する場合の通信線のように、短時間の間、動的にプログラムを保持するもの、その場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリのように、一定時間プログラムを保持しているものも含むものとする。また上記プログラムは、前述した機能の一部を実現するためのものであっても良く、さらに前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるものであっても良い。
【０２０３】
また、本発明のシミュレーション装置は、例えば、特開２０１２−１３２８２４号公報に記載されているような発生源推定装置に用いることができる。この発生源推定装置では、複数個所に設置された観測器の位置情報および計測した濃度情報を入手し、これらの情報に基づき、予め設定されている複数の仮想放出地点の中から放出地点を推定している。このとき、仮想放出地点の一つ一つについて、その仮想放出地点が発生源であり、単位量を放出したときに、各観測器でどのような濃度が測定されるか（該文献における影響関数）を推定し、その推定結果を用いている。この濃度の測定値の推定に、仮想放出地点の設定方法によらず、上述の第３の実施形態で示した移流拡散方程式を解く連続体シミュレーション装置を用いることができる。あるいは、第４の実施形態で示したようにナビエ・ストークス方程式と移流拡散方程式とを同時に解く流体シミュレーション装置を用いることができる。それにより、各観測器で測定される濃度を精度良く推定できるので、発生源も精度よく推定することができる。なお、発生源は、ガス、液体を問わない。
【０２０４】
以上、この発明の実施形態に対して図面を参照して詳述してきたが、具体的な構成はこの実施形態に限られるものではなく、この発明の要旨を逸脱しない範囲の設計変更等も含まれる。
【０２０５】
本出願を詳細にまた特定の実施態様を参照して説明したが、本発明の精神と範囲を逸脱することなく様々な変更や修正を加えることができることは当業者にとって明らかである。
本出願は、２０１２年１０月３１日出願の日本特許出願（特願２０１２-２４１２６０）に基づくものであり、その内容はここに参照として取り込まれる。
【符号の説明】
【０２０６】
１０、１０ａ、１０ｃ…流体シミュレーション装置
１０ｂ…連続体シミュレーション装置
１０１…セル構成設定部
１０２…セル構成記憶部
１０３、１０３ａ…初期条件設定部
１０３ｂ…速度場設定部
１０４、１０４ｂ…速度場記憶部
１０４ａ…密度・速度場記憶部
１０５、１０５ｂ…速度場取得部
１０５ａ…密度・速度場取得部
１０６…速度・圧力算出部
１０６ａ、１０６ｂ、１０６ｃ…解算出部
１０７…圧力場記憶部
１０７ｂ、１０７ｃ…解記憶部
【図１】
2014069421.tif 000064
【図２】
2014069421.tif 000065
【図３】
2014069421.tif 000066
【図４】
2014069421.tif 000067
【図５】
2014069421.tif 000068
【図６】
2014069421.tif 000069
【図７】
2014069421.tif 000070
【図８】
2014069421.tif 000071
【図９】
2014069421.tif 000072
【図１０】
2014069421.tif 000073
【図１１】
2014069421.tif 000074
【図１２】
2014069421.tif 000075
【図１３】
2014069421.tif 000076
【図１４】
2014069421.tif 000077
【図１５】
2014069421.tif 000078
【図１６】
2014069421.tif 000079
【図１７】
2014069421.tif 000080
【図１８】
2014069421.tif 000081
【図１９】
2014069421.tif 000082
【図２０】
2014069421.tif 000083
【国際調査報告】
2014069421.tif 000084
2014069421.tif 000085
2014069421.tif 000086
2014069421.tif 000087