JP3735704B2 - Plasma analysis apparatus and method - Google Patents

Description

【００１６】
シミュレーションは宇宙空間を６４×６４×６４に分割された空間（格子空間）とし，その内部の粒子密度，粒子速度等の平均値を格子点におけるそれぞれの値とする。そして，格子点での値をデータとして差分化した運動方程式，マックスウェル方程式等を使用して宇宙プラズマ現象をシミュレーションした。
【００１７】
シミュレーションにおけるアルベン速度Ｖ_A と光速度ｃの比をＶ_A ／ｃ＝ω_ci／ω_pi＝１．０×１０^-4とした。ω_ciはイオンのサイクロトロン角周波数である。ω_piはイオンプラズマ振動の角周波数である。
【００１８】
熱圧力と磁気圧力の比は電子とイオンについて，β_e ＝β_i ＝１とした。
太陽風の速度はｖ_SW＝4 ．０Ｖ_A である。
【００１９】
電離層に包まれた惑星を半径Ｒのプラズマ球として表現する。太陽風はｘ軸方向に左から右に流れ，磁束密度はｙ方向である。初期状態において，太陽風の動圧と惑星プラズマの圧力はほぼつりあっている。Ｌはシミュレーションする宇宙空間の大きさを示す。プラズマ球へのイオンの供給率は０．２５６ω_piである。
【００２０】
図１５（ａ），（ｂ），（ｃ）は，惑星をプラスマ球としたモデルにおけるシミュレーション結果を示す。
【００２１】
図１５（ａ）はｙ＝Ｌ／２での断面図であり，シミュレーション開始後，ω_ciｔ＝３７．５の時点でのイオン速度のｘ成分を示す。イオン速度はＶｉｘ／Ｖ_A で表している。図１３（ａ）において，惑星中心は座標（１０２．４，１０２．４）にあり，太陽風はｘ軸方向（左から右）に吹いている。図示のようにショックが形成される。惑星半径は，Ｒ＝６．４ｒ_L である（ｒ_L ＝ｖ_SW／ω_ci）である。
【００２２】
図１５（ｂ）はｙ＝Ｌ／２での断面図における磁束密度の大きさを示す。磁束密度は｜Ｂ｜／Ｂｏで表している（Ｂｏは初期条件で与えた磁束密度の大きさである）。
図１５（ｃ）はｙ＝Ｌ／２での断面図における惑星の酸素イオンの分布を示す。図１６（ａ）は ｘ＝Ｌ／２での断面図であり，シミュレーション開始後，ω_ciｔ＝３７．５の時点でのシミュレーション結果での磁束密度の大きさを示す。 図１６（ａ）において，惑星は実線の円で示され，ショック面は惑星を囲むぼんやりした円で示される。
図１６（ａ）のように，シミュレーションではショックの惑星からの高さはｚ軸の正方向で大きく，ｚ軸の負方向で小さい（ショック面は惑星に対してｚ軸の正方向にずれている。
【００２３】
これは，実際の観測結果と異なり，実際の観測結果はショックの惑星からの高さは，図１６（ｃ）に示すようにｚ軸の正方向で小さく，ｚ軸の負方向で大きい（ショック面は惑星に対してｚ軸の負方向にずれている）ことが観測されている。
【００２４】
図１６（ｂ）は ｘ＝３Ｌ／４での断面図であり，シミュレーション開始後，ω_ciｔ＝３７．５の時点での磁束密度の大きさを示す。
【００２５】
上記の太陽風と惑星プラズマの相互作用のモデルでは，惑星の高層で生じるイオン反応（光電離や荷電交換等）が扱われていなかった。そのためショック面と惑星との位置関係のシミュレーション結果が必ずしも正しいものではなかった。
【００２６】
本発明は，イオンの化学反応が生じるプラズマ系においても正確に三次元でシミュレーションできるプラズマ解析装置および方法を提供することを目的とする。
【００２７】
【課題を解決するための手段】
図１は，本発明の基本構成を示す。
【００２８】
図１において，１は初期条件保持手段であって，入力手段６により入力された初期条件を入力するものである。２はプラズマの境界条件保持手段であって，入力手段６により入力されたプラズマ境界条件を保持するものである。３はイオン化学反応条件保持手段であって，入力手段６により入力されたプラズマの化学反応についての条件を保持するものである。４はプラズマ解析手段であって，粒子の運動方程式，マックスウェル方程式，初期条件，境界条件，イオン化学反応条件およびプラズマ状態を表すデータ（イオン生成率等のデータ）によりプラズマ解析をするものである。５は解析結果出力手段であって，解析結果を出力するものである。６は入力手段である。７はプラズマの状態データ保持手段である。
【００２９】
図１の本発明の基本構成の動作を説明する。
【００３０】
入力手段６により，初期条件，プラズマの境界条件，プラズマに生じる化学反応の条件を入力する。また，プラズマのイオンサイクロトロン周波数，イオン生成率等のプラズマ状態データを入力し，状態データ保持手段７に保持する。プラズマ解析手段４は，初期条件，境界条件，イオン化学反応条件およびプラズマの状態データによりプラズマに生じる状態を解析する。出力手段５は求められたプラズマ状態を出力する。
【００３１】
化学反応条件は，例えば，太陽風と惑星の相互作用の場合，太陽風のイオン（Ｈ⁺ ）と惑星イオン（Ｏ⁺ ）とした時、
Ｈ⁺ ＋Ｏ→Ｈ＋Ｏ⁺ （惑星プラズマ）
とする等である。Ｏ，Ｈはイオン化していない原子なので，解析においてプラズマから除かれる。このような化学反応条件の場合，イオン反応のため，電流変化となってプラズマの運動などに影響するものである。後述するが，太陽風と惑星イオンとの相互作用のシミュレーションにおいて，本発明のシミュレーション装置および方法によれば，惑星とショックの高さの関係は，前述の本発明に先行する方法では，得られなかったのであるが，図８に示すように実際の観測結果と一致するようになった。
【００３２】
本発明によれば，プラズマの初期条件，境界条件だけでなく，イオン反応条件もシミュレーションの条件にしてプラズマを解析できるので，イオンの化学反応が生じるプラズマ系においても正確に三次元でシミュレーションできるようになる。
【００３３】
なお，本発明のプラズマ解析装置および方法は，宇宙プラズマに適用するだけでなく，実験室等のプラズマにも適用できるものである。
【００３４】
【発明の実施の形態】
本発明を実施するためのハイブリッドコードについて説明する。
以下において，Ｂ，Ｅ，ｖ，ｊは全て３次元ベクトルであり，それぞれベクトルＢ，ベクトルＥ，ベクトルｖ，ベクトルｊ等で表す。ベクトルｘは位置ベクトルである。Ｂは磁束密度，Ｅは電界，ｖはイオン速度，ｊは電流をそれぞれ表す。
【００３５】
ハイブリッドコードにより解析する場合の基礎方程式は次のとおりである。
【００３６】
イオンの運動方程式
【００３７】
【数２】

