TOP > 国内特許検索 > 数値計算方法、プログラムおよび記録媒体 > 明細書

明細書 :数値計算方法、プログラムおよび記録媒体

発行国 日本国特許庁(JP)
公報種別 特許公報(B2)
特許番号 特許第4304277号 (P4304277)
公開番号 特開2006-318355 (P2006-318355A)
登録日 平成21年5月15日(2009.5.15)
発行日 平成21年7月29日(2009.7.29)
公開日 平成18年11月24日(2006.11.24)
発明の名称または考案の名称 数値計算方法、プログラムおよび記録媒体
国際特許分類 G06F  17/13        (2006.01)
FI G06F 17/13
請求項の数または発明の数 7
全頁数 16
出願番号 特願2005-142404 (P2005-142404)
出願日 平成17年5月16日(2005.5.16)
審査請求日 平成19年7月31日(2007.7.31)
特許権者または実用新案権者 【識別番号】304021417
【氏名又は名称】国立大学法人東京工業大学
発明者または考案者 【氏名】青木 尊之
【氏名】今井 陽介
個別代理人の代理人 【識別番号】100064414、【弁理士】、【氏名又は名称】磯野 道造
【識別番号】100111545、【弁理士】、【氏名又は名称】多田 悦夫
審査官 【審査官】須田 勝巳
参考文献・文献 特開2002-006917(JP,A)
特開平02-210577(JP,A)
吉田 浩 外2名,ルンゲ・クッタ時間積分を用いた局所補間微分オペレータ(IDO)法による双曲型方程式解法の計算精度と数値安定性の向上,電子情報通信学会論文誌,社団法人電子情報通信学会,2003年 3月 1日,第J86-A巻,第3号,p223-231
調査した分野 G06F 17/13
特許請求の範囲 【請求項1】
空間上の複数の定義点それぞれにおける1つ以上の物理的な変数の値を、所定の物理現象を表す偏微分方程式を含む支配方程式と、前記支配方程式で使用する前記変数の値を修正するための所定の関数であって前記定義点ごとに自身の定義点以外の2点以上の前記変数の値に基づいて係数が決定される補間関数と、を用いて、所定の微小時間単位で更新する数値計算装置による数値計算方法であって、
前記数値計算装置は、記憶部と処理部とを有し、
前記記憶部は、前記支配方程式と前記補間関数とを記憶しており、
前記処理部は、前記複数の定義点のうちの特定点について前記変数の値を更新するとき、
前記支配方程式を用いて前記変数に関する第1の微分係数を算出し、
前記補間関数を用いて前記変数に関する第2の微分係数を算出し、
前記第の微分係数と前記第の微分係数との加重平均によって前記変数に関する第3の微分係数を算出し、
前記第3の微分係数と前記支配方程式とに基づいて演算を行うことで、前記変数の値を前記所定の微小時間後の値に更新する
ことを特徴とする数値計算方法。
【請求項2】
前記複数の定義点は、すべての変数やその微分係数が同じ格子上に定義されるコロケート格子点である
ことを特徴とする請求項1に記載の数値計算方法。
【請求項3】
前記補間関数は、エルミート補間関数であることを特徴とする請求項1または請求項2に記載の数値計算方法。
【請求項4】
前記補間関数は、3次関数を分子に持つ有理関数であることを特徴とする請求項1または請求項2に記載の数値計算方法。
【請求項5】
前記加重平均は、前記第の微分係数と前記第2の微分係数の比が/3対/3あるいはそれに近い値であることを特徴とする請求項1または請求項2に記載の数値計算方法。
【請求項6】
請求項1から請求項5のいずれか一項に記載の数値計算方法をコンピュータに実行させるためのプログラム。
【請求項7】
請求項6に記載のプログラムを記録した記録媒体。
発明の詳細な説明 【技術分野】
【0001】
本発明は、偏微分方程式の演算に使用する数値計算方法に関する。
【背景技術】
【0002】
近年、コンピュータを用いて自然現象を再現する取り組みが、多く行われるようになってきている。たとえば、流体力学などの数値シミュレーションがその例で、その対象物は、液体、気体、電磁波など様々である。
そして、そのような数値シミュレーションを行う場合には、2以上の変数(以下、「多変数」という)に関する偏微分方程式を、離散化された空間および時間に関して解く、という手法が用いられることが多い。なお、変数には、たとえば、速度、圧力、温度、密度などが挙げられる。
【0003】
多変数偏微分方程式における数値解法(シミュレーション)では、変数間のカップリングについて、高い時間的安定(以下、単に「安定」という)性が要求される。変数間のカップリングの安定性が低いと、短時間で計算が破綻してしまうからである。
そこで、従来、流体解析におけるSMAC(Simplified Marker and Cell)法や電磁波解析におけるFDTD(Finite Difference Time Domain)法においては、空間を離散化するときに、変数ごとの定義点をずらして設定するスタッガード格子が多く用いられていた。スタッガード格子は、変数間のカップリングの安定性が高いためである。
【0004】
また、それぞれの変数の値だけでなく、それらの空間勾配(微分係数)なども利用する計算手法であるマルチモーメントスキームにおいても、スタッガード格子が多く用いられてきた。しかし、スタッガード格子は、変数間のカップリングの安定性が高い反面、計算精度や並列計算アルゴリズムの容易さなどの点で、欠点があった。
【0005】
一方、すべての変数や空間勾配を同じ格子上に定義するコロケート格子は、スタッガード格子に比べて、計算精度や並列計算アルゴリズムの容易さなどの点で優れている。
たとえば、マルチモーメントスキームにおける計算精度は、スタッガード格子の2次に対し、コロケート格子は4次である。
また、コロケート格子は、スタッガード格子よりも、今後の発展が期待される計算手法であるAMR(Adaptive Mesh Refinement:適合格子細分化)法やCut-Cell法との組み合わせが容易であるという利点も有する。
【発明の開示】
【発明が解決しようとする課題】
【0006】
しかしながら、コロケート格子は、変数間のカップリングの安定性が低く、計算が短時間で破綻してしまうことが多いため、特殊な条件下でしか使用することができず、適用範囲が非常に限定されてしまう、という問題があった。
【0007】
たとえば、マルチモーメントスキームであるIDO(Interpolated Differential Operator:局所補間微分オペレータ)法では、間隔がΔxである3点、j-1点、j点、j+1点に関する中心補間関数F(x)を、j点を中心として、次の式(1)のような5次関数として表わすことができる。
なお、fj-1,fj,fj+1は、それぞれ、j-1点、j点、j+1点における変数の値を表わし、fx,j-1, fx,j, fx,j+1は、それぞれ、j-1点、j点、j+1点における変数の微分係数を表わす。また、xは座標を表わす。
【0008】
F(x)=A55+A44+A33+A22+fx,jx+fj ・・・(1)
また、F(x)の微分式であるFx(x)は、式(2)のように表わされる。
x(x)=5A54+4A43+3A32+2A2x+fx,j ・・・(2)
拘束条件は、F(Δx)=fj+1,F(-Δx)=fj-1,Fx(Δx)=fx,j+1, Fx(-Δx)=fx,j-1であり、これらにより、A2~A5の値が決まる。
【0009】
そして、この条件下で、支配方程式(演算に使う偏微分方程式など)中の空間微分項は、式(1)、式(2)、および、式(2)をさらに微分した式などにより、∂f/∂x=Fx(0)=fx,j,∂2f/∂x2=Fxx(0)=2A2,・・・のように求めることができる。なお、Fxx(x)は、Fx(x)を微分した式を表わす。
【0010】
ここで、多変数型の偏微分方程式である次の式(3)、式(4)を考える。なお、fとuは変数、tは時間、xは座標を表わす。
【数1】
JP0004304277B2_000002t.gif

