JP4304277B2 - 数値計算方法、プログラムおよび記録媒体 - Google Patents

Description

本発明は、偏微分方程式の演算に使用する数値計算方法に関する。

近年、コンピュータを用いて自然現象を再現する取り組みが、多く行われるようになってきている。たとえば、流体力学などの数値シミュレーションがその例で、その対象物は、液体、気体、電磁波など様々である。
そして、そのような数値シミュレーションを行う場合には、２以上の変数（以下、「多変数」という）に関する偏微分方程式を、離散化された空間および時間に関して解く、という手法が用いられることが多い。なお、変数には、たとえば、速度、圧力、温度、密度などが挙げられる。

多変数偏微分方程式における数値解法（シミュレーション）では、変数間のカップリングについて、高い時間的安定（以下、単に「安定」という）性が要求される。変数間のカップリングの安定性が低いと、短時間で計算が破綻してしまうからである。
そこで、従来、流体解析におけるＳＭＡＣ（Simplified Marker and Cell）法や電磁波解析におけるＦＤＴＤ（Finite Difference Time Domain）法においては、空間を離散化するときに、変数ごとの定義点をずらして設定するスタッガード格子が多く用いられていた。スタッガード格子は、変数間のカップリングの安定性が高いためである。

また、それぞれの変数の値だけでなく、それらの空間勾配（微分係数）なども利用する計算手法であるマルチモーメントスキームにおいても、スタッガード格子が多く用いられてきた。しかし、スタッガード格子は、変数間のカップリングの安定性が高い反面、計算精度や並列計算アルゴリズムの容易さなどの点で、欠点があった。

一方、すべての変数や空間勾配を同じ格子上に定義するコロケート格子は、スタッガード格子に比べて、計算精度や並列計算アルゴリズムの容易さなどの点で優れている。
たとえば、マルチモーメントスキームにおける計算精度は、スタッガード格子の２次に対し、コロケート格子は４次である。
また、コロケート格子は、スタッガード格子よりも、今後の発展が期待される計算手法であるＡＭＲ（Adaptive Mesh Refinement：適合格子細分化）法やCut-Cell法との組み合わせが容易であるという利点も有する。

しかしながら、コロケート格子は、変数間のカップリングの安定性が低く、計算が短時間で破綻してしまうことが多いため、特殊な条件下でしか使用することができず、適用範囲が非常に限定されてしまう、という問題があった。

たとえば、マルチモーメントスキームであるＩＤＯ（Interpolated Differential Operator：局所補間微分オペレータ）法では、間隔がΔｘである３点、ｊ−１点、ｊ点、ｊ＋１点に関する中心補間関数Ｆ（ｘ）を、ｊ点を中心として、次の式（１）のような５次関数として表わすことができる。
なお、ｆ_j-1,ｆ_j,ｆ_j+1は、それぞれ、ｊ−１点、ｊ点、ｊ＋１点における変数の値を表わし、ｆ_x,j-1, ｆ_x,j, ｆ_x,j+1は、それぞれ、ｊ−１点、ｊ点、ｊ＋１点における変数の微分係数を表わす。また、ｘは座標を表わす。

Ｆ（ｘ）＝Ａ₅ｘ⁵＋Ａ₄ｘ⁴＋Ａ₃ｘ³＋Ａ₂ｘ²＋ｆ_x,jｘ＋ｆ_j ・・・（１）
また、Ｆ（ｘ）の微分式であるＦ_x（ｘ）は、式（２）のように表わされる。
Ｆ_x（ｘ）＝５Ａ₅ｘ⁴＋４Ａ₄ｘ³＋３Ａ₃ｘ²＋２Ａ₂ｘ＋ｆ_x,j ・・・（２）
拘束条件は、Ｆ（Δｘ）＝ｆ_j+1,Ｆ（−Δｘ）＝ｆ_j-1,Ｆ_x（Δｘ）＝ｆ_x,j+1, Ｆ_x（−Δｘ）＝ｆ_x,j-1であり、これらにより、Ａ₂〜Ａ₅の値が決まる。