【００３８】
イオンの速度の式
【００３９】
【数３】

【００４０】
電子流体（質量なし）の運動方程式
【００４１】
【数４】

【００４２】
ただし，電子温度Ｔ_e は一定とする。
【００４３】
マックスウェル方程式
【００４４】
【数５】

【００４５】
【数６】

【００４６】
準中性条件
ｎ_e ＝ｎ_i ＝ｎ （６）
本実施の形態では電気伝導度は無限大と仮定する（η＝０）。
【００４７】
シミュレーションでは上記の式を差分化して使用する。
【００４８】
上記のマックスウェルの方程式（４）を差分化すると次のようになる。
【００４９】
【数７】

【００５０】
ベクトルＢ_k+1 から式（５）を用いてベクトルｊ_k+1 を得る。
【００５１】
式（１），（２）を差分化すると次のようになる。
【００５２】
【数８】

【００５３】
【数９】

【００５４】
時刻ｋの磁束密度は時刻ｋ＋１の磁束密度と時刻ｋ−１の磁束密度とにより次のように表せる。
【００５５】
【数１０】

【００５６】
時刻ｋ＋１の位置ベクトルは次のように表せる。
【００５７】
【数１１】

【００５８】
イオン速度と密度を格子点毎に足しあわせることにより格子点の速度ベクトルｖ_i と密度ｎ_i が得られる。電子流体の速度ベクトルｖ_e は電流を表す次の式で求まる。
【００５９】
【数１２】

【００６０】
電界ベクトルＥは電子速度ベクトルｖ_e と電子流体の運動方程式と準中性条件の式とにより求めることができる。
【００６１】
格子点におけるイオンの密度，速度は，格子点から近いものから大きく，遠いものは寄与が小さいようにする。
【００６２】
具体的には，図４のような形状関数を使用する。
【００６３】
密度は
【００６４】
【数１３】

【００６５】
となる。
【００６６】
速度平均は
【００６７】
【数１４】

【００６８】
で求まる。Ｎはシミュレーション対象の宇宙空間の全粒子数である。
【００６９】
図２は本発明の装置構成の実施の形態を示す。
【００７０】
図２は，イオンを粒子とし電子を流体として扱うハイブリットコードにより解析するための構成を示す。
【００７１】
１１は入力する初期条件を表す（磁束密度ベクトルＢ₀ ，Ｂ₁ ，電界ベクトルＥ₀ ，粒子密度ｎ₀ ，粒子速度ベクトルｖ_i0）。１２は入力手段である。１３は演算手段を表す。１４は演算結果の時刻ｋ＋１における電界ベクトルＥ_k+1 ，磁束密度ベクトルＢ_k+1 ，電流密度ベクトルＪ_k+1 ，粒子密度ｎ_k+1 ，粒子の位置ベクトルｘ_k+1 ，粒子の速度ベクトルｖ_k+1 である。１５は出力手段である。２０はデータ保持手段である。２１は初期条件保持部である。２２は境界条件保持部である。１７は状態データであって，イオンの生成率等のパラメータを表す。３１はハイブリッドコードにより解析する場合の演算式を保持する演算式保持部であり，上記の差分化方程式，その他演算に必要な式を保持するものである。３２は差分化イオン運動方程式である。３３は差分化電子運動方程式である。３４は差分化マックスウェル方程式である。３５は電子流体の速度の式である。３６は粒子中性条件であって，ｎ_i ＝ｎ_e ＝ｎを保持するものである。３６は粒子密度の平均，粒子速度の平均値を求める式である。３７はその他演算に必要な式である。
【００７２】
本発明のプラズマ解析装置を使用して，以下のモデルにより太陽風と惑星（火星，金星等）の相互作用をシミュレーションする場合について説明する。プラズマの生成率，プラズマ球の半径，初期の磁束密度は先行シミュレーション（前述のプラズマ塊モデル）の場合と同じである。
【００７３】
（１） 半径Ｒの惑星を表すプラズマ球の外側（１．０Ｒと１．３Ｒの昼側）で，光電離があるとする。
【００７４】
Ｏ＋ｈν→Ｏ⁺ ＋ｅ^-
（２） 同じ領域に荷電交換があるとする。
【００７５】
Ｈ⁺ ＋Ｏ→Ｈ＋Ｏ⁺
中性化した分子は計算に無関係になるので，（１）の場合はＯ⁺ がある割合で生成することになる。太陽風のＨ⁺ は変化させないものとする。生成する位置は，上記の領域内で一様に分布するように乱数を用いて決める。（２）の場合は，ある割合で太陽風が消え，その分Ｏ⁺ が生成することに等価である。即ち，１．０Ｒと１．３Ｒの間に達した太陽風のイオン（Ｈ⁺ ）のうち，あらかじめ決めておいた割合を消滅させる。消滅した位置にＯ⁺ を生成させる。この結果は粒子の平均化の処理（後述する）によりベクトルｖ_i とｎの変化にあらわれ，全体の計算に影響する。（１），（２）において，イオン速度は，温度が惑星イオンの温度と同じ速度分布になるように乱数で決める。速度分布は次の式（１３）で表されるマクスウェル分布となるようにする。
【００７６】
【数１５】