【0011】
コロケート格子上では、前記補間関数(式(1)および式(2))を用いて、次の式(5)、式(6)のような空間離散化がなされる。なお、∂f/∂t|j,∂u/∂t|jは、それぞれ、補間関数などから求めたj点での微分係数を表わす。
【数2】
JP0004304277B2_000003t.gif
そして、式(5)、式(6)からわかるように、離散化式中にそれぞれj点の変数しか現われないため、変数間のカップリングの安定性が悪い(低い)。
【0012】
そこで、本発明は、前記問題点に鑑みてなされたものであり、コロケート格子を用いた偏微分方程式の数値解法において、変数間のカップリングの安定性が高い数値計算方法を提供することを目的とする。
【課題を解決するための手段】
【0013】
前記課題を解決するために、本発明に係る数値計算方法は、各定義点上で変数に関する数値計算を行う数値計算装置による数値計算方法であって、前記数値計算装置は、記憶部と処理部を有し、前記記憶部は、前記各定義点のうち3点以上における各変数の値と、当該3点以上の各変数の値にそれぞれ対応する第1の微分係数と、前記各変数の値と前記3点以上のうち特定点を含まない定義点における前記第1の微分係数により決定される補間関数と、を記憶し、前記処理部は、前記補間関数によって前記特定点における第2の微分係数を算出し、前記第2の微分係数と前記特定点における第1の微分係数との加重平均によって前記特定点における第3の微分係数を算出することを特徴とする。
また、当該数値計算方法をコンピュータに実行させるためのプログラムを作成し、さらに、そのプログラムを記録媒体に記録することができる。
【発明の効果】
【0014】
本発明によれば、コロケート格子を用いた偏微分方程式の数値解法において、変数間のカップリングの安定性を高めることができる。
【発明を実施するための最良の形態】
【0015】
以下、本発明の実施形態に係る数値計算方法について、適宜図面を参照しながら説明する。
まず、図1を参照しながら、数値計算方法を行う数値計算装置の構成について説明する。図1は、数値計算装置の構成図である。
【0016】
数値計算装置10は、入力部1、記憶部(記録媒体)2、メモリ3、処理部4および出力部5を備えて構成されるコンピュータ装置であり、たとえば、パソコン、ワークステーション、スーパーコンピュータなどにより実現することができる。
なお、処理部4は、ここでは1つとしているが、複数設けることにより並列計算を行うようにしてもよい。
【0017】
入力部1は、数値計算装置10の使用者(以下、単に「使用者」という)が各種情報を入力する手段であり、たとえば、キーボード、マウスなどである。
記憶部2は、各種情報を記憶するものであり、たとえば、ROM(Read Only Memory)、ハードディスクなどである。記憶部2は、ここでは、流体力学の数値シミュレーションを行うために、数値計算を行う定義点である格子点、物理量(変数、微分係数など)の初期条件、支配方程式(自然現象を表わす偏微分方程式など)、補間関数、各種プログラムなどを記憶している。
【0018】
メモリ3は、処理部4が使用する一時記憶手段であり、たとえば、RAM(Random Access Memory)などである。
処理部4は、記憶部2に記憶された各種情報(プログラム、データなど)などに基づいて演算処理を行うものであり、たとえば、CPU(Central Processing Unit)である。
出力部5は、処理部4による演算処理の結果を出力するものであり、たとえば、ディスプレイなどである。
【0019】
次に、図2を参照しながら、数値シミュレーションの処理について説明する(適宜図1参照)。図2は、流体力学の数値シミュレーションの処理の大まかな流れを表わしたフローチャートである。
【0020】
まず、使用者が入力部1によって数値シミュレーションの開始に関する操作を行うと、処理部4は、記憶部2に記憶された情報に基づいて、演算を行う定義点である格子点、各物理量の初期条件などを設定する(ステップS1)。
格子点は、たとえば、2次元の場合は32×32、3次元の場合は16×16×8、といったように設定される。また、各物理量の初期条件は、それぞれの格子点における変数の値、および、それらの値の空間勾配(微分係数)などが設定される。変数としては、たとえば、速度、圧力、温度、密度などが挙げられる。
【0021】
次に、処理部4は、記憶部2に記憶された情報に基づいて、時間と空間に関する偏微分方程式などからなる支配方程式、各物理量の値を修正するための補間関数などを設定する(ステップS2)。
支配方程式としては、たとえば、圧縮性流体の運動方程式などが挙げられ、詳細については後記する。また、補間関数としては、3次関数、分母に関数を有する有理関数などが挙げられ、詳細については後記する。
【0022】
続いて、処理部4は、ステップS1およびステップS2で設定した各種情報を用い、時刻tにおける各格子点での各物理量の値に基づいて、時刻(t+1)における各格子点での各物理量の値を算出する、という処理をt=0,1,2,・・・n-1に関して行う(ステップS3~ステップS5)。
そして、この中のステップS4における処理が本発明のポイントであるが、その詳細については、図3とともに後で説明する。
【0023】
次に、処理部4は、ステップS3~ステップS5で算出した各物理量の値を、出力部5に出力する(ステップS6)。以上で、数値シミュレーションが完了する。
【0024】
続いて、図3および図4を参照しながら、図2のステップS4における処理について説明する(適宜図1参照)。図3は、各格子点と変数、1階微分係数(以下、単に「微分係数」という)の関係を表わした図である。
図3に示すように、各格子点のうちの3点、j-1点、j点(特定点)、j+1点が、間隔Δxで並んでいる。
【0025】
j-1点、j点、j+1点において、座標はそれぞれ-Δx,0,Δx、変数(速度)uの値はそれぞれuj-1,uj,uj+1、変数(速度)uの微分係数はそれぞれux,j-1,ux,j,ux,j+1、変数(圧力)pの値はそれぞれpj-1,pj,pj+1、また、変数(圧力)pの微分係数はそれぞれpx,j-1,px,j,px,j+1である。
【0026】
図4は、図2のステップS4における処理の内容の詳細を表わしたフローチャートである。なお、ステップS41~ステップS46において扱う各物理量はすべて時刻tのものであり、ステップS47において算出される各物理量の値が時刻t+1のものである。
【0027】
ここで、変数(速度)uの補間関数U(x)を次の式(7)のような、エルミート補間関数の1つである3次関数として設定する。なお、エルミート補間関数とは、定義点上の変数の値とそれらの微分係数から決定する補間関数の総称であり、ここでの補間関数として適している。
【0028】
U(x)=C33+C22+C1x+C0 ・・・(7)
また、U(x)の微分式であるUx(x)は、式(8)として表わされる。
x(x)=3C32+2C2x+C1 ・・・(8)
【0029】
拘束条件は、U(Δx)=uj+1,U(-Δx)=uj-1,Ux(Δx)=ux,j+1,Ux(-Δx)=ux,j-1であり、これらにより、C0~C3の値が決まる。なお、ここでは、j点における微分係数を拘束条件としない、すなわち、Ux(0)=ux,jを使用しないのが特徴である。もし、Ux(0)=ux,jを使用すると、ステップS42で取得するC1がux,jと等しくなってしまい、ステップS43で加重平均をとる意味がなくなってしまうからである。
【0030】
図4に戻って、以上の条件の下、処理部4は、まず、本来の微分係数ux,j(第1の微分係数)を取得する(ステップS41)。
次に、処理部4は、補間関数(式(7)および式(8))から求めた微分係数C1(第2の微分係数)を取得する(ステップS42)。つまり、j点における微分係数は、次の式(9)により、C1と求まるのである。
∂f/∂x=Ux(0)=C1 ・・・(9)
【0031】
続いて、処理部4は、近似微分係数∂u/∂x|j(第3の微分係数)を、次の式(10)により算出する(ステップS43)。なお、∂u/∂x|jは、補間関数などから求められたj点での微分係数を表わす。
【数3】
JP0004304277B2_000004t.gif