そして、この条件下で、支配方程式（演算に使う偏微分方程式など）中の空間微分項は、式（１）、式（２）、および、式（２）をさらに微分した式などにより、∂ｆ/∂ｘ＝Ｆ_x（０）＝ｆ_x,j,∂²ｆ/∂ｘ²＝Ｆ_xx（０）＝２Ａ₂,・・・のように求めることができる。なお、Ｆ_xx（ｘ）は、Ｆ_x（ｘ）を微分した式を表わす。

ここで、多変数型の偏微分方程式である次の式（３）、式（４）を考える。なお、ｆとｕは変数、ｔは時間、ｘは座標を表わす。

コロケート格子上では、前記補間関数（式（１）および式（２））を用いて、次の式（５）、式（６）のような空間離散化がなされる。なお、∂ｆ/∂ｔ｜_j,∂ｕ/∂ｔ｜_jは、それぞれ、補間関数などから求めたｊ点での微分係数を表わす。
そして、式（５）、式（６）からわかるように、離散化式中にそれぞれｊ点の変数しか現われないため、変数間のカップリングの安定性が悪い（低い）。

そこで、本発明は、前記問題点に鑑みてなされたものであり、コロケート格子を用いた偏微分方程式の数値解法において、変数間のカップリングの安定性が高い数値計算方法を提供することを目的とする。

前記課題を解決するために、本発明に係る数値計算方法は、各定義点上で変数に関する数値計算を行う数値計算装置による数値計算方法であって、前記数値計算装置は、記憶部と処理部を有し、前記記憶部は、前記各定義点のうち３点以上における各変数の値と、当該３点以上の各変数の値にそれぞれ対応する第１の微分係数と、前記各変数の値と前記３点以上のうち特定点を含まない定義点における前記第１の微分係数により決定される補間関数と、を記憶し、前記処理部は、前記補間関数によって前記特定点における第２の微分係数を算出し、前記第２の微分係数と前記特定点における第１の微分係数との加重平均によって前記特定点における第３の微分係数を算出することを特徴とする。
また、当該数値計算方法をコンピュータに実行させるためのプログラムを作成し、さらに、そのプログラムを記録媒体に記録することができる。

本発明によれば、コロケート格子を用いた偏微分方程式の数値解法において、変数間のカップリングの安定性を高めることができる。

以下、本発明の実施形態に係る数値計算方法について、適宜図面を参照しながら説明する。
まず、図１を参照しながら、数値計算方法を行う数値計算装置の構成について説明する。図１は、数値計算装置の構成図である。

数値計算装置１０は、入力部１、記憶部（記録媒体）２、メモリ３、処理部４および出力部５を備えて構成されるコンピュータ装置であり、たとえば、パソコン、ワークステーション、スーパーコンピュータなどにより実現することができる。
なお、処理部４は、ここでは１つとしているが、複数設けることにより並列計算を行うようにしてもよい。

入力部１は、数値計算装置１０の使用者（以下、単に「使用者」という）が各種情報を入力する手段であり、たとえば、キーボード、マウスなどである。
記憶部２は、各種情報を記憶するものであり、たとえば、ＲＯＭ（Read Only Memory）、ハードディスクなどである。記憶部２は、ここでは、流体力学の数値シミュレーションを行うために、数値計算を行う定義点である格子点、物理量（変数、微分係数など）の初期条件、支配方程式（自然現象を表わす偏微分方程式など）、補間関数、各種プログラムなどを記憶している。

メモリ３は、処理部４が使用する一時記憶手段であり、たとえば、ＲＡＭ（Random Access Memory）などである。
処理部４は、記憶部２に記憶された各種情報（プログラム、データなど）などに基づいて演算処理を行うものであり、たとえば、ＣＰＵ（Central Processing Unit）である。
出力部５は、処理部４による演算処理の結果を出力するものであり、たとえば、ディスプレイなどである。

次に、図２を参照しながら、数値シミュレーションの処理について説明する（適宜図１参照）。図２は、流体力学の数値シミュレーションの処理の大まかな流れを表わしたフローチャートである。