【００７７】
Ａは次になるように定める。
【００７８】
【数１６】

【００７９】
（１／２）ｍｖ_th ² ＝ｋ_B Ｔである。ｎは粒子密度，ｍはイオン質量，ｋ_B はボルツマン定数，Ｔは温度である。
【００８０】
図３は図２の構成の動作のフローチャートである。
図３の動作を説明する。
Ｓ１で磁束密度，電界，粒子速度，密度の初期条件ベクトルＢ₀ ，ベクトル Ｂ₁ ，ベクトルＥ₁ ，ベクトルｖ_i0，ｎを入力する。Ｂ₀ の初期条件は，図１４（ｂ）と同じである。
【００８１】
初期条件は，太陽風および惑星のイオン速度がマクスウェル分布になるように与える。太陽風温度と惑星プラズマの温度は同じにする。
【００８２】
Ｓ２で時刻ｋの電界ベクトルＥ_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 ）×（ベクトルＢ）の方向（−ｚ方向））で低いという結果がでていたが，実際に金星を調べた結果によると，電場の方向でショック高度は大きいことが分かっている。惑星をプラズマ塊と見做した場合にも，シミュレーションでは実際と異なった結果が得られていたが，荷電交換を考慮してシミュレーションを行なった場合には，実際の観測結果とシミュレーション結果が一致する。光電離のみ考慮した場合にはシミュレーションと実際の観測結果は一致しなかった。荷電交換がある場合には，太陽風イオンが消え，Ｏ⁺ が生成され，余分の電流が生成される。その結果，強い磁気バリアが作られる。そのため，この磁気バリアがショックの高さの違いを生じるものと考えられる。
【０１０７】
【発明の効果】
上記のように，本発明によれば，プラズマ現象の解析にイオン反応を与えることが容易にできるようになり，正確なプラズマ現象のシミュレーションが可能になる。
【０１０８】
特に，本発明を宇宙プラズマの解析装置として使用した場合には，惑星の前面に生じる衝撃波の電界方向の非対象性を再現することが可能になり，宇宙プラズマを正確にシミュレーションすることができる。
【図面の簡単な説明】
【図１】本発明の基本構成を示す図である。
【図２】本発明の装置構成の実施の形態を示す図である。
【図３】本発明の実施の形態のフローチャートを示す図である。
【図４】形状関数の例を示す図である。
【図５】宇宙プラズマのモデルの例を示す図である。
【図６】本発明のシミュレーション結果を示す図である。
【図７】本発明のシミュレーション結果を示す図である。
【図８】本発明のシミュレーション結果を示す図である。
【図９】本発明のシミュレーション結果を示す図である。
【図１０】本発明のシミュレーション結果を示す図である。
【図１１】本発明のシミュレーション結果を示す図である。
【図１２】プラズマのモデルの例を示す図である。
【図１３】プラズマ塊モデルに基づくシミュレーション結果の例を示す図である。
【図１４】プラズマ現象の解析に使用するパラメータおよび初期条件の例を示す図である。
【図１５】プラズマ塊モデルに基づくシミュレーション結果を示す図である。
【図１６】プラズマ塊モデルに基づくシミュレーション結果を示す図である。
【符号の説明】
１：初期条件保持手段
２：プラズマ境界条件保持手段
３：イオン化学反応条件保持手段
４：プラズマ解析手段
５：解析結果出力手段
６：入力手段
７：状態データ保持手段[0001]
BACKGROUND OF THE INVENTION
The present invention relates to a plasma analysis apparatus and an analysis method. In particular, the plasma density and velocity distributions and electromagnetic fields can be simulated in three dimensions even in systems where chemical reactions of ions occur, based on the equations of motion and electromagnetic field Maxwell equations described using space plasma as a simulation model. It is.
[0002]
The present invention can be used in space development related industries, plasma engineering fields, and the like.
[0003]
[Prior art]
In computer simulations of plasma phenomena such as the interaction between the solar wind and planetary plasmas, hybrid codes that use ions as particles and electrons as fluids have been used. Using a model that regards the planet as a non-magnetic solid sphere, the plasma interaction between the solar wind and the planet can be simulated. An example of applying the hybrid code to space plasma is explained in the following paper.
[0004]
(1) Leroy, M .; M.M. , D.D. Winske, C.I. C. Goodrich, C.I. S. Wu, and K.W. Papadopoulos, “The structure of perennial bow shocks”, Journal of Geophysical Research, vol. 87, pp. 5081-5094 (1982).
[0005]
FIG. 12A shows a model for analyzing the plasma interaction between the solar wind and the planet, regarding the planet as a non-magnetic solid.
[0006]
The planet 110 is regarded as a non-magnetic solid sphere, and the solar wind, which is a plasma flow, is irradiated. Solar wind plasma ions irradiating the planet are reflected on the planet surface, and some ions enter the planet and disappear inside.
[0007]
The simulation results of this model using a hybrid code are explained in detail in the following papers.
[0008]
(2) Shimazu, H .; , "Three-dimensional hybrid simulation of magnetoplasma round around an obscure", Earth, Planets and Space, vol. 51, pp 383-393 (1999).
[0009]
[Problems to be solved by the invention]
Description of Prior Art When the ionosphere is present in the upper layer of the planet, it is actually inappropriate to consider the planet as a solid and simulate the interaction with the solar wind as shown in Fig. 12 (a). . There is an ionosphere above the actual planet, and the planet can be regarded as a plasma mass. Therefore, the inventor of the present application tried to simulate the planet as a plasma mass (“3D hybrid simulation of the interaction between the solar wind and the Martian Venus-type planet”) Kaikan, November 9-12, 1999, presentation number B12-P054).
[0010]
FIG. 12B shows a model in which the planet is treated as a plasma mass of a sphere having a radius R and the solar wind irradiates the plasma mass.
[0011]
In FIG. 12B, reference numeral 110 denotes a planet, which is a plasma mass. Solar wind ions (H ⁺ ) irradiating the planet may enter the planet and mix with the ions (O ⁺ ) in the planet.
[0012]
FIG. 13 schematically shows the results of computer simulation of the plasma density, velocity distribution, and electromagnetic field by solving the equation of motion and the maximum equation of ions using the hybrid code for this model.
[0013]
In FIG. 13, reference numeral 110 denotes a planet, which is a plasma mass. 111 is a shock wave (bow shock, hereinafter referred to as shock).
As shown in FIG. 13, the solar wind is a supersonic plasma flow, and when this plasma flow enters the magnetosphere of the planet, a shock 111 is generated in the front.
[0014]
FIG. 14 (a) shows parameters for analysis using a plasma mass model. FIG. 14 (b) shows a model that gives initial conditions and boundary conditions.
In FIG. 14B, the magnetic flux density is given by the following equation.
[0015]
[Expression 1]