【0032】
そして、数値シミュレーションの条件に応じて、この式(10)のαに適切な値を代入することによって、数値シミュレーションの安定性を高めることができる。
【0033】
たとえば、多くの条件下では、α=2/3を代入することによって、数値シミュレーションの安定性を高めることができることが、実証されている。その実証例(応用例)については後記するが、ここでは、α=2/3が適していることが多いことの理由として考えられる数学的根拠について説明する。
【0034】
まず、物理量U(x)に対するj点周りのテーラー展開の式を、次の式(11)に示す。
【数4】
JP0004304277B2_000005t.gif

【0035】
そして、U(x+Δx)=uj+1,U(x)=uj,U(x-Δx)=uj-1という拘束条件の下で、式(11)を使って2次中心差分による離散近似を行うと式(12)のようになり、その式(12)を整理すると式(13)が得られる。
【数5】
JP0004304277B2_000006t.gif

【0036】
また、空間勾配Ux(x)に対するj点周りのテーラー展開の式を、次の式(14)に示す。
【数6】
JP0004304277B2_000007t.gif

【0037】
次に、Ux(x+Δx)=ux,j+1,Ux(x)=ux,j,Ux(x-Δx)=ux,j-1という拘束条件の下で、式(14)を使って前進差分、後退差分による離散近似を行い、式を整理すると次の式(15)に示すように、3階微分項を得ることができる。
【数7】
JP0004304277B2_000008t.gif