まず、使用者が入力部１によって数値シミュレーションの開始に関する操作を行うと、処理部４は、記憶部２に記憶された情報に基づいて、演算を行う定義点である格子点、各物理量の初期条件などを設定する（ステップＳ１）。
格子点は、たとえば、２次元の場合は３２×３２、３次元の場合は１６×１６×８、といったように設定される。また、各物理量の初期条件は、それぞれの格子点における変数の値、および、それらの値の空間勾配（微分係数）などが設定される。変数としては、たとえば、速度、圧力、温度、密度などが挙げられる。

次に、処理部４は、記憶部２に記憶された情報に基づいて、時間と空間に関する偏微分方程式などからなる支配方程式、各物理量の値を修正するための補間関数などを設定する（ステップＳ２）。
支配方程式としては、たとえば、圧縮性流体の運動方程式などが挙げられ、詳細については後記する。また、補間関数としては、３次関数、分母に関数を有する有理関数などが挙げられ、詳細については後記する。

続いて、処理部４は、ステップＳ１およびステップＳ２で設定した各種情報を用い、時刻ｔにおける各格子点での各物理量の値に基づいて、時刻（ｔ＋１）における各格子点での各物理量の値を算出する、という処理をｔ＝０,１,２,・・・ｎ-１に関して行う（ステップＳ３〜ステップＳ５）。
そして、この中のステップＳ４における処理が本発明のポイントであるが、その詳細については、図３とともに後で説明する。

次に、処理部４は、ステップＳ３〜ステップＳ５で算出した各物理量の値を、出力部５に出力する（ステップＳ６）。以上で、数値シミュレーションが完了する。

続いて、図３および図４を参照しながら、図２のステップＳ４における処理について説明する（適宜図１参照）。図３は、各格子点と変数、１階微分係数（以下、単に「微分係数」という）の関係を表わした図である。
図３に示すように、各格子点のうちの３点、ｊ−１点、ｊ点（特定点）、ｊ＋１点が、間隔Δｘで並んでいる。

ｊ−１点、ｊ点、ｊ＋１点において、座標はそれぞれ−Δｘ,０,Δｘ、変数（速度）ｕの値はそれぞれｕ_j-1,ｕ_j,ｕ_j+1、変数（速度）ｕの微分係数はそれぞれｕ_x,j-1,ｕ_x,j,ｕ_x,j+1、変数（圧力）ｐの値はそれぞれｐ_j-1,ｐ_j,ｐ_j+1、また、変数（圧力）ｐの微分係数はそれぞれｐ_x,j-1,ｐ_x,j,ｐ_x,j+1である。

図４は、図２のステップＳ４における処理の内容の詳細を表わしたフローチャートである。なお、ステップＳ４１〜ステップＳ４６において扱う各物理量はすべて時刻ｔのものであり、ステップＳ４７において算出される各物理量の値が時刻ｔ＋１のものである。

ここで、変数（速度）ｕの補間関数Ｕ（ｘ）を次の式（７）のような、エルミート補間関数の１つである３次関数として設定する。なお、エルミート補間関数とは、定義点上の変数の値とそれらの微分係数から決定する補間関数の総称であり、ここでの補間関数として適している。

Ｕ（ｘ）＝Ｃ₃ｘ³＋Ｃ₂ｘ²＋Ｃ₁ｘ＋Ｃ₀ ・・・（７）
また、Ｕ（ｘ）の微分式であるＵ_x（ｘ）は、式（８）として表わされる。
Ｕ_x（ｘ）＝３Ｃ₃ｘ²＋２Ｃ₂ｘ＋Ｃ₁ ・・・（８）

拘束条件は、Ｕ（Δｘ）＝ｕ_j+1,Ｕ（−Δｘ）＝ｕ_j-1,Ｕ_x（Δｘ）＝ｕ_x,j+1,Ｕ_x（−Δｘ）＝ｕ_x,j-1であり、これらにより、Ｃ₀〜Ｃ₃の値が決まる。なお、ここでは、ｊ点における微分係数を拘束条件としない、すなわち、Ｕ_x（０）＝ｕ_x,jを使用しないのが特徴である。もし、Ｕ_x（０）＝ｕ_x,jを使用すると、ステップＳ４２で取得するＣ₁がｕ_x,jと等しくなってしまい、ステップＳ４３で加重平均をとる意味がなくなってしまうからである。