[0016]
In the simulation, the outer space is divided into 64 × 64 × 64 spaces (lattice spaces), and the average values of particle density, particle velocity, etc. in the space are used as the respective values at the lattice points. Then, the space plasma phenomenon was simulated using the equations of motion, Maxwell equations, etc., which are obtained by substituting the values at the grid points as data.
[0017]
The ratio of the alben velocity V _A and the light velocity c in the simulation was V _A / c = ω _ci / ω _pi = 1.0 × 10 ⁻⁴ . ω _ci is the cyclotron angular frequency of the ion. ω _pi is the angular frequency of ion plasma oscillation.
[0018]
The ratio of thermal pressure to magnetic pressure was β _e = β _i = 1 for electrons and ions.
The speed of the solar wind is v _SW = 4. 0V _A.
[0019]
A planet wrapped in the ionosphere is expressed as a plasma sphere of radius R. The solar wind flows from left to right in the x-axis direction, and the magnetic flux density is in the y direction. In the initial state, the dynamic pressure of the solar wind and the pressure of the planetary plasma are almost balanced. L indicates the size of outer space to be simulated. The supply rate of ions to the plasma sphere is 0.256Ω _pi .
[0020]
FIGS. 15A, 15B, and 15C show simulation results in a model in which a planet is a plasma sphere.
[0021]
FIG. 15A is a cross-sectional view at y = L / 2, and shows the x component of the ion velocity at the time of ω _ci t = 37.5 after the start of the simulation. The ion velocity is expressed as Vix / _VA . In FIG. 13A, the center of the planet is at coordinates (102.4, 102.4), and the solar wind is blowing in the x-axis direction (from left to right). A shock is formed as shown. Planet radius is a is _{R = 6.4r L (r L =} v SW / ω ci).
[0022]
FIG. 15B shows the magnitude of the magnetic flux density in the cross-sectional view at y = L / 2. The magnetic flux density is represented by | B | / Bo (Bo is the magnitude of the magnetic flux density given in the initial condition).
FIG. 15C shows the distribution of planetary oxygen ions in a cross-sectional view at y = L / 2. FIG. 16A is a cross-sectional view at x = L / 2, and shows the magnitude of the magnetic flux density in the simulation result at the time of ω _ci t = 37.5 after the start of the simulation. In FIG. 16 (a), the planet is indicated by a solid circle, and the shock surface is indicated by a blurred circle surrounding the planet.
As shown in FIG. 16A, in the simulation, the height of the shock from the planet is large in the positive direction of the z axis and small in the negative direction of the z axis (the shock surface is shifted in the positive direction of the z axis with respect to the planet). Yes.
[0023]
This is different from the actual observation results, and the actual observation results show that the height of the shock from the planet is small in the positive z-axis direction and large in the negative z-axis direction (shock) as shown in FIG. It has been observed that the plane is offset in the negative z-axis direction with respect to the planet).
[0024]
FIG. 16B is a cross-sectional view at x = 3L / 4, and shows the magnitude of the magnetic flux density at the time of ω _ci t = 37.5 after starting the simulation.
[0025]
In the above model of the interaction between the solar wind and planetary plasma, ion reactions (photoionization, charge exchange, etc.) that occur in the upper layers of the planet were not handled. Therefore, the simulation results of the positional relationship between the shock surface and the planet were not always correct.
[0026]
It is an object of the present invention to provide a plasma analysis apparatus and method that can accurately perform a three-dimensional simulation even in a plasma system in which a chemical reaction of ions occurs.
[0027]
[Means for Solving the Problems]
FIG. 1 shows the basic configuration of the present invention.
[0028]
In FIG. 1, reference numeral 1 denotes an initial condition holding means for inputting the initial condition input by the input means 6. Reference numeral 2 denotes plasma boundary condition holding means for holding the plasma boundary condition input by the input means 6. Reference numeral 3 denotes ion chemical reaction condition holding means for holding conditions for the chemical reaction of the plasma input by the input means 6. 4 is a plasma analysis means for performing plasma analysis based on particle motion equations, Maxwell equations, initial conditions, boundary conditions, ion chemical reaction conditions, and data (data such as ion generation rate) representing plasma states. . Reference numeral 5 denotes analysis result output means for outputting the analysis result. 6 is an input means. Reference numeral 7 denotes plasma state data holding means.
[0029]
The operation of the basic configuration of the present invention shown in FIG. 1 will be described.
[0030]
The input means 6 inputs initial conditions, plasma boundary conditions, and conditions for chemical reactions that occur in the plasma. Further, plasma state data such as the ion cyclotron frequency and ion generation rate of the plasma is input and held in the state data holding means 7. The plasma analysis means 4 analyzes the state generated in the plasma based on the initial conditions, boundary conditions, ion chemical reaction conditions, and plasma state data. The output means 5 outputs the obtained plasma state.
[0031]
For example, in the case of the interaction between the solar wind and the planet, the chemical reaction conditions are solar wind ions (H ⁺ ) and planetary ions (O ⁺ ).
H ⁺ + O → H + O ⁺ (planet plasma)
And so on. Since O and H are non-ionized atoms, they are removed from the plasma in the analysis. In the case of such chemical reaction conditions, because of the ionic reaction, the current changes and affects the movement of the plasma. As will be described later, in the simulation of the interaction between the solar wind and planet ions, according to the simulation apparatus and method of the present invention, the relationship between the planet and the shock height cannot be obtained by the method preceding the present invention. However, as shown in FIG. 8, it became consistent with the actual observation results.
[0032]
According to the present invention, the plasma can be analyzed using not only the initial conditions and boundary conditions of the plasma but also the ion reaction conditions as the simulation conditions, so that a three-dimensional simulation can be performed accurately even in a plasma system in which an ion chemical reaction occurs. become.
[0033]
The plasma analysis apparatus and method according to the present invention can be applied not only to space plasma but also to plasma in a laboratory.
[0034]
DETAILED DESCRIPTION OF THE INVENTION
A hybrid code for carrying out the present invention will be described.
In the following, B, E, v, and j are all three-dimensional vectors, and are represented by vector B, vector E, vector v, vector j, and the like, respectively. Vector x is a position vector. B represents the magnetic flux density, E represents the electric field, v represents the ion velocity, and j represents the current.
[0035]
The basic equations for analysis using the hybrid code are as follows.
[0036]
Equation of motion of ions
[Expression 2]