【0038】
そして、式(15)を式(13)に代入すれば、次の式(16)を得ることができる。
【数8】
JP0004304277B2_000009t.gif

【0039】
また、式(10)におけるC1を、式(7)と式(8)を使って計算すると次の式(17)の通りとなる。
【数9】
JP0004304277B2_000010t.gif
さらに、式(10)に式(17)とα=2/3を代入すると次の式(18)のようになる。
【数10】
JP0004304277B2_000011t.gif

【0040】
以上の説明からわかるように、式(18)は、式(16)の右辺の左2項と一致しており、4次精度(式(16)の右辺の左から3項目以降の項におけるΔxの最低次数が4)であるので、従来の2次精度などの計算手法やその他の計算手法よりも高い精度を実現することができる。
このことが、式(10)においてα=2/3を代入することの有利性を示す数学的根拠の1つであると考えられる。
【0041】
ただし、数値シミュレーションの諸条件は千差万別なので、必ずしも常にα=2/3がベストとは限らず、計算が安定するように、必要に応じてαに2/3以外の数値(2/3に近い値、2/3よりも小さな値など)を代入すればよい。
たとえば、数値シミュレーションの対象(気体、液体、電磁波など)の波の波長が、格子間隔と比較して相対的にある程度以上長いときは、αに2/3よりも小さな数値を代入することで計算を安定させることができることもある。なぜなら、対象の波の波長が格子間隔と比較して相対的に長い(格子間隔の数倍~数百倍など)ときは、j点の物理量を算出する際に、近接格子点(j-1点、j+1点など)の物理量によって補正をする必要性が低減することもあるからである。
すなわち、本発明の特徴は、式(10)に示すように、本来の微分係数、および、補間関数から求めた微分係数、という2つの値から近似微分係数を算出することであり、α=2/3というのはその一例である。
【0042】
図4に戻って、ステップS44~ステップS46により、処理部4は変数(圧力)pに関する近似微分係数∂p/∂x|jを算出するが、ステップS44~ステップS46は、ステップS41~ステップS43においてuをpに置き換えた場合と同様であるので、説明を省略する。
【0043】
次に、処理部4は、ステップS43で算出した∂u/∂x|jとステップS46で算出した∂p/∂x|jを使って支配方程式を演算し、時刻t+1における各物理量(速度、圧力およびそれらの空間勾配)の値を算出する(ステップS47)。
たとえば、数値シミュレーションの対象物が圧縮性流体の場合、支配方程式は、次の式(19)のようになる。なお、ρは密度である。
【数11】
JP0004304277B2_000012t.gif