図４に戻って、以上の条件の下、処理部４は、まず、本来の微分係数ｕ_x,j（第１の微分係数）を取得する（ステップＳ４１）。
次に、処理部４は、補間関数（式（７）および式（８））から求めた微分係数Ｃ₁（第２の微分係数）を取得する（ステップＳ４２）。つまり、ｊ点における微分係数は、次の式（９）により、Ｃ₁と求まるのである。
∂ｆ/∂ｘ＝Ｕ_x（０）＝Ｃ₁ ・・・（９）

続いて、処理部４は、近似微分係数∂ｕ/∂ｘ｜_j（第３の微分係数）を、次の式（１０）により算出する（ステップＳ４３）。なお、∂ｕ/∂ｘ｜_jは、補間関数などから求められたｊ点での微分係数を表わす。

そして、数値シミュレーションの条件に応じて、この式（１０）のαに適切な値を代入することによって、数値シミュレーションの安定性を高めることができる。

たとえば、多くの条件下では、α＝２/３を代入することによって、数値シミュレーションの安定性を高めることができることが、実証されている。その実証例（応用例）については後記するが、ここでは、α＝２/３が適していることが多いことの理由として考えられる数学的根拠について説明する。

まず、物理量Ｕ（ｘ）に対するｊ点周りのテーラー展開の式を、次の式（１１）に示す。

そして、Ｕ（ｘ+Δｘ）＝ｕ_j+1,Ｕ（ｘ）＝ｕ_j，Ｕ（ｘ−Δｘ）＝ｕ_j-1という拘束条件の下で、式（１１）を使って２次中心差分による離散近似を行うと式（１２）のようになり、その式（１２）を整理すると式（１３）が得られる。

また、空間勾配Ｕ_x（ｘ）に対するｊ点周りのテーラー展開の式を、次の式（１４）に示す。

次に、Ｕ_x（ｘ+Δｘ）＝ｕ_x,j+1,Ｕ_x（ｘ）＝ｕ_x,j，Ｕ_x（ｘ−Δｘ）＝ｕ_x,j-1という拘束条件の下で、式（１４）を使って前進差分、後退差分による離散近似を行い、式を整理すると次の式（１５）に示すように、３階微分項を得ることができる。

そして、式（１５）を式（１３）に代入すれば、次の式（１６）を得ることができる。

また、式（１０）におけるＣ₁を、式（７）と式（８）を使って計算すると次の式（１７）の通りとなる。
さらに、式（１０）に式（１７）とα＝２／３を代入すると次の式（１８）のようになる。

以上の説明からわかるように、式（１８）は、式（１６）の右辺の左２項と一致しており、４次精度（式（１６）の右辺の左から３項目以降の項におけるΔｘの最低次数が４）であるので、従来の２次精度などの計算手法やその他の計算手法よりも高い精度を実現することができる。
このことが、式（１０）においてα＝２/３を代入することの有利性を示す数学的根拠の１つであると考えられる。

ただし、数値シミュレーションの諸条件は千差万別なので、必ずしも常にα＝２/３がベストとは限らず、計算が安定するように、必要に応じてαに２/３以外の数値（２/３に近い値、２/３よりも小さな値など）を代入すればよい。
たとえば、数値シミュレーションの対象（気体、液体、電磁波など）の波の波長が、格子間隔と比較して相対的にある程度以上長いときは、αに２/３よりも小さな数値を代入することで計算を安定させることができることもある。なぜなら、対象の波の波長が格子間隔と比較して相対的に長い（格子間隔の数倍〜数百倍など）ときは、ｊ点の物理量を算出する際に、近接格子点（ｊ−１点、ｊ＋１点など）の物理量によって補正をする必要性が低減することもあるからである。
すなわち、本発明の特徴は、式（１０）に示すように、本来の微分係数、および、補間関数から求めた微分係数、という２つの値から近似微分係数を算出することであり、α＝２/３というのはその一例である。