[0038]
Ion velocity equation
[Equation 3]

[0040]
Equation of motion of electron fluid (without mass)
[Expression 4]

[0042]
However, the electron temperature _Te is constant.
[0043]
Maxwell equation [0044]
[Equation 5]

[0045]
[Formula 6]

[0046]
Quasi-neutral condition _ne = _ni = n (6)
In this embodiment, it is assumed that the electrical conductivity is infinite (η = 0).
[0047]
In the simulation, the above formula is differentiated and used.
[0048]
The above difference of Maxwell's equation (4) is as follows.
[0049]
[Expression 7]

[0050]
A vector j _{k + 1} is obtained from the vector B _{k + 1} using equation (5).
[0051]
When formulas (1) and (2) are differentiated, they are as follows.
[0052]
[Equation 8]

[0053]
[Equation 9]

[0054]
The magnetic flux density at time k can be expressed as follows by the magnetic flux density at time k + 1 and the magnetic flux density at time k−1.
[0055]
[Expression 10]

[0056]
The position vector at time k + 1 can be expressed as follows.
[0057]
## EQU11 ##

[0058]
By adding the ion velocity and density for each lattice point, the velocity vector v _i and density n _i of the lattice point are obtained. The velocity vector v _e of the electronic fluid is obtained by the following expression representing the current.
[0059]
[Expression 12]

[0060]
Electric field vector E can be determined by the equation of motion equation and the quasi-neutral condition of the electron velocity vector v _e and e fluid.
[0061]
The density and velocity of ions at lattice points are set so as to increase from those close to the lattice point, and those far from the lattice point contribute little.
[0062]
Specifically, a shape function as shown in FIG. 4 is used.
[0063]
Density is [0064]
[Formula 13]