【0044】
この式(19)を時間積分すれば、時刻t+1における各物理量の値を算出することができる。時間積分する場合は、ルンゲ・クッタ法などを用いてもよい。また、∂u/∂x|jと∂p/∂x|jρをα=2/3を使って算出した場合は、ρも同様にして、ρを算出するための補間関数から算出された値と本来の値を元に、2/3:1/3(2:1)の加重平均によって求めた値を使うようにしてもよい。
また、式(7)の補間関数の代わりに、次の式(20)のような、分母に関数を有する有理関数を補間関数として用いてもよい。
【数12】
JP0004304277B2_000013t.gif

【0045】
なお、ここでは、分子を3次関数とし、分母を1次関数としているが、それぞれ、その他の関数を用いてもよい。また、この式(20)は、圧力に関してだけでなく、速度、その他の変数にも適用できることは言うまでもない。
このように、補間関数に式(20)のような有理関数を用いる場合も、近似微分係数∂u/∂x|jの算出の仕方は、補間関数を式(7)から式(20)に入れ替えるほかは、図4のフローチャートのステップS41~ステップS43の場合と同様である。
【0046】
数値シミュレーションの初期に、数百倍あるいは数千倍の密度比や圧力比を有する爆風や衝撃波などに関する計算の場合は、式(20)のような有理関数を用いることで、計算の安定性を高めることができる。
【0047】
次に、本実施形態の数値計算装置10による数値計算方法に関する4つの応用例について説明する。
<応用例1:浅水波方程式>
静力学大気モデルにおける力学過程は、浅水波方程式によって記述される。そして、浅水波方程式を数値シミュレーションに用いる場合には、重力ポテンシャルと流速に関して、高精度かつ安定なカップリングが必要となる。
【0048】
なお、過去に、テンパートンらが、浅水波方程式に関して、高解像度計算を行っている(文献「An efficient two-time-level semi-Lagrangian semi-implicit integration scheme」(Temperton,C.,and Stainforth,A.,1987年,「Quart.J.Roy.Meteor.Soc.(113版)1025頁~1039頁」)参照)。このテンパートンらの計算は、800×800の格子数で行われている。
【0049】
この場合、支配方程式は、回転球面上の浅水波方程式である。この浅水波方程式は、ポーラーステレオ投影面におけるカルテシアン座標(x,y)を用いて、次の式(21)~式(23)のように表わされる。
なお、tは時間、u、vはそれぞれ流速のx方向成分、y方向成分、φはジオポテンシャル高度であり、また、マップファクタをm、地球の自転の角速度をΩ、緯度をlatとし、U=u/m,V=v/m,S=m2,f=2Ωsin(lat)(コリオリパラメータ)である。
【数13】
JP0004304277B2_000014t.gif

