(57)【特許請求の範囲】
【請求項１】
入力手段と演算手段と出力手段とを備え，プラズマを解析する装置において，
イオンを粒子とし，電子を流体とみなすハイブリッドコードにより，個々のイオンの運動方程式と電子流体の運動方程式およびマックスウェル方程式とによりプラズマを解析するプラズマ解析手段を保持し，プラズマ塊と，プラズマ塊にイオンと電子流体が照射される系のプラズマ解析をするものであって，
入力されたプラズマの初期条件，境界条件および化学反応条件を保持し，
前記化学反応条件は，与えられた反応率で個々のイオンを追加する，あるいは，別種のイオンに置き換えることにより実現するものであって，
前記演算手段は，プラズマ解析手段およびプラズマの初期条件，境界条件および化学反応条件の変化に対して，プラズマ塊に照射されるイオンとプラズマ塊の外部周辺にあるイオンとの化学反応により生じるイオン，およびこれと相互作用する電子流体の，密度と速度，さらに電磁場の時間変化や空間分布を解析することを特徴とするプラズマ解析装置。
【請求項２】
プラズマ塊の境界から一定距離の範囲だけ外の領域において，前記化学反応を生じるものとして，プラズマ塊の外部および内部の物理現象を解析することを特徴とする請求項１に記載のプラズマ解析装置。
【請求項３】
前記化学反応はイオン化されていない粒子に光が照射されたことにより荷電粒子を生じる光電離反応および／または異なる種類の粒子の間で電荷の交換をともなう荷電交換反応であることを特徴とする請求項１又は２に記載のプラズマ解析装置。
【請求項４】
太陽風を外部からプラズマ塊に入射するイオンと電子流体であるとし，太陽風とプラズマ塊もしくはプラスマ塊の外部周辺でのイオンの化学反応を条件として，宇宙プラズマを解析するものであることを特徴とする請求項１，２又は３に記載のプラズマ解析装置。
【請求項５】
前記化学反応条件における光電離反応はＯ＋ｈν→Ｏ+ ＋ｅ- および荷電交換はＨ+ ＋Ｏ→Ｈ＋Ｏ+ であることを特徴とする請求項３又は４に記載のプラスマ解析装置。
【請求項６】
演算手段によりプラズマを解析する方法において，
イオンを粒子とし，電子を流体とみなすハイブリッドコードにより，個々のイオンの運動方程式と電子流体の運動方程式およびマックスウェル方程式とによりプラズマを解析するプラズマ解析手段を備え，プラズマ塊と，プラズマ塊にイオンと電子流体が照射される系のプラズマ解析をするものであって，
入力されたプラズマの初期条件，境界条件および化学反応条件を保持し，
前記化学反応条件は，与えられた反応率で個々のイオンを追加する，あるいは，別種のイオンに置き換えることにより実現するものであって，
前記プラズマ解析手段およびプラズマの初期条件，境界条件および化学反応条件の変化に対して，プラズマ塊に照射されるイオンとプラズマ塊の外部周辺にあるイオンとの化学反応により生じるイオン，およびこれと相互作用する電子流体の，密度と速度，さらに電磁場の時間変化や空間分布を解析することを特徴とするプラズマ解析方法。
【請求項７】
プラズマ塊の境界から一定距離の範囲だけ外の領域において，前記化学反応を生じるものとして，プラズマ塊の外部および内部の物理現象を解析することを特徴とする請求項６に記載のプラズマ解析方法。
【請求項８】
化学反応はイオン化されていない粒子に光が照射されたことにより荷電粒子を生じる光電離反応および／または異なる種類の粒子の間の電荷の交換をともなう荷電交換反応であることを特徴とする請求項６又は７に記載のプラスマ解析方法。
【請求項９】
太陽風を外部からプラズマ塊に入射するイオンと電子流体であるとし，太陽風とプラズマ塊もしくはプラスマ塊の外部周辺でのイオンの化学反応を条件として，宇宙プラズマを解析するものであることを特徴とする請求項６，７又は８に記載のプラズマ解析方法。
【請求項１０】
前記化学反応条件における光電離反応はＯ＋ｈν→Ｏ+ ＋ｅ- および荷電交換はＨ+ ＋Ｏ→Ｈ＋Ｏ+ であることを特徴とする請求項８又は９に記載のプラスマ解析方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
プラズマ解析装置および解析方法に関するものである。特に，宇宙空間のプラズマをシミュレーションモデルにして記述した運動方程式と電磁界のマックスウェル方程式からプラズマ密度や速度の分布と電磁界とを，イオンの化学反応が生じる系においても三次元でシミュレーションできるものである。
【０００２】
本発明は，宇宙開発関連産業，プラズマ工学分野等で利用できるものである。
【０００３】
【従来の技術】
太陽風と惑星のプラズマの相互作用等のプラズマ現象の計算機シミュレーションにおいて，イオンを粒子，電子を流体として扱うハイブリットコードが使われてきた。惑星を非磁性体の固体球とみなすモデルを使用し，太陽風と惑星とのプラズマ相互作用をシミュレーションすることができる。ハイブリットコードを宇宙プラズマに適用した例は次の論文に説明されている。
【０００４】
（１）Ｌｅｒｏｙ，Ｍ．Ｍ．，Ｄ．Ｗｉｎｓｋｅ，Ｃ．Ｃ．Ｇｏｏｄｒｉｃｈ，Ｃ．Ｓ．Ｗｕ，ａｎｄ  Ｋ．Ｐａｐａｄｏｐｏｕｌｏｓ，”Ｔｈｅ  ｓｔｒｕｃｔｕｒｅ  ｏｆ  ｐｅｒｐｅｎｄｉｃｕｌａｒ  ｂｏｗ  ｓｈｏｃｋｓ”，Ｊｏｕｒｎａｌ  ｏｆ  Ｇｅｏｐｈｙｓｉｃａｌ  Ｒｅｓｅａｒｃｈ，  ｖｏｌ．８７，ｐｐ．５０８１−５０９４（１９８２）。
【０００５】
図１２（ａ）は惑星を非磁性体の固体とみなし，太陽風と惑星とのプラズマ相互作用を解析するモデルを示す。
【０００６】
惑星１１０を非磁性体の固体球とみなし，プラズマ流である太陽風が照射されることを示している。惑星に照射する太陽風のプラズマイオンは惑星表面で反射し，一部イオンは惑星内部に侵入し，内部で消滅するモデルである。
【０００７】
このモデルについて，ハイブリッドコードによりシミュレーションした結果は次の論文等で詳細に説明されている。
【０００８】
（２）Ｓｈｉｍａｚｕ，Ｈ．，”Ｔｈｒｅｅ−ｄｉｍｅｎｓｉｏｎａｌ  ｈｙｂｒｉｄ  ｓｉｍｕｌａｔｉｏｎ  ｏｆ  ｍａｇｎｅｔｉｚｅｄ  ｐｌａｓｍａｆｌｏｗ  ａｒｏｕｎｄ  ａｎ  ｏｂｓｔａｃｌｅ”，Ｅａｒｔｈ，Ｐｌａｎｅｔｓ  ａｎｄ  Ｓｐａｃｅ，  ｖｏｌ．５１，ｐｐ３８３−３９３（１９９９）。
【０００９】
【発明が解決しようとする課題】
先行技術の説明
惑星の上層に電離層が存在する場合には，図１２（ａ）のように惑星を固体と見做して，太陽風との相互作用をシミュレーションすることは実際に合わないものである。実際の惑星の上層には電離層があり，惑星をプラズマ塊と見做すことができる。そこで，本出願の発明者は，惑星をプラズマ塊としてシミュレーションを試みた（「太陽風と火星金星型惑星との相互作用の３次元ハイブリッドシミュレーション」地球惑星圏学会第１０６回総会・講演会，仙台市民会館，１９９９年１１月９日−１２日，講演番号Ｂ１２−Ｐ０５４で発表）。
【００１０】
図１２（ｂ）は惑星を半径Ｒの球のプラズマ塊として扱い，太陽風がプラズマ塊に照射するモデルを示す。
【００１１】
図１２（ｂ）において，１１０は惑星であり，プラズマ塊である。惑星に照射する太陽風のイオン（Ｈ+ ）は惑星内部に入り，惑星内のイオン（Ｏ+ ）と混合することもある。
【００１２】
このモデルについて，ハイブリットコードを使用して，イオンの運動方程式とマックス方程式を解いて，プラズマ密度，速度分布および電磁界を計算機シミュレーションした結果を模式的に表したものを図１３に示す。
【００１３】
図１３において，１１０は惑星であり，プラズマ塊である。１１１は衝撃波（バウショックであり，以後ショックと称する）。
図１３に示されるように太陽風は超音速のプラズマ流であり，このプラズマ流が惑星の磁気圏に入る時に前面にショック１１１が生じる。
【００１４】
図１４（ａ）は，プラズマ塊モデルにより解析する時のパラメータを示す。図１４図（ｂ）は初期条件および境界条件を与えるモデルを示す。
図１４（ｂ）において，磁束密度は次式により与えられる。
【００１５】
【数１】
0003735704.tif 000002
【００１６】
シミュレーションは宇宙空間を６４×６４×６４に分割された空間（格子空間）とし，その内部の粒子密度，粒子速度等の平均値を格子点におけるそれぞれの値とする。そして，格子点での値をデータとして差分化した運動方程式，マックスウェル方程式等を使用して宇宙プラズマ現象をシミュレーションした。
【００１７】
シミュレーションにおけるアルベン速度ＶA と光速度ｃの比をＶA ／ｃ＝ωci／ωpi＝１．０×１０-4とした。ωciはイオンのサイクロトロン角周波数である。ωpiはイオンプラズマ振動の角周波数である。
【００１８】
熱圧力と磁気圧力の比は電子とイオンについて，βe ＝βi ＝１とした。
太陽風の速度はｖSW＝4 ．０ＶA である。
【００１９】
電離層に包まれた惑星を半径Ｒのプラズマ球として表現する。太陽風はｘ軸方向に左から右に流れ，磁束密度はｙ方向である。初期状態において，太陽風の動圧と惑星プラズマの圧力はほぼつりあっている。Ｌはシミュレーションする宇宙空間の大きさを示す。プラズマ球へのイオンの供給率は０．２５６ωpiである。
【００２０】
図１５（ａ），（ｂ），（ｃ）は，惑星をプラスマ球としたモデルにおけるシミュレーション結果を示す。
【００２１】
図１５（ａ）はｙ＝Ｌ／２での断面図であり，シミュレーション開始後，ωciｔ＝３７．５の時点でのイオン速度のｘ成分を示す。イオン速度はＶｉｘ／ＶA で表している。図１３（ａ）において，惑星中心は座標（１０２．４，１０２．４）にあり，太陽風はｘ軸方向（左から右）に吹いている。図示のようにショックが形成される。惑星半径は，Ｒ＝６．４ｒL である（ｒL ＝ｖSW／ωci）である。
【００２２】
図１５（ｂ）はｙ＝Ｌ／２での断面図における磁束密度の大きさを示す。磁束密度は｜Ｂ｜／Ｂｏで表している（Ｂｏは初期条件で与えた磁束密度の大きさである）。
図１５（ｃ）はｙ＝Ｌ／２での断面図における惑星の酸素イオンの分布を示す。図１６（ａ）は  ｘ＝Ｌ／２での断面図であり，シミュレーション開始後，ωciｔ＝３７．５の時点でのシミュレーション結果での磁束密度の大きさを示す。  図１６（ａ）において，惑星は実線の円で示され，ショック面は惑星を囲むぼんやりした円で示される。
図１６（ａ）のように，シミュレーションではショックの惑星からの高さはｚ軸の正方向で大きく，ｚ軸の負方向で小さい（ショック面は惑星に対してｚ軸の正方向にずれている。
【００２３】
これは，実際の観測結果と異なり，実際の観測結果はショックの惑星からの高さは，図１６（ｃ）に示すようにｚ軸の正方向で小さく，ｚ軸の負方向で大きい（ショック面は惑星に対してｚ軸の負方向にずれている）ことが観測されている。
【００２４】
図１６（ｂ）は  ｘ＝３Ｌ／４での断面図であり，シミュレーション開始後，ωciｔ＝３７．５の時点での磁束密度の大きさを示す。
【００２５】
上記の太陽風と惑星プラズマの相互作用のモデルでは，惑星の高層で生じるイオン反応（光電離や荷電交換等）が扱われていなかった。そのためショック面と惑星との位置関係のシミュレーション結果が必ずしも正しいものではなかった。
【００２６】
本発明は，イオンの化学反応が生じるプラズマ系においても正確に三次元でシミュレーションできるプラズマ解析装置および方法を提供することを目的とする。
【００２７】
【課題を解決するための手段】
図１は，本発明の基本構成を示す。
【００２８】
図１において，１は初期条件保持手段であって，入力手段６により入力された初期条件を入力するものである。２はプラズマの境界条件保持手段であって，入力手段６により入力されたプラズマ境界条件を保持するものである。３はイオン化学反応条件保持手段であって，入力手段６により入力されたプラズマの化学反応についての条件を保持するものである。４はプラズマ解析手段であって，粒子の運動方程式，マックスウェル方程式，初期条件，境界条件，イオン化学反応条件およびプラズマ状態を表すデータ（イオン生成率等のデータ）によりプラズマ解析をするものである。５は解析結果出力手段であって，解析結果を出力するものである。６は入力手段である。７はプラズマの状態データ保持手段である。
【００２９】
図１の本発明の基本構成の動作を説明する。
【００３０】
入力手段６により，初期条件，プラズマの境界条件，プラズマに生じる化学反応の条件を入力する。また，プラズマのイオンサイクロトロン周波数，イオン生成率等のプラズマ状態データを入力し，状態データ保持手段７に保持する。プラズマ解析手段４は，初期条件，境界条件，イオン化学反応条件およびプラズマの状態データによりプラズマに生じる状態を解析する。出力手段５は求められたプラズマ状態を出力する。
【００３１】
化学反応条件は，例えば，太陽風と惑星の相互作用の場合，太陽風のイオン（Ｈ+ ）と惑星イオン（Ｏ+ ）とした時、
Ｈ+ ＋Ｏ→Ｈ＋Ｏ+ （惑星プラズマ）
とする等である。Ｏ，Ｈはイオン化していない原子なので，解析においてプラズマから除かれる。このような化学反応条件の場合，イオン反応のため，電流変化となってプラズマの運動などに影響するものである。後述するが，太陽風と惑星イオンとの相互作用のシミュレーションにおいて，本発明のシミュレーション装置および方法によれば，惑星とショックの高さの関係は，前述の本発明に先行する方法では，得られなかったのであるが，図８に示すように実際の観測結果と一致するようになった。
【００３２】
本発明によれば，プラズマの初期条件，境界条件だけでなく，イオン反応条件もシミュレーションの条件にしてプラズマを解析できるので，イオンの化学反応が生じるプラズマ系においても正確に三次元でシミュレーションできるようになる。
【００３３】
なお，本発明のプラズマ解析装置および方法は，宇宙プラズマに適用するだけでなく，実験室等のプラズマにも適用できるものである。
【００３４】
【発明の実施の形態】
本発明を実施するためのハイブリッドコードについて説明する。
以下において，Ｂ，Ｅ，ｖ，ｊは全て３次元ベクトルであり，それぞれベクトルＢ，ベクトルＥ，ベクトルｖ，ベクトルｊ等で表す。ベクトルｘは位置ベクトルである。Ｂは磁束密度，Ｅは電界，ｖはイオン速度，ｊは電流をそれぞれ表す。
【００３５】
ハイブリッドコードにより解析する場合の基礎方程式は次のとおりである。
【００３６】
イオンの運動方程式
【００３７】
【数２】
0003735704.tif 000003
【００３８】
イオンの速度の式
【００３９】
【数３】
0003735704.tif 000004
【００４０】
電子流体（質量なし）の運動方程式
【００４１】
【数４】
0003735704.tif 000005
【００４２】
ただし，電子温度Ｔe は一定とする。
【００４３】
マックスウェル方程式
【００４４】
【数５】
0003735704.tif 000006
【００４５】
【数６】
0003735704.tif 000007
【００４６】
準中性条件
ｎe ＝ｎi ＝ｎ                                （６）
本実施の形態では電気伝導度は無限大と仮定する（η＝０）。
【００４７】
シミュレーションでは上記の式を差分化して使用する。
【００４８】
上記のマックスウェルの方程式（４）を差分化すると次のようになる。
【００４９】
【数７】
0003735704.tif 000008
【００５０】
ベクトルＢk+1 から式（５）を用いてベクトルｊk+1 を得る。
【００５１】
式（１），（２）を差分化すると次のようになる。
【００５２】
【数８】
0003735704.tif 000009
【００５３】
【数９】
0003735704.tif 000010
【００５４】
時刻ｋの磁束密度は時刻ｋ＋１の磁束密度と時刻ｋ−１の磁束密度とにより次のように表せる。
【００５５】
【数１０】
0003735704.tif 000011
【００５６】
時刻ｋ＋１の位置ベクトルは次のように表せる。
【００５７】
【数１１】
0003735704.tif 000012
【００５８】
イオン速度と密度を格子点毎に足しあわせることにより格子点の速度ベクトルｖi と密度ｎi が得られる。電子流体の速度ベクトルｖe は電流を表す次の式で求まる。
【００５９】
【数１２】
0003735704.tif 000013
【００６０】
電界ベクトルＥは電子速度ベクトルｖe と電子流体の運動方程式と準中性条件の式とにより求めることができる。
【００６１】
格子点におけるイオンの密度，速度は，格子点から近いものから大きく，遠いものは寄与が小さいようにする。
【００６２】
具体的には，図４のような形状関数を使用する。
【００６３】
密度は
【００６４】
【数１３】
0003735704.tif 000014
【００６５】
となる。
【００６６】
速度平均は
【００６７】
【数１４】
0003735704.tif 000015
【００６８】
で求まる。Ｎはシミュレーション対象の宇宙空間の全粒子数である。
【００６９】
図２は本発明の装置構成の実施の形態を示す。
【００７０】
図２は，イオンを粒子とし電子を流体として扱うハイブリットコードにより解析するための構成を示す。
【００７１】
１１は入力する初期条件を表す（磁束密度ベクトルＢ0 ，Ｂ1 ，電界ベクトルＥ0 ，粒子密度ｎ0 ，粒子速度ベクトルｖi0）。１２は入力手段である。１３は演算手段を表す。１４は演算結果の時刻ｋ＋１における電界ベクトルＥk+1 ，磁束密度ベクトルＢk+1 ，電流密度ベクトルＪk+1 ，粒子密度ｎk+1 ，粒子の位置ベクトルｘk+1 ，粒子の速度ベクトルｖk+1 である。１５は出力手段である。２０はデータ保持手段である。２１は初期条件保持部である。２２は境界条件保持部である。１７は状態データであって，イオンの生成率等のパラメータを表す。３１はハイブリッドコードにより解析する場合の演算式を保持する演算式保持部であり，上記の差分化方程式，その他演算に必要な式を保持するものである。３２は差分化イオン運動方程式である。３３は差分化電子運動方程式である。３４は差分化マックスウェル方程式である。３５は電子流体の速度の式である。３６は粒子中性条件であって，ｎi ＝ｎe ＝ｎを保持するものである。３６は粒子密度の平均，粒子速度の平均値を求める式である。３７はその他演算に必要な式である。
【００７２】
本発明のプラズマ解析装置を使用して，以下のモデルにより太陽風と惑星（火星，金星等）の相互作用をシミュレーションする場合について説明する。プラズマの生成率，プラズマ球の半径，初期の磁束密度は先行シミュレーション（前述のプラズマ塊モデル）の場合と同じである。
【００７３】
（１）  半径Ｒの惑星を表すプラズマ球の外側（１．０Ｒと１．３Ｒの昼側）で，光電離があるとする。
【００７４】
Ｏ＋ｈν→Ｏ+ ＋ｅ- 
（２）  同じ領域に荷電交換があるとする。
【００７５】
Ｈ+ ＋Ｏ→Ｈ＋Ｏ+ 
中性化した分子は計算に無関係になるので，（１）の場合はＯ+ がある割合で生成することになる。太陽風のＨ+ は変化させないものとする。生成する位置は，上記の領域内で一様に分布するように乱数を用いて決める。（２）の場合は，ある割合で太陽風が消え，その分Ｏ+ が生成することに等価である。即ち，１．０Ｒと１．３Ｒの間に達した太陽風のイオン（Ｈ+ ）のうち，あらかじめ決めておいた割合を消滅させる。消滅した位置にＯ+ を生成させる。この結果は粒子の平均化の処理（後述する）によりベクトルｖi とｎの変化にあらわれ，全体の計算に影響する。（１），（２）において，イオン速度は，温度が惑星イオンの温度と同じ速度分布になるように乱数で決める。速度分布は次の式（１３）で表されるマクスウェル分布となるようにする。
【００７６】
【数１５】
0003735704.tif 000016
【００７７】
Ａは次になるように定める。
【００７８】
【数１６】
0003735704.tif 000017
【００７９】
（１／２）ｍｖth2 ＝ｋB Ｔである。ｎは粒子密度，ｍはイオン質量，ｋB はボルツマン定数，Ｔは温度である。
【００８０】
図３は図２の構成の動作のフローチャートである。
図３の動作を説明する。
Ｓ１で磁束密度，電界，粒子速度，密度の初期条件ベクトルＢ0 ，ベクトル  Ｂ1 ，ベクトルＥ1 ，ベクトルｖi0，ｎを入力する。Ｂ0 の初期条件は，図１４（ｂ）と同じである。
【００８１】
初期条件は，太陽風および惑星のイオン速度がマクスウェル分布になるように与える。太陽風温度と惑星プラズマの温度は同じにする。
【００８２】
Ｓ２で時刻ｋの電界ベクトルＥk とＢk-1 とによりマクスウェルの方程式に従ってＢk+1 を求める。次にＳ３でマクスウェル方程式により時刻ｋ＋１の磁束密度ベクトルＢk+1 を使って，時刻ｋ＋１の電流密度ベクトルＪk+1 を求める。次にＳ４でイオンの粒子密度ｎi ，電子密度ｎe （＝ｎ）と速度ｖi と時刻ｋ＋１の電流密度ベクトルＪk+1 を使って, 時刻ｋ＋１の電子の速度ベクトルｖe を求める。
【００８３】
次に，Ｓ５で，時刻ｋ＋１の電子の速度ベクトルｖe ，磁束密度ベクトルＢk+1 ，電流密度ベクトルＪk+1 ，荷電密度ｎe ，ボルツマン定数ｋB ，電子温度Ｔe とにより，時刻ｋ＋１の電界ベクトルＥk+1 を求める。次に時刻をｋ＝ｋ＋２により更新して，Ｓ２以降の処理を繰り返す。
【００８４】
上記のメインルーチンに対して，境界条件，イオン反応条件，粒子密度の平均化，粒子速度の平均化はサブルーチン（Ｓ１１，Ｓ１２，Ｓ１３，Ｓ１４）で処理する。
【００８５】
Ｓ１１で太陽風とプラズマ塊の境界条件を入力する。惑星を半径Ｒのプラズマ球とする。太陽風のイオンは陽子（Ｈ+ ）であり，半径Ｒの惑星のイオンは酸素イオン（Ｏ+ ）であるとする。
【００８６】
Ｓ１２でイオン反応条件を入力する。
【００８７】
前述したように，（１）惑星を含むプラズマ球の外側（１．０Ｒと１．３Ｒの昼側）で，光電離（Ｏ＋ｈν→Ｏ+ ＋ｅ- ）があるとする。本実施の形態では生成率は０．２５６ωpiとした。
【００８８】
（３）  同じ領域に荷電交換（Ｈ+ ＋Ｏ→Ｈ＋Ｏ+ ）があるとする。本実施の形態では生成率は０．２１２ωpiとした。
【００８９】
Ｓ１３で，粒子密度平均を求める。粒子密度の平均は，宇宙空間の全粒子を元に図４の形状関数をもとに求める。
【００９０】
図４の形状関数の横軸の−１と＋１はシミュレーションの対象にする宇宙空間を６４×６４×６４に分割した格子の１つの面に対応するので，格子点での平均粒子密度はほぼ格子内の粒子密度の平均になる。
【００９１】
Ｓ１４で粒子速度の平均を求める。粒子密度と同様に図４の形状関数を使用して空間の全粒子により平均を求める。
【００９２】
Ｓ１３，Ｓ１４の結果は，Ｓ４とＳ５で呼び出されて使用される。
【００９３】
図５は上記のモデルの説明図である。
【００９４】
１１０は惑星であって，半径Ｒのプラズマ球である。イオンは酸素イオン（Ｏ+ ）である。半径１．０Ｒと１．３Ｒの領域で，光電離反応（Ｏ＋ｈν→Ｏ+ ＋ｅ- ）および荷電交換（Ｈ+ ＋Ｏ→Ｈ＋Ｏ+ ）がある。
【００９５】
図６（ａ），図６（ｂ）はシミュレーション結果を示す。
図６（ａ）と図６（ｂ）は，それぞれ，時刻３７．５／ωｃｉの光電離のみ，および荷電交換がある場合の，ｙ＝Ｌ／2 面における電子速度のｘ成分を示す。惑星中心は（１０２．４，１０２．４）の位置にある。
【００９６】
図６（ａ）（ｂ）に示されるように，バウショックが生成される。
【００９７】
図７（ａ）と図７（ｂ）は，それぞれ，光電離のみがある場合と，荷電交換がある場合に，ｙ＝Ｌ／２面における，磁束密度を示す。磁気障壁が図示のように形成される。図７（ａ）は（ｂ）を比較する時，磁気障壁の強さは，荷電交換が含まれる時の方が大きいことが示されている。
【００９８】
図８（ａ）は光電離のみがある場合のｘ＝Ｌ／２面の磁場強度を示す。図８（ｂ）は荷電交換がある場合のｘ＝Ｌ／２面の磁場強度を示す。
【００９９】
惑星中心から測定されたショックの高さは，図８（ａ）において，＋ｚ方向において惑星中心に対して，他の側より７．９８パーセント大きい。非対称性の方向は，光電離も荷電交換もない前回のシミュレーションと一致する。しかし，荷電交換が含まれる図８（ｂ）の場合は，ショック高度の惑星中心に対する非対称性は図８（ａ）の場合と反対になり，実際の観測結果に一致する。惑星からのショックの高さは，電場方向（−ｚ方向）において反対方向側より７．７９％大きい。
【０１００】
ショック高度が異なる原因を確かめるために，イオンの速度分布を考察した。
【０１０１】
図９（ａ），（ｂ）は，ｙ＝ｚ＝Ｌ／２ラインに沿ってのイオン速度の位相空間分布（ｖz −ｘ）を示す。図９（ａ）は光電離のみがある場合であり，図９（ｂ）は荷電交換がある場合を示す。図９（ａ）と（ｂ）の違いは，〔ｘ／（ｃ／ωpi）＝７０付近〕のエッジ部分である。荷電交換がある場合には，そこで惑星イオンは−ｚ方向に加速される。
【０１０２】
この加速に対する理由として，一般化されたオームの法則におけるホール項（ｊ×Ｂ項，ここにｊは電流であり，ホール電流と呼ばれる）を考える。図１０は，ｊ×Ｂのｚ成分を示す（ｘ成分およびｙ成分は，ｚ成分より小さい）。図１０（ａ）と図１０（ｂ）を比較すると，その大きさは，荷電交換がある時の方が大きい。Ｊ×Ｂのｚ成分は，主にｊx とＢy により生成される（ｘ軸は太陽風の流れの方向であり，ｙ軸は周囲の磁場方向である）。電流ｊｘは，荷電交換が生じたことによりイオンと電子の間の速度の差により生成される。
【０１０３】
イオンは，ホール項に起因する電場により影響され，そして−ｚ方向の側に加速される。図１１（ａ），（ｂ）は，ｘ＝８０ｃ／ωｐｉ，ｙ＝Ｌ／２，Ｚ＝８０ｃ／ωｐｉにおける速度空間（ｖz −ｖx ）の分布を示す。図１１（ａ）は光電離のみがある場合であり，図１１（ｂ）は荷電交換がある場合である。荷電交換がある場合には，イオンが電場により散乱されることを示している。
【０１０４】
上記は荷電交換がＨ+ ＋Ｏ→Ｈ＋Ｏ+ の場合のシミュレーションであるが，さらに，次の荷電交換反応がある場合についてもシミュレーションを行なった。
【０１０５】
Ｈ+ ＋Ｈ→Ｈ＋Ｈ+ 
その結果のショック高度の非対称の方向は，酸素の場合と同様であることが示された。従って，ホール電流がショックの非対称を生成するものであることが確かめられ，それは実際の観測結果と一致している。
【０１０６】
上記のように，前回のシミュレーションでは衝撃波（ショック）の高さが電場方向（−（ベクトルｖe ）×（ベクトルＢ）の方向（−ｚ方向））で低いという結果がでていたが，実際に金星を調べた結果によると，電場の方向でショック高度は大きいことが分かっている。惑星をプラズマ塊と見做した場合にも，シミュレーションでは実際と異なった結果が得られていたが，荷電交換を考慮してシミュレーションを行なった場合には，実際の観測結果とシミュレーション結果が一致する。光電離のみ考慮した場合にはシミュレーションと実際の観測結果は一致しなかった。荷電交換がある場合には，太陽風イオンが消え，Ｏ+ が生成され，余分の電流が生成される。その結果，強い磁気バリアが作られる。そのため，この磁気バリアがショックの高さの違いを生じるものと考えられる。
【０１０７】
【発明の効果】
上記のように，本発明によれば，プラズマ現象の解析にイオン反応を与えることが容易にできるようになり，正確なプラズマ現象のシミュレーションが可能になる。
【０１０８】
特に，本発明を宇宙プラズマの解析装置として使用した場合には，惑星の前面に生じる衝撃波の電界方向の非対象性を再現することが可能になり，宇宙プラズマを正確にシミュレーションすることができる。
【図面の簡単な説明】
【図１】本発明の基本構成を示す図である。
【図２】本発明の装置構成の実施の形態を示す図である。
【図３】本発明の実施の形態のフローチャートを示す図である。
【図４】形状関数の例を示す図である。
【図５】宇宙プラズマのモデルの例を示す図である。
【図６】本発明のシミュレーション結果を示す図である。
【図７】本発明のシミュレーション結果を示す図である。
【図８】本発明のシミュレーション結果を示す図である。
【図９】本発明のシミュレーション結果を示す図である。
【図１０】本発明のシミュレーション結果を示す図である。
【図１１】本発明のシミュレーション結果を示す図である。
【図１２】プラズマのモデルの例を示す図である。
【図１３】プラズマ塊モデルに基づくシミュレーション結果の例を示す図である。
【図１４】プラズマ現象の解析に使用するパラメータおよび初期条件の例を示す図である。
【図１５】プラズマ塊モデルに基づくシミュレーション結果を示す図である。
【図１６】プラズマ塊モデルに基づくシミュレーション結果を示す図である。
【符号の説明】
１：初期条件保持手段
２：プラズマ境界条件保持手段
３：イオン化学反応条件保持手段
４：プラズマ解析手段
５：解析結果出力手段
６：入力手段
７：状態データ保持手段
【図１】
0003735704.tif 000018
【図２】
0003735704.tif 000019
【図３】
0003735704.tif 000020
【図４】
0003735704.tif 000021
【図５】
0003735704.tif 000022
【図６】
0003735704.tif 000023
【図７】
0003735704.tif 000024
【図８】
0003735704.tif 000025
【図９】
0003735704.tif 000026
【図１０】
0003735704.tif 000027
【図１１】
0003735704.tif 000028
【図１２】
0003735704.tif 000029
【図１３】
0003735704.tif 000030
【図１４】
0003735704.tif 000031
【図１５】
0003735704.tif 000032
【図１６】
0003735704.tif 000033