数値計算方法、プログラムおよび記録媒体
【課題】保存形ＩＤＯ法を含めたＩＤＯ法に適用することで計算精度を向上させることができる数値計算方法を提供すること。
【解決手段】本発明は、流体に関する数値計算を行う局所補間微分オペレータ法を実行する数値計算装置による数値計算方法である。数値計算装置は、記憶部と演算部を有し、記憶部は、時間積分を行う変数の値と、その変数に関する補間関数とを記憶する。そして、演算部は、偏微分方程式に関して、時間軸上の固定点以外の流動点による影響も考慮した方程式解法であるラグランジュ法を移流項に適用するセミラグランジュ移流を、前記補間関数を使って前記偏微分方程式の移流項に対して行うことで前記時間積分を行う。
【解決手段】本発明は、流体に関する数値計算を行う局所補間微分オペレータ法を実行する数値計算装置による数値計算方法である。数値計算装置は、記憶部と演算部を有し、記憶部は、時間積分を行う変数の値と、その変数に関する補間関数とを記憶する。そして、演算部は、偏微分方程式に関して、時間軸上の固定点以外の流動点による影響も考慮した方程式解法であるラグランジュ法を移流項に適用するセミラグランジュ移流を、前記補間関数を使って前記偏微分方程式の移流項に対して行うことで前記時間積分を行う。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、偏微分方程式の演算に使用する数値計算方法、プログラムおよび記録媒体に関する。
【背景技術】
【０００２】
近年、流体（液体や気体など）、電磁波、熱などに関する数値計算（シミュレーション）についての研究や開発が活発になってきている。
そのような数値計算を行う場合、流体などの変数（速度、圧力、温度、密度など）に関する偏微分方程式を、離散化された空間および時間に関して解く、という手法が用いられるのが一般的である。なぜなら、偏微分方程式の数学的な解析解（厳密解）を求めることは現実的に不可能だからである。したがって、高精度の数値計算方法を開発することが重要となる。
【０００３】
そして、より高精度な計算を行うため、変数の値だけでなく変数の微分係数（空間勾配）の値も従属変数とするマルチモーメント手法であるＣＩＰ法（Cubic-Interpolated Pseudo-Particle Scheme）やＩＤＯ法（Interpolated Differential Operator Scheme：局所補間微分オペレータ法）が開発された。
【０００４】
ＩＤＯ法は、物理的に保存しなければならない質量、運動量、エネルギーを厳密に保存する保存形ＩＤＯ法と、高精度数値計算の範囲内で物理的な保存を満足する非保存形ＩＤＯ法とに大別することができる。保存形ＩＤＯ法は、たとえば、変数の区間積分値（以下、単に「積分値」という。）を従属変数とすることで、その変数に関する物理量を保存することができる。
【０００５】
そして、一般的に、計算の精度は保存形ＩＤＯ法のほうがよい。しかし、保存形ＩＤＯ法は値や積分値の定義点が異なるためにプログラムが複雑になるのに対し、非保存形ＩＤＯ法は値と微分値が同じ格子点上に定義されるため、プログラムの構造が簡単になり境界条件も取り扱い易くなる。つまり、保存形ＩＤＯ法と非保存形ＩＤＯ法は、それぞれ一長一短を有し、場面や状況に応じて選択および使用される。
【０００６】
そして、本出願人は、特許文献１において、非保存形ＩＤＯ法における変数間の安定化カップリングの技術について提案した。この技術によれば、非保存形ＩＤＯ法における数値計算の精度を向上させることができる。
【特許文献１】特開２００６−３１８３５５号公報
【発明の開示】
【発明が解決しようとする課題】
【０００７】
しかしながら、特許文献１の技術は、非保存形ＩＤＯ法に使用するものであり、保存形ＩＤＯ法に使用することはできない。つまり、保存形ＩＤＯ法の計算精度の向上には、特許文献１の技術とは別の技術が必要となる。
そこで、本発明は、前記問題点に鑑みてなされたものであり、保存形ＩＤＯ法を含めたＩＤＯ法に適用することで計算精度を向上させることができる数値計算方法、プログラムおよび記録媒体を提供することを目的とする。
【課題を解決するための手段】
【０００８】
前記課題を解決するために、本発明は、複数の定義点における１つ以上の変数とその微分係数を含んだ複数の値を用いて偏微分方程式に関する演算を行い、変数のうち少なくともいずれか１つの変数を時間積分することにより、流体に関する数値計算を行う局所補間微分オペレータ法を実行する数値計算装置による数値計算方法である。
数値計算装置は、記憶部と演算部を有し、記憶部は、時間積分を行う変数の値と、その変数に関する補間関数とを記憶する。
そして、演算部は、偏微分方程式に関して、時間軸上の固定点以外の流動点による影響も考慮した方程式解法であるラグランジュ法を移流項に適用するセミラグランジュ移流を、補間関数を使って偏微分方程式の移流項に対して行うことで時間積分を行う。
また、当該数値計算方法をコンピュータに実行させるためのプログラムを作成し、さらに、そのプログラムを記録媒体に記録することができる。
【発明の効果】
【０００９】
本発明によれば、保存形ＩＤＯ法を含めたＩＤＯ法において、計算精度を向上させることができる。
【発明を実施するための最良の形態】
【００１０】
以下、本発明の各実施形態について、適宜図面を参照しながら説明する。なお、本実施形態において、方程式の解法のうち、時間軸上の固定点に基づく解法をオイラー法といい、時間軸上の固定点以外の流動点による影響も考慮した解法をラグランジュ法という。また、偏微分方程式において、移流項にラグランジュ法を適用し、非移流項にオイラー法を適用する解法のことをセミラグランジュ（「セミラグ」と略記することもある。）法という。さらに、移流項に対してラグランジュ法を適用することを、「セミラグランジュ移流する（を行う）」と呼ぶ。ここで、移流項とは、ある物理量の局所変化を考える際に、対象物（流体など）自体の移動速度ベクトルとその物理量（またはその物理量の空間勾配）の空間勾配の内積に負の符号をつけた項のことを指す。また、移流項以外の項を非移流項という。
【００１１】
図１および図２は各実施形態の説明に共通であり、まず、図１を参照しながら、数値計算方法を行う数値計算装置の構成について説明する。図１は、数値計算装置の構成図である。
数値計算装置１０は、入力部１、記憶部２、メモリ３、演算部４および出力部５を備えて構成されるコンピュータ装置であり、たとえば、パソコン、ワークステーション、スーパーコンピュータなどにより実現することができる。
なお、演算部４は、ここでは１つとしているが、複数設けることにより並列計算を行うようにしてもよい。
【００１２】
入力部１は、数値計算装置１０の使用者（以下、単に「使用者」という）が各種情報を入力する手段であり、たとえば、キーボード、マウスなどである。
記憶部２は、各種情報を記憶するものであり、たとえば、ＲＯＭ（Read Only Memory）、ハードディスクなどである。記憶部２は、ここでは、流体の数値計算（シミュレーション）を行うために、数値計算を行う場所を示す定義点である格子点、物理量（変数）やその微分係数の初期条件、支配方程式（自然現象を表わす偏微分方程式など）、補間関数、各種プログラムなどを記憶している。
【００１３】
メモリ３は、演算部４が使用する一時記憶手段であり、たとえば、ＲＡＭ（Random Access Memory）などである。
演算部４は、記憶部２に記憶された各種情報（プログラム、データなど）などに基づいて演算（処理）を行うものであり、たとえば、ＣＰＵ（Central Processing Unit）である。
出力部５は、演算部４による演算の結果を出力するものであり、たとえば、ディスプレイなどである。
【００１４】
次に、図２を参照しながら、数値計算の演算について説明する（適宜図１参照）。図２は、流体の数値計算の演算の大まかな流れを表わしたフローチャートである。
【００１５】
まず、使用者が入力部１によって数値計算の開始に関する操作を行うと、演算部４は、記憶部２に記憶された情報に基づいて、演算を行う場所を示す定義点である格子点、各物理量の初期条件などを設定する（ステップＳ１）。
格子点は、たとえば、１次元の場合は「１２８」、２次元の場合は「３２×３２」、３次元の場合は「１６×１６×８」、といったように設定される。また、各物理量の初期条件としては、それぞれの格子点における変数の値、および、それらの値の微分係数（空間勾配）などが設定される。変数としては、たとえば、流体の速度、圧力、温度、密度などが挙げられる。
【００１６】
次に、演算部４は、記憶部２に記憶された情報に基づいて、時間と空間に関する偏微分方程式などからなる支配方程式、各物理量の値を修正するための補間関数などを設定する（ステップＳ２）。
支配方程式としては、たとえば、圧縮性流体の運動方程式、流体の質量保存方程式などが挙げられ、詳細については後記する。また、補間関数としては、２次関数、３次関数、分母に関数を有する有理関数などが挙げられ、詳細については後記する。
【００１７】
続いて、演算部４は、ステップＳ１およびステップＳ２で設定した各種情報を用い、時刻「ｔ」における各格子点での各物理量の値に基づいて、時刻「ｔ＋１」における各格子点での各物理量の値を算出する、という演算をｔ＝０,１,２,・・・ｎに関して行う（ステップＳ３〜ステップＳ５）。
そして、このステップＳ４における演算が本発明のポイントであるが、その詳細については後記する。
【００１８】
次に、演算部４は、ステップＳ３〜ステップＳ５で算出した各物理量の値を、出力部５に出力する（ステップＳ６）。以上で、数値計算が完了する。
【００１９】
（第１実施形態）
続いて、本発明の第１実施形態について説明する。第１実施形態は、保存形ＩＤＯ法を用いた数値計算の例である。第１実施形態について説明する前に、その理解を助けるため、図９、図１０Ａおよび図１０Ｂを参照しながら、保存形ＩＤＯ法を用いた数値計算の比較例（従来技術）について説明する。図９は、比較例の保存形ＩＤＯ法における変数イメージを示した図である。図１０Ａおよび図１０Ｂは、比較例の保存形ＩＤＯ法におけるルンゲクッタ法による時間更新（時間積分）を示した図である。
【００２０】
なお、第１実施形態において、移流項における微分値の算出には風上補間（２つの格子点のうち風下側（流れの下流側）の格子点を中心にして行う補間）を用い、非移流項における微分値の算出には中心補間（３つの格子点のうちの真ん中の格子点を中心にして行う補間）を用いる。
【００２１】
比較例では、支配方程式として、拡散項を有する一次元質量保存の偏微分方程式である式（１）を使用する。なお、流体に関して、ρは密度、ｕは速度、ｔは時間、ｘは座標を表す変数である。また、κ（カッパ）は拡散係数（物質の拡散の度合を示す係数）であり、あらかじめ与えられる。なお、速度ｕを求めるための式が別に必要であるが、その記載を省略し、速度ｕはその式によってすでに求められているものとする。
【数１】
【００２２】
そして、実際に数値計算を行うことができるように、式（１）を変形して式（２）とする。なお、式（２）において「i」が付与された各変数は、複数の格子点のうち「i」番目の格子点上の変数であることを意味する。また、∂ρ/∂ｘ｜iは補間関数を用いた風上補間から求められる値である（詳細は後記）。さらに、∂ｕ/∂ｘ｜iおよび∂2ρ/∂ｘ2｜iは補間関数を用いた中心補間から求められる値である（詳細は後記）。なお、本実施形態において、右側に「｜i」が付与された項は、「座標「i」における、補間関数を用いて算出される値」であることを意味する。
【数２】
【００２３】
さらに、ρの積分値を式（３）と定義し、その積分値に関する方程式を式（４）と定義する。なお、式（３）において記載されている「i+1/2」は、積分区間の両端（「i」と「i+1」）の中間地点を表している。また、式（４）において、「i+1」が付与された各変数は「i+1」番目の格子点上の変数であることを意味する。
【数３】
【００２４】
また、∂ρ/∂ｘ｜iの風上補間に用いられる補間関数を式（５）と定義する。式（６）〜（８）は、その拘束条件である。これらの式（５）〜（８）から、補間係数ａ，ｂおよびｃが式（９）〜（１１）のように求められる。また、式（１２）に示したように、ρiを微分したρx,i（＝∂ρ/∂ｘ｜i）はｂと等しくなる。
【数４】
【００２５】
このような条件の下で、図１０Ａおよび図１０Ｂに示したように、ルンゲクッタ法による時間更新（時間積分）を行う。ここでの目的は、時刻「ｎ」における密度（ρｎｉ）と各係数（ｋ1,i〜ｋ4,i）から時刻「ｎ＋１」における密度（ρn＋1ｉ）を求めること（式（１０１５）参照）と、時刻「ｎ」における密度の積分値（ρ-X,nｉ+1/2）と各係数（ｍ1,i+1/2〜ｍ4,i+1/2）から時刻「ｎ＋１」における密度（ρ-X,n+1ｉ+1/2）を求めること（式（１０１６）参照）である。式（１００１）〜式（１０１４）は、各係数（ｋ1,i〜ｋ4,i，ｍ1,i+1/2〜ｍ4,i+1/2）を求めるための式である。
【００２６】
まず、式（１００１）と式（１００２）によって、ｋ1,i（以下、「ｋ1」とも表記する。ｋ2,i〜ｋ4,iについても同様）とｍ1,i+1/2（以下、「ｍ1」とも表記する。ｍ2,i+1/2〜ｍ4,i+1/2についても同様）を求める。ここで、ｕi，ｕi+1，ρi，ρi+1の値は、時刻「ｎ」における値としてすでに求められているものを使用することができる。なお、移流項である「−ｕi∂ρ/∂ｘ｜i」における「∂ρ/∂ｘ｜i」は、前記したように風上補間によって求めることができる。
【００２７】
また、∂ｕ/∂ｘ｜iはｕの補間関数から求めることができる。ここでは、たとえば、ｕの補間関数として４次関数を使用し、中心補間によって∂ｕ/∂ｘ｜iを求めることができる。その場合、３つの格子点の座標を「i-1」，「i」および「i+1」とし、補間関数をｆ（Ｘ）、その積分値（前記式（８）と同様）をｆＸ（Ｘ）すると、その補間関数は、ｆ（i-1），ｆ（i），ｆ（i+1），ｆＸ（i-1/2）およびｆＸ（i+1/2）の拘束条件から決定することができる。∂2ρ/∂ｘ2｜iおよび非移流項中の∂ρ/∂ｘ｜iについても、同様にして中心補間から求めることができる。
【００２８】
Δｔは、時刻「ｎ」と時刻「ｎ＋１」の差の時間長（時間刻み幅）である。なお、この式（１００１）と式（１００２）によって演算が行われる段階では、すべての格子点に関して演算が行われ、その後、次の演算に進む。以降についても同様である。
【００２９】
次に、式（１００３）と式（１００４）によって、ρ*ｉ（ρn+1ｉの暫定値）とρ-X,*ｉ+1/2（ρ-X,n+1ｉ+1/2の暫定値）を求める。ここで、ρi，ρ-Xｉ+1/2は、時刻「ｎ」における値としてすでに求められているものを使用することができる。また、ｋ1とｍ1は、それぞれ式（１００１）と式（１００２）で求められた値である。
【００３０】
そして、式（１００５）と式（１００６）によって、ｋ2,iとｍ2,i+1/2を求める。ここで、ρ*ｉ（ρ*ｉ+1も同様）は、式（１００３）で求められた値である。また、移流項における∂ρ*/∂ｘ｜i（ρ*x,i）は、式（１２）からｂであることがわかり、前記した式（１０）の右辺に、式（１００３）と式（１００４）から求められたρ*ｉ，ρ*ｉ-1，ρ-X,*ｉ-1/2を代入することにより求められる。また、∂2ρ*/∂ｘ2｜iおよび非移流項中の∂ρ*/∂ｘ｜iについては、中心補間から求めることができる。
【００３１】
以下、同様にして、式（１００７）〜式（１０１４）を使うことにより、ｋ1,i〜ｋ4,i，ｍ1,i+1/2〜ｍ4,i+1/2の値をすべて求めることができる。そして、前記したように、式（１０１５）と式（１０１６）を使って、ρn＋1ｉとρ-X,n+1ｉ+1/2を求めることができる。なお、式（１００７）と式（１０１１）においてρ*ｉは新たに更新され、式（１００８）と式（１０１２）においてρ-X,*ｉ+1/2は新たに更新される。
【００３２】
次に、図３〜図６を参照しながら、本発明の第１実施形態の部分セミラグ保存形ＩＤＯ法（セミラグランジュ法を部分的に適用した保存形ＩＤＯ法）について説明する。
第１実施形態では、計算に、前記した比較例の場合と同じ式（１）〜式（１１）を使用するので、重複説明を適宜省略する。なお、第１実施形態は、式（２）における右辺の移流項である第１項の「−ｕi∂ρ/∂ｘ｜i」に関してセミラグランジュ移流を行う点で、比較例における計算と異なっている。
【００３３】
第１実施形態では、式（５）に対して、Ｘ＝ｘi−ｕiΔｔを代入することで、Δｔの分だけセミラグランジュ移流したρの値であるρSLi,Δtを算出することができる（式（１３）参照）。なお、式（１３）において、Ｆni（Ｘ），ａnおよびｂnは、時刻「ｎ」におけるそれぞれの値であることを意味している。また、式（１３）におけるΔｔをΔｔ／２と置き換えることで、Δｔ／２の分だけセミラグランジュ移流したρの値であるρSLi,Δt/2を算出することができる。
【数５】
【００３４】
ここで、図３は、第１実施形態の部分セミラグ保存形ＩＤＯ法における変数イメージを示した図である。図３に示したように、このρSLi,Δtの値を用いて、時刻「ｎ＋１」における密度（ρn＋1ｉ）を求めることができる（詳細は後記）。
【００３５】
次に、図４Ａ〜図４Ｃを参照しながら、第１実施形態における数値計算装置の演算（処理）について説明する（適宜他図参照）。図４Ａ〜図４Ｃは、第１実施形態における数値計算装置の演算を示したフローチャートである。数値計算装置１０の演算部４は、記憶部２とメモリ３を使用しながら各演算を行う。そして、この演算の目的は、比較例の場合と同様、時刻「ｎ」における密度（ρｉ）やその積分値（ρ-Xｉ+1/2（ρ-X,nｉ+1/2））などに基づいて、時刻「ｎ＋１」における密度（ρn+1ｉ）と積分値（ρ-X,n+1ｉ+1/2）を求めることである（ステップＳ４１０参照）。つまり、比較例の場合と同様、時刻「ｎ」におけるρi，ρ-Xｉ+1/2，ｕi，ｕi+1，ρi，ρi+1などの値はすでに求められているものとして、以下の演算を行う。また、∂ｕ/∂ｘ｜i，∂2ρ/∂ｘ2｜i，∂ρ/∂ｘ｜iなどは、前記した比較例の場合と同様、中心補間によって求められる。
【００３６】
演算部４は、式（２）を、右辺の移流項である第１項の「−ｕi∂ρ/∂ｘ｜i」を削除して（省いて）使用し（式（４０１１）参照）、また、式（４）はそのまま使用する（式（４０１２）参照）（ステップＳ４０１）。
次に、演算部４は、式（１３）を用いて、セミラグランジュ移流させた値であるρSLi,ΔtおよびρSLi,Δt/2を取得する（ステップＳ４０２）。
【００３７】
続いて、演算部４は、式（４０３１）と式（４０３２）によって、ｋ1,i（ｋ1）とｍ1,i+1/2（ｍ1）を求める。ここで、∂ｕ/∂ｘ｜iはｕの補間関数（説明を省略）から求めることができる。Δｔは、時刻「ｎ」と時刻「ｎ＋１」の差の時間長である。
【００３８】
次に、演算部４は、式（４０４１）と式（４０４２）によって、ρ*ｉ（ρn+1ｉの暫定値）とρ-X,*ｉ+1/2（ρ-X,n+1ｉ+1/2の暫定値）を取得する（ステップＳ４０４）。ここで、ｋ1とｍ1はステップＳ４０３で取得された値であり、ρSLi,Δt/2はステップＳ４０２で取得された値である。
【００３９】
続いて、演算部４は、式（４０５１）と式（４０５２）によって、ｋ2,iとｍ2,i+1/2を取得する（ステップＳ４０５）。ここで、ρ*ｉ（ρ*ｉ+1も同様）は、ステップＳ４０４で取得された値である。また、∂ｕ/∂ｘ｜i，∂2ρ*/∂ｘ2｜i，∂ρ*/∂ｘ｜iなどは、中心補間によって求められる。
【００４０】
以下、演算部４は、同様にして、ステップＳ４０６〜ステップＳ４０９の演算を行うことによって、ｋ1,i〜ｋ4,i，ｍ1,i+1/2〜ｍ4,i+1/2のすべての値を取得することができる。なお、ステップＳ４０６およびステップＳ４０８において、ρ*ｉおよびρ-X,*ｉ+1/2はそれぞれ新たに更新される。
そして、演算部４は、ステップＳ４０９までの演算で得た各値と式（４１０１）と式（４１０２）によって、ρn＋1ｉとρ-X,n+1ｉ+1/2を取得する（ステップＳ４１０）。
このようにして、第１実施形態の数値計算装置１０によれば、保存形ＩＤＯ法における偏微分方程式に対して部分的に（２つのうち片方の方程式に）セミラグランジュ法を適用することで、時間更新を行うことができる。
【００４１】
次に、図５および図６を参照しながら、第１実施形態の数値計算装置による数値計算方法の計算結果例について説明する。図５は、第１実施形態の数値計算装置による数値計算方法の計算結果例を示すグラフである。ここでは、実際に爆発実験を行い、大気圧に対する加圧（量）に関して、実測値（実験データ）、第１実施形態による計算結果、および、比較例による計算結果の三者を比較した。
【００４２】
図５に示すグラフにおいて、横軸は爆源からの距離（メートル）、縦軸は大気圧をゼロとした場合の加圧（パスカル）を表している。符号５１の破線は、実験で得られた各距離の地点におけるピーク加圧の実測値を示す。符号５１の破線を見れば、爆発による衝撃波が少しずつ弱まりながら遠方に伝搬していることがわかる。
【００４３】
符号５２の実線は、第１実施形態の部分セミラグ保存形ＩＤＯ法により計算したある瞬間の衝撃波プロファイルである。そのピーク値が符号５１の破線とほぼ一致しており、計算結果の精度が高いことがわかる。なお、この衝撃波プロファイルは、時間の経過とともに、その形をほぼ維持したまま右方に平行移動することになる。
一方、符号５３の二点鎖線は、比較例の非保存形ＩＤＯ法により計算したある瞬間の衝撃波プロファイルを示す。そのピーク値が符号５１の破線と大きく離れており、計算結果の精度が第１実施形態の部分セミラグ保存形ＩＤＯ法の場合よりも低いことがわかる。
【００４４】
同じ実験に関する別のグラフについて、図６を参照しながら説明する。図６は、第１実施形態の数値計算装置による数値計算方法の計算結果例を示す両対数グラフである。ここでも、図５の場合と同様に、加圧に関して、実測値（実験データ）、第１実施形態による計算結果、および、比較例による計算結果の三者を比較した。
【００４５】
図６に示す両対数グラフにおいて、横軸は爆源からの距離（メートル）、縦軸は大気圧をゼロとした場合のピーク加圧（パスカル）を表している。符号６１の破線は、実験で得られた各距離の地点におけるピーク加圧の実測値を示す。符号６１の破線を見れば、爆発による衝撃波が弱まりながら遠方に伝搬していることがわかる。
【００４６】
符号６２の実線は、第１実施形態の部分セミラグ保存形ＩＤＯ法により計算したピーク加圧に関する衝撃波プロファイルを示す。符号６２の実線は、符号６１の破線とほぼ一致しており、計算結果の精度が高いことがわかる。
一方、符号６３の二点鎖線は、比較例の非保存形ＩＤＯ法により計算したピーク加圧に関する衝撃波プロファイルである。符号６３の二点鎖線は、符号６１の破線と大きく離れており、計算結果の精度が第１実施形態の部分セミラグ保存形ＩＤＯ法の場合よりも低いことがわかる。
【００４７】
つまり、図５および図６からわかるように、第１実施形態の部分セミラグ保存形ＩＤＯ法によれば、比較例の非保存形ＩＤＯ法の場合に比べて、爆源から近距離（たとえば１０ｍ程度）に関しても高い精度での計算結果を得ることができ、また、爆源から遠距離（たとえば２００ｍ以上）に関してはさらに高い精度での計算結果を得ることができる。
【００４８】
このような第１実施形態の部分セミラグ保存形ＩＤＯ法を応用すれば、車体周りの空気や熱の流れ、風力発電の空力解析、気象計算、津波予測など、多くの流体計算を高精度で行うことができる。また、爆発現象や、水と空気が混在する流れのような、空間的な密度変化が１０００倍以上になる流体計算が可能となる。さらに、高速流体（高レイノルズ数流れ）で頻繁に発生する乱流に対し、定義された格子間隔よりも小さな渦による影響を渦粘性項無しでシミュレーションする、いわゆる陰的なラージ・エディ・シミュレーションになっていて、特別な乱流モデルを導入すること無しに複雑な流体計算を行うことができる。
【００４９】
さらに説明すると、一般に、多くの利点を有している保存形ＩＤＯ法の時間積分は、ＣＩＰ法のようなセミラグランジュ法ではなく、オイラー法で解くためにルンゲクッタ法などが用いられる。しかし、通常のルンゲクッタ法を用いてアンダーシュート（プロファイルにおける下方への飛び出し誤差）やオーバーシュート（プロファイルにおける上方への飛び出し誤差）を抑制することは非常に難しい。そこで、ルンゲクッタ法による時間積分の各段階において、保存量となる積分値については従来の保存形ＩＤＯ法による計算を行い、保存が保証されない格子点上の変数の時間積分に関する偏微分方程式を移流項と非移流項に分割し、セミラグランジュ法を適用し、さらに望ましくは後記するような単調性を有する補間関数を使用することで、アンダーシュートやオーバーシュートを抑制し、高い計算精度を実現することができる。
【００５０】
（第２実施形態）
次に、本発明の第２実施形態について説明する。第２実施形態は、非保存形ＩＤＯ法を用いた数値計算の例である。第２実施形態について説明する前に、その理解を助けるため、図１１、図１２Ａおよび図１２Ｂを参照しながら、非保存形ＩＤＯ法を用いた数値計算の比較例（従来技術）について説明する。図１１は、比較例の非保存形ＩＤＯ法における変数イメージを示した図である。図１２Ａおよび図１２Ｂは、比較例の非保存形ＩＤＯ法におけるルンゲクッタ法による時間更新（時間積分）を示した図である。なお、以下、第１実施形態と重複する説明については適宜省略する。また、第２実施形態において、移流項における微分値の算出には風上補間を用い、非移流項における２階以上の微分値の算出には中心補間を用いる。
【００５１】
比較例では、支配方程式として、拡散項を有する一次元質量保存の偏微分方程式である式（１４）（式（１）と同様）を使用する。つまり、流体に関して、ρは密度、ｕは速度、ｔは時間、ｘは座標を表す変数である。なお、速度ｕを求めるための式が別に必要であるが、その記載を省略し、ｕ，ｕx（ｕの一階微分）の値はその式などからすでに求められているものとする。
【数６】
【００５２】
そして、実際に数値計算を行うことができるように、式（１４）を変形して式（１５）とする。
【数７】
【００５３】
さらに、ρの微分値を式（１６）と定義し、その微分値に関する方程式を式（１７）と定義する。なお、∂ρx/∂ｘ｜iは、後記する補間関数（式（１８））による風上補間によって求められる値である。また、非移流項における∂2ｕ/∂ｘ2｜iおよび∂3ρ/∂ｘ3｜iは、中心補間によって求められる。
【数８】
【００５４】
また、ρに関する補間関数を式（１８）と定義する。式（１９）〜（２２）は、その拘束条件である。これらの式（１９）〜（２２）から、補間係数ａ〜ｄが式（２３）〜（２６）のように求められる。また、式（２７）に示すように、ρiを二階微分した∂ρx/∂ｘ｜iは２ｂと等しくなる。
【数９】
【００５５】
このような条件の下で、図１２Ａおよび図１２Ｂに示したように、ルンゲクッタ法による時間更新（時間積分）を行う。ここでの目的は、時刻「ｎ」における密度（ρｎｉ）と各係数（ｋ1,i〜ｋ4,i）から時刻「ｎ＋１」における密度（ρn＋1ｉ）を求めること（式（１２１５）参照）と、時刻「ｎ」における密度の微分値（ρnx,i）と各係数（ｍ1,i〜ｍ4,i）から時刻「ｎ＋１」における密度（ρn＋1x,i）を求めること（式（１２１６）参照）である。式（１２０１）〜式（１２１４）は、各係数（ｋ1,i〜ｋ4,i，ｍ1,i〜ｍ4,i）を求めるための式である。
【００５６】
まず、式（１２０１）と式（１２０２）によって、ｋ1,i（以下、「ｋ1」とも表記する。ｋ2,i〜ｋ4,iについても同様）とｍ1,i（以下、「ｍ1」とも表記する。ｍ2,i〜ｍ4,iについても同様）を求める。ここで、ｕi，ｕx,i，ρiおよびρx,iは、時刻「ｎ」における値としてすでに求められている値を使用することができる。また、前記したように、移流項における∂ρ/∂ｘ｜iと∂ρx/∂ｘ｜iは風上補間により求めることができる。さらに、非移流項における∂2ρ/∂ｘ2｜i，∂2ｕ/∂ｘ2｜iおよび∂3ρ/∂ｘ3｜iは、中心補間によって求められる。Δｔは、時刻「ｎ」と時刻「ｎ＋１」の差の時間長である。
【００５７】
次に、式（１２０３）と式（１２０４）によって、ρ*ｉ（ρn+1ｉの暫定値）とρ*x,i（ρn+1x,iの暫定値）を求める。ここで、ρiとρx,iは、時刻「ｎ」における値としてすでに求められている値を使用することができる。また、ｋ1とｍ1は、それぞれ式（１２０１）と式（１２０２）で求められた値である。
【００５８】
そして、式（１２０５）と式（１２０６）によって、ｋ2,iとｍ2,iを求める。ここで、ρ*ｉは、式（１２０３）で求められた値である。また、∂ρ*/∂ｘ｜iは、式（１２０４）で求められたρ*x,iである。
∂ρ*x/∂ｘ｜iは、式（２７）と式（２４）から求められる値であり（式（２７）の∂ρx/∂ｘ｜iに相当）、具体的には、式（２４）において式（１２０３）と式（１２０４）から求められたρ*ｉ，ρ*ｉ-1，ρ*x,iおよびρ*x,ｉ-1を代入し、それによって得られたｂに２を乗算することで求められる。また、非移流項における∂2ρ*/∂ｘ2｜i，∂2ｕ/∂ｘ2｜iおよび∂3ρ*/∂ｘ3｜iは、中心補間によって求められる。
【００５９】
以下、同様にして、式（１２０７）〜式（１２１４）を使うことにより、ｋ1,i〜ｋ4,i，ｍ1,i〜ｍ4,iの値をすべて求めることができる。そして、前記したように、式（１２１５）と式（１２１６）を使って、ρn＋1ｉとρn＋1x,ｉを求めることができる。なお、式（１２０７）と式（１２１１）においてρ*ｉは新たに更新され、式（１２０８）と式（１２１２）においてρ*x,ｉは新たに更新される。
【００６０】
次に、図７および図８Ａ〜図８Ｃを参照しながら、本発明の第２実施形態であるセミラグ非保存形ＩＤＯ法（セミラグランジュ法を適用した非保存形ＩＤＯ法）について説明する。
第２実施形態では、計算に、前記した比較例の場合と同じ式（１４）〜式（２６）を使用するので、重複説明を適宜省略する。なお、第２実施形態は、式（１５）における右辺の移流項である第１項の「−ｕiρx,i」と、式（１７）における右辺の移流項である第１項の「−ｕi∂ρx/∂ｘ｜i」とに関してセミラグランジュ移流を行う点で、比較例における計算と異なっている。
【００６１】
第２実施形態では、式（１８）に対して、Ｘ＝ｘi−ｕiΔｔを代入することで、Δｔの分だけセミラグランジュ移流したρの値であるρSLi,Δtを算出することができる（式（２８）参照）。なお、Ｆni（Ｘ），ａn，ｂn，ｃnおよびｄnは、時刻「ｎ」におけるそれぞれの値であることを意味している。また、式（２８）におけるΔｔをΔｔ／２と置き換えることで、Δｔ／２の分だけセミラグランジュ移流したρの値であるρSLi,Δt/2を算出することができる。
【数１０】
【００６２】
また、式（１８）を微分した式（説明を省略）に対して、Ｘ＝ｘi−ｕiΔｔを代入することで、Δｔの分だけセミラグランジュ移流したρx,iの値であるρSLx,i,Δtを算出することができる（式（２９）参照）。また、式（２９）におけるΔｔをΔｔ／２と置き換えることで、Δｔ／２の分だけセミラグランジュ移流したρx,iの値であるρSLx,i,Δt/2を算出することができる。
【数１１】
【００６３】
ここで、図７は、第２実施形態のセミラグ非保存形ＩＤＯ法における変数イメージを示した図である。図７に示したように、このρSLi,Δt（ρSLi,Δt/2）およびρSLx,i,Δt（ρSLx,i,Δt/2）の値を用いて、時刻「ｎ＋１」における密度（ρn＋1ｉ）を求めることができる（詳細は後記）。
【００６４】
次に、図８Ａ〜図８Ｃを参照しながら、第２実施形態における数値計算装置の演算（処理）について説明する（適宜他図参照）。図８Ａ〜図８Ｃは、第２実施形態における数値計算装置の演算を示したフローチャートである。数値計算装置１０の演算部４は、記憶部２とメモリ３を使用しながら各演算を行う。そして、この演算の目的は、比較例の場合と同様、時刻「ｎ」における密度（ρｉ）やその微分値（ρx,ｉ）などに基づいて、時刻「ｎ＋１」における密度（ρn+1ｉ）と微分値（ρn+1x,ｉ）を求めることである（ステップＳ４３０参照）。つまり、比較例の場合と同様、時刻「ｎ」におけるρi，ρx,i，ｕx,iなどの値はすでに求められているものとして以下の演算を行う。なお、非移流項における∂2ρ/∂ｘ2｜i，∂2ｕ/∂ｘ2｜iおよび∂3ρ/∂ｘ3｜iは、中心補間によって求められる。
【００６５】
演算部４は、式（１５）について右辺の移流項である第１項の「−ｕiρx,i」を削除して（省いて）使用し（式（４２１１）参照）、また、式（１７）について右辺の移流項である第１項の「−ｕi∂ρx/∂ｘ｜i」を削除して（省いて）使用する（式（４２１２）参照）（ステップＳ４２１）。
次に、演算部４は、式（２８）および式（２９）を用いて、セミラグランジュ移流させた値であるρSLi,Δt，ρSLi,Δt/2，ρSLx,i,ΔtおよびρSLx,i,Δt/2を取得する（ステップＳ４２２）。
【００６６】
続いて、演算部４は、式（４２３１）と式（４２３２）によって、ｋ1,i（ｋ1）とｍ1,i（ｍ1）を求める。ここで、∂2ρ/∂ｘ2｜i，∂2ｕ/∂ｘ2｜iおよび∂3ρ/∂ｘ3｜iは、中心補間によって求められる。Δｔは、時刻「ｎ」と時刻「ｎ＋１」の差の時間長である。
【００６７】
次に、演算部４は、式（４２４１）と式（４２４２）によって、ρ*ｉ（ρn+1ｉの暫定値）とρ*x,i（ρn+1x,iの暫定値）を取得する（ステップＳ４２４）。ここで、ｋ1とｍ1はステップＳ４２３で取得された値であり、ρSLi,Δt/2およびρSLx,i,Δt/2はステップＳ４２２で取得された値である。
【００６８】
続いて、演算部４は、式（４２５１）と式（４２５２）によって、ｋ2,iとｍ2,iを取得する（ステップＳ４２５）。ここで、ρ*ｉおよびρ*x,ｉは、ステップＳ４２４で取得された値である。また、∂2ρ*/∂ｘ2｜i，∂2ｕ/∂ｘ2｜iおよび∂3ρ*/∂ｘ3｜iは、中心補間によって求められる。
【００６９】
以下、演算部４は、同様にして、ステップＳ４２６〜ステップＳ４２９の演算を行うことによって、ｋ1,i〜ｋ4,i，ｍ1,i〜ｍ4,iのすべての値を取得することができる。なお、ステップＳ４２６およびステップＳ４２８において、ρ*ｉおよびρ*x,ｉはそれぞれ新たに更新される。
そして、演算部４は、ステップＳ４２９までの演算で得た各値と式（４３０１）と式（４３０２）によって、ρn＋1ｉとρn+1x,iを取得する（ステップＳ４３０）。
【００７０】
このようにして、第２実施形態の数値計算装置１０によれば、非保存形ＩＤＯ法における偏微分方程式に対してセミラグランジュ法を適用することで、時間更新を行うことができる。また、その詳細な効果についての説明を省略するが、図５および図６を参照して説明した第１実施形態の場合と同様の効果を得ることができる。
【００７１】
（変形例）
第１実施形態では、補間関数として、式（５）に示した式を使用した場合について説明したが、それ以外に、式（３０）に示す分母に関数を有する有理関数を使用してもよい。β，ａおよびｃは、拘束条件により決定される係数である。
【数１２】
【００７２】
また、第２実施形態では、補間関数として、式（１８）に示した式を使用した場合について説明したが、それ以外に、式（３１）に示す分母に関数を有する有理関数を使用してもよい。ａ，ｂ，ｃおよびｄは、拘束条件により決定される係数である。
【数１３】
そして、αとＢは、数値計算の諸条件により、式（３１）が単調性を有するように決定すればよい。たとえば、α＝１とし、Ｂを式（３２）および式（３３）により決定すればよい。
【数１４】
【００７３】
このような式（３０）や式（３１）に示すように、単調性を有する（単調増加または単調減少する）補間関数を使用すれば、密度比や圧力比などが極端に大きな（たとえば１０００倍以上）現象について数値計算を行う場合でも、その計算精度を高く維持することができる。つまり、密度比や圧力比などが極端に大きな現象に関する数値計算を行う場合、変数の数値の大きな場面での誤差の影響が大きくなる。そのため、単調性を有する補間関数を使用することで、急激な空間勾配のプロファイルに対してアンダーシュートやオーバーシュートをなくし、計算を安定させ、高い精度を維持することができるようになる。
【００７４】
図１３の（ａ）は、前記した２つの比較例においてκ（拡散係数）＝０の場合の計算結果のプロファイルの例を示した図である。図１３の（ａ）に示すように、比較例による計算結果（ドット）は、解析解（実線）に対してアンダーシュートやオーバーシュートを起こすことがあった。つまり、比較例では、密度比や圧力比などが極端に大きな現象についての数値計算が破綻することがあった。
【００７５】
図１３の（ｂ）は、本実施形態（第１実施形態および第２実施形態に前記変形例の補間関数を適用したもの）においてκ（拡散係数）＝０の場合の計算結果のプロファイルの例を示した図である。図１３の（ｂ）に示すように、本実施形態において上記のような有理関数を用いた場合の計算結果（ドット）は、解析解（実線）に対してアンダーシュートやオーバーシュートを起こすことがなくなる。つまり、本実施形態では、密度比や圧力比などが極端に大きな現象についての数値計算が安定するという効果を奏する。
なお、単調性を有する補間関数として、式（３０）および式（３１）を示したが、これらに限定されるものではない。
【００７６】
また、当業者であれば、前記した各実施形態とその変形例にかかる数値計算方法をコンピュータに実行させるためのプログラムを作成し、さらに、そのプログラムをＣＤ（Compact Disc）やＤＶＤ（Digital Versatile Disk）などの記録媒体に記録することができる。
【００７７】
以上で本実施形態の説明を終えるが、本発明の態様はこれらに限定されるものではない。
たとえば、本実施形態では、説明を簡潔にするため、一次元空間の偏微分方程式について説明したが、２次元以上の空間の偏微分方程式を使用してもかまわない。
また、流体の変数として、密度を中心に説明したが、速度、圧力、温度などの別の変数に関しても同様に、本発明のセミラグランジュ移流による数値計算を適用することができる。
【００７８】
さらに、時間積分のルンゲクッタ法については、４段で行う手法について説明したが、２段や３段など別の手法で行うようにしてもよい。
その他、数値計算装置や各フローチャートなどの具体的な構成について、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。
【図面の簡単な説明】
【００７９】
【図１】数値計算装置の構成図である。
【図２】流体の数値計算の演算の大まかな流れを表わしたフローチャートである。
【図３】第１実施形態の部分セミラグ保存形ＩＤＯ法における変数イメージを示した図である。
【図４Ａ】第１実施形態における数値計算装置の演算を示したフローチャートである。
【図４Ｂ】第１実施形態における数値計算装置の演算を示したフローチャートである。
【図４Ｃ】第１実施形態における数値計算装置の演算を示したフローチャートである。
【図５】第１実施形態の数値計算装置による数値計算方法の計算結果例を示すグラフである。
【図６】第１実施形態の数値計算装置による数値計算方法の計算結果例を示す両対数グラフである。
【図７】第２実施形態のセミラグ非保存形ＩＤＯ法における変数イメージを示した図である。
【図８Ａ】第２実施形態における数値計算装置の演算を示したフローチャートである。
【図８Ｂ】第２実施形態における数値計算装置の演算を示したフローチャートである。
【図８Ｃ】第２実施形態における数値計算装置の演算を示したフローチャートである。
【図９】比較例の保存形ＩＤＯ法における変数イメージを示した図である。
【図１０Ａ】比較例の保存形ＩＤＯ法におけるルンゲクッタ法による時間更新（時間積分）を示した図である。
【図１０Ｂ】比較例の保存形ＩＤＯ法におけるルンゲクッタ法による時間更新（時間積分）を示した図である。
【図１１】比較例の非保存形ＩＤＯ法における変数イメージを示した図である。
【図１２Ａ】比較例の非保存形ＩＤＯ法におけるルンゲクッタ法による時間更新（時間積分）を示した図である。
【図１２Ｂ】比較例の非保存形ＩＤＯ法におけるルンゲクッタ法による時間更新（時間積分）を示した図である。
【図１３】（ａ）は、比較例による計算結果のプロファイルの例を示した図であり、（ｂ）は、本実施形態による計算結果のプロファイルの例を示した図である。
【符号の説明】
【００８０】
１ 入力部
２ 記憶部
３ メモリ
４ 演算部
５ 出力部
１０ 数値計算装置
【技術分野】
【０００１】
本発明は、偏微分方程式の演算に使用する数値計算方法、プログラムおよび記録媒体に関する。
【背景技術】
【０００２】
近年、流体（液体や気体など）、電磁波、熱などに関する数値計算（シミュレーション）についての研究や開発が活発になってきている。
そのような数値計算を行う場合、流体などの変数（速度、圧力、温度、密度など）に関する偏微分方程式を、離散化された空間および時間に関して解く、という手法が用いられるのが一般的である。なぜなら、偏微分方程式の数学的な解析解（厳密解）を求めることは現実的に不可能だからである。したがって、高精度の数値計算方法を開発することが重要となる。
【０００３】
そして、より高精度な計算を行うため、変数の値だけでなく変数の微分係数（空間勾配）の値も従属変数とするマルチモーメント手法であるＣＩＰ法（Cubic-Interpolated Pseudo-Particle Scheme）やＩＤＯ法（Interpolated Differential Operator Scheme：局所補間微分オペレータ法）が開発された。
【０００４】
ＩＤＯ法は、物理的に保存しなければならない質量、運動量、エネルギーを厳密に保存する保存形ＩＤＯ法と、高精度数値計算の範囲内で物理的な保存を満足する非保存形ＩＤＯ法とに大別することができる。保存形ＩＤＯ法は、たとえば、変数の区間積分値（以下、単に「積分値」という。）を従属変数とすることで、その変数に関する物理量を保存することができる。
【０００５】
そして、一般的に、計算の精度は保存形ＩＤＯ法のほうがよい。しかし、保存形ＩＤＯ法は値や積分値の定義点が異なるためにプログラムが複雑になるのに対し、非保存形ＩＤＯ法は値と微分値が同じ格子点上に定義されるため、プログラムの構造が簡単になり境界条件も取り扱い易くなる。つまり、保存形ＩＤＯ法と非保存形ＩＤＯ法は、それぞれ一長一短を有し、場面や状況に応じて選択および使用される。
【０００６】
そして、本出願人は、特許文献１において、非保存形ＩＤＯ法における変数間の安定化カップリングの技術について提案した。この技術によれば、非保存形ＩＤＯ法における数値計算の精度を向上させることができる。
【特許文献１】特開２００６−３１８３５５号公報
【発明の開示】
【発明が解決しようとする課題】
【０００７】
しかしながら、特許文献１の技術は、非保存形ＩＤＯ法に使用するものであり、保存形ＩＤＯ法に使用することはできない。つまり、保存形ＩＤＯ法の計算精度の向上には、特許文献１の技術とは別の技術が必要となる。
そこで、本発明は、前記問題点に鑑みてなされたものであり、保存形ＩＤＯ法を含めたＩＤＯ法に適用することで計算精度を向上させることができる数値計算方法、プログラムおよび記録媒体を提供することを目的とする。
【課題を解決するための手段】
【０００８】
前記課題を解決するために、本発明は、複数の定義点における１つ以上の変数とその微分係数を含んだ複数の値を用いて偏微分方程式に関する演算を行い、変数のうち少なくともいずれか１つの変数を時間積分することにより、流体に関する数値計算を行う局所補間微分オペレータ法を実行する数値計算装置による数値計算方法である。
数値計算装置は、記憶部と演算部を有し、記憶部は、時間積分を行う変数の値と、その変数に関する補間関数とを記憶する。
そして、演算部は、偏微分方程式に関して、時間軸上の固定点以外の流動点による影響も考慮した方程式解法であるラグランジュ法を移流項に適用するセミラグランジュ移流を、補間関数を使って偏微分方程式の移流項に対して行うことで時間積分を行う。
また、当該数値計算方法をコンピュータに実行させるためのプログラムを作成し、さらに、そのプログラムを記録媒体に記録することができる。
【発明の効果】
【０００９】
本発明によれば、保存形ＩＤＯ法を含めたＩＤＯ法において、計算精度を向上させることができる。
【発明を実施するための最良の形態】
【００１０】
以下、本発明の各実施形態について、適宜図面を参照しながら説明する。なお、本実施形態において、方程式の解法のうち、時間軸上の固定点に基づく解法をオイラー法といい、時間軸上の固定点以外の流動点による影響も考慮した解法をラグランジュ法という。また、偏微分方程式において、移流項にラグランジュ法を適用し、非移流項にオイラー法を適用する解法のことをセミラグランジュ（「セミラグ」と略記することもある。）法という。さらに、移流項に対してラグランジュ法を適用することを、「セミラグランジュ移流する（を行う）」と呼ぶ。ここで、移流項とは、ある物理量の局所変化を考える際に、対象物（流体など）自体の移動速度ベクトルとその物理量（またはその物理量の空間勾配）の空間勾配の内積に負の符号をつけた項のことを指す。また、移流項以外の項を非移流項という。
【００１１】
図１および図２は各実施形態の説明に共通であり、まず、図１を参照しながら、数値計算方法を行う数値計算装置の構成について説明する。図１は、数値計算装置の構成図である。
数値計算装置１０は、入力部１、記憶部２、メモリ３、演算部４および出力部５を備えて構成されるコンピュータ装置であり、たとえば、パソコン、ワークステーション、スーパーコンピュータなどにより実現することができる。
なお、演算部４は、ここでは１つとしているが、複数設けることにより並列計算を行うようにしてもよい。
【００１２】
入力部１は、数値計算装置１０の使用者（以下、単に「使用者」という）が各種情報を入力する手段であり、たとえば、キーボード、マウスなどである。
記憶部２は、各種情報を記憶するものであり、たとえば、ＲＯＭ（Read Only Memory）、ハードディスクなどである。記憶部２は、ここでは、流体の数値計算（シミュレーション）を行うために、数値計算を行う場所を示す定義点である格子点、物理量（変数）やその微分係数の初期条件、支配方程式（自然現象を表わす偏微分方程式など）、補間関数、各種プログラムなどを記憶している。
【００１３】
メモリ３は、演算部４が使用する一時記憶手段であり、たとえば、ＲＡＭ（Random Access Memory）などである。
演算部４は、記憶部２に記憶された各種情報（プログラム、データなど）などに基づいて演算（処理）を行うものであり、たとえば、ＣＰＵ（Central Processing Unit）である。
出力部５は、演算部４による演算の結果を出力するものであり、たとえば、ディスプレイなどである。
【００１４】
次に、図２を参照しながら、数値計算の演算について説明する（適宜図１参照）。図２は、流体の数値計算の演算の大まかな流れを表わしたフローチャートである。
【００１５】
まず、使用者が入力部１によって数値計算の開始に関する操作を行うと、演算部４は、記憶部２に記憶された情報に基づいて、演算を行う場所を示す定義点である格子点、各物理量の初期条件などを設定する（ステップＳ１）。
格子点は、たとえば、１次元の場合は「１２８」、２次元の場合は「３２×３２」、３次元の場合は「１６×１６×８」、といったように設定される。また、各物理量の初期条件としては、それぞれの格子点における変数の値、および、それらの値の微分係数（空間勾配）などが設定される。変数としては、たとえば、流体の速度、圧力、温度、密度などが挙げられる。
【００１６】
次に、演算部４は、記憶部２に記憶された情報に基づいて、時間と空間に関する偏微分方程式などからなる支配方程式、各物理量の値を修正するための補間関数などを設定する（ステップＳ２）。
支配方程式としては、たとえば、圧縮性流体の運動方程式、流体の質量保存方程式などが挙げられ、詳細については後記する。また、補間関数としては、２次関数、３次関数、分母に関数を有する有理関数などが挙げられ、詳細については後記する。
【００１７】
続いて、演算部４は、ステップＳ１およびステップＳ２で設定した各種情報を用い、時刻「ｔ」における各格子点での各物理量の値に基づいて、時刻「ｔ＋１」における各格子点での各物理量の値を算出する、という演算をｔ＝０,１,２,・・・ｎに関して行う（ステップＳ３〜ステップＳ５）。
そして、このステップＳ４における演算が本発明のポイントであるが、その詳細については後記する。
【００１８】
次に、演算部４は、ステップＳ３〜ステップＳ５で算出した各物理量の値を、出力部５に出力する（ステップＳ６）。以上で、数値計算が完了する。
【００１９】
（第１実施形態）
続いて、本発明の第１実施形態について説明する。第１実施形態は、保存形ＩＤＯ法を用いた数値計算の例である。第１実施形態について説明する前に、その理解を助けるため、図９、図１０Ａおよび図１０Ｂを参照しながら、保存形ＩＤＯ法を用いた数値計算の比較例（従来技術）について説明する。図９は、比較例の保存形ＩＤＯ法における変数イメージを示した図である。図１０Ａおよび図１０Ｂは、比較例の保存形ＩＤＯ法におけるルンゲクッタ法による時間更新（時間積分）を示した図である。
【００２０】
なお、第１実施形態において、移流項における微分値の算出には風上補間（２つの格子点のうち風下側（流れの下流側）の格子点を中心にして行う補間）を用い、非移流項における微分値の算出には中心補間（３つの格子点のうちの真ん中の格子点を中心にして行う補間）を用いる。
【００２１】
比較例では、支配方程式として、拡散項を有する一次元質量保存の偏微分方程式である式（１）を使用する。なお、流体に関して、ρは密度、ｕは速度、ｔは時間、ｘは座標を表す変数である。また、κ（カッパ）は拡散係数（物質の拡散の度合を示す係数）であり、あらかじめ与えられる。なお、速度ｕを求めるための式が別に必要であるが、その記載を省略し、速度ｕはその式によってすでに求められているものとする。
【数１】
【００２２】
そして、実際に数値計算を行うことができるように、式（１）を変形して式（２）とする。なお、式（２）において「i」が付与された各変数は、複数の格子点のうち「i」番目の格子点上の変数であることを意味する。また、∂ρ/∂ｘ｜iは補間関数を用いた風上補間から求められる値である（詳細は後記）。さらに、∂ｕ/∂ｘ｜iおよび∂2ρ/∂ｘ2｜iは補間関数を用いた中心補間から求められる値である（詳細は後記）。なお、本実施形態において、右側に「｜i」が付与された項は、「座標「i」における、補間関数を用いて算出される値」であることを意味する。
【数２】
【００２３】
さらに、ρの積分値を式（３）と定義し、その積分値に関する方程式を式（４）と定義する。なお、式（３）において記載されている「i+1/2」は、積分区間の両端（「i」と「i+1」）の中間地点を表している。また、式（４）において、「i+1」が付与された各変数は「i+1」番目の格子点上の変数であることを意味する。
【数３】
【００２４】
また、∂ρ/∂ｘ｜iの風上補間に用いられる補間関数を式（５）と定義する。式（６）〜（８）は、その拘束条件である。これらの式（５）〜（８）から、補間係数ａ，ｂおよびｃが式（９）〜（１１）のように求められる。また、式（１２）に示したように、ρiを微分したρx,i（＝∂ρ/∂ｘ｜i）はｂと等しくなる。
【数４】
【００２５】
このような条件の下で、図１０Ａおよび図１０Ｂに示したように、ルンゲクッタ法による時間更新（時間積分）を行う。ここでの目的は、時刻「ｎ」における密度（ρｎｉ）と各係数（ｋ1,i〜ｋ4,i）から時刻「ｎ＋１」における密度（ρn＋1ｉ）を求めること（式（１０１５）参照）と、時刻「ｎ」における密度の積分値（ρ-X,nｉ+1/2）と各係数（ｍ1,i+1/2〜ｍ4,i+1/2）から時刻「ｎ＋１」における密度（ρ-X,n+1ｉ+1/2）を求めること（式（１０１６）参照）である。式（１００１）〜式（１０１４）は、各係数（ｋ1,i〜ｋ4,i，ｍ1,i+1/2〜ｍ4,i+1/2）を求めるための式である。
【００２６】
まず、式（１００１）と式（１００２）によって、ｋ1,i（以下、「ｋ1」とも表記する。ｋ2,i〜ｋ4,iについても同様）とｍ1,i+1/2（以下、「ｍ1」とも表記する。ｍ2,i+1/2〜ｍ4,i+1/2についても同様）を求める。ここで、ｕi，ｕi+1，ρi，ρi+1の値は、時刻「ｎ」における値としてすでに求められているものを使用することができる。なお、移流項である「−ｕi∂ρ/∂ｘ｜i」における「∂ρ/∂ｘ｜i」は、前記したように風上補間によって求めることができる。
【００２７】
また、∂ｕ/∂ｘ｜iはｕの補間関数から求めることができる。ここでは、たとえば、ｕの補間関数として４次関数を使用し、中心補間によって∂ｕ/∂ｘ｜iを求めることができる。その場合、３つの格子点の座標を「i-1」，「i」および「i+1」とし、補間関数をｆ（Ｘ）、その積分値（前記式（８）と同様）をｆＸ（Ｘ）すると、その補間関数は、ｆ（i-1），ｆ（i），ｆ（i+1），ｆＸ（i-1/2）およびｆＸ（i+1/2）の拘束条件から決定することができる。∂2ρ/∂ｘ2｜iおよび非移流項中の∂ρ/∂ｘ｜iについても、同様にして中心補間から求めることができる。
【００２８】
Δｔは、時刻「ｎ」と時刻「ｎ＋１」の差の時間長（時間刻み幅）である。なお、この式（１００１）と式（１００２）によって演算が行われる段階では、すべての格子点に関して演算が行われ、その後、次の演算に進む。以降についても同様である。
【００２９】
次に、式（１００３）と式（１００４）によって、ρ*ｉ（ρn+1ｉの暫定値）とρ-X,*ｉ+1/2（ρ-X,n+1ｉ+1/2の暫定値）を求める。ここで、ρi，ρ-Xｉ+1/2は、時刻「ｎ」における値としてすでに求められているものを使用することができる。また、ｋ1とｍ1は、それぞれ式（１００１）と式（１００２）で求められた値である。
【００３０】
そして、式（１００５）と式（１００６）によって、ｋ2,iとｍ2,i+1/2を求める。ここで、ρ*ｉ（ρ*ｉ+1も同様）は、式（１００３）で求められた値である。また、移流項における∂ρ*/∂ｘ｜i（ρ*x,i）は、式（１２）からｂであることがわかり、前記した式（１０）の右辺に、式（１００３）と式（１００４）から求められたρ*ｉ，ρ*ｉ-1，ρ-X,*ｉ-1/2を代入することにより求められる。また、∂2ρ*/∂ｘ2｜iおよび非移流項中の∂ρ*/∂ｘ｜iについては、中心補間から求めることができる。
【００３１】
以下、同様にして、式（１００７）〜式（１０１４）を使うことにより、ｋ1,i〜ｋ4,i，ｍ1,i+1/2〜ｍ4,i+1/2の値をすべて求めることができる。そして、前記したように、式（１０１５）と式（１０１６）を使って、ρn＋1ｉとρ-X,n+1ｉ+1/2を求めることができる。なお、式（１００７）と式（１０１１）においてρ*ｉは新たに更新され、式（１００８）と式（１０１２）においてρ-X,*ｉ+1/2は新たに更新される。
【００３２】
次に、図３〜図６を参照しながら、本発明の第１実施形態の部分セミラグ保存形ＩＤＯ法（セミラグランジュ法を部分的に適用した保存形ＩＤＯ法）について説明する。
第１実施形態では、計算に、前記した比較例の場合と同じ式（１）〜式（１１）を使用するので、重複説明を適宜省略する。なお、第１実施形態は、式（２）における右辺の移流項である第１項の「−ｕi∂ρ/∂ｘ｜i」に関してセミラグランジュ移流を行う点で、比較例における計算と異なっている。
【００３３】
第１実施形態では、式（５）に対して、Ｘ＝ｘi−ｕiΔｔを代入することで、Δｔの分だけセミラグランジュ移流したρの値であるρSLi,Δtを算出することができる（式（１３）参照）。なお、式（１３）において、Ｆni（Ｘ），ａnおよびｂnは、時刻「ｎ」におけるそれぞれの値であることを意味している。また、式（１３）におけるΔｔをΔｔ／２と置き換えることで、Δｔ／２の分だけセミラグランジュ移流したρの値であるρSLi,Δt/2を算出することができる。
【数５】
【００３４】
ここで、図３は、第１実施形態の部分セミラグ保存形ＩＤＯ法における変数イメージを示した図である。図３に示したように、このρSLi,Δtの値を用いて、時刻「ｎ＋１」における密度（ρn＋1ｉ）を求めることができる（詳細は後記）。
【００３５】
次に、図４Ａ〜図４Ｃを参照しながら、第１実施形態における数値計算装置の演算（処理）について説明する（適宜他図参照）。図４Ａ〜図４Ｃは、第１実施形態における数値計算装置の演算を示したフローチャートである。数値計算装置１０の演算部４は、記憶部２とメモリ３を使用しながら各演算を行う。そして、この演算の目的は、比較例の場合と同様、時刻「ｎ」における密度（ρｉ）やその積分値（ρ-Xｉ+1/2（ρ-X,nｉ+1/2））などに基づいて、時刻「ｎ＋１」における密度（ρn+1ｉ）と積分値（ρ-X,n+1ｉ+1/2）を求めることである（ステップＳ４１０参照）。つまり、比較例の場合と同様、時刻「ｎ」におけるρi，ρ-Xｉ+1/2，ｕi，ｕi+1，ρi，ρi+1などの値はすでに求められているものとして、以下の演算を行う。また、∂ｕ/∂ｘ｜i，∂2ρ/∂ｘ2｜i，∂ρ/∂ｘ｜iなどは、前記した比較例の場合と同様、中心補間によって求められる。
【００３６】
演算部４は、式（２）を、右辺の移流項である第１項の「−ｕi∂ρ/∂ｘ｜i」を削除して（省いて）使用し（式（４０１１）参照）、また、式（４）はそのまま使用する（式（４０１２）参照）（ステップＳ４０１）。
次に、演算部４は、式（１３）を用いて、セミラグランジュ移流させた値であるρSLi,ΔtおよびρSLi,Δt/2を取得する（ステップＳ４０２）。
【００３７】
続いて、演算部４は、式（４０３１）と式（４０３２）によって、ｋ1,i（ｋ1）とｍ1,i+1/2（ｍ1）を求める。ここで、∂ｕ/∂ｘ｜iはｕの補間関数（説明を省略）から求めることができる。Δｔは、時刻「ｎ」と時刻「ｎ＋１」の差の時間長である。
【００３８】
次に、演算部４は、式（４０４１）と式（４０４２）によって、ρ*ｉ（ρn+1ｉの暫定値）とρ-X,*ｉ+1/2（ρ-X,n+1ｉ+1/2の暫定値）を取得する（ステップＳ４０４）。ここで、ｋ1とｍ1はステップＳ４０３で取得された値であり、ρSLi,Δt/2はステップＳ４０２で取得された値である。
【００３９】
続いて、演算部４は、式（４０５１）と式（４０５２）によって、ｋ2,iとｍ2,i+1/2を取得する（ステップＳ４０５）。ここで、ρ*ｉ（ρ*ｉ+1も同様）は、ステップＳ４０４で取得された値である。また、∂ｕ/∂ｘ｜i，∂2ρ*/∂ｘ2｜i，∂ρ*/∂ｘ｜iなどは、中心補間によって求められる。
【００４０】
以下、演算部４は、同様にして、ステップＳ４０６〜ステップＳ４０９の演算を行うことによって、ｋ1,i〜ｋ4,i，ｍ1,i+1/2〜ｍ4,i+1/2のすべての値を取得することができる。なお、ステップＳ４０６およびステップＳ４０８において、ρ*ｉおよびρ-X,*ｉ+1/2はそれぞれ新たに更新される。
そして、演算部４は、ステップＳ４０９までの演算で得た各値と式（４１０１）と式（４１０２）によって、ρn＋1ｉとρ-X,n+1ｉ+1/2を取得する（ステップＳ４１０）。
このようにして、第１実施形態の数値計算装置１０によれば、保存形ＩＤＯ法における偏微分方程式に対して部分的に（２つのうち片方の方程式に）セミラグランジュ法を適用することで、時間更新を行うことができる。
【００４１】
次に、図５および図６を参照しながら、第１実施形態の数値計算装置による数値計算方法の計算結果例について説明する。図５は、第１実施形態の数値計算装置による数値計算方法の計算結果例を示すグラフである。ここでは、実際に爆発実験を行い、大気圧に対する加圧（量）に関して、実測値（実験データ）、第１実施形態による計算結果、および、比較例による計算結果の三者を比較した。
【００４２】
図５に示すグラフにおいて、横軸は爆源からの距離（メートル）、縦軸は大気圧をゼロとした場合の加圧（パスカル）を表している。符号５１の破線は、実験で得られた各距離の地点におけるピーク加圧の実測値を示す。符号５１の破線を見れば、爆発による衝撃波が少しずつ弱まりながら遠方に伝搬していることがわかる。
【００４３】
符号５２の実線は、第１実施形態の部分セミラグ保存形ＩＤＯ法により計算したある瞬間の衝撃波プロファイルである。そのピーク値が符号５１の破線とほぼ一致しており、計算結果の精度が高いことがわかる。なお、この衝撃波プロファイルは、時間の経過とともに、その形をほぼ維持したまま右方に平行移動することになる。
一方、符号５３の二点鎖線は、比較例の非保存形ＩＤＯ法により計算したある瞬間の衝撃波プロファイルを示す。そのピーク値が符号５１の破線と大きく離れており、計算結果の精度が第１実施形態の部分セミラグ保存形ＩＤＯ法の場合よりも低いことがわかる。
【００４４】
同じ実験に関する別のグラフについて、図６を参照しながら説明する。図６は、第１実施形態の数値計算装置による数値計算方法の計算結果例を示す両対数グラフである。ここでも、図５の場合と同様に、加圧に関して、実測値（実験データ）、第１実施形態による計算結果、および、比較例による計算結果の三者を比較した。
【００４５】
図６に示す両対数グラフにおいて、横軸は爆源からの距離（メートル）、縦軸は大気圧をゼロとした場合のピーク加圧（パスカル）を表している。符号６１の破線は、実験で得られた各距離の地点におけるピーク加圧の実測値を示す。符号６１の破線を見れば、爆発による衝撃波が弱まりながら遠方に伝搬していることがわかる。
【００４６】
符号６２の実線は、第１実施形態の部分セミラグ保存形ＩＤＯ法により計算したピーク加圧に関する衝撃波プロファイルを示す。符号６２の実線は、符号６１の破線とほぼ一致しており、計算結果の精度が高いことがわかる。
一方、符号６３の二点鎖線は、比較例の非保存形ＩＤＯ法により計算したピーク加圧に関する衝撃波プロファイルである。符号６３の二点鎖線は、符号６１の破線と大きく離れており、計算結果の精度が第１実施形態の部分セミラグ保存形ＩＤＯ法の場合よりも低いことがわかる。
【００４７】
つまり、図５および図６からわかるように、第１実施形態の部分セミラグ保存形ＩＤＯ法によれば、比較例の非保存形ＩＤＯ法の場合に比べて、爆源から近距離（たとえば１０ｍ程度）に関しても高い精度での計算結果を得ることができ、また、爆源から遠距離（たとえば２００ｍ以上）に関してはさらに高い精度での計算結果を得ることができる。
【００４８】
このような第１実施形態の部分セミラグ保存形ＩＤＯ法を応用すれば、車体周りの空気や熱の流れ、風力発電の空力解析、気象計算、津波予測など、多くの流体計算を高精度で行うことができる。また、爆発現象や、水と空気が混在する流れのような、空間的な密度変化が１０００倍以上になる流体計算が可能となる。さらに、高速流体（高レイノルズ数流れ）で頻繁に発生する乱流に対し、定義された格子間隔よりも小さな渦による影響を渦粘性項無しでシミュレーションする、いわゆる陰的なラージ・エディ・シミュレーションになっていて、特別な乱流モデルを導入すること無しに複雑な流体計算を行うことができる。
【００４９】
さらに説明すると、一般に、多くの利点を有している保存形ＩＤＯ法の時間積分は、ＣＩＰ法のようなセミラグランジュ法ではなく、オイラー法で解くためにルンゲクッタ法などが用いられる。しかし、通常のルンゲクッタ法を用いてアンダーシュート（プロファイルにおける下方への飛び出し誤差）やオーバーシュート（プロファイルにおける上方への飛び出し誤差）を抑制することは非常に難しい。そこで、ルンゲクッタ法による時間積分の各段階において、保存量となる積分値については従来の保存形ＩＤＯ法による計算を行い、保存が保証されない格子点上の変数の時間積分に関する偏微分方程式を移流項と非移流項に分割し、セミラグランジュ法を適用し、さらに望ましくは後記するような単調性を有する補間関数を使用することで、アンダーシュートやオーバーシュートを抑制し、高い計算精度を実現することができる。
【００５０】
（第２実施形態）
次に、本発明の第２実施形態について説明する。第２実施形態は、非保存形ＩＤＯ法を用いた数値計算の例である。第２実施形態について説明する前に、その理解を助けるため、図１１、図１２Ａおよび図１２Ｂを参照しながら、非保存形ＩＤＯ法を用いた数値計算の比較例（従来技術）について説明する。図１１は、比較例の非保存形ＩＤＯ法における変数イメージを示した図である。図１２Ａおよび図１２Ｂは、比較例の非保存形ＩＤＯ法におけるルンゲクッタ法による時間更新（時間積分）を示した図である。なお、以下、第１実施形態と重複する説明については適宜省略する。また、第２実施形態において、移流項における微分値の算出には風上補間を用い、非移流項における２階以上の微分値の算出には中心補間を用いる。
【００５１】
比較例では、支配方程式として、拡散項を有する一次元質量保存の偏微分方程式である式（１４）（式（１）と同様）を使用する。つまり、流体に関して、ρは密度、ｕは速度、ｔは時間、ｘは座標を表す変数である。なお、速度ｕを求めるための式が別に必要であるが、その記載を省略し、ｕ，ｕx（ｕの一階微分）の値はその式などからすでに求められているものとする。
【数６】
【００５２】
そして、実際に数値計算を行うことができるように、式（１４）を変形して式（１５）とする。
【数７】
【００５３】
さらに、ρの微分値を式（１６）と定義し、その微分値に関する方程式を式（１７）と定義する。なお、∂ρx/∂ｘ｜iは、後記する補間関数（式（１８））による風上補間によって求められる値である。また、非移流項における∂2ｕ/∂ｘ2｜iおよび∂3ρ/∂ｘ3｜iは、中心補間によって求められる。
【数８】
【００５４】
また、ρに関する補間関数を式（１８）と定義する。式（１９）〜（２２）は、その拘束条件である。これらの式（１９）〜（２２）から、補間係数ａ〜ｄが式（２３）〜（２６）のように求められる。また、式（２７）に示すように、ρiを二階微分した∂ρx/∂ｘ｜iは２ｂと等しくなる。
【数９】
【００５５】
このような条件の下で、図１２Ａおよび図１２Ｂに示したように、ルンゲクッタ法による時間更新（時間積分）を行う。ここでの目的は、時刻「ｎ」における密度（ρｎｉ）と各係数（ｋ1,i〜ｋ4,i）から時刻「ｎ＋１」における密度（ρn＋1ｉ）を求めること（式（１２１５）参照）と、時刻「ｎ」における密度の微分値（ρnx,i）と各係数（ｍ1,i〜ｍ4,i）から時刻「ｎ＋１」における密度（ρn＋1x,i）を求めること（式（１２１６）参照）である。式（１２０１）〜式（１２１４）は、各係数（ｋ1,i〜ｋ4,i，ｍ1,i〜ｍ4,i）を求めるための式である。
【００５６】
まず、式（１２０１）と式（１２０２）によって、ｋ1,i（以下、「ｋ1」とも表記する。ｋ2,i〜ｋ4,iについても同様）とｍ1,i（以下、「ｍ1」とも表記する。ｍ2,i〜ｍ4,iについても同様）を求める。ここで、ｕi，ｕx,i，ρiおよびρx,iは、時刻「ｎ」における値としてすでに求められている値を使用することができる。また、前記したように、移流項における∂ρ/∂ｘ｜iと∂ρx/∂ｘ｜iは風上補間により求めることができる。さらに、非移流項における∂2ρ/∂ｘ2｜i，∂2ｕ/∂ｘ2｜iおよび∂3ρ/∂ｘ3｜iは、中心補間によって求められる。Δｔは、時刻「ｎ」と時刻「ｎ＋１」の差の時間長である。
【００５７】
次に、式（１２０３）と式（１２０４）によって、ρ*ｉ（ρn+1ｉの暫定値）とρ*x,i（ρn+1x,iの暫定値）を求める。ここで、ρiとρx,iは、時刻「ｎ」における値としてすでに求められている値を使用することができる。また、ｋ1とｍ1は、それぞれ式（１２０１）と式（１２０２）で求められた値である。
【００５８】
そして、式（１２０５）と式（１２０６）によって、ｋ2,iとｍ2,iを求める。ここで、ρ*ｉは、式（１２０３）で求められた値である。また、∂ρ*/∂ｘ｜iは、式（１２０４）で求められたρ*x,iである。
∂ρ*x/∂ｘ｜iは、式（２７）と式（２４）から求められる値であり（式（２７）の∂ρx/∂ｘ｜iに相当）、具体的には、式（２４）において式（１２０３）と式（１２０４）から求められたρ*ｉ，ρ*ｉ-1，ρ*x,iおよびρ*x,ｉ-1を代入し、それによって得られたｂに２を乗算することで求められる。また、非移流項における∂2ρ*/∂ｘ2｜i，∂2ｕ/∂ｘ2｜iおよび∂3ρ*/∂ｘ3｜iは、中心補間によって求められる。
【００５９】
以下、同様にして、式（１２０７）〜式（１２１４）を使うことにより、ｋ1,i〜ｋ4,i，ｍ1,i〜ｍ4,iの値をすべて求めることができる。そして、前記したように、式（１２１５）と式（１２１６）を使って、ρn＋1ｉとρn＋1x,ｉを求めることができる。なお、式（１２０７）と式（１２１１）においてρ*ｉは新たに更新され、式（１２０８）と式（１２１２）においてρ*x,ｉは新たに更新される。
【００６０】
次に、図７および図８Ａ〜図８Ｃを参照しながら、本発明の第２実施形態であるセミラグ非保存形ＩＤＯ法（セミラグランジュ法を適用した非保存形ＩＤＯ法）について説明する。
第２実施形態では、計算に、前記した比較例の場合と同じ式（１４）〜式（２６）を使用するので、重複説明を適宜省略する。なお、第２実施形態は、式（１５）における右辺の移流項である第１項の「−ｕiρx,i」と、式（１７）における右辺の移流項である第１項の「−ｕi∂ρx/∂ｘ｜i」とに関してセミラグランジュ移流を行う点で、比較例における計算と異なっている。
【００６１】
第２実施形態では、式（１８）に対して、Ｘ＝ｘi−ｕiΔｔを代入することで、Δｔの分だけセミラグランジュ移流したρの値であるρSLi,Δtを算出することができる（式（２８）参照）。なお、Ｆni（Ｘ），ａn，ｂn，ｃnおよびｄnは、時刻「ｎ」におけるそれぞれの値であることを意味している。また、式（２８）におけるΔｔをΔｔ／２と置き換えることで、Δｔ／２の分だけセミラグランジュ移流したρの値であるρSLi,Δt/2を算出することができる。
【数１０】
【００６２】
また、式（１８）を微分した式（説明を省略）に対して、Ｘ＝ｘi−ｕiΔｔを代入することで、Δｔの分だけセミラグランジュ移流したρx,iの値であるρSLx,i,Δtを算出することができる（式（２９）参照）。また、式（２９）におけるΔｔをΔｔ／２と置き換えることで、Δｔ／２の分だけセミラグランジュ移流したρx,iの値であるρSLx,i,Δt/2を算出することができる。
【数１１】
【００６３】
ここで、図７は、第２実施形態のセミラグ非保存形ＩＤＯ法における変数イメージを示した図である。図７に示したように、このρSLi,Δt（ρSLi,Δt/2）およびρSLx,i,Δt（ρSLx,i,Δt/2）の値を用いて、時刻「ｎ＋１」における密度（ρn＋1ｉ）を求めることができる（詳細は後記）。
【００６４】
次に、図８Ａ〜図８Ｃを参照しながら、第２実施形態における数値計算装置の演算（処理）について説明する（適宜他図参照）。図８Ａ〜図８Ｃは、第２実施形態における数値計算装置の演算を示したフローチャートである。数値計算装置１０の演算部４は、記憶部２とメモリ３を使用しながら各演算を行う。そして、この演算の目的は、比較例の場合と同様、時刻「ｎ」における密度（ρｉ）やその微分値（ρx,ｉ）などに基づいて、時刻「ｎ＋１」における密度（ρn+1ｉ）と微分値（ρn+1x,ｉ）を求めることである（ステップＳ４３０参照）。つまり、比較例の場合と同様、時刻「ｎ」におけるρi，ρx,i，ｕx,iなどの値はすでに求められているものとして以下の演算を行う。なお、非移流項における∂2ρ/∂ｘ2｜i，∂2ｕ/∂ｘ2｜iおよび∂3ρ/∂ｘ3｜iは、中心補間によって求められる。
【００６５】
演算部４は、式（１５）について右辺の移流項である第１項の「−ｕiρx,i」を削除して（省いて）使用し（式（４２１１）参照）、また、式（１７）について右辺の移流項である第１項の「−ｕi∂ρx/∂ｘ｜i」を削除して（省いて）使用する（式（４２１２）参照）（ステップＳ４２１）。
次に、演算部４は、式（２８）および式（２９）を用いて、セミラグランジュ移流させた値であるρSLi,Δt，ρSLi,Δt/2，ρSLx,i,ΔtおよびρSLx,i,Δt/2を取得する（ステップＳ４２２）。
【００６６】
続いて、演算部４は、式（４２３１）と式（４２３２）によって、ｋ1,i（ｋ1）とｍ1,i（ｍ1）を求める。ここで、∂2ρ/∂ｘ2｜i，∂2ｕ/∂ｘ2｜iおよび∂3ρ/∂ｘ3｜iは、中心補間によって求められる。Δｔは、時刻「ｎ」と時刻「ｎ＋１」の差の時間長である。
【００６７】
次に、演算部４は、式（４２４１）と式（４２４２）によって、ρ*ｉ（ρn+1ｉの暫定値）とρ*x,i（ρn+1x,iの暫定値）を取得する（ステップＳ４２４）。ここで、ｋ1とｍ1はステップＳ４２３で取得された値であり、ρSLi,Δt/2およびρSLx,i,Δt/2はステップＳ４２２で取得された値である。
【００６８】
続いて、演算部４は、式（４２５１）と式（４２５２）によって、ｋ2,iとｍ2,iを取得する（ステップＳ４２５）。ここで、ρ*ｉおよびρ*x,ｉは、ステップＳ４２４で取得された値である。また、∂2ρ*/∂ｘ2｜i，∂2ｕ/∂ｘ2｜iおよび∂3ρ*/∂ｘ3｜iは、中心補間によって求められる。
【００６９】
以下、演算部４は、同様にして、ステップＳ４２６〜ステップＳ４２９の演算を行うことによって、ｋ1,i〜ｋ4,i，ｍ1,i〜ｍ4,iのすべての値を取得することができる。なお、ステップＳ４２６およびステップＳ４２８において、ρ*ｉおよびρ*x,ｉはそれぞれ新たに更新される。
そして、演算部４は、ステップＳ４２９までの演算で得た各値と式（４３０１）と式（４３０２）によって、ρn＋1ｉとρn+1x,iを取得する（ステップＳ４３０）。
【００７０】
このようにして、第２実施形態の数値計算装置１０によれば、非保存形ＩＤＯ法における偏微分方程式に対してセミラグランジュ法を適用することで、時間更新を行うことができる。また、その詳細な効果についての説明を省略するが、図５および図６を参照して説明した第１実施形態の場合と同様の効果を得ることができる。
【００７１】
（変形例）
第１実施形態では、補間関数として、式（５）に示した式を使用した場合について説明したが、それ以外に、式（３０）に示す分母に関数を有する有理関数を使用してもよい。β，ａおよびｃは、拘束条件により決定される係数である。
【数１２】
【００７２】
また、第２実施形態では、補間関数として、式（１８）に示した式を使用した場合について説明したが、それ以外に、式（３１）に示す分母に関数を有する有理関数を使用してもよい。ａ，ｂ，ｃおよびｄは、拘束条件により決定される係数である。
【数１３】
そして、αとＢは、数値計算の諸条件により、式（３１）が単調性を有するように決定すればよい。たとえば、α＝１とし、Ｂを式（３２）および式（３３）により決定すればよい。
【数１４】
【００７３】
このような式（３０）や式（３１）に示すように、単調性を有する（単調増加または単調減少する）補間関数を使用すれば、密度比や圧力比などが極端に大きな（たとえば１０００倍以上）現象について数値計算を行う場合でも、その計算精度を高く維持することができる。つまり、密度比や圧力比などが極端に大きな現象に関する数値計算を行う場合、変数の数値の大きな場面での誤差の影響が大きくなる。そのため、単調性を有する補間関数を使用することで、急激な空間勾配のプロファイルに対してアンダーシュートやオーバーシュートをなくし、計算を安定させ、高い精度を維持することができるようになる。
【００７４】
図１３の（ａ）は、前記した２つの比較例においてκ（拡散係数）＝０の場合の計算結果のプロファイルの例を示した図である。図１３の（ａ）に示すように、比較例による計算結果（ドット）は、解析解（実線）に対してアンダーシュートやオーバーシュートを起こすことがあった。つまり、比較例では、密度比や圧力比などが極端に大きな現象についての数値計算が破綻することがあった。
【００７５】
図１３の（ｂ）は、本実施形態（第１実施形態および第２実施形態に前記変形例の補間関数を適用したもの）においてκ（拡散係数）＝０の場合の計算結果のプロファイルの例を示した図である。図１３の（ｂ）に示すように、本実施形態において上記のような有理関数を用いた場合の計算結果（ドット）は、解析解（実線）に対してアンダーシュートやオーバーシュートを起こすことがなくなる。つまり、本実施形態では、密度比や圧力比などが極端に大きな現象についての数値計算が安定するという効果を奏する。
なお、単調性を有する補間関数として、式（３０）および式（３１）を示したが、これらに限定されるものではない。
【００７６】
また、当業者であれば、前記した各実施形態とその変形例にかかる数値計算方法をコンピュータに実行させるためのプログラムを作成し、さらに、そのプログラムをＣＤ（Compact Disc）やＤＶＤ（Digital Versatile Disk）などの記録媒体に記録することができる。
【００７７】
以上で本実施形態の説明を終えるが、本発明の態様はこれらに限定されるものではない。
たとえば、本実施形態では、説明を簡潔にするため、一次元空間の偏微分方程式について説明したが、２次元以上の空間の偏微分方程式を使用してもかまわない。
また、流体の変数として、密度を中心に説明したが、速度、圧力、温度などの別の変数に関しても同様に、本発明のセミラグランジュ移流による数値計算を適用することができる。
【００７８】
さらに、時間積分のルンゲクッタ法については、４段で行う手法について説明したが、２段や３段など別の手法で行うようにしてもよい。
その他、数値計算装置や各フローチャートなどの具体的な構成について、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。
【図面の簡単な説明】
【００７９】
【図１】数値計算装置の構成図である。
【図２】流体の数値計算の演算の大まかな流れを表わしたフローチャートである。
【図３】第１実施形態の部分セミラグ保存形ＩＤＯ法における変数イメージを示した図である。
【図４Ａ】第１実施形態における数値計算装置の演算を示したフローチャートである。
【図４Ｂ】第１実施形態における数値計算装置の演算を示したフローチャートである。
【図４Ｃ】第１実施形態における数値計算装置の演算を示したフローチャートである。
【図５】第１実施形態の数値計算装置による数値計算方法の計算結果例を示すグラフである。
【図６】第１実施形態の数値計算装置による数値計算方法の計算結果例を示す両対数グラフである。
【図７】第２実施形態のセミラグ非保存形ＩＤＯ法における変数イメージを示した図である。
【図８Ａ】第２実施形態における数値計算装置の演算を示したフローチャートである。
【図８Ｂ】第２実施形態における数値計算装置の演算を示したフローチャートである。
【図８Ｃ】第２実施形態における数値計算装置の演算を示したフローチャートである。
【図９】比較例の保存形ＩＤＯ法における変数イメージを示した図である。
【図１０Ａ】比較例の保存形ＩＤＯ法におけるルンゲクッタ法による時間更新（時間積分）を示した図である。
【図１０Ｂ】比較例の保存形ＩＤＯ法におけるルンゲクッタ法による時間更新（時間積分）を示した図である。
【図１１】比較例の非保存形ＩＤＯ法における変数イメージを示した図である。
【図１２Ａ】比較例の非保存形ＩＤＯ法におけるルンゲクッタ法による時間更新（時間積分）を示した図である。
【図１２Ｂ】比較例の非保存形ＩＤＯ法におけるルンゲクッタ法による時間更新（時間積分）を示した図である。
【図１３】（ａ）は、比較例による計算結果のプロファイルの例を示した図であり、（ｂ）は、本実施形態による計算結果のプロファイルの例を示した図である。
【符号の説明】
【００８０】
１ 入力部
２ 記憶部
３ メモリ
４ 演算部
５ 出力部
１０ 数値計算装置
【特許請求の範囲】
【請求項１】
複数の定義点における１つ以上の変数とその微分係数を含んだ複数の値を用いて偏微分方程式に関する演算を行い、前記変数のうち少なくともいずれか１つの変数を時間積分することにより、流体に関する数値計算を行う局所補間微分オペレータ法を実行する数値計算装置による数値計算方法であって、
前記数値計算装置は、記憶部と演算部を有し、
前記記憶部は、前記時間積分を行う変数の値と、その変数に関する補間関数と、を記憶し、
前記演算部は、前記偏微分方程式に関して、時間軸上の固定点以外の流動点による影響も考慮した方程式解法であるラグランジュ法を移流項に適用するセミラグランジュ移流を、前記補間関数を使って前記偏微分方程式の移流項に対して行うことで前記時間積分を行う
ことを特徴とする数値計算方法。
【請求項２】
前記局所補間微分オペレータ法は、前記時間積分を行う変数の区間積分値を保存する保存形局所補間微分オペレータ法であり、
前記演算部は、前記区間積分値を演算する場合には前記セミラグランジュ移流を行わず、前記時間積分を行う変数の値を演算する場合には前記セミラグランジュ移流を行うことで、前記時間積分を行う
ことを特徴とする請求項１に記載の数値計算方法。
【請求項３】
前記演算部は、ルンゲクッタ法によって前記時間積分を行う
ことを特徴とする請求項１または請求項２に記載の数値計算方法。
【請求項４】
前記補間関数は、単調増加または単調減少する有理関数である
ことを特徴とする請求項１から請求項３のいずれか一項に記載の数値計算方法。
【請求項５】
請求項１から請求項４のいずれか一項に記載の数値計算方法をコンピュータに実行させるためのプログラム。
【請求項６】
請求項５に記載のプログラムを記録した記録媒体。
【請求項１】
複数の定義点における１つ以上の変数とその微分係数を含んだ複数の値を用いて偏微分方程式に関する演算を行い、前記変数のうち少なくともいずれか１つの変数を時間積分することにより、流体に関する数値計算を行う局所補間微分オペレータ法を実行する数値計算装置による数値計算方法であって、
前記数値計算装置は、記憶部と演算部を有し、
前記記憶部は、前記時間積分を行う変数の値と、その変数に関する補間関数と、を記憶し、
前記演算部は、前記偏微分方程式に関して、時間軸上の固定点以外の流動点による影響も考慮した方程式解法であるラグランジュ法を移流項に適用するセミラグランジュ移流を、前記補間関数を使って前記偏微分方程式の移流項に対して行うことで前記時間積分を行う
ことを特徴とする数値計算方法。
【請求項２】
前記局所補間微分オペレータ法は、前記時間積分を行う変数の区間積分値を保存する保存形局所補間微分オペレータ法であり、
前記演算部は、前記区間積分値を演算する場合には前記セミラグランジュ移流を行わず、前記時間積分を行う変数の値を演算する場合には前記セミラグランジュ移流を行うことで、前記時間積分を行う
ことを特徴とする請求項１に記載の数値計算方法。
【請求項３】
前記演算部は、ルンゲクッタ法によって前記時間積分を行う
ことを特徴とする請求項１または請求項２に記載の数値計算方法。
【請求項４】
前記補間関数は、単調増加または単調減少する有理関数である
ことを特徴とする請求項１から請求項３のいずれか一項に記載の数値計算方法。
【請求項５】
請求項１から請求項４のいずれか一項に記載の数値計算方法をコンピュータに実行させるためのプログラム。
【請求項６】
請求項５に記載のプログラムを記録した記録媒体。
【図１】
【図２】
【図３】
【図４Ａ】
【図４Ｂ】
【図４Ｃ】
【図５】
【図６】
【図７】
【図８Ａ】
【図８Ｂ】
【図８Ｃ】
【図９】
【図１０Ａ】
【図１０Ｂ】
【図１１】
【図１２Ａ】
【図１２Ｂ】
【図１３】
【図２】
【図３】
【図４Ａ】
【図４Ｂ】
【図４Ｃ】
【図５】
【図６】
【図７】
【図８Ａ】
【図８Ｂ】
【図８Ｃ】
【図９】
【図１０Ａ】
【図１０Ｂ】
【図１１】
【図１２Ａ】
【図１２Ｂ】
【図１３】
【公開番号】特開２００８−１９７９４９（Ｐ２００８−１９７９４９Ａ）
【公開日】平成２０年８月２８日（２００８．８．２８）
【国際特許分類】
【出願番号】特願２００７−３２９６７（Ｐ２００７−３２９６７）
【出願日】平成１９年２月１４日（２００７．２．１４）
【出願人】（３０４０２１４１７）国立大学法人東京工業大学 (1,821)
【Ｆターム（参考）】
【公開日】平成２０年８月２８日（２００８．８．２８）
【国際特許分類】
【出願日】平成１９年２月１４日（２００７．２．１４）
【出願人】（３０４０２１４１７）国立大学法人東京工業大学 (1,821)
【Ｆターム（参考）】
[ Back to top ]