図４に戻って、ステップＳ４４〜ステップＳ４６により、処理部４は変数（圧力）ｐに関する近似微分係数∂ｐ/∂ｘ｜_jを算出するが、ステップＳ４４〜ステップＳ４６は、ステップＳ４１〜ステップＳ４３においてｕをｐに置き換えた場合と同様であるので、説明を省略する。

次に、処理部４は、ステップＳ４３で算出した∂ｕ/∂ｘ｜_jとステップＳ４６で算出した∂ｐ/∂ｘ｜_jを使って支配方程式を演算し、時刻ｔ+１における各物理量（速度、圧力およびそれらの空間勾配）の値を算出する（ステップＳ４７）。
たとえば、数値シミュレーションの対象物が圧縮性流体の場合、支配方程式は、次の式（１９）のようになる。なお、ρは密度である。

この式（１９）を時間積分すれば、時刻ｔ+１における各物理量の値を算出することができる。時間積分する場合は、ルンゲ・クッタ法などを用いてもよい。また、∂ｕ/∂ｘ｜_jと∂ｐ/∂ｘ｜_jρをα＝２/３を使って算出した場合は、ρも同様にして、ρを算出するための補間関数から算出された値と本来の値を元に、２／３：１／３（２：１）の加重平均によって求めた値を使うようにしてもよい。
また、式（７）の補間関数の代わりに、次の式（２０）のような、分母に関数を有する有理関数を補間関数として用いてもよい。

なお、ここでは、分子を３次関数とし、分母を１次関数としているが、それぞれ、その他の関数を用いてもよい。また、この式（２０）は、圧力に関してだけでなく、速度、その他の変数にも適用できることは言うまでもない。
このように、補間関数に式（２０）のような有理関数を用いる場合も、近似微分係数∂ｕ/∂ｘ｜_jの算出の仕方は、補間関数を式（７）から式（２０）に入れ替えるほかは、図４のフローチャートのステップＳ４１〜ステップＳ４３の場合と同様である。

数値シミュレーションの初期に、数百倍あるいは数千倍の密度比や圧力比を有する爆風や衝撃波などに関する計算の場合は、式（２０）のような有理関数を用いることで、計算の安定性を高めることができる。

次に、本実施形態の数値計算装置１０による数値計算方法に関する４つの応用例について説明する。
＜応用例１：浅水波方程式＞
静力学大気モデルにおける力学過程は、浅水波方程式によって記述される。そして、浅水波方程式を数値シミュレーションに用いる場合には、重力ポテンシャルと流速に関して、高精度かつ安定なカップリングが必要となる。

なお、過去に、テンパートンらが、浅水波方程式に関して、高解像度計算を行っている（文献「An efficient two-time-level semi-Lagrangian semi-implicit integration scheme」（Temperton,C.,and Stainforth,A.,１９８７年,「Quart.J.Roy.Meteor.Soc.（１１３版）1025頁〜1039頁」）参照）。このテンパートンらの計算は、８００×８００の格子数で行われている。

この場合、支配方程式は、回転球面上の浅水波方程式である。この浅水波方程式は、ポーラーステレオ投影面におけるカルテシアン座標（ｘ,ｙ）を用いて、次の式（２１）〜式（２３）のように表わされる。
なお、ｔは時間、ｕ、ｖはそれぞれ流速のｘ方向成分、ｙ方向成分、φはジオポテンシャル高度であり、また、マップファクタをｍ、地球の自転の角速度をΩ、緯度をlatとし、Ｕ＝ｕ/ｍ，Ｖ＝ｖ/ｍ，Ｓ＝ｍ²，ｆ＝２Ωsin（lat）（コリオリパラメータ）である。

計算対象領域は、北極を中心とした正方領域であり、マップファクタの基準緯度を北緯６０度としている。また、境界条件として、赤道近傍で固定壁条件を与える。初期条件として、1984年2月28日12時の気圧５００ｍｂに対する解析に基づいて流速分布、ジオポテンシャル高度分布を与え、５分間隔で２日間分の積分を行う。

このような条件下において、図５（ａ）は、変数間のカップリングを行わない従来の計算手法で計算した場合のジオポテンシャル高度の誤差を表わした図、図５（ｂ）は、本実施形態の数値計算方法により計算した場合のジオポテンシャル高度の誤差を表わした図であり、格子点数はともに１００×１００である。