[0065]
It becomes.
[0066]
Speed average is [0067]
[Expression 14]

[0068]
It is obtained by N is the total number of particles in the space to be simulated.
[0069]
FIG. 2 shows an embodiment of the apparatus configuration of the present invention.
[0070]
FIG. 2 shows a configuration for analyzing with a hybrid code that handles ions as particles and electrons as fluids.
[0071]
11 represents initial conditions to be input (magnetic flux density vectors B ₀ and B ₁ , electric field vector E ₀ , particle density n ₀ , particle velocity vector v _i0 ). Reference numeral 12 denotes input means. Reference numeral 13 denotes a calculation means. 14 is the electric field vector E _{k + 1} , magnetic flux density vector B _{k + 1} , current density vector J _{k + 1} , particle density n _{k + 1} , particle position vector x _{k + 1} , particle velocity The vector v _{k + 1} . Reference numeral 15 denotes output means. Reference numeral 20 denotes data holding means. Reference numeral 21 denotes an initial condition holding unit. Reference numeral 22 denotes a boundary condition holding unit. Reference numeral 17 denotes state data representing parameters such as ion generation rate. Reference numeral 31 denotes an arithmetic expression holding unit that holds an arithmetic expression in the case of analysis using a hybrid code, and holds the above-described difference equation and other expressions necessary for the calculation. 32 is a differential ion motion equation. Reference numeral 33 denotes a difference electron motion equation. Reference numeral 34 denotes a differentiated Maxwell equation. 35 is an equation for the velocity of the electronic fluid. Reference numeral 36 denotes a particle neutral condition, which holds n _i = n _e = n. Reference numeral 36 denotes an equation for obtaining an average particle density and an average particle velocity. 37 is an expression required for other calculations.
[0072]
The case where the interaction between the solar wind and planets (Mars, Venus, etc.) is simulated using the following model using the plasma analysis apparatus of the present invention will be described. The plasma generation rate, plasma sphere radius, and initial magnetic flux density are the same as in the previous simulation (the plasma lump model described above).
[0073]
(1) Suppose that there is photoionization outside the plasma sphere representing the planet of radius R (the dayside of 1.0R and 1.3R).
[0074]
O + hν → O ⁺ + e ⁻
(2) Assume that there is charge exchange in the same region.
[0075]
H ⁺ + O → H + O ⁺
Since neutralized molecules are irrelevant to calculation, in the case of (1), O ⁺ is generated at a certain ratio. The solar wind H ⁺ is not changed. The positions to be generated are determined using random numbers so that they are uniformly distributed in the above-described region. In the case of (2), the solar wind disappears at a certain rate, which is equivalent to the generation of O ⁺ . That is, a predetermined ratio of the solar wind ions (H ⁺ ) reaching between 1.0R and 1.3R is eliminated. O ⁺ is generated at the disappeared position. This result appears in the change of the vectors v _i and n by the particle averaging process (described later) and affects the overall calculation. In (1) and (2), the ion velocity is determined by random numbers so that the temperature has the same velocity distribution as the planetary ion temperature. The velocity distribution is set to be a Maxwell distribution represented by the following equation (13).
[0076]
[Expression 15]

[0077]
A is determined as follows.
[0078]
[Expression 16]

