(57)【要約】    【課題】  移動する固体相を含む非圧縮性流体の計算を簡単に行い得る数値計算方法、及びプログラムを提供する。    【解決手段】  数値流体力学における非圧縮性流体の数値計算方法において、固体相を表現するパラメータＨを０≦Ｈ≦１と定義し、該パラメータＨを含む項を圧力ポアソン方程式に導入する。これにより、移動する固体相が流体相に与える効果を当該圧力方程式に反映させ、当該圧力ポアソン方程式を解くことで流体相と固体相を区別せずに一括して計算を行うことができるようにする。【選択図】      図５
2007025818.tif 000002
【特許請求の範囲】
【請求項１】
移動する固体相を含む非圧縮性流体の数値計算方法において、固体相を表現するパラメータを設定する工程と、該パラメータを含んだ圧力方程式を解く工程とを有することを特徴とする数値計算方法。
【請求項２】
前記パラメータは固体相の体積比率であることを特徴とする請求項１に記載の数値計算方法。
【請求項３】
前記パラメータをＨとすると、０≦Ｈ≦１であることを特徴とする請求項１の数値計算方法。
【請求項４】
前記パラメータをＨとすると、０＜Ｈ≦１であることを特徴とする請求項１の数値計算方法。
【請求項５】
前記流体相は、複数の流体成分により構成され、前記圧力方程式は、
当該複数の流体成分のそれぞれに係る前記パラメータの総和を含んでいることを特徴とする請求項１に記載の数値計算方法。
【請求項６】
計算対象の物理空間を或る間隔で離散化するための格子を生成する第１の工程と、前記第１の工程で生成された格子内の単位セルの固体相を移動させる第２の工程と、前記単位セルでの流体相の仮の流速を求める第３の工程と、前記第３の工程で求めた仮の流速を用いて圧力方程式を解いて圧力を求める第４の工程と、前記第４の工程で求められた圧力を用いて前記仮の流速を修正し、真の流速を求める第５の工程とを有する非圧縮性流体の数値計算方法において、
前記第４の工程は、前記単位セル中の固体相の体積比率を表すパラメータを含んだ圧力方程式を解くことを特徴とする数値計算方法。
【請求項７】
請求項１乃至６に記載の数値計算方法を実行するコンピュータ読み取り可能なプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、数値流体力学における非圧縮性流体の数値計算方法、及びプログラムに関する。
【背景技術】
【０００２】
一般的に、或る物理体系に対して数値計算を行う場合、その物理体系を離散格子として差分化している。その物理体系が非圧縮性流体場である場合の数値計算手法としては、流速を段階的に求めるフラクショナルステップ法、ＭＡＣ法系統のＳＭＡＣ（Ｓｉｍｐｌｉｆｉｅｄ  Ｍａｒｋｅｒ  ａｎｄ  Ｃｅｌｌ）法などがよく用いられる。
【０００３】
これらフラクショナルステップ法、ＳＭＡＣ法では、図９に示すように、「まず、現在時刻における流速の移流を行い、仮の流速を求めておく（９ａ）。次に、仮の流速を用いて圧力ポアソン方程式を解き、圧力を求める（９ｂ）。そして、求められた圧力を用いて仮の流速を修正して真の流速を求め、次の時刻における流速を求める（９ｃ）。」という３段階の処理を行って求めるべき流体場を計算する。
【０００４】
以下、この計算手法を具体的に説明するが、ここでは、粘性のない非圧縮流体の場合について説明する。
【０００５】
粘性のない流体の基礎方程式は、質量保存式（連続の式）
【０００６】
【数４６】
2007025818.tif 000003
【０００７】
と、運動方程式（Ｎａｖｉｅｒ−Ｓｔｏｋｅｓ方程式）
【０００８】
【数４７】
2007025818.tif 000004
【０００９】
であるが、非圧縮性流体の場合は、
【００１０】
【数４８】
2007025818.tif 000005
【００１１】
が成り立つため、この数式４８を用いると数式４６と数式４７は、それぞれ、
【００１２】
【数４９】
2007025818.tif 000006
【００１３】
【数５０】
2007025818.tif 000007
【００１４】
となる。ここで、uは流速、pは圧力、ρは密度であり、この数式４９と数式５０が粘性のない非圧縮性流体の基礎方程式である。
【００１５】
この基礎方程式に基づいて流体場を計算する場合、フラクショナルステップ法やＳＭＡＣ法では、数式４９の流速と数式５０の圧力については、陰的に評価、すなわち、新しい時刻（次の時刻）で評価する。一方、数式５０の流速については、陽的に評価、すなわち、現在の時刻で評価する。この評価法則に則って数式４９と数式５０を離散化すると、数式４９、数式５０は、それぞれ、
【００１６】
【数５１】
2007025818.tif 000008
【００１７】
【数５２】
2007025818.tif 000009
【００１８】
となる。
さらに、数式５２を
【００１９】
【数５３】
2007025818.tif 000010
【００２０】
と、
【００２１】
【数５４】
2007025818.tif 000011
【００２２】
という２式に分解する。ここで、添字「*」は仮の値であることを示し、u*は仮の流速を示している。また、添字「ｎ」は、現在時刻での値、添字「ｎ＋１」は、新しい時刻での値であることを示している（以下、同様）。
【００２３】
数式５１と数式５４から新しい時刻における流速un+1を消去すると、結局、以下の圧力ポアソン方程式
【００２４】
【数５５】
2007025818.tif 000012
【００２５】
が得られる。この数式５５の圧力ポアソン方程式を解いて新しい時刻における圧力pn+1が求められれば、数式５４より、新しい時刻における流速un+1が求まることとなる。このようなフラクショナルステップ法、ＳＭＡＣ法における数値解法のアルゴリズムを図１０に示す。
【００２６】
ところで、計算対象の物理体系に移動を伴う何らかの固体相（障害物等）が存在する場合は、固体相を除く流体相の領域に対してのみ上記の計算を行い、流体相と固体相の境界に対しては或る速度を別に指定する方法が一般的である。すなわち、固体相が或る速度で移動する場合には、その固体相の移動に応じた速度を境界速度として設定することになる。つまり、流体相が満たすべき条件は、
【００２７】
【数５６】
2007025818.tif 000013
【００２８】
と、
【００２９】
【数５７】
2007025818.tif 000014
【００３０】
である。換言すれば、流体相では、数式５６を適用し、固体相と流体相の境界部分では、数式５７を適用するというように、場合分けをして計算を行なう。
【００３１】
上記の従来法を用いた問題例として図１１に示す物理体系を想定する。この物理体系では、２次元空間に固体相（壁）１１ａ，１１ｂが存在し、それら固体相１１ａ，１１ｂの間に流体相１１ｃが存在している。また、流体相１１ｃの左右の境界は自由境界面１１ｅとし、同種の流体が自由に流入・流出できるものとする。また、壁１１ｂは固定され、壁１１ａは時間と共に一定の速度で下方向に移動するものとする。この物理体系では、結局、時間と共に、壁１１ａ，１１ｂよって流体相１１ｃの一部が左右の自由境界面１１ｅから流出することとなる。
【００３２】
このような問題を従来法で計算しようとすると、図１２のように、流体相１２bと壁１２aとの境界１２cに対し、壁１２aの移動速度１２dと同じ速度を各時間ステップで設定することになる。1ステップの計算が終わると、計算対象の物理体系は、図１３に示すように、壁１３aの移動によって前の時刻の境界１３dが境界１３cに移動し、この移動に伴って流体相１３bの一部は計算対象の物理体系の外部に流出した状態となる。
【００３３】
この他の方法としては、Ｐａｎの研究（非特許文献１参照）や、梶島らの研究（非特許文献２参照）がある。これらの研究では、固体相が流体に与える影響を上記の数式４９のＮａｖｉｅｒ−Ｓｔｏｋｅｓ方程式の外力項として取り込んで計算を行っている。
【非特許文献１】Ｙ．Ｐａｎ，Ｓ．Ｂａｎｅｒｊｅｅ，Ｐｈｙｓ．Ｆｌｕｉｄｓ  ９（１２），３７８６−３８０７（１９９７）
【非特許文献２】梶島、瀧口、浜崎、三宅、日本機械学会論文集（Ｂ編）、６６巻６４７号、１７３４−１７４１（２０００）
【発明の開示】
【発明が解決しようとする課題】
【００３４】
しかしながら、従来法では、上記の通り、流体相と固体相に分けて計算を行う必要があり、さらに、流体相と固体相の境界に対しては或る流速を別に設定するという特別な処理を行う必要がある。
【００３５】
また、流体境界に対して或る流速を設定することは、多くの場合において簡単ではなく、固体相がより複雑な形状を持つ場合は、流速指定の処理はより困難なものとなる。例えば図１１〜図１３に示した例では、移動固体相(上部の壁１２ａ）は、鉛直方向に移動速度１２dで移動するため、流体境界に設定すべき速度の値は単純に移動速度１２ｄと等しい値で良かった。しかし、これは非常に単純化した例であり、より一般的な問題では上記の例よりもはるかに複雑な挙動を扱う。例えば図１４のように、プロペラ１４ａが流体相１４ｂの中を或る速度で回転するような問題を考えた場合、プロペラ１４ａと流体相１４ｂの境界を含むセルにおいて、流体境界に速度を設定する処理は、非常に複雑になる。
【００３６】
また、従来法は、流体相が１相のみで構成されている場合を想定したものが多く、工学上重要な多相流への適用に関しては課題や未知な点が多い。
【００３７】
本発明は、このような背景の下になされたもので、移動する固体相を含む非圧縮性流体の計算を簡単に行い得る数値計算方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【００３８】
上記目的を達成するため、本発明は、移動する固体相を含む非圧縮性流体の数値計算方法において、固体相を表現するパラメータを設定する工程と、該パラメータを含んだ圧力方程式を解く工程とを有することを特徴とする。
【発明の効果】
【００３９】
本発明によれば、非圧縮性流体に関して、固体相を表現するパラメータを含む項を圧力方程式に導入する。これにより、移動する固体相が流体相に与える効果を当該圧力方程式に反映させているので、流体相と固体相を区別せずに一括して計算を行うことができる。よって、移動する固体相を含む非圧縮性流体の計算を簡単に行い得る数値計算手法、及びプログラムを提供することが可能となる。
【発明を実施するための最良の形態】
【００４０】
本発明の実施の形態による数値計算方法は、図１に示すように、入力部１ａ、ＣＰＵ１ｂ、出力部１ｃ、ＲＯＭ１ｄ、ＲＡＭ１ｅ、補助記憶装置１ｆを有する一般的な計算機上で実現されるものである。ここで、入力部１ａは各種データやパラメータ、命令などを入力し、ＣＰＵ１ｂは、プログラムを実行し、各部を制御する。また、出力部１ｃは、実行結果やメッセージなどを出力する。ＲＯＭ１ｄは、固定データやプログラムを記憶しており、ＲＡＭ１ｅは、ワークメモリとして利用され、入力データや中間結果を記憶する。また補助記憶装置１ｆには、後述する第１〜第４の実施の形態に係る数値計算方法を用いて流体のシュミレーション等を行なうＪａｖａ（登録商標）プログラム等のアプリケーションプログラムや実行結果のデータなどが不揮発に格納されている。
【００４１】
本発明の実施の形態による数値計算方法では、計算対象の物理体系の空間を、例えば図２に示すように離散化するための格子２aを生成し、複数のセルを形成する。この物理体系には、固体相２bと流体相２cが存在する。なお、流体相２cは1相のみで構成されるものとする。このとき、固体相２bと流体相２cの境界が含まれる単位セル２d（もしくはその拡大図２e）に着目した際、単位セル２d（又は２e）中の流体相２fの体積比率をα（０≦α≦１）とする。
【００４２】
このとき、単位セル２ｄ（又は２e）における固体相２gの体積比率は１−αであるが、固体相の形状を表現するパラメータＨは、Ｈ＝α、すなわち０≦Ｈ≦１と定義と定義する。ただし、このパラメータＨは、単位セル中の固体相が１００％になった場合はＨ＝０、０％になった場合はＨ＝１となるように設定される。
【００４３】
ここで、図３のように、前述の図１１と同じような問題を考えることにする。図３のように、固体相としての壁３ａ，３ｂが物理体系に存在し、それら壁３ａ，３ｂの間に流体相３ｃが存在するものとする。格子３ｅに壁３ａと流体相３ｃの境界３dが含まれているとき、格子３ｅの中の或る単位セル３ｆに注目してみる。このとき、３ｇに示す流体相の体積比率（斜線部）をαと定義する。上部の壁３ａが時間と共に下向きに移動すると考えると、その現象は図４に示すように、流体相の体積比率αの時間変化として表現されるが、先にパラメータＨをＨ＝αと定義してあるため、この現象はパラメータＨの時間変化としても表現できる。
【００４４】
単位セル３ｆの状態は、最初は、４ａのように（Ｈ＝０．８５）であったとする。或る時間ステップ（Δｔ）経過した後は、単位セル３ｆの状態は、４ｂのように変化する。このとき、パラメータＨの値は、０．８５から０．５５に変化する。さらにΔｔ時間経過後には、単位セル３ｆの状態は、４ｃのように変化する。このときパラメータＨの値は、０．５５から０．３５に変化する。
【００４５】
以下では、上記のように、単位セル中のパラメータＨの変化で固体相の移動現象を表現する数値解法を示すことにする。
【００４６】
[第１の実施の形態]
上記のような定義を行った物理体系における基礎方程式は、質量保存式
【００４７】
【数１】
2007025818.tif 000015
【００４８】
と、運動量保存式
【００４９】
【数２】
2007025818.tif 000016
【００５０】
である。ここで、流体相は１相のみで成り立っているとする。
流体相が非圧縮性流体の場合には、数式１と数式２は、それぞれ
【００５１】
【数３】
2007025818.tif 000017
【００５２】
【数４】
2007025818.tif 000018
【００５３】
と変形することができる。ここで、数式３の流速uを陰的に評価して離散化すると、
【００５４】
【数５】
2007025818.tif 000019
【００５５】
となる。ここで、新しい時刻においてはＨn+1＝αn+1が成り立たなければならないので、数式５は、
【００５６】
【数６】
2007025818.tif 000020
【００５７】
となる。また、数式４の第２項の移流項を陽的に評価し、圧力を陰的に評価して
【００５８】
【数７】
2007025818.tif 000021
【００５９】
【数８】
2007025818.tif 000022
【００６０】
と２式に分解して離散化する。
【００６１】
数式６と数式８から新しい時刻における流速un+1を消去して整理すると、結局、圧力ポアソン方程式
【００６２】
【数９】
2007025818.tif 000023
【００６３】
を得ることができる。数式９の圧力ポアソン方程式において、仮の流速u*は、数式７より求められ、パラメータＨn+1は外部より与えられる（固体相の移動法則は外部から設定される）。このため、数式９の圧力ポアソン方程式を数値的に解くことにより、新しい時刻における圧力pn+1を算出することができる。この新しい時刻における圧力pn+1を数式８に代入することで、新しい時刻における流速un+1を算出することができる。また、その新しい時刻における流速un+1を数式５に代入することで、新しい時刻における流体相の体積比率αn+1を算出することができる。
【００６４】
ここで、数式９の圧力ポアソン方程式に着目すると、右辺第２項の(Ｈn+1−αn)/（Δt）は、移動する固体相が流体相に及ぼす影響を表す項である。そこで、この項を数式９のように圧力ポアソン方程式に導入することで、従来の流体境界への複雑な速度指定を行う必要がなくなり、移動する固体相の影響を考慮した数値計算をより簡単に行うことが可能となる。
【００６５】
すなわち、従来は、移動固体相の影響を考慮するために、流体相では数式５６、固体相と流体相の境界部分では数式５７を適用するというように場合分けをして計算を行っていたが、この方法は多くの場合処理が複雑になりがちであった。
【００６６】
これに対し、第１の実施の形態では、数式９のように、固体相を表現するパラメータＨを移動固体相が流体相に与える影響を示す項として圧力方程式に導入している。このため、この圧力方程式を解くことで、流体境界への速度指定と言う特殊かつ複雑な処理を行わずに流体相と固体相の計算を一括して行うことができるようになる。これにより、流体相と固体相の区別をする必要がなくなり、移動する固体相を含む流体計算を従来法よりも簡単に行なうことが可能となる。
【００６７】
なお、上記の実施の形態において、計算対象の物理体系によっては、固体相と流体相の境界において、想定される流速よりも非常に大きい流速値が算出されてしまう場合があることが判っている。その場合は、完全な固体相領域においてもＨ＝０と定義せず、微小な正の実数ε（ε＜＜１）を導入し、Ｈ＝ε、すなわち０＜Ｈ≦１と定義して計算を行うことにより、固体相と流体相境界においても適切な流速を算出できるようになる。
【００６８】
このような第１の実施の形態における数値計算のアルゴリズムを、図５に示す。
【００６９】
[第２の実施の形態]
次に、流体相がｎ相の多成分から成り立っている場合の数値解法を説明する。今、流体相の中の成分ｃ(ｃ＝1,2,…,n)に着目すると、成分ｃに対する非圧縮性流体の基礎方程式は、
【００７０】
【数１０】
2007025818.tif 000024
【００７１】
【数１１】
2007025818.tif 000025
【００７２】
となる。なお、αc、uc、ρcはそれぞれ成分cの体積比率、流速、密度である。また、圧力pは全ての成分において等しいと仮定した。
【００７３】
ここで、数式１０の成分ｃの流速ucを陰的に評価して離散化すると、
【００７４】
【数１２】
2007025818.tif 000026
【００７５】
となる。ここで、数式１２を全ての成分に対して考え、それらの和をとると、
【００７６】
【数１３】
2007025818.tif 000027
【００７７】
となる。
【００７８】
新しい時刻においては、
【００７９】
【数１４】
2007025818.tif 000028
【００８０】
が成り立たなければならないので、数式１３は、
【００８１】
【数１５】
2007025818.tif 000029
【００８２】
となる。また、数式１１の第２項の移流項を陽的に評価し、圧力を陰的に評価して、
【００８３】
【数１６】
2007025818.tif 000030
【００８４】
【数１７】
2007025818.tif 000031
【００８５】
と２式に分解して離散化する。
【００８６】
数式１５と数式１７から下記の数式１８で示される新しい時刻における流速
【００８７】
【数１８】
2007025818.tif 000032
【００８８】
を消去して整理すると、結局、圧力ポアソン方程式
【００８９】
【数１９】
2007025818.tif 000033
【００９０】
を得ることができる。この数式１９の圧力ポアソン方程式において、下記の数式２０で示される仮の流速
【００９１】
【数２０】
2007025818.tif 000034
【００９２】
は、数式１６より求められ、新しい時刻におけるパラメータHn+1は外部より与えられる（固体相の移動法則は外部から設定される）ため、数式１９の圧力ポアソン方程式を数値的に解くことにより、新しい時刻における圧力pn+1を算出することができる。この新しい時刻における圧力pn+1を数式１７に代入することで、上記の数式１８で示される新しい時刻における流速を算出することができ、また、その新しい時刻における流速を数式１２に代入することで、下記の数式２１で示される新しい時刻における流体相の体積比率
【００９３】
【数２１】
2007025818.tif 000035
【００９４】
を算出することができる。
【００９５】
なお、上記実施例において、計算対象の物理体系によっては、固体相と流体相の境界において想定される流速よりも非常に大きい流速値が算出されてしまう場合があることが判っている。その場合は、完全な固体相領域においてもＨ＝０と定義せず、微小な正の実数ε（ε＜＜１）を導入し、Ｈ＝ε、すなわち０＜Ｈ≦１と定義して計算を行うことにより、固体相と流体相境界においても適切な流速を算出できるようになる。第２の実施の形態における数値計算のアルゴリズムを、図６に示す。
【００９６】
このように、流体相が多成分より成る物理体系においても、１成分の流体相における数値解法を自然に拡張する形で計算を行うことができる。
【００９７】
すなわち、流体相が多成分から成る場合、流体境界に対して速度を指定する従来法では、固体相が多成分流体相へ及ぼす影響を考慮して計算する方法が明確ではなかった。しかし、第２の実施の形態では、流体相が多成分から成る場合においても１成分からの自然な拡張による数式１９の圧力ポアソン方程式を解くことで、流体相と固体相の計算を一括して行うことができるようになる。これにより、流体相と固体相の区別をする必要がなくなり、固体相が多成分流体相へ及ぼす影響を考慮した形で、移動する固体相を含む多成分流体計算を従来法よりも簡単に行なうことが可能となる。
【００９８】
[第３の実施の形態]
第１,第２の実施の形態と同様に、第３,第４の実施の形態に係る数値計算方法も、図１に示したような構成の計算機上で実現されるものである。また、第１,第２の実施の形態と同様に、第３,第４の実施の形態に係る数値計算方法でも、例えば図２に示すように計算対象の物理体系の空間を離散化して格子２aを生成する。さらに、第１,第２の実施の形態と同様に、第３,第４の実施の形態に係る数値計算方法でも、図３、図４に示した２次元スクイーズ問題を考えている。
【００９９】
ただし、第１,第２の実施の形態では、流体相の体積比率を「α」で表記していたのに対し、第３,第４の実施の形態では、流体相の体積比率を「ｈ」で表記している点で相違する（図３参照）。また、第１,第２の実施の形態では、固体相の形状を表現するパラメータの記号を「Ｈ」としていたが、第３，第４の実施の形態では、固体相の形状を表現するパラメータの記号を「ｈ」としている点で相違する（図４参照）。なお、第３,第４の実施の形態においても第１,第２の実施の形態と同様に、固体相の形状を表現するパラメータｈを０≦ｈ≦１と定義する。このパラメータｈは、単位セル中の固体相が１００％になった場合はｈ＝０、０％になった場合はｈ＝１となるように設定される。
【０１００】
従って、図３に示した上部の壁３aが時間と共に下向きに移動する場合を考えると、その現象は、パラメータｈの時間変化として表現することができ、図４に示した例では、単位セルにおけるパラメータｈは、「０．８５→０．５５→０．３５」と変化している。
【０１０１】
以下、同様に、単位セル中のパラメータｈの変化で固体相の移動現象を表現することにする。
【０１０２】
上記のような表記・定義を行った物理体系における基礎方程式は、質量保存式
【０１０３】
【数２２】
2007025818.tif 000036
【０１０４】
と、運動量保存式
【０１０５】
【数２３】
2007025818.tif 000037
【０１０６】
である。
【０１０７】
以下、本第３の実施の形態では、「非保存形式」による記述に基づいた数値解法を示す。
【０１０８】
流体相が非圧縮性流体の場合には、数式２２の質量保存式を下記の数式２４のように変形することができる。
【０１０９】
【数２４】
2007025818.tif 000038
【０１１０】
また、流体相が非圧縮性流体の場合には、数式２３の運動量保存式も下記の数式２５のように変形することができる。
【０１１１】
【数２５】
2007025818.tif 000039
【０１１２】
ここで、数式２４、数式２５の左辺第１項に出てくる微分（Ｄ/Ｄt）はLagrange微分と呼ばれるものであり、
【０１１３】
【数２６】
2007025818.tif 000040
【０１１４】
である。
ここで、数式２４を離散化すると、
【０１１５】
【数２７】
2007025818.tif 000041
【０１１６】
となる。ただし、パラメータｈは、
【０１１７】
【数２８】
2007025818.tif 000042
【０１１８】
と定義する。ここで、前述のように、添字「*」は仮の値であることを示し、添字「ｎ」は、現在時刻に係る値、添字「ｎ＋１」は、新しい時刻に係る値であることを示している。すなわち、数式２７の流速un+1は、陰的に新しい時刻で評価したものである。
【０１１９】
また、数式２５を離散化すると、
【０１２０】
【数２９】
2007025818.tif 000043
【０１２１】
となる（圧力は新しい時刻で評価）。ただし、仮の流速u*は、
【０１２２】
【数３０】
2007025818.tif 000044
【０１２３】
と定義する。
【０１２４】
数式２７と数式２９から新しい時刻におけるun+1を消去して整理すると、結局、圧力ポアソン方程式
【０１２５】
【数３１】
2007025818.tif 000045
【０１２６】
を得ることができる。この数式３１の圧力ポアソン方程式において、h*とu*はそれぞれ数式２８と数式３０より求められ、新しい時刻におけるパラメータhn+1は外部より与えられる（固体相の移動法則は外部から設定される）。このため、数式３１の圧力ポアソン方程式を数値的に解くことにより、新しい時刻における圧力pn+1を算出することができる。この新しい時刻における圧力pn+1を数式２９に代入することで、最終的に求めたい新しい時刻における流速un+1を算出することができる。
【０１２７】
ここで数式３１の圧力ポアソン方程式に着目すると、右辺第２項の（hn+1−h*）/（h*Δt2）は、移動固体相が流体相に及ぼす影響を表す項である。この項を数式３１のように圧力ポアソン方程式に導入することで、従来の流体境界への複雑な速度指定を行う必要がなくなり、移動固体相が流体相に及ぼす影響を考慮した数値計算をより簡単に行うことが可能となる。第３の実施の形態における数値計算のアルゴリズムを、図７に示す。
【０１２８】
[第４の実施の形態]
次に、第３の実施の形態の変形例として、保存形式による記述に基づいた数値解法を示す。
【０１２９】
数式２２に対して、Ｍ＝hρと定義すると、この数式２２は、
【０１３０】
【数３２】
2007025818.tif 000046
【０１３１】
と書くことができる。また、数式２３に対して、Ｊ＝Ｍu＝hρuと定義すると、この数式２３は、
【０１３２】
【数３３】
2007025818.tif 000047
【０１３３】
と書くことができる。
【０１３４】
ここで、数式３２を離散化して整理すると、
【０１３５】
【数３４】
2007025818.tif 000048
【０１３６】
となる。ここで、添字の意味は前述の通りであり、数式３４における流速un+1は、陽的に現在の時刻で評価したものである。
【０１３７】
また、数式３３に対して離散化を行うと、
【０１３８】
【数３５】
2007025818.tif 000049
【０１３９】
となる。なお、数式３５において、圧力は新しい時刻で評価している。この数式３５は、
【０１４０】
【数３６】
2007025818.tif 000050
【０１４１】
と定義すると、
【０１４２】
【数３７】
2007025818.tif 000051
【０１４３】
【数３８】
2007025818.tif 000052
【０１４４】
と変形することができる。
【０１４５】
ここで、Ｊ＝Ｍuであることと、数式３４のＭn+1はＭnを流速unで保存移流させたものである、すなわち
【０１４６】
【数３９】
2007025818.tif 000053
【０１４７】
であることを考えると、数式３６と数式３８は、それぞれ
【０１４８】
【数４０】
2007025818.tif 000054
【０１４９】
【数４１】
2007025818.tif 000055
【０１５０】
と書くことができる。
【０１５１】
また、流体が非圧縮性であることにより、数式２２は、
【０１５２】
【数４２】
2007025818.tif 000056
【０１５３】
となるから、これを離散化すると
【０１５４】
【数４３】
2007025818.tif 000057
【０１５５】
が得られる（流速は新しい時刻で評価）。
【０１５６】
数式４１と数式４３から新しい時刻における流速un+1を消去して整理すると、結局、圧力ポアソン方程式
【０１５７】
【数４４】
2007025818.tif 000058
【０１５８】
を得ることができる。
【０１５９】
この数式４４の圧力ポアソン方程式において、Ｍn+1、
【０１６０】
【数４５】
2007025818.tif 000059
【０１６１】
は、それぞれ、数式３４、数式４０より求められ、hn+1は外部より与えられるため、数式４４の圧力ポアソン方程式を数値的に解くことにより、新しい時刻における圧力pn+1を算出することができる。この新しい時刻における圧力pn+1を数式４１に代入することで、最終的に求めたい新しい時刻における流速un+1を算出することができる。
【０１６２】
ここで、数式４４の圧力ポアソン方程式には、各々の項に固体相の体積比率ｈの効果が含まれているため、従来の流体境界への複雑な速度指定を行う必要がなくなり、移動固体相の影響を考慮した数値計算をより簡単に行うことが可能となる。この第４の実施の形態における数値計算のアルゴリズムを、図８に示す。
【０１６３】
このように、第３,第４の実施の形態では、計算対象の物理空間を離散化した場合の各単位セルに含まれる固体相の体積比率を０以上１以下の実数で表す。そして、各単位セルに含まれる固体相が移動する際に流体相に及ぼす影響を示す項として、当該固体相の体積比率を含んだ項を圧力ポアソン方程式に取り込むことにより、移動固体相を含む流体計算を従来よりも簡単かつ柔軟に行うことができる。
【０１６４】
すなわち、従来は移動固体相の影響を考慮するために、流体部分では数式５６、固体相と流体相の境界部分では数式５７を適用するというように、場合分けをして計算を行っていたが、この方法は処理が複雑になりがちであった。
【０１６５】
しかし、第３、又は第４の実施の形態では、数式３１、又は数式４４の圧力ポアソン方程式を解くことで、流体境界への速度指定と言う特殊かつ複雑な処理や数式適用の場合分けを行わずに流体相と固体相の計算を一括して行うことができるようになる。これにより、流体相と固体相の区別をする必要がなくなり、従来法よりも移動固体相を含む流体計算を簡単に行なうことが可能となる。
【０１６６】
なお、本発明の目的は、実施形態の機能を実現するソフトウェアのプログラムコードを記録した記憶媒体を、システム或いは装置に供給し、そのシステム或いは装置のコンピュータ（またはＣＰＵやＭＰＵ等）が記憶媒体に格納されたプログラムコードを読み出して実行することによっても達成される。
【０１６７】
この場合、記憶媒体から読み出されたプログラムコード自体が前述した実施の形態の機能を実現することになり、そのプログラムコード及び該プログラムコードを記憶した記憶媒体は本発明を構成することになる。
【０１６８】
又、プログラムコードを供給するための記憶媒体としては、例えば、フロッピー（登録商標）ディスク、ハードディスク、光磁気ディスク、ＣＤ−ＲＯＭ、ＣＤ−Ｒ、ＣＤ−ＲＷ、ＤＶＤ−ＲＯＭ、ＤＶＤ−ＲＡＭ、ＤＶＤ−ＲＷ、ＤＶＤ＋ＲＷ、磁気テープ、不揮発性のメモリカード、ＲＯＭ等を用いることができる。または、プログラムコードをネットワークを介してダウンロードしてもよい。
【０１６９】
また、コンピュータが読み出したプログラムコードを実行することにより、上記実施の形態の機能が実現されるだけでなく、そのプログラムコードの指示に基づき、コンピュータ上で稼動しているＯＳ（オペレーティングシステム）等が実際の処理の一部または全部を行い、その処理によって前述した実施形態の機能が実現される場合も含まれる。
【０１７０】
更に、記憶媒体から読み出されたプログラムコードが、コンピュータに挿入された機能拡張ボードやコンピュータに接続された機能拡張ユニットに備わるメモリに書き込まれた後、そのプログラムコードの指示に基づき、その機能拡張ボードや機能拡張ユニットに備わるＣＰＵ等が実際の処理の一部または全部を行い、その処理によって前述した実施形態の機能が実現される場合も含まれる。
【図面の簡単な説明】
【０１７１】
【図１】本発明の実施の形態による数値計算方法を実現する計算機の概略構成を示すブロック図である。
【図２】流体相と固体相を含んだ離散格子（単位セル）を示す図である。
【図３】２次元スクイーズ問題を例示した図である（パラメータ説明用）。
【図４】単位セル中の固体相の体積比率を示すパラメータの時間変化を説明するための図である。
【図５】本発明の第１の実施の形態に係る数値解法のアルゴリズムを示すフローチャートである。
【図６】本発明の第２の実施の形態に係る数値解法のアルゴリズムを示すフローチャートである。
【図７】本発明の第３の実施の形態に係る数値解法のアルゴリズムを示すフローチャートである。
【図８】本発明の第４の実施の形態に係る数値解法のアルゴリズムを示すフローチャートである。
【図９】従来のフラクショナルステップ法、ＳＭＡＣ法の基本的な解法を示すフローチャートである。
【図１０】従来のフラクショナルステップ法、ＳＭＡＣ法に係る数値解法のアルゴリズムを示すフローチャートである。
【図１１】２次元スクイーズ問題例を示す図である（初期状態）。
【図１２】上記２次元スクイーズ問題例における従来の流体境界への速度指定を説明するための図である。
【図１３】上記２次元スクイーズ問題例での固体相移動後の状態を示す図である。
【図１４】より一般的で複雑な問題例を示す図である。
【符号の説明】
【０１７２】
１ａ…入力部
１ｂ…ＣＰＵ
１ｄ…ＲＯＭ
1ｅ…ＲＡＭ
１ｆ…補助記憶装置
Ｈ，ｈ…パラメータ（固体相の体積比率）
【図１】
2007025818.tif 000060
【図２】
2007025818.tif 000061
【図３】
2007025818.tif 000062
【図４】
2007025818.tif 000063
【図５】
2007025818.tif 000064
【図６】
2007025818.tif 000065
【図７】
2007025818.tif 000066
【図８】
2007025818.tif 000067
【図９】
2007025818.tif 000068
【図１０】
2007025818.tif 000069
【図１１】
2007025818.tif 000070
【図１２】
2007025818.tif 000071
【図１３】
2007025818.tif 000072
【図１４】
2007025818.tif 000073