それぞれの図において、図の中心は北極を表わし、曲線５０１および曲線５１１は赤道を表わし、曲線５０２などおよび曲線５１２などはテンパートンらの計算によって得られたジオポテンシャル高度が等しい地点を表わすものである。また、符号５０３および符号５１３は、ジオポテンシャル高度の誤差の大きさと色との関係を表わしたものである。

この２つの図を比較してわかるように、本実施形態の数値計算方法による計算結果は、従来の計算手法よりもはるかに精度が高く、また、格子点数が８×８倍のテンパートンらの計算とほぼ同じ計算結果が得られている。

＜応用例２：乱流計算＞
乱流の計算には、各物理量をフーリエ展開し、波数空間上で計算を行うスペクトル法が主に用いられており、最も精度の高い計算手法として知られている。
一方、本実施形態の数値計算方法は、次の式（２４）〜式（２６）に示す支配方程式を使用する。式（２４）と式（２５）は、２次元非圧縮性ナビエ・ストークス方程式であり、式（２６）は、連続式である。
なお、ｕ、ｖはそれぞれ流速のｘ方向成分、ｙ方向成分、ｘとｙは座標、ｔは時間、Ｐは圧力、Ｒｅはレイノルズ数を表わす。

そして、格子点数が５１２×５１２の正方領域に、所定の乱流プロファイル（ここでは、渦度ω＝∂ｖ/∂ｘ−∂ｕ/∂ｙ）を与える。このとき、周期境界条件を用いて式（２４）〜式（２６）を解けば、図６のような結果が得られる。

図６は、散逸（Dissipation）とＥＴＴ（Eddy Turnover Time）との関係を示した図である。散逸とは、流体統計量の１つで、単位時間当たりに運動エネルギーが熱エネルギーに変わる量の大きさを表わすものであり、散逸をΦ、エンストロフィをε、レイノルズ数をＲｅとすると、Φ＝２ε/Ｒｅである。なお、渦度をωとすると、ε＝ω²/２である。
また、ＥＴＴは、無次元化した時間であり、初期における積分長をＬ、速度の標準偏差をＵ_r.m.sとすると、１ＥＴＴ＝Ｌ/Ｕ_r.m.sと表わされる。

図６に示すように、スペクトル法による計算結果（曲線６０１：Spectral）と、本実施形態の数値計算方法による計算結果（図中のドット：ＩＤＯ-ＳＣ）が、ほぼ一致している。
乱流現象は時間的、空間的に不規則であるという特徴を有するが、本実施形態の数値計算方法によれば、散逸だけでなく、速度の標準偏差、歪度、扁平度、エンストロフィなど、他の統計量に関しても、スペクトル法と同程度の、高い安定性の計算を実現することができる。

＜応用例３：２次元キャビティフロー問題＞
キャビティフローの問題は、非圧縮性流体計算のベンチマークとして広く知られており、従来では、ギアらによる多重格子法の計算結果が最も高精度とされている（文献「High-Resolution for incompressible flow using the Navier-Stokes equations and a multitigrid method」（Ghia,U.,Ghia,K.N.,and Shin,C.T.,１９８２年,「J.Comput.Phys.（４８版）387頁〜411頁」）参照）。

キャビティフローの問題は、圧力と流速のカップリングが重要であるため、従来のＩＤＯ法では、流速を正しく評価できなかった。図７（ａ）は、格子点数４０×４０で、従来のＩＤＯ法による計算を行った場合の結果を示した図である。
なお、従来のＩＤＯ法による計算も、本実施形態の数値計算方法による計算も、使用する支配方程式は、前記した式（２４）〜式（２６）である。