【0050】
計算対象領域は、北極を中心とした正方領域であり、マップファクタの基準緯度を北緯60度としている。また、境界条件として、赤道近傍で固定壁条件を与える。初期条件として、1984年2月28日12時の気圧500mbに対する解析に基づいて流速分布、ジオポテンシャル高度分布を与え、5分間隔で2日間分の積分を行う。
【0051】
このような条件下において、図5(a)は、変数間のカップリングを行わない従来の計算手法で計算した場合のジオポテンシャル高度の誤差を表わした図、図5(b)は、本実施形態の数値計算方法により計算した場合のジオポテンシャル高度の誤差を表わした図であり、格子点数はともに100×100である。
【0052】
それぞれの図において、図の中心は北極を表わし、曲線501および曲線511は赤道を表わし、曲線502などおよび曲線512などはテンパートンらの計算によって得られたジオポテンシャル高度が等しい地点を表わすものである。また、符号503および符号513は、ジオポテンシャル高度の誤差の大きさと色との関係を表わしたものである。
【0053】
この2つの図を比較してわかるように、本実施形態の数値計算方法による計算結果は、従来の計算手法よりもはるかに精度が高く、また、格子点数が8×8倍のテンパートンらの計算とほぼ同じ計算結果が得られている。
【0054】
<応用例2:乱流計算>
乱流の計算には、各物理量をフーリエ展開し、波数空間上で計算を行うスペクトル法が主に用いられており、最も精度の高い計算手法として知られている。
一方、本実施形態の数値計算方法は、次の式(24)~式(26)に示す支配方程式を使用する。式(24)と式(25)は、2次元非圧縮性ナビエ・ストークス方程式であり、式(26)は、連続式である。
なお、u、vはそれぞれ流速のx方向成分、y方向成分、xとyは座標、tは時間、Pは圧力、Reはレイノルズ数を表わす。
【数14】
JP0004304277B2_000015t.gif