[0079]
(1/2) mv _th ² = k _B T. n is the particle density, m is the ion mass, k _B is the Boltzmann constant, and T is the temperature.
[0080]
FIG. 3 is a flowchart of the operation of the configuration of FIG.
The operation of FIG. 3 will be described.
In S1, the initial condition vector B ₀ , vector B ₁ , vector E ₁ , vector v _i0 , n of magnetic flux density, electric field, particle velocity, and density are input. The initial condition for B ₀ is the same as in FIG.
[0081]
Initial conditions are given so that the solar wind and planetary ion velocities have a Maxwell distribution. The solar wind temperature and planetary plasma temperature should be the same.
[0082]
In S2, B _{k + 1} is obtained from the electric field vectors E _k and B _k-1 at time k according to Maxwell's equations. Next, in S3, the current density vector J _{k + 1} at time k + 1 is obtained using the magnetic flux density vector B _{k + 1} at time k + 1 according to the Maxwell equation. Then the S4 ion particle density n _i, with the current density vector J _{k + 1} of the electron density n _e (= n) between the speed v _i and time k + _1, obtains the time k + 1 of the electron velocity vector v _e.
[0083]
Next, in S5, the time k + 1 of the electron velocity vector v _e, the magnetic flux density vector B _{k + 1,} the current density vector J _{k + 1,} charge density n _e, Boltzmann constant k _B, by an electronic temperature T _e, the time An electric field vector E _{k + 1} of _{k + 1} is obtained. Next, the time is updated with k = k + 2, and the processes after S2 are repeated.
[0084]
In contrast to the above main routine, boundary conditions, ion reaction conditions, particle density averaging, and particle velocity averaging are processed in a subroutine (S11, S12, S13, S14).
[0085]
In S11, the boundary condition between the solar wind and the plasma mass is input. Let the planet be a plasma sphere of radius R. The solar wind ions are protons (H ⁺ ), and the planetary ions with radius R are oxygen ions (O ⁺ ).
[0086]
In S12, ion reaction conditions are input.
[0087]
As described above, (1) it is assumed that there is photoionization (O + hν → O ⁺ + e ⁻ ) outside the plasma sphere including the planet (the dayside of 1.0R and 1.3R). In this embodiment, the generation rate is 0.256 ω _pi .
[0088]
(3) It is assumed that there is charge exchange (H ⁺ + O → H + O ⁺ ) in the same region. In this embodiment, the generation rate is 0.212 ω _pi .
[0089]
In S13, the average particle density is obtained. The average particle density is obtained based on the shape function shown in FIG. 4 based on all particles in outer space.
[0090]
Since -1 and +1 on the horizontal axis of the shape function in FIG. 4 correspond to one surface of the lattice obtained by dividing the cosmic space to be simulated into 64 × 64 × 64, the average particle density at the lattice point is approximately lattice. It becomes the average of the particle density inside.
[0091]
In S14, the average particle velocity is obtained. As with the particle density, the shape function of FIG. 4 is used to determine the average for all particles in the space.
[0092]
The results of S13 and S14 are called up and used in S4 and S5.
[0093]
FIG. 5 is an explanatory diagram of the above model.
[0094]
Reference numeral 110 denotes a planet, which is a plasma sphere having a radius R. The ions are oxygen ions (O ⁺ ). There are photoionization reactions (O + hv → O ⁺ + e ⁻ ) and charge exchange (H ⁺ + O → H + O ⁺ ) in the radius 1.0R and 1.3R regions.
[0095]
FIG. 6A and FIG. 6B show simulation results.
FIGS. 6 (a) and 6 (b) show the x component of the electron velocity on the y = L / 2 plane when only the photoionization at time 37.5 / ωci and there is charge exchange, respectively. The center of the planet is at the position (102.4, 102.4).
[0096]
As shown in FIGS. 6A and 6B, a bow shock is generated.
[0097]
FIG. 7A and FIG. 7B show the magnetic flux density on the y = L / 2 plane when there is only photoionization and when there is charge exchange, respectively. A magnetic barrier is formed as shown. FIG. 7 (a) shows that when comparing (b), the strength of the magnetic barrier is greater when charge exchange is involved.
[0098]
FIG. 8A shows the magnetic field strength of the x = L / 2 plane when there is only photoionization. FIG. 8B shows the magnetic field strength of the x = L / 2 plane when there is charge exchange.
[0099]
The height of the shock measured from the center of the planet is 7.98 percent greater than the other side relative to the center of the planet in the + z direction in FIG. 8 (a). The direction of asymmetry is consistent with the previous simulation with no photoionization or charge exchange. However, in the case of FIG. 8B in which charge exchange is included, the asymmetry of the shock height with respect to the planet center is opposite to that in FIG. 8A and matches the actual observation result. The height of the shock from the planet is 7.79% greater in the electric field direction (-z direction) than the opposite direction side.
[0100]
In order to confirm the cause of the difference in shock height, the velocity distribution of ions was considered.
[0101]
FIGS. 9A and 9B show the phase space distribution (v _z −x) of ion velocity along the y = z = L / 2 line. FIG. 9A shows the case where only photoionization is present, and FIG. 9B shows the case where there is charge exchange. The difference between FIGS. 9A and 9B is the edge portion of [near x / (c / ω _pi ) = 70]. In the presence of charge exchange, planetary ions are then accelerated in the -z direction.
[0102]
As a reason for this acceleration, consider the hole term in the generalized Ohm's law (j × B term, where j is a current and is called a hole current). FIG. 10 shows the z component of j × B (the x component and the y component are smaller than the z component). Comparing FIG. 10 (a) and FIG. 10 (b), the size is larger when there is charge exchange. Z component of J × B is mainly generated by the j _x and B _y (x-axis is the direction of the solar wind flow, the y-axis is the magnetic field direction of the circumference). Current jx is generated by the difference in velocity between ions and electrons due to charge exchange occurring.
[0103]
The ions are affected by the electric field due to the Hall term and are accelerated to the -z side. FIGS. 11A and 11B show the distribution of velocity space (v _z −v _x ) at x = 80c / ωpi, y = L / 2, and Z = 80c / ωpi. FIG. 11A shows a case where there is only photoionization, and FIG. 11B shows a case where there is charge exchange. When there is charge exchange, it indicates that ions are scattered by the electric field.
[0104]
The above is a simulation in the case where the charge exchange is H ⁺ + O → H + O ⁺ , but the simulation was also performed in the case where there was a next charge exchange reaction.
[0105]
H ⁺ + H → H + H ⁺
The result shows that the direction of asymmetrical shock height is the same as that of oxygen. Therefore, it is confirmed that the Hall current generates a shock asymmetry, which is consistent with the actual observation.
[0106]
As described above, in the previous simulation, the shock wave (shock) height was low in the electric field direction (− (vector v _e ) × (vector B) direction (−z direction)). According to the results of investigating Venus, the shock height is large in the direction of the electric field. When the planet was regarded as a plasma mass, the simulation showed a different result from the actual one. However, when the simulation was performed in consideration of charge exchange, the actual observation and the simulation result coincided. . When only photoionization was considered, the simulation and the actual observation result did not agree. When there is a charge exchange, solar wind ions disappear, O ⁺ is generated, and extra current is generated. As a result, a strong magnetic barrier is created. Therefore, this magnetic barrier is considered to cause a difference in shock height.
[0107]
【The invention's effect】
As described above, according to the present invention, an ion reaction can be easily given to the analysis of a plasma phenomenon, and an accurate simulation of the plasma phenomenon becomes possible.
[0108]
In particular, when the present invention is used as a space plasma analysis device, it becomes possible to reproduce the non-objectivity of the electric field direction of the shock wave generated in front of the planet, and the space plasma can be accurately simulated.
[Brief description of the drawings]
FIG. 1 is a diagram showing a basic configuration of the present invention.
FIG. 2 is a diagram showing an embodiment of a device configuration of the present invention.
FIG. 3 is a flowchart of an embodiment of the present invention.
FIG. 4 is a diagram illustrating an example of a shape function.
FIG. 5 is a diagram showing an example of a space plasma model.
FIG. 6 is a diagram showing a simulation result of the present invention.
FIG. 7 is a diagram showing a simulation result of the present invention.
FIG. 8 is a diagram showing a simulation result of the present invention.
FIG. 9 is a diagram showing a simulation result of the present invention.
FIG. 10 is a diagram showing a simulation result of the present invention.
FIG. 11 is a diagram showing a simulation result of the present invention.
FIG. 12 is a diagram showing an example of a plasma model.
FIG. 13 is a diagram illustrating an example of a simulation result based on a plasma mass model.
FIG. 14 is a diagram showing an example of parameters and initial conditions used for analysis of a plasma phenomenon.
FIG. 15 is a diagram showing a simulation result based on a plasma mass model.
FIG. 16 is a diagram showing a simulation result based on a plasma mass model.
[Explanation of symbols]
1: initial condition holding means 2: plasma boundary condition holding means 3: ion chemical reaction condition holding means 4: plasma analysis means 5: analysis result output means 6: input means 7: state data holding means

