移流方程式を数値的に解く方法、並びに、ブラソフ方程式及びマクスウェル方程式を数値的に解く方法、並びに、プログラム。
【課題】数値拡散が少ない高精度の移流方程式の数値解法を提供すること。
【解決手段】離散化された空間上で定義された物理量と速度場とから、この物理量の移流を求める移流方程式を数値計算によって解く。記憶手段には、物理量ｆ（ｘ、ｔ0）及びその０からｍｍａｘ次までのモーメント積分量Ｍm（ｘ、ｔ0）を記憶する。空間プロファイル算出ステップでは、ｆ（ｘ、ｔn）及びＭm（ｘ、ｔn）の束縛条件を課して、ｆ（ｘ、ｔn）及びＭm（ｘ、ｔn）の空間プロファイルＦ（ｘ、ｔn）及びＧm（ｘ、ｔn）を決定する。移流項計算ステップでは、ｆ（ｘ、ｔ*）＝Ｆ（ｘ−ｕ×ｄｔ、ｔn）及びＭm（ｘ、ｔ*）＝Ｇm（ｘ−ｕ×ｄｔ、ｔn）により中間時間ステップでの物理量ｆ及びＭmを求める（ここでｕは移流速度ベクトル）。非移流項加算ステップでは、非移流項をｆ（ｘ、ｔ*）及びＭm（ｘ、ｔ*）から算出し、これを用いて次の時間ステップでの物理量ｆ（ｘ、ｔn+1）及びモーメント積分量Ｍm（ｘ、ｔn+1）を求める。
【解決手段】離散化された空間上で定義された物理量と速度場とから、この物理量の移流を求める移流方程式を数値計算によって解く。記憶手段には、物理量ｆ（ｘ、ｔ0）及びその０からｍｍａｘ次までのモーメント積分量Ｍm（ｘ、ｔ0）を記憶する。空間プロファイル算出ステップでは、ｆ（ｘ、ｔn）及びＭm（ｘ、ｔn）の束縛条件を課して、ｆ（ｘ、ｔn）及びＭm（ｘ、ｔn）の空間プロファイルＦ（ｘ、ｔn）及びＧm（ｘ、ｔn）を決定する。移流項計算ステップでは、ｆ（ｘ、ｔ*）＝Ｆ（ｘ−ｕ×ｄｔ、ｔn）及びＭm（ｘ、ｔ*）＝Ｇm（ｘ−ｕ×ｄｔ、ｔn）により中間時間ステップでの物理量ｆ及びＭmを求める（ここでｕは移流速度ベクトル）。非移流項加算ステップでは、非移流項をｆ（ｘ、ｔ*）及びＭm（ｘ、ｔ*）から算出し、これを用いて次の時間ステップでの物理量ｆ（ｘ、ｔn+1）及びモーメント積分量Ｍm（ｘ、ｔn+1）を求める。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、物理現象を表現する基本的な方程式である移流方程式を数値的に解く方法に関する。更には、移流方程式を数値的に解くためのプログラムに関する。
【背景技術】
【０００２】
移流方程式とは、ある物理量がある速度場によって移動（伝播）する現象を記述した方程式である。この移流現象は流体運動を初め、多くの物理現象に関わるため、現象を記述する基本的かつ重要な方程式であるといえる。移流方程式には解析解が存在する一方、空間上に離散化されたグリッド上で移流方程式を数値的に解く場合は、一転して多くの困難を伴うことが古くから知られている。その重要性かつ困難性から、移流方程式を数値的に解く手法を開発することは、学術面のみならず、産業上においても重要な課題として残されている。
【０００３】
移流方程式を数値的に解く場合、移流させる物理量をコンピュータのメモリー空間上に離散化されたグリッド上に記憶するため、グリッド間の情報は失われる。そこで、記憶してある量を用いてグリッド間のプロファイルを再構築する必要がある。古典的な１次風上差分法は、注目するグリッド上での移流速度から、風上方向にあるグリッド間の情報を線形補間することによって、物理量のプロファイルを再構築する手法である。この手法は数値的に安定であり、古典的手法の代表的なものである。図１０に１次風上差分法の計算例を示す。横軸が実空間、縦軸が物理量を表し、「＋」が計算結果で、破線が解析解である。図でわかるように、本計算手法は数値粘性が大きいため、時間と共にプロファイルが減衰してしまい、長時間計算には適さないことが知られている。一方、近年開発された手法として、ＣＩＰ（Ｃｏｎｓｔｒａｉｎｅｄ Ｉｎｔｅｒｐｏｌａｔｉｏｎ Ｐｒｏｆｉｌｅ）法が挙げられる（非特許文献１）。図１１に１次風上差分法及びＣＩＰ法の概念を示す。横軸が実空間、黒丸がグリッドを表す。グリッド上に離散化された物理量を数値的に移流させる際（図１１ａ）、１次風上差分法では物理量のみを移流させようとするため（図１１ｂ）、プロファイルが減衰してしまう（図１１ｃ）。ＣＩＰ法は物理量のみならず、物理量の空間プロファイルの勾配（微分値）の時間発展をも同時に解き進める手法である。図１１ｄで示されているように、風上方向にあるグリッド間のプロファイルを再構築する際に、物理量と勾配の情報も使うため、高精度かつ安定に解く手法として近年注目を浴びてきている。ＣＩＰ法の派生である保存保証型ＣＩＰ法（ＣＩＰ−ＣＳＬ２法、非特許文献２）は、物理量とそのグリッド間の積分量を記憶する。グリッド上の物理量とグリッド間の積分量の情報を使い、グリッド間のプロファイルを高精度で再構築する移流方程式解法である。この解法は、ＣＩＰ法と同程度の精度を達成し、またＣＩＰ法に対して質量保存が保証されるメリットがある。しかし、不連続面において若干の数値振動が生じる点は、極めて高精度な移流方程式の解を要求される分野、例えば、プラズマ物理学の分野では依然大きな問題点として残されている。
【０００４】
前記プラズマ物理学の分野では、上に挙げた保存性、無振動性を特徴とする手法が近年提案されている。非特許文献３では、プラズマ物理現象を記述するブラソフ方程式を解くことを目的として、ＰＦＣ（Ｐｏｓｉｔｉｖｅ ａｎｄ Ｆｌｕｘ Ｃｏｎｓｅｒｖａｔｉｖｅ）法が開発されている。本手法は移流方程式を保存形で解くことにより、保存則を満たし、不連続面近傍では流束に制限をかけることにより、数値振動を少なくする手法である。さらに近年では、流束制限手法を工夫することにより、無振動性を保証する手法も開発された（非特許文献４）。以上の手法は、極めて高い精度が要求されるプラズマ物理分野における移流方程式である、ブラソフ方程式の解法として知られている。しかしながら、数値拡散の抑制については未だ不十分であり、数値拡散のない（極めて少ない）移流方程式解法が求められている。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】タカシ・ヤベ（Takashi Yabe）、フェン・ショウ（Feng Xiao）、及びタカユキ・ウツミ（Takayuki Utsumi）、"ザ コンストレインド インターポレーション プロファイル メソッド フォー マルチフェーズ アナリシス（The Constrained Interpolation Profile Method for Multiphase Analysis）"、ジャーナル オブ コンピュテーショナル フィジクス（Journal of Computational Physics）、オランダ王国、アムステルダム、２００１年５月、第１６０巻、ｐ．５５６−５９３
【非特許文献２】タカシ・ヤベ（Takashi Yabe）、リョウタロウ・タナカ（Ryotaro Tanaka）、タカシ・ナカムラ（Takashi Nakamura）、及びフェン・ショウ（Feng Xiao）、マンスリー ウェザー レビュー（Monthly Weather Review）、アメリカ合衆国、ボストン、２００１年２月、第１２９巻、ｐ．３３２−３４４
【非特許文献３】フランシス・フィルべ（Francis Filbet）、エリック・ソネンドラッカー（Eric Sonnendrucker）、及びピエール・ベルランド（Pierre Bertrand）、"コンザバティブ ニューメリカル スキームズ フォー ザ ブラソフ エクエーション（Conservative Numerical Schemes for the Vlasov Equation）"、ジャーナル オブ コンピュテーショナル フィジクス（Journal of Computational Physics）、オランダ王国、アムステルダム、２００１年９月、第１７２巻、ｐ．１６６−１８７
【非特許文献４】タカユキ・ウメダ（Takayuki Umeda）、ケンタロウ・トガノ（Kentaro Togano）、及びタツキ・オギノ（Tatsuki Ogino）、"トゥーディメンショナル フルエレクトロマグネティック ブラソフ コード ウィズ コンサバティブ スキーム アンド イッツ アプリケーション トゥ マグネティック リコネクション（Two-dimensional full-electromagnetic Vlasov code with conservative scheme and its application to magnetic reconnection）", コンピューター フィジクスコミュニケーションズ（Computer Physics Communications）、オランダ王国、アムステルダム、２００９年３月、第１８０巻、ｐ．３６５−３７４
【発明の概要】
【発明が解決しようとする課題】
【０００６】
上述した従来の解法では、そのいずれも、高い精度を要求される分野、例えば、プラズマ物理学への適用という点では、プロファイル再構築精度は不十分である。例えば図６左列及び図７左列では、ＣＩＰ−ＣＳＬ２法を用いた１次元移流方程式及び２次元移流方程式の数値計算結果を示している。プロファイルがなまってしまい、解析解から大きく外れてしまっている。このように、数値振動や数値拡散が無視できなくなるといった欠点があった。これは、多次元問題（図７左列）で特に顕著である。このように、既存の解法を遥かに上回る高精度の移流方程式解法が必要とされている。このような数値解法及び当該解法を実現するためのプログラムを提供することが本発明の目的である。
【課題を解決するための手段】
【０００７】
数値精度の低さが原因で発生する数値振動や数値拡散は、言いかえれば物理量の情報エントロピーが増大してしまっていることに相当する。独立変数ｘに対して物理量ｆ（ｘ）を定義すると、一般的に情報エントロピーは∫ｆ（ｘ）×ｂ（ｆ（ｘ））ｄｘと、情報量の関数ｂを用いて表現される。そこでｂをｘについてテイラー展開すると、情報エントロピーはｆ（ｘ）のｍ次モーメント積分量Ｍm（ｘ）＝∫ｘmｆ（ｘ）ｄｘ／ｍ！（ここで、ｍは０からｍｍａｘまでの整数を表す。）の線形結合ΣＭm（ｘ）（ｍ＝０、・・・、ｍｍａｘ）で表現されることがわかる。よって、数値分散や数値拡散を抑制するためには、情報エントロピーの保存、つまり、各モーメント積分量が保存するように移流方程式を解くことが本質的である。例えば、熱平衡状態を表現するガウス（正規）分布の情報エントロピーは２次のモーメント積分量Ｍ2（ｘ）＝∫ｘ2ｆ（ｘ）ｄｘ／２で表現される。前記ｍｍａｘの選択は任意であるが、ガウス分布は最も基本的な物理量のプロファイルであることから、少なくとも２次までのモーメント積分量を保存させることが本質である。１次元移流方程式解法においては、ｍｍａｘ＝０は前記ＣＩＰ−ＣＳＬ２法に一致する。本発明に対応する実施例では、前記情報エントロピーの考察に基づき、ｍｍａｘ＝２としている。
【０００８】
以上の考察に基づき、本発明は、物理量とそのｍｍａｘ次までのモーメント積分量を独立に記憶し、更新することを特徴とする移流方程式解法であり、以下に説明する解法をマルチモーメント移流法（Ｍｕｌｔｉ−Ｍｏｍｅｎｔ Ａｄｖｅｃｔｉｏｎ法、以下ＭＭＡ法）と呼ぶ。
【０００９】
本解法は、初期時刻ｔ0での空間上の各グリッドにおける物理量ｆ（ｘ、ｔ0）（ここで、ベクトルｘは空間上での位置ベクトル）と、その物理量ｆのｍｍａｘ次までのモーメント積分量Ｍm（ｘ、ｔ0）と、最大時間ステップ数ｎｍａｘと、を記憶手段に記憶する初期値記憶ステップ（図１に示すＳ１０１）と、空間のグリッド上で定義された移流速度ベクトルｕから、上流点及び上流方向のグリッドを決める上流点算出ステップ（Ｓ１０２）と、前記物理量ｆ（ｘ、ｔn）及び前記各モーメント積分量Ｍm（ｘ、ｔn）の束縛条件を課して、物理量ｆ（ｘ、ｔn）の空間プロファイルＦ（ｘ、ｔn）（ここで、ｔnは第ｎ時間ステップでの時刻を表し、ｎは０以上の整数である。）及び各モーメント積分量Ｍm（ｘ、ｔn）の空間プロファイルＧm（ｘ、ｔn）と、を求める空間プロファイル算出ステップ（Ｓ１０３）と、前記空間プロファイルに基づく次式ｆ（ｘ、ｔ*）＝Ｆ（ｘ−ｕ×ｄｔ、ｔn）及びＭm（ｘ、ｔ*）＝Ｇm（ｘ−ｕ×ｄｔ、ｔn）に従って（ここで、ｄｔ＝ｔn+1−ｔnは時間刻み幅、t*は移流後の中間時間ステップである。）、中間時間ステップでの物理量ｆ（ｘ、ｔ*）及び各モーメント積分量Ｍm（ｘ、ｔ*）を算出する移流項計算ステップ（Ｓ１０４）と、前記物理量ｆ（ｘ、ｔ*）及び各モーメント積分量Ｍm（ｘ、ｔ*）に基づき、物理量ｆの非移流項Ｓ（ｘ、ｔ*）及び各モーメント積分量Ｍmの非移流項Ｓm（ｘ、ｔ*）を算出し、これら物理量ｆの非移流項Ｓ（ｘ、ｔ*）及び各モーメント積分量Ｍmの非移流項Ｓm（ｘ、ｔ*）を前記物理量ｆ（ｘ、ｔ*）及び各モーメント積分量Ｍm（ｘ、ｔ*）に加算して、次の時間ステップに対応する物理量ｆ（ｘ、ｔn+1）及び各モーメント積分量Ｍm（ｘ、ｔn+1）を算出する非移流項加算ステップ（Ｓ１０５）と、計算結果を出力する出力ステップ（Ｓ１０６）と、から構成されている。本発明により、数値拡散と数値振動が極めて少ない移流方程式の解を得ることが出来る。
【００１０】
プラズマ物理分野での移流方程式であるブラソフ方程式及びマクスウェル方程式の解法は、初期時刻ｔ0での位相空間上におけるプラズマ分布関数ｆ（ｘ、ｖ、ｔ0）（ここでベクトルｘは実空間、ベクトルｖは速度空間での位置ベクトルを表す。）と、その各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ0）と、実空間上で定義された電場ベクトルＥ（ｘ、ｔ0）と、磁場ベクトルＢ（ｘ、ｔ0）と、最大時間ステップ数ｎｍａｘと、を記憶手段に記憶する初期値記憶ステップ（図２に示すＳ２０１）と、プラズマ分布関数ｆ（ｘ、ｖ、ｔn）と、速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔn）と、から実空間上における移流更新を行い、第一中間時間ステップでのプラズマ分布関数ｆ（ｘ、ｖ、ｔ*）と、速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ*）と、を算出する実空間時間更新ステップ（Ｓ２０２）と、０次の速度空間モーメント積分量Ｍ0（ｘ、ｖ、ｔ*）と、１次の速度空間モーメント積分量Ｍ1（ｘ、ｖ、ｔ*）と、から、電荷密度と、電流密度と、を求める電荷・電流密度算出ステップ（Ｓ２０３）と、マクスウェル方程式を解くことにより、前記電荷・電流密度と、電場ベクトルＥ（ｘ、ｔn）と、磁場ベクトルＢ（ｘ、ｔn）と、から次の時間ステップに対応する電場ベクトルＥ（ｘ、ｔn+1）と、磁場ベクトルＢ（ｘ、ｔn+1）と、を算出する電磁場更新ステップ（Ｓ２０４）と、前記時間更新後の電場ベクトルＥ（ｘ、ｔn+1）と、磁場ベクトルＢ（ｘ、ｔn+1）と、から速度空間上における移流速度ベクトルである加速度場ベクトルＡ（ｘ、ｖ、ｔn+1）を算出する加速度場算出ステップ（Ｓ２０５）と、速度空間のグリッド上で定義された加速度ベクトルＡ（ｘ、ｖ、ｔn+1）から、速度空間上の上流点及び上流方向のグリッドを決める速度空間上流点算出ステップ（Ｓ２０６）と、速度空間上でプラズマ分布関数ｆ（ｘ、ｖ、ｔ*）及び前記各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ*）の束縛条件を課して、プラズマ分布関数ｆ（ｘ、ｖ、ｔ*）の速度空間プロファイルＦ（ｘ、ｖ、ｔ*）及び各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ*）の速度空間プロファイルＧm（ｘ、ｖ、ｔ*）を求める速度空間プロファイル算出ステップ（Ｓ２０７）と、前記空間プロファイルに基づく次式ｆ（ｘ、ｖ、ｔ**）＝Ｆ（ｘ、ｖ−Ａ×ｄｔ、ｔ*）及びＭm（ｘ、ｖ、ｔ**）＝Ｇm（ｘ、ｖ−Ａ×ｄｔ、ｔ*）に従って、第二中間時間ステップでのプラズマ分布関数ｆ（ｘ、ｖ、ｔ**）及び各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ**）を算出する速度空間移流項計算ステップ（Ｓ２０８）と、前記プラズマ分布関数ｆ（ｘ、ｖ、ｔ**）及び各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ**）に基づき、プラズマ分布関数ｆの非移流項Ｓ（ｘ、ｖ、ｔ**）及び各速度空間モーメント積分量Ｍmの非移流項Ｓm（ｘ、ｖ、ｔ**）を算出し、これらプラズマ分布関数ｆの非移流項Ｓ（ｘ、ｖ、ｔ**）及び各速度空間モーメント積分量Ｍmの非移流項Ｓm（ｘ、ｖ、ｔ**）を前記プラズマ分布関数ｆ（ｘ、ｖ、ｔ**）及び各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ**）に加算して、次の時間ステップに対応するプラズマ分布関数ｆ（ｘ、ｖ、ｔn+1）及び各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔn+1）を算出する速度空間非移流項加算ステップ（Ｓ２０９）と、計算結果を出力する出力ステップ（Ｓ２１０）と、から構成されている。本発明により、従来の手法に比べて極めて高い信頼性を持つブラソフ方程式及びマクスウェル方程式の解を得ることができる。
【図面の簡単な説明】
【００１１】
【図１】本実施の形態における移流方程式解法装置における処理手順を示すフローチャートである。
【図２】本実施の形態におけるブラソフ方程式及びマクスウェル方程式解法装置における処理手順を示すフローチャートである。
【図３】本実施の形態における移流方程式及びブラソフ・マクスウェル方程式解法装置のハードウェア構成を示すブロック図である。
【図４】１次元ＭＭＡ法を説明する模式図である。
【図５】２次元ＭＭＡ法を説明する模式図である。
【図６】１次元移流方程式の数値計算結果を示すグラフである。
【図７】２次元移流方程式（剛体回転）の数値計算結果を示すグラフである。
【図８】ＭＭＡ法のブラソフ方程式及びマクスウェル方程式への適用例を示すグラフである。
【図９】ＭＭＡ法のブラソフ方程式及びマクスウェル方程式への適用例を示すグラフである。
【図１０】１次風上差分法の計算例を示すグラフである。
【図１１】１次風上差分法及びＣＩＰ法の概念を示す説明図である。
【発明を実施するための形態】
【００１２】
添付図を参照しながら、以下に、本発明の実施形態を詳細に説明する。
まず、図３を用いて、移流方程式解法装置１のハードウェア構成について説明する。移流方程式解析装置１は、ＣＰＵ（図３装置３０１）と、ＲＡＭ（３０２）と、ＨＤＤ（ハードディスクドライブ、３０３）と、マウス（３０４）と、キーボード（３０５）と、ディスプレイ（３０６）と、を備えている。また、各構成部はバス（３００）によってそれぞれ接続されている。
【００１３】
ここで、ＣＰＵ（３０１）は移流方程式解析装置１の全体の制御を司る。ＲＡＭ（３０２）はＣＰＵ（３０１）のワークエリアとして使用される。初期値記憶ステップ（図１に示すＳ１０１と図２に示すＳ２０１）における初期値の保持や値の時間更新（図１に示すＳ１０２からＳ１０５と図２に示すＳ２０２からＳ２０９）はＲＡＭ（３０２）を介して行われる。ＨＤＤ（３０３）はＣＰＵ（３０１）の制御にしたがって、データの読み込み・書き込みに使用される。出力ステップ（図１に示すＳ１０６と図２に示すＳ２１０）において、計算結果の保存はＨＤＤ（３０３）を介して行われる。マウス（３０４）はカーソルの移動や範囲選択、ウィンドウのサイズ変更などに使用される。キーボード（３０５）は、文字・数字などの入力のためのキーを備え、移流方程式解法プログラムはキーボードを介したコマンド入力によってＣＰＵ（３０１）により実行される。ディスプレイ（３０６）は、液晶ディスプレイなどが相当し、出力結果の表示に使用される。
【００１４】
続いて、図１を用いて、ＣＰＵ（３０１）が移流方程式解法プログラムを実行することにより実現される移流方程式を数値的に解く処理手順について説明する。処理手順は、初期値記憶ステップ（Ｓ１０１）と、上流点算出ステップ（Ｓ１０２）と、空間プロファイル算出ステップ（Ｓ１０３）と、移流項計算ステップ（Ｓ１０４）と、非移流項加算ステップ（Ｓ１０５）と、出力ステップ（Ｓ１０６）と、から構成されている。
【００１５】
初期値記憶ステップ（Ｓ１０１）では、空間グリッド上で定義された物理量ｆ（ｘ、ｔ0）のデータ群と、各グリッド間の０次からｍｍａｘ次までのモーメント積分量Ｍm（ｘ、ｔ0）（ｍ＝０、１、２、・・・、ｍｍａｘ−１、ｍｍａｘ）のデータ群と、最大時間ステップ数ｎｍａｘと、をマウス（３０４）及びキーボード（３０５）によりユーザーが入力し、ＲＡＭ（３０２）に記憶する。
【００１６】
ｎ＝１時間ステップでは、
上流点算出ステップ（Ｓ１０２）では、空間上で定義された速度場ｕから、上流点及び上流方向のグリッドを算出する。
【００１７】
空間プロファイル算出ステップ（Ｓ１０３）では、まず、注目するグリッドと、前記上流点算出ステップ（Ｓ１０２）で算出した上流方向のグリッドの間の、物理量ｆ（ｘ、ｔ0）と、各モーメント積分量Ｍm（ｘ，ｔ0）と、の空間プロファイルを、多項式で表現する。多項式を決定する各項の未知係数は、注目するグリッドと、その上流方向のグリッドと、で記憶されている物理量ｆ（ｘ、ｔ0）と、注目するグリッドとその上流方向のグリッドの間で記憶されている各モーメント積分量Ｍm（ｘ、ｔ0）と、の束縛条件から決定される。これにより、ｔ＝０におけるｆ（ｘ、ｔ0）と、Ｍm（ｘ、ｔ0）と、の空間プロファイルＦ（ｘ、ｔ0）と、Ｇm（ｘ、ｔ0）と、を算出する。
【００１８】
移流項計算ステップ（Ｓ１０４）では、時間ステップの中間値であるｆ（ｘ、ｔ*）と、Ｍm（ｘ、ｔ*）と、をｆ（ｘ、ｔ*）＝ Ｆ（ｘ−ｕ×ｄｔ、ｔ0）と、Ｍm（ｘ、ｔ*）＝ Ｇm（ｘ−ｕ×ｄｔ、ｔ0）と、により時間更新する（ｄｔ＝ｔ1−ｔ0は数値計算の時間刻みを表す）。
【００１９】
非移流項加算ステップ（Ｓ１０５）では、物理量ｆの非移流項Ｓ（ｘ、ｔ*）と、各モーメント積分量Ｍmの非移流項Ｓm（ｘ、ｔ*）とを、ｆ（ｘ、ｔ*）と、Ｍm（ｘ、ｔ*）と、から算出する。Ｓ（ｘ、ｔ*）と、Ｓm（ｘ、ｔ*）と、をｆ（ｘ、ｔ*）と、Ｍm（ｘ、ｔ*）と、に加算し、ｆ（ｘ、ｔ1）と、Ｍm（ｘ、ｔ1）と、を算出する。
【００２０】
時間ステップ数を更新し、ｎ＝２とする。
時間ステップｎ＝２以降では、１ステップ前の値（ｎ＝２の場合、ｎ＝１での値）を基にして、前記上流点算出ステップ（Ｓ１０２）から非移流項加算ステップ（Ｓ１０６）までを、ｎ＝ｎｍａｘになるまで繰り返す。
【００２１】
出力ステップ（Ｓ１０７）は、計算結果ｆ（ｘ、ｖ、ｔnmax）をＨＤＤ（３０３）に出力する。
続いて、図２を用いてＣＰＵ（３０１）がプラズマ物理現象を記述する方程式であるブラソフ方程式及びマクスウェル方程式解法プログラムを実行することにより実現される、ブラソフ方程式及びマクスウェル方程式を数値的に解く処理手順について説明する。ブラソフ方程式は、
【００２２】
【数１】
【００２３】
で表される、粒子種ｓの分布関数ｆsの位相空間上での移流を記述する線形移流方程式である。ここで、ベクトルｘは実空間での位置ベクトル、ベクトルｖは速度空間での位置ベクトル（速度場ベクトル）、ｑsは粒子種ｓの電荷、ｍsは粒子種ｓの質量、ベクトルＥ、ベクトルＢはそれぞれ電場ベクトル、磁場ベクトルを表す。上式は、次のように、実空間と速度空間での移流方程式に分けてそれぞれ計算を行う。
【００２４】
【数２】
【００２５】
但しベクトルＡは速度空間上における移流速度ベクトルである加速度場ベクトルであり、次式で計算される。
【００２６】
【数３】
【００２７】
実空間における移流は例えばＣＩＰ−ＣＳＬ２法などの汎用的な手法を用いて時間更新する。電場ベクトル、磁場ベクトルが与えられれば、前記ＭＭＡ法により速度空間上における移流方程式を解くことができる。電場ベクトル、磁場ベクトルの時間発展は別途、マクスウェル方程式
【００２８】
【数４】
【００２９】
によって記述される。ここで、ρe、ベクトルｊはそれぞれ電荷密度、電流密度ベクトルを表し、それぞれ、
【００３０】
【数５】
【００３１】
で定義される。右辺に現れる量は、ＭＭＡ法では独立に記憶・更新される０次及び１次のモーメント積分量である。この様にして得られた電荷密度と電流密度ベクトルを使ってマクスウェル方程式を解くことにより、電磁場の時間発展を追うことができる。
【００３２】
図２を用いて、本実施の形態にかかるブラソフ方程式及びマクスウェル方程式解法装置の処理手順について説明する。処理手順は、初期値記憶ステップ（Ｓ２０１）と、実空間時間更新ステップ（Ｓ２０２）と、電荷・電流密度算出ステップ（Ｓ２０３）と、電磁場更新ステップ（Ｓ２０４）と、加速度場算出ステップ（Ｓ２０５）と、速度空間上流点算出ステップ（Ｓ２０６）と、速度空間プロファイル算出ステップ（Ｓ２０８）と、速度空間移流項計算ステップ（Ｓ２０８）と、速度空間非移流項加算ステップ（Ｓ２０９）と、出力ステップ（Ｓ２１０）と、から構成されている。これらは、ＣＰＵ（３０１）によって実行される。
【００３３】
初期値記憶ステップ（Ｓ２０１）では、位相空間グリッド上で定義された分布関数ｆ（ｘ、ｖ、ｔ0）と、速度空間のグリッド間の０次からｍｍａｘ次までのモーメント積分量Ｍm（ｘ、ｖ、ｔ0）（ｍ＝０，１，２、・・・、ｍｍａｘ−１、ｍｍａｘ）と、実空間グリッド上に定義された速度場ベクトルｖと、実空間グリッド上に定義された電場ベクトルＥ（ｘ、ｔ0）と、磁場ベクトルＢ（ｘ、ｔ0）と、速度空間上における移流速度である加速度場ベクトルＡ（ｘ、ｖ、ｔ0）と、最大時間ステップ数ｎｍａｘと、をマウス（３０４）及びキーボード（３０５）によりユーザーが入力し、ＲＡＭ（３０２）に記憶する。
【００３４】
ｎ＝１時間ステップでは、
実空間時間更新ステップ（Ｓ２０２）では、速度ベクトルｖから、ｆ（ｘ、ｖ、ｔ0）と、Ｍm（ｘ、ｖ、ｔ0）と、の実空間上における移流を、例えば、ＣＩＰ−ＣＳＬ２法などの汎用的な手法を用いて時間更新し、第一中間時間ステップでの値ｆ（ｘ、ｖ、ｔ*）と、Ｍm（ｘ、ｖ、ｔ*）と、を得る。
【００３５】
電荷・電流密度算出ステップ（Ｓ２０３）では、更新した０次モーメントＭ0（ｘ、ｖ、ｔ*）と、１次モーメントＭ1（ｘ、ｖ、ｔ*）と、から電荷密度ρe（ｘ、ｔ*）と、電流密度ベクトルｊ（ｘ、ｔ*）と、を算出する。
【００３６】
電磁場更新ステップ（Ｓ２０４）では、算出した電荷密度ρe（ｘ、ｔ*）と、電流密度ベクトルｊ（ｘ、ｔ*）と、からマクスウェル方程式を解き、電場ベクトルＥ（ｘ、ｔ0）と、磁場ベクトルＢ（ｘ、ｔ0）と、から、時間更新後の電場ベクトルＥ（ｘ、ｔ1）と、磁場ベクトルＢ（ｘ、ｔ1）と、を算出する。
【００３７】
加速度場算出ステップ（Ｓ２０５）では、前記時間更新後の電場ベクトルＥ（ｘ、ｔ1）と、磁場ベクトルＢ（ｘ、ｔ1）と、から前記加速度場ベクトルＡ（ｘ、ｖ、ｔ1）を算出する。
【００３８】
速度空間上流点算出ステップ（Ｓ２０６）では、前記加速度場ベクトル算出ステップ（Ｓ２０５）で算出した加速度場ベクトルＡ（ｘ、ｖ、ｔ1）から、速度空間における上流点及び上流方向のグリッドを算出する。
【００３９】
速度空間プロファイル算出ステップ（Ｓ２０７）では、まず、注目するグリッドと、前記上流点算出ステップ（Ｓ２０６）で算出した上流方向のグリッドの間の、物理量ｆ（ｘ、ｖ、ｔ*）と、各モーメント積分量Ｍm（ｘ、ｖ、ｔ*）と、の速度空間プロファイルを、多項式で表現する。多項式を決定する各項の未知係数は、注目するグリッドと、その上流方向のグリッドと、で記憶されている物理量ｆ（ｘ、ｖ、ｔ*）と、注目するグリッドとその上流方向のグリッドの間で記憶されている各モーメント積分量Ｍm（ｘ、ｖ、ｔ*）と、の束縛条件から決定される。これにより、ｔ＝０におけるｆ（ｘ、ｖ、ｔ*）と、Ｍm（ｘ、ｖ、ｔ*）と、の速度空間プロファイルＦ（ｘ、ｖ、ｔ*）と、Ｇm（ｘ、ｖ、ｔ*）と、を算出する。
【００４０】
速度空間移流項計算ステップ（Ｓ２０８）では、前記加速度場算出ステップ（Ｓ２０５）で得られた加速度場ベクトルＡ（ｘ、ｖ、ｔ1）より、第二中間時間ステップの値であるｆ（ｘ、ｖ、ｔ**）と、Ｍm（ｘ、ｖ、ｔ**）と、をｆ（ｘ、ｖ、ｔ**）＝ Ｆ（ｘ、ｖ−Ａ×ｄｔ、ｔ*）と、Ｍm（ｘ、ｖ、ｔ**）＝ Ｇm（ｘ、ｖ−Ａ×ｄｔ、ｔ*）と、により時間更新する。
【００４１】
速度空間非移流項加算ステップ（Ｓ２０９）では、物理量ｆの非移流項Ｓ（ｘ、ｖ、ｔ**）と、各モーメント積分量Ｍmの非移流項Ｓm（ｘ、ｖ、ｔ**）とを、ｆ（ｘ、ｖ、ｔ**）と、Ｍm（ｘ、ｖ、ｔ**）と、から算出する。Ｓ（ｘ、ｖ、ｔ**）と、Ｓm（ｘ、ｖ、ｔ**）と、をｆ（ｘ，ｖ，ｔ**）と、Ｍm（ｘ，ｖ，ｔ**）と、に加算することにより、ｆ（ｘ、ｖ、ｔ1）と、Ｍm（ｘ、ｖ、ｔ1）と、を算出する。
【００４２】
時間ステップ数を更新し、ｎ＝２とする。
時間ステップｎ＝２以降では、１ステップ前の値（ｎ＝２の場合、ｎ＝１での値）を基にして、前記実空間時間更新ステップ（Ｓ２０２）から、前記速度空間非移流項加算ステップ（Ｓ２０９）まで、をｎ＝ｎｍａｘになるまで繰り返す。
【００４３】
出力ステップ（Ｓ２１０）は、計算結果であるｆ（ｘ、ｖ、ｔnmax）と、電場ベクトルＥ（ｘ、ｔnmax）と、磁場ベクトルＢ（ｘ、ｔnmax）と、をＨＤＤ（３０３）に出力する。
【００４４】
この発明により、ブラソフ方程式とマクスウェル方程式を数値的に解くことが出来る。
以下では、移流方程式解法の定式化・実装について説明する。
＜１次元ＭＭＡ法＞
１次元ＭＭＡ法は、物理量とその０次、１次、２次モーメント積分量を用いて、１次元移流方程式を数値的に解く解法である。扱う基礎方程式は保存型の移流方程式で、次の通りである。
【００４５】
【数６】
【００４６】
ここで、時刻ｔ、空間座標ｘに対してｆが物理量、ｕが速度場、
【００４７】
【数７】
【００４８】
はｆのｍ次モーメントの空間積分量（モーメント積分量）である。式（１．２）−（１．３）は式（１．１）から厳密に導かれる。
コンピュータは、図４で示すように、グリッド上の物理量ｆi及びグリッド間のモーメント積分量
【００４９】
【数８】
【００５０】
の、計４つの変数をＲＡＭ（３０２）に記憶する（初期値記憶ステップＳ１０１）。また、コンピュータは、グリッド上で速度場ｕi及びその空間微分量（∂ｕ／∂ｘ）iを別途求め、ＲＡＭ（３０２）に記憶する。
【００５１】
物理量ｆとモーメント積分量Ｍmのグリッド間の空間プロファイルを、下記の区分化された多項式で近似する。
【００５２】
【数９】
【００５３】
ここで、
【００５４】
【数１０】
【００５５】
記憶量であるグリッド上の物理量及びグリッド間のモーメント積分量を束縛条件として用いて、前記多項式の係数Ｃkiを決定する。束縛条件は、
【００５６】
【数１１】
【００５７】
で与えられる。ここで、
【００５８】
【数１２】
【００５９】
は、ｉ番目のグリッドから見た上流方向のグリッド及びセルの番号であり、
【００６０】
【数１３】
【００６１】
は、ｉ番目の座標ｘiから見た上流点までの距離であり、Δｔは数値計算の時間刻みであり、ｘi＋ξiが上流点の位置である（上流点算出ステップＳ１０２）。以上より、係数Ｃkiが次の様に求まる。
【００６２】
【数１４】
【００６３】
ここで、Δｘ = ｘiup−ｘiである。
このようにして、上記係数を要素とする物理量及びモーメント積分量のグリッド間の関数を空間プロファイル（１．４）、（１．５）として決定する（空間プロファイル算出ステップＳ１０３）。続いて、決定された空間プロファイルを用いて、グリッド上の物理量及びグリッド間のモーメント積分量を、以下に述べる方法を用いて時間更新する。基礎方程式（１．１）−（１．３）は、移流項（左辺第２項）と非移流項（右辺）と、を含む。まず移流項計算は、セミラグランジュ法と呼ばれる、上流点の値を速度場に沿って移流する方法を用いる（移流項計算ステップＳ１０４）。
【００６４】
【数１５】
【００６５】
ここで、左上添え字は時間ステップを表わし、アスタリスクは、該当する量が移流項計算後の中間値であることを示す。式（１．１３）、（１．１４）は、０次、１次、２次モーメント積分量の保存を保証する。次にグリッド上の物理量及びグリッド間のモーメント積分量の非移流項を、移流項計算後のグリッド上の物理量及びグリッド間のモーメント積分量を用いて、以下に述べる方法で算出する。
【００６６】
【数１６】
【００６７】
このようにして決定された非移流項を用いて、グリッド上の物理量及びグリッド間のモーメント積分量に非移流項を加算する（非移流項加算ステップＳ１０５）。
【００６８】
【数１７】
【００６９】
以上のステップを経て、ｎ時間ステップのグリッド上の物理量 nｆi及びグリッド間のモーメント積分量 nＭmi+1/2から、（ｎ＋１）時間ステップのグリッド上の物理量（１．１５）及びグリッド間のモーメント積分量（１．１３）、（１．１６）が得られる。
＜２次元ＭＭＡ法＞
２次元ＭＭＡ法は、物理量とその０次モーメント積分量、各次元方向の１次、２次モーメント積分量を用いて、２次元移流方程式を数値的に解く解法である。扱う基礎方程式は保存型の移流方程式で、次の通りである。
【００７０】
【数１８】
【００７１】
ここで、時刻ｔ、空間座標ｘ、ｙに対してｆが物理量、ｕ、ｖがｘ、ｙ方向の速度場である。
【００７２】
【数１９】
【００７３】
はそれぞれ、ｆの０次モーメントの空間積分量、ｘ方向のｍ次モーメントの空間積分量、ｙ方向のｍ次モーメントの空間積分量（モーメント積分量）である。今後、都合上、Ｍx0=Ｍy0=Ｍ0として表記することがある。式（２．２）−（２．４）は式（２．１）から厳密に導かれる。
【００７４】
コンピュータは、図５で示すように、グリッド上の物理量ｆij及びグリッド間のモーメント積分量
【００７５】
【数２０】
【００７６】
の、計６つの変数をＲＡＭ（３０２）に記憶する（初期値記憶ステップＳ１０１）。また、コンピュータは、グリッド上でｘ、ｙ方向の速度場ｕij、 ｖij及びその空間微分量（∂ｕ／∂ｘ）ij、（∂ｖ／∂ｙ）ijを別途求め、ＲＡＭ（３０２）に記憶する。
【００７７】
物理量ｆとモーメント積分量Ｍxm、Ｍymのグリッド間のプロファイルを、下記の区分化された多項式で近似する。
【００７８】
【数２１】
【００７９】
今後Ｇx0＝Ｇy0＝Ｇ0として表記することがある。
記憶量であるグリッド上の物理量とグリッド間のモーメント積分量を束縛条件として用いて、前記多項式の係数Ｃlkijを決定する。束縛条件は、
【００８０】
【数２２】
【００８１】
で与えられる。ここで、
【００８２】
【数２３】
【００８３】
は、（ｉ，ｊ）番目のグリッドから見た上流ｘ方向のグリッド及びセルの番号と上流ｙ方向のグリッド及びセルの番号であり、
【００８４】
【数２４】
【００８５】
は、（ｉ，ｊ）番目の座標（ｘi，ｙj）から見た上流点までのｘ方向及びｙ方向の距離であり、（ｘi＋ξij，ｙj＋ηij）が上流点の位置である（上流点算出ステップＳ１０２）。以上より、係数Ｃlkijが次の様に求まる。
【００８６】
【数２５】
【００８７】
ここで、Δｘ = ｘiup−ｘi、Δｙ = ｙjup−ｙjである。
このようにして、上記係数を要素とする物理量及びモーメント積分量のグリッド間の関数を空間プロファイル（２．８）−（２．１０）として決定する（空間プロファイル算出ステップＳ１０３）。続いて、決定された空間プロファイルを用いて、グリッド上の物理量及びグリッド間のモーメント積分量を、以下に述べる方法を用いて時間更新する。基礎方程式（２．１）−（２．４）は、移流項（左辺第２項、３項）と非移流項（右辺）と、を含む。まずグリッド上の物理量の時間更新は、セミラグランジュ法による移流項計算（移流項計算ステップＳ１０４）と、非移流項計算（非移流項加算ステップＳ１０５）と、をルンゲ・クッタ法を用いて行う。
【００８８】
【数２６】
【００８９】
ここで、ｋmaxはルンゲ・クッタ法の次数を表す。式（２．２２）、（２．２３）は特性曲線方程式であり、この式から、（ｉ，ｊ）番目の座標（ｘi，ｙj）から見た上流点のルンゲ・クッタ法のサブステップ毎の座標（ｘi＋ kξij，ｙj＋ kηij）が決定される。
【００９０】
kＳijはルンゲ・クッタ法のサブステップ毎のｆの非移流項であり、 kαはルンゲ・クッタ法の重み係数である。式（２．２４）がｆの時間更新式であり、右辺第１項がセミラグランジュ法による移流項計算（移流項計算ステップＳ１０４）、右辺第２項が非移流項の加算（非移流項加算ステップＳ１０５）である。この式から、ルンゲ・クッタ法のサブステップ毎のグリッド上の物理量 kｆijが求まり、ｋ＝ｋmaxのとき、kmaxｆij＝n+1ｆijであり、（ｎ＋１）時間ステップのグリッド上の物理量n+1ｆijが得られる。
【００９１】
次にグリッド間のモーメント積分量の時間更新方法について述べる。多次元移流問題では、グリッド間のモーメント積分量のセミラグランジュ法による移流計算は非常に煩雑になるため、有限体積法による移流項計算（移流項計算ステップＳ１０４）と、非移流項加算（非移流項加算ステップＳ１０５）と、を以下で行う。式（２．２）−（２．４）を次の保存形式で表現する。
【００９２】
【数２７】
【００９３】
ここで式（２．２５）−（２．２７）の右辺第１項から第４項が移流項であり、物理量の線積分量で表現されている。そこで先に決定された物理量の空間プロファイル（２．８）を用いて、線積分量を構築する。
【００９４】
【数２８】
【００９５】
これらを用いて式（２．２５）−（２．２７）の移流項の計算を行う（移流項計算ステップＳ１０４）。式（２．２５）−（２．２７）は保存型で書かれているため、０次、１次、２次モーメント積分量の保存が保証されている。続いて、ｘ、ｙ方向の１次及び２次モーメント積分量の非移流項である式（２．２６）、（２．２７）の右辺第５項を
【００９６】
【数２９】
【００９７】
により算出し、式（２．２６）、（２．２７）に加算する（非移流項加算ステップＳ１０５）。以上より、式（２．２５）−（２．２７）を有限体積法で計算することができ、結果、（ｎ＋１）時間ステップのグリッド間のモーメント積分量 n+1Ｍxm、 n+1Ｍymが得られる。但し計算を安定に行うために、式（２．２５）−（２．２７）は、ルンゲ・クッタ法を用いて時間更新する必要がある。
＜計算例＞
前記ＭＭＡ法による移流方程式の解法例を以下に示す。
【００９８】
図６上段はガウス分布形状の１次元線形移流計算結果であり、左よりＣＩＰ−ＣＳＬ２法、ＭＭＡ法の結果を表す。図中で、横軸が１次元空間で、「＋」は計算結果を、破線は解析解を表す。下段は解析解からの誤差を表す。ＣＩＰ−ＣＳＬ２法では、移流と共にプロファイルがなまってしまい、数値拡散の影響が大きく出ている。また、裾野のあたりでは数値振動が出ていることがわかる。一方、ＭＭＡ法はプロファイルを正確に再現しており、誤差の大きさもＣＩＰ−ＣＳＬ２法に比べて３桁以上小さい。
【００９９】
図７は、ガウス分布形状の２次元剛体回転問題の数値計算結果である。軸が２次元空間、背景及び等高線が物理量を表す。プロファイルを最上段のグラフで示されている矢印の方向に回転し続けた。左の列がＣＩＰ−ＣＳＬ２法の、右の列がＭＭＡ法の計算結果である。それぞれの使用メモリー量を合わせるために、グリッド数を変えている。上から順に、初期状態、３０回転後、３００回転後の結果を表す。２次元剛体回転問題では、数値拡散の影響が顕著に現れる。ＣＩＰ−ＣＳＬ２法による剛体回転問題結果では、３０回転後においてすでにプロファイルの数値拡散が見られ、３００回転後では原形をとどめていない。一方、ＭＭＡ法は、３００回転後においても初期の形状を保ったまま、剛体回転運動を数値的に記述することができている。
【０１００】
図８は、ＭＭＡ法のブラソフ方程式への適用例である。初期に互いに異なる向きの電子ビームを与えた。（ａ）から（ｄ）の順に、その位相空間分布の時間発展を示す。横軸が１次元実空間、縦軸が１次元速度空間、背景及び等高線が位相空間分布を表す。プラズマ物理分野でよく知られる２流体不安定が励起されており、その特徴的な構造である、位相空間上における電子ホールの形成が再現できている。
【０１０１】
図９は、ＭＭＡ法のブラソフ方程式への適用例である。イオン、電子が存在する中での、磁場に垂直方向に伝播するプラズマ波動である、バーンシュタイン波の分散関係である。横軸が波数、縦軸が周波数、背景が電場強度を表す。本計算は、プラズマの旋回運動（剛体回転）を長時間計算する必要があるため、従来のブラソフシミュレーション手法では再現することができなかった。左側がＭＭＡ法で求めたブラソフシミュレーション結果で、右側が粒子法の一種である電磁粒子（ＰＩＣ）コードで計算した結果である。両者共に同じ結果が得られ、ＭＭＡ法による計算が正しいことがわかる。また、ブラソフシミュレーションは電磁粒子計算に比べて数値ノイズが極めて少ないことが特徴である。
【産業上の利用可能性】
【０１０２】
移流方程式は物質の移動を記述する基本的な方程式であるため、その数値的解法は、数値流体力学（ＣＦＤ）や、宇宙・核融合分野におけるプラズマ現象への適用のみならず、大気流体現象、燃焼流体現象、コンピュータ画像処理など、幅広い分野に応用可能である。
【符号の説明】
【０１０３】
１…移流方程式解法装置、３００…バス、３０１…ＣＰＵ、３０２…ＲＡＭ、３０３…ＨＤＤ、３０４…マウス、３０５…キーボード、３０６…ディスプレイ
【技術分野】
【０００１】
本発明は、物理現象を表現する基本的な方程式である移流方程式を数値的に解く方法に関する。更には、移流方程式を数値的に解くためのプログラムに関する。
【背景技術】
【０００２】
移流方程式とは、ある物理量がある速度場によって移動（伝播）する現象を記述した方程式である。この移流現象は流体運動を初め、多くの物理現象に関わるため、現象を記述する基本的かつ重要な方程式であるといえる。移流方程式には解析解が存在する一方、空間上に離散化されたグリッド上で移流方程式を数値的に解く場合は、一転して多くの困難を伴うことが古くから知られている。その重要性かつ困難性から、移流方程式を数値的に解く手法を開発することは、学術面のみならず、産業上においても重要な課題として残されている。
【０００３】
移流方程式を数値的に解く場合、移流させる物理量をコンピュータのメモリー空間上に離散化されたグリッド上に記憶するため、グリッド間の情報は失われる。そこで、記憶してある量を用いてグリッド間のプロファイルを再構築する必要がある。古典的な１次風上差分法は、注目するグリッド上での移流速度から、風上方向にあるグリッド間の情報を線形補間することによって、物理量のプロファイルを再構築する手法である。この手法は数値的に安定であり、古典的手法の代表的なものである。図１０に１次風上差分法の計算例を示す。横軸が実空間、縦軸が物理量を表し、「＋」が計算結果で、破線が解析解である。図でわかるように、本計算手法は数値粘性が大きいため、時間と共にプロファイルが減衰してしまい、長時間計算には適さないことが知られている。一方、近年開発された手法として、ＣＩＰ（Ｃｏｎｓｔｒａｉｎｅｄ Ｉｎｔｅｒｐｏｌａｔｉｏｎ Ｐｒｏｆｉｌｅ）法が挙げられる（非特許文献１）。図１１に１次風上差分法及びＣＩＰ法の概念を示す。横軸が実空間、黒丸がグリッドを表す。グリッド上に離散化された物理量を数値的に移流させる際（図１１ａ）、１次風上差分法では物理量のみを移流させようとするため（図１１ｂ）、プロファイルが減衰してしまう（図１１ｃ）。ＣＩＰ法は物理量のみならず、物理量の空間プロファイルの勾配（微分値）の時間発展をも同時に解き進める手法である。図１１ｄで示されているように、風上方向にあるグリッド間のプロファイルを再構築する際に、物理量と勾配の情報も使うため、高精度かつ安定に解く手法として近年注目を浴びてきている。ＣＩＰ法の派生である保存保証型ＣＩＰ法（ＣＩＰ−ＣＳＬ２法、非特許文献２）は、物理量とそのグリッド間の積分量を記憶する。グリッド上の物理量とグリッド間の積分量の情報を使い、グリッド間のプロファイルを高精度で再構築する移流方程式解法である。この解法は、ＣＩＰ法と同程度の精度を達成し、またＣＩＰ法に対して質量保存が保証されるメリットがある。しかし、不連続面において若干の数値振動が生じる点は、極めて高精度な移流方程式の解を要求される分野、例えば、プラズマ物理学の分野では依然大きな問題点として残されている。
【０００４】
前記プラズマ物理学の分野では、上に挙げた保存性、無振動性を特徴とする手法が近年提案されている。非特許文献３では、プラズマ物理現象を記述するブラソフ方程式を解くことを目的として、ＰＦＣ（Ｐｏｓｉｔｉｖｅ ａｎｄ Ｆｌｕｘ Ｃｏｎｓｅｒｖａｔｉｖｅ）法が開発されている。本手法は移流方程式を保存形で解くことにより、保存則を満たし、不連続面近傍では流束に制限をかけることにより、数値振動を少なくする手法である。さらに近年では、流束制限手法を工夫することにより、無振動性を保証する手法も開発された（非特許文献４）。以上の手法は、極めて高い精度が要求されるプラズマ物理分野における移流方程式である、ブラソフ方程式の解法として知られている。しかしながら、数値拡散の抑制については未だ不十分であり、数値拡散のない（極めて少ない）移流方程式解法が求められている。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】タカシ・ヤベ（Takashi Yabe）、フェン・ショウ（Feng Xiao）、及びタカユキ・ウツミ（Takayuki Utsumi）、"ザ コンストレインド インターポレーション プロファイル メソッド フォー マルチフェーズ アナリシス（The Constrained Interpolation Profile Method for Multiphase Analysis）"、ジャーナル オブ コンピュテーショナル フィジクス（Journal of Computational Physics）、オランダ王国、アムステルダム、２００１年５月、第１６０巻、ｐ．５５６−５９３
【非特許文献２】タカシ・ヤベ（Takashi Yabe）、リョウタロウ・タナカ（Ryotaro Tanaka）、タカシ・ナカムラ（Takashi Nakamura）、及びフェン・ショウ（Feng Xiao）、マンスリー ウェザー レビュー（Monthly Weather Review）、アメリカ合衆国、ボストン、２００１年２月、第１２９巻、ｐ．３３２−３４４
【非特許文献３】フランシス・フィルべ（Francis Filbet）、エリック・ソネンドラッカー（Eric Sonnendrucker）、及びピエール・ベルランド（Pierre Bertrand）、"コンザバティブ ニューメリカル スキームズ フォー ザ ブラソフ エクエーション（Conservative Numerical Schemes for the Vlasov Equation）"、ジャーナル オブ コンピュテーショナル フィジクス（Journal of Computational Physics）、オランダ王国、アムステルダム、２００１年９月、第１７２巻、ｐ．１６６−１８７
【非特許文献４】タカユキ・ウメダ（Takayuki Umeda）、ケンタロウ・トガノ（Kentaro Togano）、及びタツキ・オギノ（Tatsuki Ogino）、"トゥーディメンショナル フルエレクトロマグネティック ブラソフ コード ウィズ コンサバティブ スキーム アンド イッツ アプリケーション トゥ マグネティック リコネクション（Two-dimensional full-electromagnetic Vlasov code with conservative scheme and its application to magnetic reconnection）", コンピューター フィジクスコミュニケーションズ（Computer Physics Communications）、オランダ王国、アムステルダム、２００９年３月、第１８０巻、ｐ．３６５−３７４
【発明の概要】
【発明が解決しようとする課題】
【０００６】
上述した従来の解法では、そのいずれも、高い精度を要求される分野、例えば、プラズマ物理学への適用という点では、プロファイル再構築精度は不十分である。例えば図６左列及び図７左列では、ＣＩＰ−ＣＳＬ２法を用いた１次元移流方程式及び２次元移流方程式の数値計算結果を示している。プロファイルがなまってしまい、解析解から大きく外れてしまっている。このように、数値振動や数値拡散が無視できなくなるといった欠点があった。これは、多次元問題（図７左列）で特に顕著である。このように、既存の解法を遥かに上回る高精度の移流方程式解法が必要とされている。このような数値解法及び当該解法を実現するためのプログラムを提供することが本発明の目的である。
【課題を解決するための手段】
【０００７】
数値精度の低さが原因で発生する数値振動や数値拡散は、言いかえれば物理量の情報エントロピーが増大してしまっていることに相当する。独立変数ｘに対して物理量ｆ（ｘ）を定義すると、一般的に情報エントロピーは∫ｆ（ｘ）×ｂ（ｆ（ｘ））ｄｘと、情報量の関数ｂを用いて表現される。そこでｂをｘについてテイラー展開すると、情報エントロピーはｆ（ｘ）のｍ次モーメント積分量Ｍm（ｘ）＝∫ｘmｆ（ｘ）ｄｘ／ｍ！（ここで、ｍは０からｍｍａｘまでの整数を表す。）の線形結合ΣＭm（ｘ）（ｍ＝０、・・・、ｍｍａｘ）で表現されることがわかる。よって、数値分散や数値拡散を抑制するためには、情報エントロピーの保存、つまり、各モーメント積分量が保存するように移流方程式を解くことが本質的である。例えば、熱平衡状態を表現するガウス（正規）分布の情報エントロピーは２次のモーメント積分量Ｍ2（ｘ）＝∫ｘ2ｆ（ｘ）ｄｘ／２で表現される。前記ｍｍａｘの選択は任意であるが、ガウス分布は最も基本的な物理量のプロファイルであることから、少なくとも２次までのモーメント積分量を保存させることが本質である。１次元移流方程式解法においては、ｍｍａｘ＝０は前記ＣＩＰ−ＣＳＬ２法に一致する。本発明に対応する実施例では、前記情報エントロピーの考察に基づき、ｍｍａｘ＝２としている。
【０００８】
以上の考察に基づき、本発明は、物理量とそのｍｍａｘ次までのモーメント積分量を独立に記憶し、更新することを特徴とする移流方程式解法であり、以下に説明する解法をマルチモーメント移流法（Ｍｕｌｔｉ−Ｍｏｍｅｎｔ Ａｄｖｅｃｔｉｏｎ法、以下ＭＭＡ法）と呼ぶ。
【０００９】
本解法は、初期時刻ｔ0での空間上の各グリッドにおける物理量ｆ（ｘ、ｔ0）（ここで、ベクトルｘは空間上での位置ベクトル）と、その物理量ｆのｍｍａｘ次までのモーメント積分量Ｍm（ｘ、ｔ0）と、最大時間ステップ数ｎｍａｘと、を記憶手段に記憶する初期値記憶ステップ（図１に示すＳ１０１）と、空間のグリッド上で定義された移流速度ベクトルｕから、上流点及び上流方向のグリッドを決める上流点算出ステップ（Ｓ１０２）と、前記物理量ｆ（ｘ、ｔn）及び前記各モーメント積分量Ｍm（ｘ、ｔn）の束縛条件を課して、物理量ｆ（ｘ、ｔn）の空間プロファイルＦ（ｘ、ｔn）（ここで、ｔnは第ｎ時間ステップでの時刻を表し、ｎは０以上の整数である。）及び各モーメント積分量Ｍm（ｘ、ｔn）の空間プロファイルＧm（ｘ、ｔn）と、を求める空間プロファイル算出ステップ（Ｓ１０３）と、前記空間プロファイルに基づく次式ｆ（ｘ、ｔ*）＝Ｆ（ｘ−ｕ×ｄｔ、ｔn）及びＭm（ｘ、ｔ*）＝Ｇm（ｘ−ｕ×ｄｔ、ｔn）に従って（ここで、ｄｔ＝ｔn+1−ｔnは時間刻み幅、t*は移流後の中間時間ステップである。）、中間時間ステップでの物理量ｆ（ｘ、ｔ*）及び各モーメント積分量Ｍm（ｘ、ｔ*）を算出する移流項計算ステップ（Ｓ１０４）と、前記物理量ｆ（ｘ、ｔ*）及び各モーメント積分量Ｍm（ｘ、ｔ*）に基づき、物理量ｆの非移流項Ｓ（ｘ、ｔ*）及び各モーメント積分量Ｍmの非移流項Ｓm（ｘ、ｔ*）を算出し、これら物理量ｆの非移流項Ｓ（ｘ、ｔ*）及び各モーメント積分量Ｍmの非移流項Ｓm（ｘ、ｔ*）を前記物理量ｆ（ｘ、ｔ*）及び各モーメント積分量Ｍm（ｘ、ｔ*）に加算して、次の時間ステップに対応する物理量ｆ（ｘ、ｔn+1）及び各モーメント積分量Ｍm（ｘ、ｔn+1）を算出する非移流項加算ステップ（Ｓ１０５）と、計算結果を出力する出力ステップ（Ｓ１０６）と、から構成されている。本発明により、数値拡散と数値振動が極めて少ない移流方程式の解を得ることが出来る。
【００１０】
プラズマ物理分野での移流方程式であるブラソフ方程式及びマクスウェル方程式の解法は、初期時刻ｔ0での位相空間上におけるプラズマ分布関数ｆ（ｘ、ｖ、ｔ0）（ここでベクトルｘは実空間、ベクトルｖは速度空間での位置ベクトルを表す。）と、その各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ0）と、実空間上で定義された電場ベクトルＥ（ｘ、ｔ0）と、磁場ベクトルＢ（ｘ、ｔ0）と、最大時間ステップ数ｎｍａｘと、を記憶手段に記憶する初期値記憶ステップ（図２に示すＳ２０１）と、プラズマ分布関数ｆ（ｘ、ｖ、ｔn）と、速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔn）と、から実空間上における移流更新を行い、第一中間時間ステップでのプラズマ分布関数ｆ（ｘ、ｖ、ｔ*）と、速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ*）と、を算出する実空間時間更新ステップ（Ｓ２０２）と、０次の速度空間モーメント積分量Ｍ0（ｘ、ｖ、ｔ*）と、１次の速度空間モーメント積分量Ｍ1（ｘ、ｖ、ｔ*）と、から、電荷密度と、電流密度と、を求める電荷・電流密度算出ステップ（Ｓ２０３）と、マクスウェル方程式を解くことにより、前記電荷・電流密度と、電場ベクトルＥ（ｘ、ｔn）と、磁場ベクトルＢ（ｘ、ｔn）と、から次の時間ステップに対応する電場ベクトルＥ（ｘ、ｔn+1）と、磁場ベクトルＢ（ｘ、ｔn+1）と、を算出する電磁場更新ステップ（Ｓ２０４）と、前記時間更新後の電場ベクトルＥ（ｘ、ｔn+1）と、磁場ベクトルＢ（ｘ、ｔn+1）と、から速度空間上における移流速度ベクトルである加速度場ベクトルＡ（ｘ、ｖ、ｔn+1）を算出する加速度場算出ステップ（Ｓ２０５）と、速度空間のグリッド上で定義された加速度ベクトルＡ（ｘ、ｖ、ｔn+1）から、速度空間上の上流点及び上流方向のグリッドを決める速度空間上流点算出ステップ（Ｓ２０６）と、速度空間上でプラズマ分布関数ｆ（ｘ、ｖ、ｔ*）及び前記各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ*）の束縛条件を課して、プラズマ分布関数ｆ（ｘ、ｖ、ｔ*）の速度空間プロファイルＦ（ｘ、ｖ、ｔ*）及び各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ*）の速度空間プロファイルＧm（ｘ、ｖ、ｔ*）を求める速度空間プロファイル算出ステップ（Ｓ２０７）と、前記空間プロファイルに基づく次式ｆ（ｘ、ｖ、ｔ**）＝Ｆ（ｘ、ｖ−Ａ×ｄｔ、ｔ*）及びＭm（ｘ、ｖ、ｔ**）＝Ｇm（ｘ、ｖ−Ａ×ｄｔ、ｔ*）に従って、第二中間時間ステップでのプラズマ分布関数ｆ（ｘ、ｖ、ｔ**）及び各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ**）を算出する速度空間移流項計算ステップ（Ｓ２０８）と、前記プラズマ分布関数ｆ（ｘ、ｖ、ｔ**）及び各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ**）に基づき、プラズマ分布関数ｆの非移流項Ｓ（ｘ、ｖ、ｔ**）及び各速度空間モーメント積分量Ｍmの非移流項Ｓm（ｘ、ｖ、ｔ**）を算出し、これらプラズマ分布関数ｆの非移流項Ｓ（ｘ、ｖ、ｔ**）及び各速度空間モーメント積分量Ｍmの非移流項Ｓm（ｘ、ｖ、ｔ**）を前記プラズマ分布関数ｆ（ｘ、ｖ、ｔ**）及び各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ**）に加算して、次の時間ステップに対応するプラズマ分布関数ｆ（ｘ、ｖ、ｔn+1）及び各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔn+1）を算出する速度空間非移流項加算ステップ（Ｓ２０９）と、計算結果を出力する出力ステップ（Ｓ２１０）と、から構成されている。本発明により、従来の手法に比べて極めて高い信頼性を持つブラソフ方程式及びマクスウェル方程式の解を得ることができる。
【図面の簡単な説明】
【００１１】
【図１】本実施の形態における移流方程式解法装置における処理手順を示すフローチャートである。
【図２】本実施の形態におけるブラソフ方程式及びマクスウェル方程式解法装置における処理手順を示すフローチャートである。
【図３】本実施の形態における移流方程式及びブラソフ・マクスウェル方程式解法装置のハードウェア構成を示すブロック図である。
【図４】１次元ＭＭＡ法を説明する模式図である。
【図５】２次元ＭＭＡ法を説明する模式図である。
【図６】１次元移流方程式の数値計算結果を示すグラフである。
【図７】２次元移流方程式（剛体回転）の数値計算結果を示すグラフである。
【図８】ＭＭＡ法のブラソフ方程式及びマクスウェル方程式への適用例を示すグラフである。
【図９】ＭＭＡ法のブラソフ方程式及びマクスウェル方程式への適用例を示すグラフである。
【図１０】１次風上差分法の計算例を示すグラフである。
【図１１】１次風上差分法及びＣＩＰ法の概念を示す説明図である。
【発明を実施するための形態】
【００１２】
添付図を参照しながら、以下に、本発明の実施形態を詳細に説明する。
まず、図３を用いて、移流方程式解法装置１のハードウェア構成について説明する。移流方程式解析装置１は、ＣＰＵ（図３装置３０１）と、ＲＡＭ（３０２）と、ＨＤＤ（ハードディスクドライブ、３０３）と、マウス（３０４）と、キーボード（３０５）と、ディスプレイ（３０６）と、を備えている。また、各構成部はバス（３００）によってそれぞれ接続されている。
【００１３】
ここで、ＣＰＵ（３０１）は移流方程式解析装置１の全体の制御を司る。ＲＡＭ（３０２）はＣＰＵ（３０１）のワークエリアとして使用される。初期値記憶ステップ（図１に示すＳ１０１と図２に示すＳ２０１）における初期値の保持や値の時間更新（図１に示すＳ１０２からＳ１０５と図２に示すＳ２０２からＳ２０９）はＲＡＭ（３０２）を介して行われる。ＨＤＤ（３０３）はＣＰＵ（３０１）の制御にしたがって、データの読み込み・書き込みに使用される。出力ステップ（図１に示すＳ１０６と図２に示すＳ２１０）において、計算結果の保存はＨＤＤ（３０３）を介して行われる。マウス（３０４）はカーソルの移動や範囲選択、ウィンドウのサイズ変更などに使用される。キーボード（３０５）は、文字・数字などの入力のためのキーを備え、移流方程式解法プログラムはキーボードを介したコマンド入力によってＣＰＵ（３０１）により実行される。ディスプレイ（３０６）は、液晶ディスプレイなどが相当し、出力結果の表示に使用される。
【００１４】
続いて、図１を用いて、ＣＰＵ（３０１）が移流方程式解法プログラムを実行することにより実現される移流方程式を数値的に解く処理手順について説明する。処理手順は、初期値記憶ステップ（Ｓ１０１）と、上流点算出ステップ（Ｓ１０２）と、空間プロファイル算出ステップ（Ｓ１０３）と、移流項計算ステップ（Ｓ１０４）と、非移流項加算ステップ（Ｓ１０５）と、出力ステップ（Ｓ１０６）と、から構成されている。
【００１５】
初期値記憶ステップ（Ｓ１０１）では、空間グリッド上で定義された物理量ｆ（ｘ、ｔ0）のデータ群と、各グリッド間の０次からｍｍａｘ次までのモーメント積分量Ｍm（ｘ、ｔ0）（ｍ＝０、１、２、・・・、ｍｍａｘ−１、ｍｍａｘ）のデータ群と、最大時間ステップ数ｎｍａｘと、をマウス（３０４）及びキーボード（３０５）によりユーザーが入力し、ＲＡＭ（３０２）に記憶する。
【００１６】
ｎ＝１時間ステップでは、
上流点算出ステップ（Ｓ１０２）では、空間上で定義された速度場ｕから、上流点及び上流方向のグリッドを算出する。
【００１７】
空間プロファイル算出ステップ（Ｓ１０３）では、まず、注目するグリッドと、前記上流点算出ステップ（Ｓ１０２）で算出した上流方向のグリッドの間の、物理量ｆ（ｘ、ｔ0）と、各モーメント積分量Ｍm（ｘ，ｔ0）と、の空間プロファイルを、多項式で表現する。多項式を決定する各項の未知係数は、注目するグリッドと、その上流方向のグリッドと、で記憶されている物理量ｆ（ｘ、ｔ0）と、注目するグリッドとその上流方向のグリッドの間で記憶されている各モーメント積分量Ｍm（ｘ、ｔ0）と、の束縛条件から決定される。これにより、ｔ＝０におけるｆ（ｘ、ｔ0）と、Ｍm（ｘ、ｔ0）と、の空間プロファイルＦ（ｘ、ｔ0）と、Ｇm（ｘ、ｔ0）と、を算出する。
【００１８】
移流項計算ステップ（Ｓ１０４）では、時間ステップの中間値であるｆ（ｘ、ｔ*）と、Ｍm（ｘ、ｔ*）と、をｆ（ｘ、ｔ*）＝ Ｆ（ｘ−ｕ×ｄｔ、ｔ0）と、Ｍm（ｘ、ｔ*）＝ Ｇm（ｘ−ｕ×ｄｔ、ｔ0）と、により時間更新する（ｄｔ＝ｔ1−ｔ0は数値計算の時間刻みを表す）。
【００１９】
非移流項加算ステップ（Ｓ１０５）では、物理量ｆの非移流項Ｓ（ｘ、ｔ*）と、各モーメント積分量Ｍmの非移流項Ｓm（ｘ、ｔ*）とを、ｆ（ｘ、ｔ*）と、Ｍm（ｘ、ｔ*）と、から算出する。Ｓ（ｘ、ｔ*）と、Ｓm（ｘ、ｔ*）と、をｆ（ｘ、ｔ*）と、Ｍm（ｘ、ｔ*）と、に加算し、ｆ（ｘ、ｔ1）と、Ｍm（ｘ、ｔ1）と、を算出する。
【００２０】
時間ステップ数を更新し、ｎ＝２とする。
時間ステップｎ＝２以降では、１ステップ前の値（ｎ＝２の場合、ｎ＝１での値）を基にして、前記上流点算出ステップ（Ｓ１０２）から非移流項加算ステップ（Ｓ１０６）までを、ｎ＝ｎｍａｘになるまで繰り返す。
【００２１】
出力ステップ（Ｓ１０７）は、計算結果ｆ（ｘ、ｖ、ｔnmax）をＨＤＤ（３０３）に出力する。
続いて、図２を用いてＣＰＵ（３０１）がプラズマ物理現象を記述する方程式であるブラソフ方程式及びマクスウェル方程式解法プログラムを実行することにより実現される、ブラソフ方程式及びマクスウェル方程式を数値的に解く処理手順について説明する。ブラソフ方程式は、
【００２２】
【数１】
【００２３】
で表される、粒子種ｓの分布関数ｆsの位相空間上での移流を記述する線形移流方程式である。ここで、ベクトルｘは実空間での位置ベクトル、ベクトルｖは速度空間での位置ベクトル（速度場ベクトル）、ｑsは粒子種ｓの電荷、ｍsは粒子種ｓの質量、ベクトルＥ、ベクトルＢはそれぞれ電場ベクトル、磁場ベクトルを表す。上式は、次のように、実空間と速度空間での移流方程式に分けてそれぞれ計算を行う。
【００２４】
【数２】
【００２５】
但しベクトルＡは速度空間上における移流速度ベクトルである加速度場ベクトルであり、次式で計算される。
【００２６】
【数３】
【００２７】
実空間における移流は例えばＣＩＰ−ＣＳＬ２法などの汎用的な手法を用いて時間更新する。電場ベクトル、磁場ベクトルが与えられれば、前記ＭＭＡ法により速度空間上における移流方程式を解くことができる。電場ベクトル、磁場ベクトルの時間発展は別途、マクスウェル方程式
【００２８】
【数４】
【００２９】
によって記述される。ここで、ρe、ベクトルｊはそれぞれ電荷密度、電流密度ベクトルを表し、それぞれ、
【００３０】
【数５】
【００３１】
で定義される。右辺に現れる量は、ＭＭＡ法では独立に記憶・更新される０次及び１次のモーメント積分量である。この様にして得られた電荷密度と電流密度ベクトルを使ってマクスウェル方程式を解くことにより、電磁場の時間発展を追うことができる。
【００３２】
図２を用いて、本実施の形態にかかるブラソフ方程式及びマクスウェル方程式解法装置の処理手順について説明する。処理手順は、初期値記憶ステップ（Ｓ２０１）と、実空間時間更新ステップ（Ｓ２０２）と、電荷・電流密度算出ステップ（Ｓ２０３）と、電磁場更新ステップ（Ｓ２０４）と、加速度場算出ステップ（Ｓ２０５）と、速度空間上流点算出ステップ（Ｓ２０６）と、速度空間プロファイル算出ステップ（Ｓ２０８）と、速度空間移流項計算ステップ（Ｓ２０８）と、速度空間非移流項加算ステップ（Ｓ２０９）と、出力ステップ（Ｓ２１０）と、から構成されている。これらは、ＣＰＵ（３０１）によって実行される。
【００３３】
初期値記憶ステップ（Ｓ２０１）では、位相空間グリッド上で定義された分布関数ｆ（ｘ、ｖ、ｔ0）と、速度空間のグリッド間の０次からｍｍａｘ次までのモーメント積分量Ｍm（ｘ、ｖ、ｔ0）（ｍ＝０，１，２、・・・、ｍｍａｘ−１、ｍｍａｘ）と、実空間グリッド上に定義された速度場ベクトルｖと、実空間グリッド上に定義された電場ベクトルＥ（ｘ、ｔ0）と、磁場ベクトルＢ（ｘ、ｔ0）と、速度空間上における移流速度である加速度場ベクトルＡ（ｘ、ｖ、ｔ0）と、最大時間ステップ数ｎｍａｘと、をマウス（３０４）及びキーボード（３０５）によりユーザーが入力し、ＲＡＭ（３０２）に記憶する。
【００３４】
ｎ＝１時間ステップでは、
実空間時間更新ステップ（Ｓ２０２）では、速度ベクトルｖから、ｆ（ｘ、ｖ、ｔ0）と、Ｍm（ｘ、ｖ、ｔ0）と、の実空間上における移流を、例えば、ＣＩＰ−ＣＳＬ２法などの汎用的な手法を用いて時間更新し、第一中間時間ステップでの値ｆ（ｘ、ｖ、ｔ*）と、Ｍm（ｘ、ｖ、ｔ*）と、を得る。
【００３５】
電荷・電流密度算出ステップ（Ｓ２０３）では、更新した０次モーメントＭ0（ｘ、ｖ、ｔ*）と、１次モーメントＭ1（ｘ、ｖ、ｔ*）と、から電荷密度ρe（ｘ、ｔ*）と、電流密度ベクトルｊ（ｘ、ｔ*）と、を算出する。
【００３６】
電磁場更新ステップ（Ｓ２０４）では、算出した電荷密度ρe（ｘ、ｔ*）と、電流密度ベクトルｊ（ｘ、ｔ*）と、からマクスウェル方程式を解き、電場ベクトルＥ（ｘ、ｔ0）と、磁場ベクトルＢ（ｘ、ｔ0）と、から、時間更新後の電場ベクトルＥ（ｘ、ｔ1）と、磁場ベクトルＢ（ｘ、ｔ1）と、を算出する。
【００３７】
加速度場算出ステップ（Ｓ２０５）では、前記時間更新後の電場ベクトルＥ（ｘ、ｔ1）と、磁場ベクトルＢ（ｘ、ｔ1）と、から前記加速度場ベクトルＡ（ｘ、ｖ、ｔ1）を算出する。
【００３８】
速度空間上流点算出ステップ（Ｓ２０６）では、前記加速度場ベクトル算出ステップ（Ｓ２０５）で算出した加速度場ベクトルＡ（ｘ、ｖ、ｔ1）から、速度空間における上流点及び上流方向のグリッドを算出する。
【００３９】
速度空間プロファイル算出ステップ（Ｓ２０７）では、まず、注目するグリッドと、前記上流点算出ステップ（Ｓ２０６）で算出した上流方向のグリッドの間の、物理量ｆ（ｘ、ｖ、ｔ*）と、各モーメント積分量Ｍm（ｘ、ｖ、ｔ*）と、の速度空間プロファイルを、多項式で表現する。多項式を決定する各項の未知係数は、注目するグリッドと、その上流方向のグリッドと、で記憶されている物理量ｆ（ｘ、ｖ、ｔ*）と、注目するグリッドとその上流方向のグリッドの間で記憶されている各モーメント積分量Ｍm（ｘ、ｖ、ｔ*）と、の束縛条件から決定される。これにより、ｔ＝０におけるｆ（ｘ、ｖ、ｔ*）と、Ｍm（ｘ、ｖ、ｔ*）と、の速度空間プロファイルＦ（ｘ、ｖ、ｔ*）と、Ｇm（ｘ、ｖ、ｔ*）と、を算出する。
【００４０】
速度空間移流項計算ステップ（Ｓ２０８）では、前記加速度場算出ステップ（Ｓ２０５）で得られた加速度場ベクトルＡ（ｘ、ｖ、ｔ1）より、第二中間時間ステップの値であるｆ（ｘ、ｖ、ｔ**）と、Ｍm（ｘ、ｖ、ｔ**）と、をｆ（ｘ、ｖ、ｔ**）＝ Ｆ（ｘ、ｖ−Ａ×ｄｔ、ｔ*）と、Ｍm（ｘ、ｖ、ｔ**）＝ Ｇm（ｘ、ｖ−Ａ×ｄｔ、ｔ*）と、により時間更新する。
【００４１】
速度空間非移流項加算ステップ（Ｓ２０９）では、物理量ｆの非移流項Ｓ（ｘ、ｖ、ｔ**）と、各モーメント積分量Ｍmの非移流項Ｓm（ｘ、ｖ、ｔ**）とを、ｆ（ｘ、ｖ、ｔ**）と、Ｍm（ｘ、ｖ、ｔ**）と、から算出する。Ｓ（ｘ、ｖ、ｔ**）と、Ｓm（ｘ、ｖ、ｔ**）と、をｆ（ｘ，ｖ，ｔ**）と、Ｍm（ｘ，ｖ，ｔ**）と、に加算することにより、ｆ（ｘ、ｖ、ｔ1）と、Ｍm（ｘ、ｖ、ｔ1）と、を算出する。
【００４２】
時間ステップ数を更新し、ｎ＝２とする。
時間ステップｎ＝２以降では、１ステップ前の値（ｎ＝２の場合、ｎ＝１での値）を基にして、前記実空間時間更新ステップ（Ｓ２０２）から、前記速度空間非移流項加算ステップ（Ｓ２０９）まで、をｎ＝ｎｍａｘになるまで繰り返す。
【００４３】
出力ステップ（Ｓ２１０）は、計算結果であるｆ（ｘ、ｖ、ｔnmax）と、電場ベクトルＥ（ｘ、ｔnmax）と、磁場ベクトルＢ（ｘ、ｔnmax）と、をＨＤＤ（３０３）に出力する。
【００４４】
この発明により、ブラソフ方程式とマクスウェル方程式を数値的に解くことが出来る。
以下では、移流方程式解法の定式化・実装について説明する。
＜１次元ＭＭＡ法＞
１次元ＭＭＡ法は、物理量とその０次、１次、２次モーメント積分量を用いて、１次元移流方程式を数値的に解く解法である。扱う基礎方程式は保存型の移流方程式で、次の通りである。
【００４５】
【数６】
【００４６】
ここで、時刻ｔ、空間座標ｘに対してｆが物理量、ｕが速度場、
【００４７】
【数７】
【００４８】
はｆのｍ次モーメントの空間積分量（モーメント積分量）である。式（１．２）−（１．３）は式（１．１）から厳密に導かれる。
コンピュータは、図４で示すように、グリッド上の物理量ｆi及びグリッド間のモーメント積分量
【００４９】
【数８】
【００５０】
の、計４つの変数をＲＡＭ（３０２）に記憶する（初期値記憶ステップＳ１０１）。また、コンピュータは、グリッド上で速度場ｕi及びその空間微分量（∂ｕ／∂ｘ）iを別途求め、ＲＡＭ（３０２）に記憶する。
【００５１】
物理量ｆとモーメント積分量Ｍmのグリッド間の空間プロファイルを、下記の区分化された多項式で近似する。
【００５２】
【数９】
【００５３】
ここで、
【００５４】
【数１０】
【００５５】
記憶量であるグリッド上の物理量及びグリッド間のモーメント積分量を束縛条件として用いて、前記多項式の係数Ｃkiを決定する。束縛条件は、
【００５６】
【数１１】
【００５７】
で与えられる。ここで、
【００５８】
【数１２】
【００５９】
は、ｉ番目のグリッドから見た上流方向のグリッド及びセルの番号であり、
【００６０】
【数１３】
【００６１】
は、ｉ番目の座標ｘiから見た上流点までの距離であり、Δｔは数値計算の時間刻みであり、ｘi＋ξiが上流点の位置である（上流点算出ステップＳ１０２）。以上より、係数Ｃkiが次の様に求まる。
【００６２】
【数１４】
【００６３】
ここで、Δｘ = ｘiup−ｘiである。
このようにして、上記係数を要素とする物理量及びモーメント積分量のグリッド間の関数を空間プロファイル（１．４）、（１．５）として決定する（空間プロファイル算出ステップＳ１０３）。続いて、決定された空間プロファイルを用いて、グリッド上の物理量及びグリッド間のモーメント積分量を、以下に述べる方法を用いて時間更新する。基礎方程式（１．１）−（１．３）は、移流項（左辺第２項）と非移流項（右辺）と、を含む。まず移流項計算は、セミラグランジュ法と呼ばれる、上流点の値を速度場に沿って移流する方法を用いる（移流項計算ステップＳ１０４）。
【００６４】
【数１５】
【００６５】
ここで、左上添え字は時間ステップを表わし、アスタリスクは、該当する量が移流項計算後の中間値であることを示す。式（１．１３）、（１．１４）は、０次、１次、２次モーメント積分量の保存を保証する。次にグリッド上の物理量及びグリッド間のモーメント積分量の非移流項を、移流項計算後のグリッド上の物理量及びグリッド間のモーメント積分量を用いて、以下に述べる方法で算出する。
【００６６】
【数１６】
【００６７】
このようにして決定された非移流項を用いて、グリッド上の物理量及びグリッド間のモーメント積分量に非移流項を加算する（非移流項加算ステップＳ１０５）。
【００６８】
【数１７】
【００６９】
以上のステップを経て、ｎ時間ステップのグリッド上の物理量 nｆi及びグリッド間のモーメント積分量 nＭmi+1/2から、（ｎ＋１）時間ステップのグリッド上の物理量（１．１５）及びグリッド間のモーメント積分量（１．１３）、（１．１６）が得られる。
＜２次元ＭＭＡ法＞
２次元ＭＭＡ法は、物理量とその０次モーメント積分量、各次元方向の１次、２次モーメント積分量を用いて、２次元移流方程式を数値的に解く解法である。扱う基礎方程式は保存型の移流方程式で、次の通りである。
【００７０】
【数１８】
【００７１】
ここで、時刻ｔ、空間座標ｘ、ｙに対してｆが物理量、ｕ、ｖがｘ、ｙ方向の速度場である。
【００７２】
【数１９】
【００７３】
はそれぞれ、ｆの０次モーメントの空間積分量、ｘ方向のｍ次モーメントの空間積分量、ｙ方向のｍ次モーメントの空間積分量（モーメント積分量）である。今後、都合上、Ｍx0=Ｍy0=Ｍ0として表記することがある。式（２．２）−（２．４）は式（２．１）から厳密に導かれる。
【００７４】
コンピュータは、図５で示すように、グリッド上の物理量ｆij及びグリッド間のモーメント積分量
【００７５】
【数２０】
【００７６】
の、計６つの変数をＲＡＭ（３０２）に記憶する（初期値記憶ステップＳ１０１）。また、コンピュータは、グリッド上でｘ、ｙ方向の速度場ｕij、 ｖij及びその空間微分量（∂ｕ／∂ｘ）ij、（∂ｖ／∂ｙ）ijを別途求め、ＲＡＭ（３０２）に記憶する。
【００７７】
物理量ｆとモーメント積分量Ｍxm、Ｍymのグリッド間のプロファイルを、下記の区分化された多項式で近似する。
【００７８】
【数２１】
【００７９】
今後Ｇx0＝Ｇy0＝Ｇ0として表記することがある。
記憶量であるグリッド上の物理量とグリッド間のモーメント積分量を束縛条件として用いて、前記多項式の係数Ｃlkijを決定する。束縛条件は、
【００８０】
【数２２】
【００８１】
で与えられる。ここで、
【００８２】
【数２３】
【００８３】
は、（ｉ，ｊ）番目のグリッドから見た上流ｘ方向のグリッド及びセルの番号と上流ｙ方向のグリッド及びセルの番号であり、
【００８４】
【数２４】
【００８５】
は、（ｉ，ｊ）番目の座標（ｘi，ｙj）から見た上流点までのｘ方向及びｙ方向の距離であり、（ｘi＋ξij，ｙj＋ηij）が上流点の位置である（上流点算出ステップＳ１０２）。以上より、係数Ｃlkijが次の様に求まる。
【００８６】
【数２５】
【００８７】
ここで、Δｘ = ｘiup−ｘi、Δｙ = ｙjup−ｙjである。
このようにして、上記係数を要素とする物理量及びモーメント積分量のグリッド間の関数を空間プロファイル（２．８）−（２．１０）として決定する（空間プロファイル算出ステップＳ１０３）。続いて、決定された空間プロファイルを用いて、グリッド上の物理量及びグリッド間のモーメント積分量を、以下に述べる方法を用いて時間更新する。基礎方程式（２．１）−（２．４）は、移流項（左辺第２項、３項）と非移流項（右辺）と、を含む。まずグリッド上の物理量の時間更新は、セミラグランジュ法による移流項計算（移流項計算ステップＳ１０４）と、非移流項計算（非移流項加算ステップＳ１０５）と、をルンゲ・クッタ法を用いて行う。
【００８８】
【数２６】
【００８９】
ここで、ｋmaxはルンゲ・クッタ法の次数を表す。式（２．２２）、（２．２３）は特性曲線方程式であり、この式から、（ｉ，ｊ）番目の座標（ｘi，ｙj）から見た上流点のルンゲ・クッタ法のサブステップ毎の座標（ｘi＋ kξij，ｙj＋ kηij）が決定される。
【００９０】
kＳijはルンゲ・クッタ法のサブステップ毎のｆの非移流項であり、 kαはルンゲ・クッタ法の重み係数である。式（２．２４）がｆの時間更新式であり、右辺第１項がセミラグランジュ法による移流項計算（移流項計算ステップＳ１０４）、右辺第２項が非移流項の加算（非移流項加算ステップＳ１０５）である。この式から、ルンゲ・クッタ法のサブステップ毎のグリッド上の物理量 kｆijが求まり、ｋ＝ｋmaxのとき、kmaxｆij＝n+1ｆijであり、（ｎ＋１）時間ステップのグリッド上の物理量n+1ｆijが得られる。
【００９１】
次にグリッド間のモーメント積分量の時間更新方法について述べる。多次元移流問題では、グリッド間のモーメント積分量のセミラグランジュ法による移流計算は非常に煩雑になるため、有限体積法による移流項計算（移流項計算ステップＳ１０４）と、非移流項加算（非移流項加算ステップＳ１０５）と、を以下で行う。式（２．２）−（２．４）を次の保存形式で表現する。
【００９２】
【数２７】
【００９３】
ここで式（２．２５）−（２．２７）の右辺第１項から第４項が移流項であり、物理量の線積分量で表現されている。そこで先に決定された物理量の空間プロファイル（２．８）を用いて、線積分量を構築する。
【００９４】
【数２８】
【００９５】
これらを用いて式（２．２５）−（２．２７）の移流項の計算を行う（移流項計算ステップＳ１０４）。式（２．２５）−（２．２７）は保存型で書かれているため、０次、１次、２次モーメント積分量の保存が保証されている。続いて、ｘ、ｙ方向の１次及び２次モーメント積分量の非移流項である式（２．２６）、（２．２７）の右辺第５項を
【００９６】
【数２９】
【００９７】
により算出し、式（２．２６）、（２．２７）に加算する（非移流項加算ステップＳ１０５）。以上より、式（２．２５）−（２．２７）を有限体積法で計算することができ、結果、（ｎ＋１）時間ステップのグリッド間のモーメント積分量 n+1Ｍxm、 n+1Ｍymが得られる。但し計算を安定に行うために、式（２．２５）−（２．２７）は、ルンゲ・クッタ法を用いて時間更新する必要がある。
＜計算例＞
前記ＭＭＡ法による移流方程式の解法例を以下に示す。
【００９８】
図６上段はガウス分布形状の１次元線形移流計算結果であり、左よりＣＩＰ−ＣＳＬ２法、ＭＭＡ法の結果を表す。図中で、横軸が１次元空間で、「＋」は計算結果を、破線は解析解を表す。下段は解析解からの誤差を表す。ＣＩＰ−ＣＳＬ２法では、移流と共にプロファイルがなまってしまい、数値拡散の影響が大きく出ている。また、裾野のあたりでは数値振動が出ていることがわかる。一方、ＭＭＡ法はプロファイルを正確に再現しており、誤差の大きさもＣＩＰ−ＣＳＬ２法に比べて３桁以上小さい。
【００９９】
図７は、ガウス分布形状の２次元剛体回転問題の数値計算結果である。軸が２次元空間、背景及び等高線が物理量を表す。プロファイルを最上段のグラフで示されている矢印の方向に回転し続けた。左の列がＣＩＰ−ＣＳＬ２法の、右の列がＭＭＡ法の計算結果である。それぞれの使用メモリー量を合わせるために、グリッド数を変えている。上から順に、初期状態、３０回転後、３００回転後の結果を表す。２次元剛体回転問題では、数値拡散の影響が顕著に現れる。ＣＩＰ−ＣＳＬ２法による剛体回転問題結果では、３０回転後においてすでにプロファイルの数値拡散が見られ、３００回転後では原形をとどめていない。一方、ＭＭＡ法は、３００回転後においても初期の形状を保ったまま、剛体回転運動を数値的に記述することができている。
【０１００】
図８は、ＭＭＡ法のブラソフ方程式への適用例である。初期に互いに異なる向きの電子ビームを与えた。（ａ）から（ｄ）の順に、その位相空間分布の時間発展を示す。横軸が１次元実空間、縦軸が１次元速度空間、背景及び等高線が位相空間分布を表す。プラズマ物理分野でよく知られる２流体不安定が励起されており、その特徴的な構造である、位相空間上における電子ホールの形成が再現できている。
【０１０１】
図９は、ＭＭＡ法のブラソフ方程式への適用例である。イオン、電子が存在する中での、磁場に垂直方向に伝播するプラズマ波動である、バーンシュタイン波の分散関係である。横軸が波数、縦軸が周波数、背景が電場強度を表す。本計算は、プラズマの旋回運動（剛体回転）を長時間計算する必要があるため、従来のブラソフシミュレーション手法では再現することができなかった。左側がＭＭＡ法で求めたブラソフシミュレーション結果で、右側が粒子法の一種である電磁粒子（ＰＩＣ）コードで計算した結果である。両者共に同じ結果が得られ、ＭＭＡ法による計算が正しいことがわかる。また、ブラソフシミュレーションは電磁粒子計算に比べて数値ノイズが極めて少ないことが特徴である。
【産業上の利用可能性】
【０１０２】
移流方程式は物質の移動を記述する基本的な方程式であるため、その数値的解法は、数値流体力学（ＣＦＤ）や、宇宙・核融合分野におけるプラズマ現象への適用のみならず、大気流体現象、燃焼流体現象、コンピュータ画像処理など、幅広い分野に応用可能である。
【符号の説明】
【０１０３】
１…移流方程式解法装置、３００…バス、３０１…ＣＰＵ、３０２…ＲＡＭ、３０３…ＨＤＤ、３０４…マウス、３０５…キーボード、３０６…ディスプレイ
【特許請求の範囲】
【請求項１】
実空間と、その空間上で定義された物理量と、前記空間上で定義された移流速度場と、から前記物理量の移流を求める移流方程式を数値計算によって解く方法であって、
初期時刻ｔ0での空間上の各グリッドにおける物理量ｆ（ｘ、ｔ0）（ここで、ｘは空間上での位置ベクトル）と、その物理量ｆの複数のモーメント積分量Ｍm（ｘ、ｔ0）（ここで、ｍは０からｍｍａｘまでの整数で、モーメントの次数を表す。）と、最大時間ステップ数ｎｍａｘと、を記憶手段に記憶する初期値記憶ステップと、
空間のグリッド上で定義された移流速度ベクトルｕから、上流点及び上流方向のグリッドを決める上流点算出ステップと、
前記物理量ｆ（ｘ、ｔn）及び前記各モーメント積分量Ｍm（ｘ、ｔn）の束縛条件を課して、物理量ｆ（ｘ、ｔn）の空間プロファイルＦ（ｘ、ｔn）（ここで、ｔnは第ｎ時間ステップでの時刻を表し、ｎは０以上の整数である。）及び各モーメント量Ｍm（ｘ、ｔn）の空間プロファイルＧm（ｘ、ｔn）を求める空間プロファイル算出ステップと、
前記空間プロファイルに基づく次式
ｆ（ｘ、ｔ*）＝Ｆ（ｘ−ｕ×ｄｔ、ｔn）及び
Ｍm（ｘ、ｔ*）＝Ｇm（ｘ−ｕ×ｄｔ、ｔn）
に従って（ここで、ｄｔ＝ｔn+1−ｔnは時間刻み幅、t*は移流後の中間時間ステップである。）、中間時間ステップでの物理量ｆ（ｘ、ｔ*）及び各モーメント積分量Ｍm（ｘ、ｔ*）を算出する移流項計算ステップと、
前記物理量ｆ（ｘ、ｔ*）及び各モーメント積分量Ｍm（ｘ、ｔ*）に基づき、物理量ｆの非移流項Ｓ（ｘ、ｔ*）及び各モーメント積分量Ｍmの非移流項Ｓm（ｘ、ｔ*）を算出し、これら物理量ｆの非移流項Ｓ（ｘ、ｔ*）及び各モーメント積分量Ｍmの非移流項Ｓm（ｘ、ｔ*）を前記物理量ｆ（ｘ、ｔ*）及び各モーメント積分量Ｍm（ｘ、ｔ*）に加算して、次の時間ステップに対応する物理量ｆ（ｘ、ｔn+1）及び各モーメント積分量Ｍm（ｘ、ｔn+1）を算出する非移流項加算ステップと、
計算結果を出力する出力ステップと、
を含むことを特徴とする移流方程式を数値的に解く方法。
【請求項２】
位相空間と、その空間上で定義された物理量であるプラズマ分布関数ｆ（ｘ、ｖ、ｔ）（ここでベクトルｘは実空間での位置ベクトル、ベクトルｖは速度空間での位置ベクトル（速度場ベクトル））と、その実空間上で定義された物理量である電場ベクトルＥ（ｘ、ｔ）と、磁場ベクトルＢ（ｘ、ｔ）と、前記速度場ベクトルｖと、前記電場ベクトルＥ（ｘ、ｔ）と、磁場ベクトルＢ（ｘ、ｔ）と、から算出される、速度空間上での移流速度である加速度場ベクトルＡ（ｘ、ｖ、ｔ）と、から前記プラズマ分布関数の移流を求める移流方程式であるブラソフ方程式と、マクスウェル方程式と、を数値計算によって解く方法において、
初期時刻ｔ0におけるプラズマ分布関数ｆ（ｘ、ｖ、ｔ0）と、その速度空間における複数のモーメント積分量Ｍm（ｘ、ｖ、ｔ0）（ここで、ｍは０からｍｍａｘまでの整数で、モーメントの次数を表す。）と、電場ベクトルＥ（ｘ、ｔ0）と、磁場ベクトルＢ（ｘ、ｔ0）と、最大時間ステップ数ｎｍａｘと、を記憶手段に記憶する初期値記憶ステップと、
プラズマ分布関数ｆと、速度空間モーメント積分量Ｍmと、の実空間上における時間更新を行い、第一中間時間ステップにおけるプラズマ分布関数ｆ（ｘ、ｖ、ｔ*）と、速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ*）と、を算出する実空間時間更新ステップと、
０次の速度空間モーメント積分量Ｍ0（ｘ、ｖ、ｔ*）と、１次の速度空間モーメント積分量Ｍ1（ｘ、ｖ、ｔ*）と、から、電荷密度と、電流密度と、を求める電荷・電流密度算出ステップと、
前記電荷・電流密度算出ステップで算出した、電荷密度と、電流密度と、からマクスウェル方程式を解き、前記電場ベクトルＥ（ｘ、ｔn）と、磁場ベクトルＢ（ｘ、ｔn）と、から次の時間ステップに対応する電場ベクトルＥ（ｘ、ｔn+1）と、磁場ベクトルＢ（ｘ、ｔn+1）と、を算出する電磁場更新ステップと、
前記時間更新後の電場ベクトルＥ（ｘ、ｔn+1）と、磁場ベクトルＢ（ｘ、ｔn+1）と、から前記加速度場ベクトルＡ（ｘ、ｖ、ｔn+1）を算出する加速度場算出ステップと、
速度空間上における移流速度である前記加速度ベクトルＡ（ｘ、ｖ、ｔn+1）から、速度空間における上流点及び上流方向のグリッドを算出する速度空間上流点算出ステップと、
前記プラズマ分布関数ｆ（ｘ、ｖ、ｔ*）及び前記各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ*）の束縛条件を課して、プラズマ分布関数ｆ（ｘ、ｖ、ｔ*）の速度空間プロファイルＦ（ｘ、ｖ、ｔ*）及び各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ*）の速度空間プロファイルＧm（ｘ、ｖ、ｔ*）を求める速度空間プロファイル算出ステップと、
前記空間プロファイルに基づく次式
ｆ（ｘ、ｖ、ｔ**）＝Ｆ（ｘ、ｖ−Ａ×ｄｔ、ｔ*）及び
Ｍm（ｘ、ｖ、ｔ**）＝Ｇm（ｘ、ｖ−Ａ×ｄｔ、ｔ*）
に従って、第二中間時間ステップでのプラズマ分布関数ｆ（ｘ、ｖ、ｔ**）及び各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ**）を算出する速度空間移流項計算ステップと、
前記プラズマ分布関数ｆ（ｘ、ｖ、ｔ**）及び各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ**）に基づき、プラズマ分布関数ｆの非移流項Ｓ（ｘ、ｖ、ｔ**）及び各速度空間モーメント積分量Ｍmの非移流項Ｓm（ｘ、ｖ、ｔ**）を算出し、これらプラズマ分布関数ｆの非移流項Ｓ（ｘ、ｖ、ｔ**）及び各速度空間モーメント積分量Ｍmの非移流項Ｓm（ｘ、ｖ、ｔ**）を前記プラズマ分布関数ｆ（ｘ、ｖ、ｔ**）及び各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ**）に加算して、次の時間ステップに対応するプラズマ分布関数ｆ（ｘ、ｖ、ｔn+1）及び各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔn+1）を算出する速度空間非移流項加算ステップと、
計算結果を出力する出力ステップと、
を含むことを特徴とする移流方程式を解く方法。
【請求項３】
前記電磁場更新ステップでは、実空間上での時間更新後の前記０次の速度空間モーメント積分量Ｍ0（ｘ、ｖ、ｔ*）と、１次の速度空間モーメント積分量Ｍ1（ｘ、ｖ、ｔ*）と、から、電荷密度と、電流密度と、を求めることを特徴とする請求項２記載の移流方程式及びマクスウェル方程式を数値的に解く方法。
【請求項４】
コンピュータに、請求項１に記載の移流方程式を数値的に解く方法に含まれる各ステップを実行させるためのプログラム。
【請求項５】
コンピュータに、請求項２又は請求項３に記載の移流方程式及びマクスウェル方程式を数値的に解く方法に含まれる各ステップを実行させるためのプログラム。
【請求項１】
実空間と、その空間上で定義された物理量と、前記空間上で定義された移流速度場と、から前記物理量の移流を求める移流方程式を数値計算によって解く方法であって、
初期時刻ｔ0での空間上の各グリッドにおける物理量ｆ（ｘ、ｔ0）（ここで、ｘは空間上での位置ベクトル）と、その物理量ｆの複数のモーメント積分量Ｍm（ｘ、ｔ0）（ここで、ｍは０からｍｍａｘまでの整数で、モーメントの次数を表す。）と、最大時間ステップ数ｎｍａｘと、を記憶手段に記憶する初期値記憶ステップと、
空間のグリッド上で定義された移流速度ベクトルｕから、上流点及び上流方向のグリッドを決める上流点算出ステップと、
前記物理量ｆ（ｘ、ｔn）及び前記各モーメント積分量Ｍm（ｘ、ｔn）の束縛条件を課して、物理量ｆ（ｘ、ｔn）の空間プロファイルＦ（ｘ、ｔn）（ここで、ｔnは第ｎ時間ステップでの時刻を表し、ｎは０以上の整数である。）及び各モーメント量Ｍm（ｘ、ｔn）の空間プロファイルＧm（ｘ、ｔn）を求める空間プロファイル算出ステップと、
前記空間プロファイルに基づく次式
ｆ（ｘ、ｔ*）＝Ｆ（ｘ−ｕ×ｄｔ、ｔn）及び
Ｍm（ｘ、ｔ*）＝Ｇm（ｘ−ｕ×ｄｔ、ｔn）
に従って（ここで、ｄｔ＝ｔn+1−ｔnは時間刻み幅、t*は移流後の中間時間ステップである。）、中間時間ステップでの物理量ｆ（ｘ、ｔ*）及び各モーメント積分量Ｍm（ｘ、ｔ*）を算出する移流項計算ステップと、
前記物理量ｆ（ｘ、ｔ*）及び各モーメント積分量Ｍm（ｘ、ｔ*）に基づき、物理量ｆの非移流項Ｓ（ｘ、ｔ*）及び各モーメント積分量Ｍmの非移流項Ｓm（ｘ、ｔ*）を算出し、これら物理量ｆの非移流項Ｓ（ｘ、ｔ*）及び各モーメント積分量Ｍmの非移流項Ｓm（ｘ、ｔ*）を前記物理量ｆ（ｘ、ｔ*）及び各モーメント積分量Ｍm（ｘ、ｔ*）に加算して、次の時間ステップに対応する物理量ｆ（ｘ、ｔn+1）及び各モーメント積分量Ｍm（ｘ、ｔn+1）を算出する非移流項加算ステップと、
計算結果を出力する出力ステップと、
を含むことを特徴とする移流方程式を数値的に解く方法。
【請求項２】
位相空間と、その空間上で定義された物理量であるプラズマ分布関数ｆ（ｘ、ｖ、ｔ）（ここでベクトルｘは実空間での位置ベクトル、ベクトルｖは速度空間での位置ベクトル（速度場ベクトル））と、その実空間上で定義された物理量である電場ベクトルＥ（ｘ、ｔ）と、磁場ベクトルＢ（ｘ、ｔ）と、前記速度場ベクトルｖと、前記電場ベクトルＥ（ｘ、ｔ）と、磁場ベクトルＢ（ｘ、ｔ）と、から算出される、速度空間上での移流速度である加速度場ベクトルＡ（ｘ、ｖ、ｔ）と、から前記プラズマ分布関数の移流を求める移流方程式であるブラソフ方程式と、マクスウェル方程式と、を数値計算によって解く方法において、
初期時刻ｔ0におけるプラズマ分布関数ｆ（ｘ、ｖ、ｔ0）と、その速度空間における複数のモーメント積分量Ｍm（ｘ、ｖ、ｔ0）（ここで、ｍは０からｍｍａｘまでの整数で、モーメントの次数を表す。）と、電場ベクトルＥ（ｘ、ｔ0）と、磁場ベクトルＢ（ｘ、ｔ0）と、最大時間ステップ数ｎｍａｘと、を記憶手段に記憶する初期値記憶ステップと、
プラズマ分布関数ｆと、速度空間モーメント積分量Ｍmと、の実空間上における時間更新を行い、第一中間時間ステップにおけるプラズマ分布関数ｆ（ｘ、ｖ、ｔ*）と、速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ*）と、を算出する実空間時間更新ステップと、
０次の速度空間モーメント積分量Ｍ0（ｘ、ｖ、ｔ*）と、１次の速度空間モーメント積分量Ｍ1（ｘ、ｖ、ｔ*）と、から、電荷密度と、電流密度と、を求める電荷・電流密度算出ステップと、
前記電荷・電流密度算出ステップで算出した、電荷密度と、電流密度と、からマクスウェル方程式を解き、前記電場ベクトルＥ（ｘ、ｔn）と、磁場ベクトルＢ（ｘ、ｔn）と、から次の時間ステップに対応する電場ベクトルＥ（ｘ、ｔn+1）と、磁場ベクトルＢ（ｘ、ｔn+1）と、を算出する電磁場更新ステップと、
前記時間更新後の電場ベクトルＥ（ｘ、ｔn+1）と、磁場ベクトルＢ（ｘ、ｔn+1）と、から前記加速度場ベクトルＡ（ｘ、ｖ、ｔn+1）を算出する加速度場算出ステップと、
速度空間上における移流速度である前記加速度ベクトルＡ（ｘ、ｖ、ｔn+1）から、速度空間における上流点及び上流方向のグリッドを算出する速度空間上流点算出ステップと、
前記プラズマ分布関数ｆ（ｘ、ｖ、ｔ*）及び前記各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ*）の束縛条件を課して、プラズマ分布関数ｆ（ｘ、ｖ、ｔ*）の速度空間プロファイルＦ（ｘ、ｖ、ｔ*）及び各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ*）の速度空間プロファイルＧm（ｘ、ｖ、ｔ*）を求める速度空間プロファイル算出ステップと、
前記空間プロファイルに基づく次式
ｆ（ｘ、ｖ、ｔ**）＝Ｆ（ｘ、ｖ−Ａ×ｄｔ、ｔ*）及び
Ｍm（ｘ、ｖ、ｔ**）＝Ｇm（ｘ、ｖ−Ａ×ｄｔ、ｔ*）
に従って、第二中間時間ステップでのプラズマ分布関数ｆ（ｘ、ｖ、ｔ**）及び各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ**）を算出する速度空間移流項計算ステップと、
前記プラズマ分布関数ｆ（ｘ、ｖ、ｔ**）及び各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ**）に基づき、プラズマ分布関数ｆの非移流項Ｓ（ｘ、ｖ、ｔ**）及び各速度空間モーメント積分量Ｍmの非移流項Ｓm（ｘ、ｖ、ｔ**）を算出し、これらプラズマ分布関数ｆの非移流項Ｓ（ｘ、ｖ、ｔ**）及び各速度空間モーメント積分量Ｍmの非移流項Ｓm（ｘ、ｖ、ｔ**）を前記プラズマ分布関数ｆ（ｘ、ｖ、ｔ**）及び各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔ**）に加算して、次の時間ステップに対応するプラズマ分布関数ｆ（ｘ、ｖ、ｔn+1）及び各速度空間モーメント積分量Ｍm（ｘ、ｖ、ｔn+1）を算出する速度空間非移流項加算ステップと、
計算結果を出力する出力ステップと、
を含むことを特徴とする移流方程式を解く方法。
【請求項３】
前記電磁場更新ステップでは、実空間上での時間更新後の前記０次の速度空間モーメント積分量Ｍ0（ｘ、ｖ、ｔ*）と、１次の速度空間モーメント積分量Ｍ1（ｘ、ｖ、ｔ*）と、から、電荷密度と、電流密度と、を求めることを特徴とする請求項２記載の移流方程式及びマクスウェル方程式を数値的に解く方法。
【請求項４】
コンピュータに、請求項１に記載の移流方程式を数値的に解く方法に含まれる各ステップを実行させるためのプログラム。
【請求項５】
コンピュータに、請求項２又は請求項３に記載の移流方程式及びマクスウェル方程式を数値的に解く方法に含まれる各ステップを実行させるためのプログラム。
【図１】
【図２】
【図３】
【図４】
【図１０】
【図１１】
【図５】
【図６】
【図７】
【図８】
【図９】
【図２】
【図３】
【図４】
【図１０】
【図１１】
【図５】
【図６】
【図７】
【図８】
【図９】
【公開番号】特開２０１１−１５９０３２（Ｐ２０１１−１５９０３２Ａ）
【公開日】平成２３年８月１８日（２０１１．８．１８）
【国際特許分類】
【出願番号】特願２０１０−１９０２０（Ｐ２０１０−１９０２０）
【出願日】平成２２年１月２９日（２０１０．１．２９）
【出願人】（５０４１３９６６２）国立大学法人名古屋大学 (996)
【Ｆターム（参考）】
【公開日】平成２３年８月１８日（２０１１．８．１８）
【国際特許分類】
【出願日】平成２２年１月２９日（２０１０．１．２９）
【出願人】（５０４１３９６６２）国立大学法人名古屋大学 (996)
【Ｆターム（参考）】
[ Back to top ]