【0055】
そして、格子点数が512×512の正方領域に、所定の乱流プロファイル(ここでは、渦度ω=∂v/∂x-∂u/∂y)を与える。このとき、周期境界条件を用いて式(24)~式(26)を解けば、図6のような結果が得られる。
【0056】
図6は、散逸(Dissipation)とETT(Eddy Turnover Time)との関係を示した図である。散逸とは、流体統計量の1つで、単位時間当たりに運動エネルギーが熱エネルギーに変わる量の大きさを表わすものであり、散逸をΦ、エンストロフィをε、レイノルズ数をReとすると、Φ=2ε/Reである。なお、渦度をωとすると、ε=ω2/2である。
また、ETTは、無次元化した時間であり、初期における積分長をL、速度の標準偏差をUr.m.sとすると、1ETT=L/Ur.m.sと表わされる。
【0057】
図6に示すように、スペクトル法による計算結果(曲線601:Spectral)と、本実施形態の数値計算方法による計算結果(図中のドット:IDO-SC)が、ほぼ一致している。
乱流現象は時間的、空間的に不規則であるという特徴を有するが、本実施形態の数値計算方法によれば、散逸だけでなく、速度の標準偏差、歪度、扁平度、エンストロフィなど、他の統計量に関しても、スペクトル法と同程度の、高い安定性の計算を実現することができる。
【0058】
<応用例3:2次元キャビティフロー問題>
キャビティフローの問題は、非圧縮性流体計算のベンチマークとして広く知られており、従来では、ギアらによる多重格子法の計算結果が最も高精度とされている(文献「High-Resolution for incompressible flow using the Navier-Stokes equations and a multitigrid method」(Ghia,U.,Ghia,K.N.,and Shin,C.T.,1982年,「J.Comput.Phys.(48版)387頁~411頁」)参照)。
【0059】
キャビティフローの問題は、圧力と流速のカップリングが重要であるため、従来のIDO法では、流速を正しく評価できなかった。図7(a)は、格子点数40×40で、従来のIDO法による計算を行った場合の結果を示した図である。
なお、従来のIDO法による計算も、本実施形態の数値計算方法による計算も、使用する支配方程式は、前記した式(24)~式(26)である。
【0060】
図7(a)において、xとyは座標を表わす。また、vは、y=0.5の地点におけるx座標に対応するx方向の流速を表わし、曲線701(IDO)の近傍のドット(Ghia et al)がギアらの計算結果(格子点数256×256)、曲線701が従来の計算手法による計算結果(格子点数40×40)である。さらに、uは、x=0.5の地点におけるy座標に対応するy方向の流速を表わし、曲線702(IDO)の近傍のドット(Ghia et al)がギアらの計算結果(格子点数256×256)、曲線702が従来の計算手法による計算結果(格子点数40×40)である。図7(a)を見れば、従来の計算結果の精度がよくないことがわかる。
【0061】
一方、図7(b)は、本実施形態の数値計算装置による計算結果を示した図である。図7(b)において、xとyは座標を表わす。
また、vは、y=0.5の地点におけるx座標に対応するx方向の流速を表わし、曲線711(IDO-SC)の近傍のドット(Ghia et al)がギアらの計算結果(格子点数256×256)、曲線711が本実施形態の数値計算方法による計算結果(格子点数40×40)である。さらに、uは、x=0.5の地点におけるy座標に対応するy方向の流速を表わし、曲線712(IDO-SC)の近傍のドット(Ghia et al)がギアらの計算結果(格子点数256×256)、曲線712が本実施形態の数値計算方法による計算結果(格子点数40×40)である。
【0062】
図7(b)を見れば、本実施形態の数値計算方法による計算結果は、40×40という非常に少ない格子点数ながらも、格子点数256×256のギアらの計算結果とほぼ一致し、精度の高いことがわかる。
【0063】
<応用例4:衝撃波問題>
初期条件として大きな圧力比や密度比が与えられる衝撃波問題において、従来のIDO法などによる計算では、圧力や密度と流速のカップリングが悪く、計算が初期から不安定になり破綻していた。
【0064】
ここでは、密閉された管に密度と圧力の異なる2つの流体が隔膜を隔てて存在し、その隔膜を瞬間的に取り除いたときに、高圧の流体が低圧の流体側に流入し、低圧の流体が急速に圧縮されて衝撃波が発生する、という1次元衝撃波管問題について数値シミュレーションを行った場合について説明する。
【0065】
この場合、支配方程式は、1次元圧縮性ナビエ・ストークス方程式と連続式であり、非保存形で表わすと、次の式(27)~式(29)のようになる。
なお、tは時間、xは座標、ρは密度、pは圧力、qは人工粘性、uは流速、eは内部エネルギーである。
【数15】
JP0004304277B2_000016t.gif