Claims (10)

In an apparatus for analyzing plasma comprising an input means, an arithmetic means, and an output means,
With a hybrid code that treats ions as particles and electrons as fluids, it maintains a plasma analysis means for analyzing plasma using the equations of motion of individual ions, the equations of motion of electron fluids, and Maxwell's equations. Analyzing plasma of a system irradiated with ions and electron fluid,
Maintains initial plasma conditions, boundary conditions, and chemical reaction conditions.
The chemical reaction conditions are realized by adding individual ions at a given reaction rate or replacing them with different kinds of ions,
The calculation means is a plasma analysis means and ions generated by a chemical reaction between ions irradiated on the plasma mass and ions outside the plasma mass in response to changes in the initial conditions, boundary conditions and chemical reaction conditions of the plasma , A plasma analyzer that analyzes the density and velocity of the electron fluid that interacts with it and the time variation and spatial distribution of the electromagnetic field .   2. The plasma analyzing apparatus according to claim 1, wherein a physical phenomenon inside and outside the plasma mass is analyzed as the chemical reaction occurring in a region outside the range of a certain distance from the boundary of the plasma mass.   The chemical reaction is a photoionization reaction in which charged particles are generated by irradiating light to non-ionized particles and / or a charge exchange reaction with charge exchange between different kinds of particles. Item 3. The plasma analysis apparatus according to Item 1 or 2.   It is assumed that the solar wind is an ion and an electron fluid incident on the plasma mass from outside, and the space plasma is analyzed under the condition of the chemical reaction between the solar wind and the ion around the plasma mass or the plasma mass. The plasma analysis apparatus according to claim 1, 2 or 3. 5. The plasma analyzing apparatus according to claim 3, wherein the photoionization reaction under the chemical reaction conditions is O + hv → O ⁺ + e ⁻ and the charge exchange is H ⁺ + O → H + O ⁺ . In a method for analyzing plasma by means of computation,
A hybrid code that considers ions as particles and electrons as fluids, and equipped with plasma analysis means for analyzing plasma using the equations of motion of individual ions, the equations of motion of electron fluids, and Maxwell's equations. And plasma analysis of the system irradiated with electron fluid,
Maintains initial plasma conditions, boundary conditions, and chemical reaction conditions.
The chemical reaction conditions are realized by adding individual ions at a given reaction rate or replacing them with different kinds of ions,
In response to changes in the plasma analysis means and the initial conditions, boundary conditions, and chemical reaction conditions of the plasma, ions generated by a chemical reaction between ions irradiated on the plasma mass and ions outside the plasma mass, and mutual A plasma analysis method characterized by analyzing the density and velocity of the acting electron fluid, and the time variation and spatial distribution of the electromagnetic field .   7. The plasma analysis method according to claim 6, wherein a physical phenomenon inside and outside the plasma mass is analyzed on the assumption that the chemical reaction occurs in a region outside a range of a certain distance from the boundary of the plasma mass.   The chemical reaction is a photoionization reaction in which charged particles are generated by irradiating light to non-ionized particles and / or a charge exchange reaction involving exchange of charges between different kinds of particles. 6. The plasma analysis method according to 6 or 7.   It is assumed that the solar wind is an ion and an electron fluid incident on the plasma mass from outside, and the space plasma is analyzed under the condition of the chemical reaction between the solar wind and the ion around the plasma mass or the plasma mass. The plasma analysis method according to claim 6, 7 or 8. 10. The plasma analysis method according to claim 8, wherein the photoionization reaction under the chemical reaction conditions is O + hv → O ⁺ + e ⁻ and the charge exchange is H ⁺ + O → H + O ⁺ .

Images