図７（ａ）において、ｘとｙは座標を表わす。また、ｖは、ｙ＝０．５の地点におけるｘ座標に対応するｘ方向の流速を表わし、曲線７０１（ＩＤＯ）の近傍のドット（Ghia et al）がギアらの計算結果（格子点数２５６×２５６）、曲線７０１が従来の計算手法による計算結果（格子点数４０×４０）である。さらに、ｕは、ｘ＝０．５の地点におけるｙ座標に対応するｙ方向の流速を表わし、曲線７０２（ＩＤＯ）の近傍のドット（Ghia et al）がギアらの計算結果（格子点数２５６×２５６）、曲線７０２が従来の計算手法による計算結果（格子点数４０×４０）である。図７（ａ）を見れば、従来の計算結果の精度がよくないことがわかる。

一方、図７（ｂ）は、本実施形態の数値計算装置による計算結果を示した図である。図７（ｂ）において、ｘとｙは座標を表わす。
また、ｖは、ｙ＝０．５の地点におけるｘ座標に対応するｘ方向の流速を表わし、曲線７１１（ＩＤＯ-ＳＣ）の近傍のドット（Ghia et al）がギアらの計算結果（格子点数２５６×２５６）、曲線７１１が本実施形態の数値計算方法による計算結果（格子点数４０×４０）である。さらに、ｕは、ｘ＝０．５の地点におけるｙ座標に対応するｙ方向の流速を表わし、曲線７１２（ＩＤＯ-ＳＣ）の近傍のドット（Ghia et al）がギアらの計算結果（格子点数２５６×２５６）、曲線７１２が本実施形態の数値計算方法による計算結果（格子点数４０×４０）である。

図７（ｂ）を見れば、本実施形態の数値計算方法による計算結果は、４０×４０という非常に少ない格子点数ながらも、格子点数２５６×２５６のギアらの計算結果とほぼ一致し、精度の高いことがわかる。

＜応用例４：衝撃波問題＞
初期条件として大きな圧力比や密度比が与えられる衝撃波問題において、従来のＩＤＯ法などによる計算では、圧力や密度と流速のカップリングが悪く、計算が初期から不安定になり破綻していた。

ここでは、密閉された管に密度と圧力の異なる２つの流体が隔膜を隔てて存在し、その隔膜を瞬間的に取り除いたときに、高圧の流体が低圧の流体側に流入し、低圧の流体が急速に圧縮されて衝撃波が発生する、という１次元衝撃波管問題について数値シミュレーションを行った場合について説明する。

この場合、支配方程式は、１次元圧縮性ナビエ・ストークス方程式と連続式であり、非保存形で表わすと、次の式（２７）〜式（２９）のようになる。
なお、ｔは時間、ｘは座標、ρは密度、ｐは圧力、ｑは人工粘性、ｕは流速、ｅは内部エネルギーである。

図８（ａ）は、圧力比１０：１の場合の圧力と座標の関係図であり、図８（ｂ）は、圧力比１００００：１の場合の圧力と座標の関係図である。
図８（ａ）において、横軸は座標、縦軸は流体の圧力（Pressure）である。ここでは、初期に、ｘ＝１．０よりも左側に圧力１の流体が存在し、ｘ＝１．０よりも右側に圧力０．１の流体が存在し、その条件下で、ｘ＝１．０の地点に設けられていた隔膜を瞬間的に取り除いてから微小時間後のある時刻における流体のそれぞれの座標における圧力を示している。

曲線８０１（Exact）は論理的に求めた解析解（厳密解）であり、ドット（ＩＤＯ-ＳＣ）で示された本実施形態の数値計算方法による計算結果は、曲線８０１とほぼ一致していることがわかる。

また、図８（ｂ）において、横軸は座標、縦軸（対数表示）は流体の圧力（Pressure）である。ここでは、初期に、ｘ＝１．０よりも左側に圧力１（１０⁰）の流体が存在し、ｘ＝１．０よりも右側に圧力１０^-4の流体が存在し、その条件下で、ｘ＝１．０の地点に設けられていた隔膜を瞬間的に取り除いてから微小時間後のある時刻における流体のそれぞれの座標における圧力を示している。