【0066】
図8(a)は、圧力比10:1の場合の圧力と座標の関係図であり、図8(b)は、圧力比10000:1の場合の圧力と座標の関係図である。
図8(a)において、横軸は座標、縦軸は流体の圧力(Pressure)である。ここでは、初期に、x=1.0よりも左側に圧力1の流体が存在し、x=1.0よりも右側に圧力0.1の流体が存在し、その条件下で、x=1.0の地点に設けられていた隔膜を瞬間的に取り除いてから微小時間後のある時刻における流体のそれぞれの座標における圧力を示している。
【0067】
曲線801(Exact)は論理的に求めた解析解(厳密解)であり、ドット(IDO-SC)で示された本実施形態の数値計算方法による計算結果は、曲線801とほぼ一致していることがわかる。
【0068】
また、図8(b)において、横軸は座標、縦軸(対数表示)は流体の圧力(Pressure)である。ここでは、初期に、x=1.0よりも左側に圧力1(100)の流体が存在し、x=1.0よりも右側に圧力10-4の流体が存在し、その条件下で、x=1.0の地点に設けられていた隔膜を瞬間的に取り除いてから微小時間後のある時刻における流体のそれぞれの座標における圧力を示している。
【0069】
曲線802(Exact)は論理的に求めた解析解(厳密解)であり、ドット(IDO-SC)で示された本実施形態の数値計算方法による計算結果は、曲線802とほぼ一致していることがわかる。
本実施形態の数値計算方法によれば、図8(a)の場合のような一般的な衝撃波管問題だけでなく、図8(b)の場合のような圧力比などがさらに大きな衝撃波管問題や壁面との相互作用を含む問題においても、安定的かつ高精度の高い計算を行うことができる。なお、密度比や圧力比などが極端に大きな場合は、補間関数として、3次関数ではなく、分母に関数を有する有理関数、たとえば、分母に1次関数、分子に3次関数を有する有理関数などを用いることが望ましい。
【0070】
以上説明したように、本実施形態の数値計算方法によれば、コロケート格子を用いた多変数偏微分方程式において、変数間の安定なカップリングを実現することができる。また、この数値計算方法をコンピュータに実行させるためのプログラムを作成し、そのプログラムを記憶部(記録媒体)2(図1参照)に記憶(記録)することができる。
【0071】
以上で実施形態の説明を終えるが、本発明の態様はこれらに限定されるものではない。
たとえば、本発明は、電磁気に関する数値シミュレーションや、レイリーテーラー不安定性の解析など、偏微分方程式を使った自然現象の再現全般に広く適用することができる。
その他、数値計算装置や各フローチャートなどの具体的な構成について、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。
【図面の簡単な説明】
【0072】
【図1】数値計算装置の構成図である。
【図2】数値シミュレーションの処理の大まかな流れを表わしたフローチャートである。
【図3】各格子点と変数、微分係数の関係を表わした図である。
【図4】図2のステップS4における処理の内容を表わしたフローチャートである。
【図5】応用例1の計算結果を表わす図であり、(a)は、従来の計算手法で計算した場合の誤差を表わした図、(b)は、本実施形態の数値計算装置により計算した場合の誤差を表わした図である。
【図6】応用例2の計算結果を表わす図であり、散逸とETTとの関係を示した図である。
【図7】応用例3の計算結果を表わす図であり、(a)は、従来のIDO法による計算結果を示した図、(b)は、本実施形態の数値計算装置による計算結果を示した図である。
【図8】応用例4の計算結果を表わす図であり、(a)は、圧力比10:1の場合の圧力と座標の関係図、(b)は、圧力比10000:1の場合の圧力と座標の関係図である。
【符号の説明】
【0073】
1 入力部
2 記憶部
3 メモリ
4 処理部
5 出力部
10 数値計算装置
図面
【図1】
0
【図2】
1
【図3】
2
【図4】
3
【図6】
4
【図7】
5
【図8】
6
【図5】
7