曲線８０２（Exact）は論理的に求めた解析解（厳密解）であり、ドット（ＩＤＯ-ＳＣ）で示された本実施形態の数値計算方法による計算結果は、曲線８０２とほぼ一致していることがわかる。
本実施形態の数値計算方法によれば、図８（ａ）の場合のような一般的な衝撃波管問題だけでなく、図８（ｂ）の場合のような圧力比などがさらに大きな衝撃波管問題や壁面との相互作用を含む問題においても、安定的かつ高精度の高い計算を行うことができる。なお、密度比や圧力比などが極端に大きな場合は、補間関数として、３次関数ではなく、分母に関数を有する有理関数、たとえば、分母に１次関数、分子に３次関数を有する有理関数などを用いることが望ましい。

以上説明したように、本実施形態の数値計算方法によれば、コロケート格子を用いた多変数偏微分方程式において、変数間の安定なカップリングを実現することができる。また、この数値計算方法をコンピュータに実行させるためのプログラムを作成し、そのプログラムを記憶部（記録媒体）２（図１参照）に記憶（記録）することができる。

以上で実施形態の説明を終えるが、本発明の態様はこれらに限定されるものではない。
たとえば、本発明は、電磁気に関する数値シミュレーションや、レイリーテーラー不安定性の解析など、偏微分方程式を使った自然現象の再現全般に広く適用することができる。
その他、数値計算装置や各フローチャートなどの具体的な構成について、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。

数値計算装置の構成図である。 数値シミュレーションの処理の大まかな流れを表わしたフローチャートである。 各格子点と変数、微分係数の関係を表わした図である。 図２のステップＳ４における処理の内容を表わしたフローチャートである。 応用例１の計算結果を表わす図であり、（ａ）は、従来の計算手法で計算した場合の誤差を表わした図、（ｂ）は、本実施形態の数値計算装置により計算した場合の誤差を表わした図である。 応用例２の計算結果を表わす図であり、散逸とＥＴＴとの関係を示した図である。 応用例３の計算結果を表わす図であり、（ａ）は、従来のＩＤＯ法による計算結果を示した図、（ｂ）は、本実施形態の数値計算装置による計算結果を示した図である。 応用例４の計算結果を表わす図であり、（ａ）は、圧力比１０：１の場合の圧力と座標の関係図、（ｂ）は、圧力比１００００：１の場合の圧力と座標の関係図である。

符号の説明

１ 入力部
２ 記憶部
３ メモリ
４ 処理部
５ 出力部
１０ 数値計算装置

Claims (7)

1. 空間上の複数の定義点それぞれにおける１つ以上の物理的な変数の値を、所定の物理現象を表す偏微分方程式を含む支配方程式と、前記支配方程式で使用する前記変数の値を修正するための所定の関数であって前記定義点ごとに自身の定義点以外の２点以上の前記変数の値に基づいて係数が決定される補間関数と、を用いて、所定の微小時間単位で更新する数値計算装置による数値計算方法であって、
   前記数値計算装置は、記憶部と処理部とを有し、
   前記記憶部は、前記支配方程式と前記補間関数とを記憶しており、
   前記処理部は、前記複数の定義点のうちの特定点について前記変数の値を更新するとき、
   前記支配方程式を用いて前記変数に関する第１の微分係数を算出し、
   前記補間関数を用いて前記変数に関する第２の微分係数を算出し、
   前記第１の微分係数と前記第２の微分係数との加重平均によって前記変数に関する第３の微分係数を算出し、
   前記第３の微分係数と前記支配方程式とに基づいて演算を行うことで、前記変数の値を前記所定の微小時間後の値に更新する
   ことを特徴とする数値計算方法。
2. 前記複数の定義点は、すべての変数やその微分係数が同じ格子上に定義されるコロケート格子点である
   ことを特徴とする請求項１に記載の数値計算方法。
3. 前記補間関数は、エルミート補間関数であることを特徴とする請求項１または請求項２に記載の数値計算方法。
4. 前記補間関数は、３次関数を分子に持つ有理関数であることを特徴とする請求項１または請求項２に記載の数値計算方法。
5. 前記加重平均は、前記第１の微分係数と前記第２の微分係数の比が１/３対２/３あるいはそれに近い値であることを特徴とする請求項１または請求項２に記載の数値計算方法。
6. 請求項１から請求項５のいずれか一項に記載の数値計算方法をコンピュータに実行させるためのプログラム。
7. 請求項６に記載のプログラムを記録した記録媒体。