JP6157941B2

JP6157941B2 - プラズマシミュレーション方法及びプラズマシミュレーションプログラム

Info

Publication number: JP6157941B2
Application number: JP2013124691A
Authority: JP
Inventors: 大宮下
Original assignee: Sumitomo Heavy Industries Ltd
Current assignee: Sumitomo Heavy Industries Ltd
Priority date: 2013-06-13
Filing date: 2013-06-13
Publication date: 2017-07-05
Anticipated expiration: 2033-06-13
Also published as: JP2015001769A

Description

本発明は、粒子モデル及び流体モデルを用いたプラズマシミュレーション方法及びプラズマシミュレーションプログラムに関する。

従来のプラズマシミュレーション方法としては、例えば下記非特許文献１に記載されたものが知られている。このプラズマシミュレーション方法では、粒子モデル及び流体モデルを用いて計算を行う、いわゆるハイブリッドＰＩＣ（Particle In Cell）法を利用し、プラズマシミュレーションを行うことが図られている。

Masaru MIYASHITA etal.,"Plasma Simulation of Reactive Plasma Deposition Equipment Using ThreeDimensional Hybrid PIC/MC method", Frontier of Applied Plasma TechnologyVol.5 No.2 July 2012

ここで、上記プラズマシミュレーション方法では、流体モデルを用いて計算を行う際、自己誘導磁場の計算時において、無限遠の境界条件を再現するために計算領域を大きくする場合がある。しかし、このように計算領域を大きくすると、例えば膨大な規模の格子について電子流体方程式を取り扱うことになることから、計算が冗長となり計算時間かかり過ぎる虞がある。

本発明は上記実情に鑑みてなされてものであり、計算コストを低減することができるプラズマシミュレーション方法及びプラズマシミュレーションプログラムを提供することを課題とする。

上記課題を解決するために、本発明に係るプラズマシミュレーション方法は、粒子モデル及び流体モデルを用いたプラズマシミュレーション方法であって、粒子モデルを用いて、イオン及び中性粒子に関する計算を行う粒子モデル計算ステップと、流体モデルを用いて、電子に関する計算を行う流体モデル計算ステップと、を備え、流体モデル計算ステップは、電子流体方程式を用いて電子温度及び電位を求める第１ステップと、自己誘導磁場を求める第２ステップと、移動度を更新する第３ステップと、を有し、第２ステップにおける計算領域は、第１及び第３ステップにおける計算領域よりも大きい。

このプラズマシミュレーション方法では、電子温度及び電位を求める場合及び移動度を更新する場合には、例えば所定領域を計算領域として電子流体方程式を計算し、自己誘導磁場を求める場合には、当該所定領域よりも大きい領域を計算領域として計算することが可能となる。すなわち、自己誘導磁場の計算領域を大きくして計算精度を維持しつつ、当該計算領域よりも小さい計算領域にて電子温度及び電位の算出並びに移動度の更新を行い、その計算を簡潔化することができる。従って、本発明によれば、計算コストを低減することが可能となる。

また、上記作用効果を好適に奏する場合として、具体的には、第１及び第３ステップにおける計算領域は、シミュレートの対象となる所定領域とされている場合がある。

また、第２ステップにおける計算領域は、第１又は第３ステップにおける計算領域を、各座標軸において等距離拡大するように設定されてなることが好ましい。この場合、シミュレーションの精度を高めることが可能となる。

また、第２ステップにおける計算領域は、第１又は第３ステップにおける計算領域と重複する領域と、第１又は第３ステップにおける計算領域と重複しない領域と、を含み、重複しない領域の計算点の密度は、重複する領域の計算点の密度よりも小さいことが好ましい。このように、重複しない領域の計算点の密度を重複する領域の計算点の密度よりも小さくすることにより、シミュレーションの精度を維持しながら、計算コストを一層低減することができる。

また、本発明に係るプラズマシミュレーションプログラムは、上記プラズマシミュレーション方法をコンピュータに実行させる。このプラズマシミュレーションプログラムにおいても、上記作用効果、すなわち、計算コストを低減することが可能となるという作用効果が奏される。

本発明によれば、計算コストを低減することが可能となる。

一実施形態に係るプラズマシミュレーション方法の対象となるＲＰＤ装置を示す概略図である。一実施形態に係るプラズマシミュレーションプログラムの構成を示す概略図である。一実施形態に係るプラズマシミュレーション方法を示すフローチャートである。図３のプラズマシミュレーション方法における自己誘導磁場計算を説明するための図である。

以下、添付図面を参照して、本発明の好適な実施形態について詳細に説明する。なお、以下の説明において同一又は相当要素には同一符号を付し、重複する説明を省略する。

図１は、一実施形態に係るプラズマシミュレーション方法の対象となるＲＰＤ装置を示す概略図である。図１に示すように、本実施形態に係るプラズマシミュレーション方法及びプラズマシミュレーションプログラムは、一例として、ＲＰＤ（Reactive Plasma Deposition）装置１をシミュレートの対象としている。そこで、まず、このＲＰＤ装置１について説明する。

ＲＰＤ装置１は、例えば、材料蒸発源Ｍａを低基板温度条件下（〜２００℃）で成膜用基板２に成膜するイオンプレーティング装置である。材料蒸発源Ｍａとしては、ＩＴＯやＺｎＯ等の透明導電材料や、ＳｉＯＮ等の絶縁封止材料が例示され、成膜用基板２としては、例えばガラス基板やプラスチック基板等の板状部材が例示される。

このＲＰＤ装置１では、真空チャンバ３において、プラズマ源４から発生したプラズマＰが、ステアリングコイル５とハースコイル６ａ及びハース磁石６ｂから構成されるプラズマビームコントローラ６とで形成される磁場によって、材料蒸発源Ｍａへ導かれる。プラズマＰ中を流れる電流が例えば１５０Ａと大きいことから、材料蒸発源Ｍａが熱蒸発され、一部がプラズマＰと相互作用を起こしてイオン化・励起される。そして、イオンはプラズマＰ中の電位勾配によって加速され、成膜粒子Ｍｂとして広がり、成膜用基板２にて成膜が行われる。

特に、ここでのＲＰＤ装置１は、プラズマＰ中を流れる電流が大きいことから、当該電流により生じる自己誘導磁場の影響が強いものとされており、この自己誘導磁場によってプラズマＰが中心軸からズレ、ねじれが発生する場合があるものとされている。

図２は、一実施形態に係るプラズマシミュレーションプログラムの構成を示す概略図である。図２に示すように、プラズマシミュレーションプログラム１０は、プラズマシミュレーション方法をコンピュータに実行させるためのプログラムであって、コンピュータが備える記録媒体５０に形成されたプログラム格納領域５０ａ内に格納されている。

このプラズマシミュレーションプログラム１０は、初期条件の入力処理を行う初期条件入力モジュール１１と、粒子モデルを用いてイオン及び中性粒子に関する計算処理を行う粒子モデル計算モジュール１２と、流体モデルを用いて電子に関する計算処理を行う流体モデル計算モジュール１３と、収束判定処理を行う収束判定モジュール１４と、を備えている。

粒子モデル計算モジュール１２は、運動方程式の計算処理を行う運動方程式計算部１２ａと、粒子境界条件の入力処理を行う粒子境界条件入力部１２ｂと、衝突計算処理を行う衝突計算部１２ｃと、を含んでいる。流体モデル計算モジュール１３は、電子流体方程式の計算処理を行う電子流体方程式計算部１３ａと、自己誘導磁場の計算処理を行う自己誘導磁場計算部１３ｂと、移動度の更新処理を行う移動度更新部１３ｃと、を含んでいる。なお、上記の各処理については、後述のプラズマシミュレーション方法の説明において詳説する。

図３は、一実施形態に係るプラズマシミュレーション方法を示すフローチャートである。本実施形態のプラズマシミュレーション方法は、粒子モデル及び流体モデルを用いて計算を行う、いわゆるハイブリッドＰＩＣ（Particle In Cell）法を利用するものである。

すなわち、本実施形態のプラズマシミュレーション方法は、速度が大きく（〜１０^６ｍ／ｓ）且つ粒子モデルでは計算コストのかかる電子について流体モデルで取り扱う一方で、速度分布を求めたいイオンや中性粒子については粒子モデルで取り扱う、電気的準中性を仮定している。

なお、電子のエネルギ分布は、１温度のマックスウェル・ボルツマン分布に近似できることが見出されることからも、流体モデルを適用することは好ましい。また、本実施形態では、３次元的にプラズマシミュレーションを実施しており、シミュレートの対象となる所定領域（求めたい領域）として真空チャンバ３内を設定している。

図３に示すように、本実施形態のプラズマシミュレーション方法では、まず、初期条件を入力する（Ｓ１）。初期条件としては、例えば、イオン及び中性粒子の初期位置及び初期速度、プラズマＰ中の初期電位、初期電子密度分布、初期磁場分布、及び、流体モデル計算における計算領域を少なくとも含んでいる。各初期値は、任意の値とすることもできるし、例えば初期磁場分布は、プラズマビームコントローラ６による磁場分布とすることもできる。

次に、粒子モデルを用いて、イオン及び中性粒子に関する計算を行う粒子モデル計算パート（粒子モデル計算ステップ）を実施する。粒子モデル計算パートでは、運動後の粒子がセルを跨がず且つ１つの粒子が１時間刻みの間に２回衝突しないように、十分に時間刻みを短く設定し、粒子の運動と衝突現象とを分離して計算させている。なお、粒子の所属するセルを簡単に計算するために、ｘｙｚ軸が直交するカルテシアン座標系を用いて数値計算を行っている（以下の数値計算において同じ）。

具体的には、下式（１）を用いて、イオン及び中性粒子の速度及び位置を計算する（Ｓ２）。下式（１）の時間積分は、時間刻みを大きくとった際にも軌道計算の誤差が小さくなるように、高次精度の古典的４次のルンゲ・クッタ法を用いている。外力としては、電場及び磁場のローレンツ力を考慮し、重力は無視している。

ｍ_α：粒子αの質量、Ｖ_α：粒子αの速度、ｑ：電化、φ：電位、
Ｂ（＝Ｂ_self＋Ｂ_ex）：磁束密度、Ｂ_self：コイル等による磁束密度、
Ｂ_ex：プラズマＰ中の電流による磁束密度。

ここでは、粒子の所属するセルの８隅にある格子点に電場・磁束密度の値を持たせ、粒子の位置へ線形補間を行っている。また、補間方法に合わせて粒子が数値計算のための力を受けないようにするために、１次の形状関数を与えて数密度分布や流速分布等の統計処理を行っている。

続いて、以下の粒子境界条件を与える（Ｓ３）。流入粒子数は、実際の蒸発レートからの推測値を用いている。なお、粒子の種類や衝突の種類は、支配的な現象及び再現したい現象を絞って選定することができる。
アルゴン，酸素分子：鏡面反射
その他の粒子：消失

続いて、粒子の衝突計算を実施する（Ｓ４）。ここでの衝突計算では、以下に示すように、少なくとも、電子−中性粒子衝突モデルによる計算と、重粒子同士衝突モデルによる計算と、を行っている。

［電子−中性粒子衝突モデル］
電気的準中性の仮定からイオン密度と電子密度は等しく、また、後述の電子流体方程式の計算結果より、あるセルにおける電子温度が求まり反応確率がわかることから、反応確率に応じてイオン化等の中性粒子の反応処理を行う。

［重粒子同士衝突モデル］
重粒子同士の衝突を計算し、運動量輸送を考慮する。本実施形態では、演算の簡略化のために剛体球衝突を仮定する。このように、重粒子同士の衝突を考慮することで、イオンが材料蒸発源Ｍａ付近の背景アルゴンガスを排除する現象を再現できる。

次に、流体モデルを用いて、電子に関する計算を行う流体モデル計算パート（流体モデル計算ステップ）を実施する。この流体モデル計算パートでは、まず、電子流体方程式を用いて、電子温度及び電位を求める（Ｓ５：第１ステップ）。具体的には、以下に示すように、運動量保存則、質量保存則及びエネルギ保存則から、電子温度及び電位を算出する。

［運動量保存則］
ＲＰＤ装置１中のプラズマＰの振動周期（〜１０^−１２ｓｅｃ）は、イオンや中性粒子の運動を計算する時間刻み（〜１０^−６ｓｅｃ）では、電子の運動が定常状態に達していると仮定できる。定常状態を仮定すると、運動量保存則から電子流束Γ_ｅを下式（２）で表すことができる。磁束密度の方向に対して複数の拡散・移流現象が考えられるが、下式（２）では、電子の流れとして、磁力線に対して平行方向及び垂直方向ともに古典拡散(μ_||)を仮定した。但し、磁力線を横切る方向（μ_⊥）に関しては、衝突周波数とラーマー周波数との比で決まる補正係数で補正している。なお、[μ_ｅ]は、後述の移動度テンソルである。

ｎ_ｅ：電子密度、Ｔ_ｅ：電子温度。

［質量保存則］
上記の運動方程式から導かれる電子流束の保存則を下式(３）式に示す。ここでは、電位φに対して線形な方程式として、有限体積法を用いて離散化して解いている。

γ_ion：単位時間・体積当たりのイオン化による電子の発生数密度。

境界条件としては、下式（４）に示すように、熱流束に壁とプラズマ電位とのボルツマンファクターを乗ずることでシースの影響を考慮する。

α：壁面上の法線ベクトルと磁束密度の成す角、
△φ：プラズマ端と壁面との電位差。

［エネルギ保存則］
エネルギ保存則の式は、下式(５)で表される。下式（５）では、電子温度Ｔ_ｅに対して線形に近似して解いており、また、一項目の対流項と二項目の拡散項との比率に応じて、風上差分及び中心差分を使い分けるハイブリッド法により離散化を評価している。なお、下式（６）に示すように、エネルギ流束の境界条件としては、流束同様に熱流束を仮定している。

続いて、自己誘導磁場を算出する（Ｓ６：第２ステップ）。自己誘導磁場の計算では、イオン及び電子の運動が定常状態に達しているとの仮定から、磁場も定常状態に達していると考えられるため、プラズマＰ中の電流が作るベクトルポテンシャル(下式（７）)の回転として自己誘導磁場Ｂ_selfを計算する（下式(８)）。なお、境界条件は、Ａ_ｘ＝Ａ_ｙ＝Ａ_ｚ＝０としている。

μ_０：真空中の誘電率。

このとき、電子流体計算全体を通して、計算を安定化させて行列の優対角性を保証するために、対角項と同符号の項を湧出し項へ移項して１反復毎に係数を修正すると共に、不足緩和係数をかけて物理量を更新する。また、各物理量に対して非現実的な解を除去するために、各物理量に最大値及び最小値を設けて制限する。なお、計算を安定化させるために、上述の安定化手法に限らず、他の安定化手法を用いてもよい。

続いて、移動度を更新する（Ｓ７：第３ステップ）。下式（９）に示すように、移動度テンソル［μ_ｅ］は、カルテシアン座標系において、磁束密度ベクトルの座標系から座標変換Λで変換される。

次に、計算結果の場における値の変動が十分に（所定値よりも）小さく、収束しているか否かを判定する（Ｓ８）。上記Ｓ８にてＮｏの場合、上記Ｓ２に戻り、イオン及び中性粒子の位置及び速度が更新されると共に、電子温度及び電位が更新される。一方、上記Ｓ８にてＹｅｓの場合、シミュレートが終了し、結果が出力される。当該出力結果には、例えば、シミュレートの対象となる所定領域における、イオン及び中性粒子の位置分布、速度分布、電位分布、電子密度分布、及び磁場分布が少なくとも含まれている。

ここで、流体モデル計算パート（上記Ｓ５〜上記Ｓ７）において、自己誘導磁場の計算時の境界条件として無限遠を再現するために、計算領域を大きくする場合がある。しかしこの場合、例えば膨大な規模の格子について電子流体方程式を取り扱うことになるために計算が冗長となり、計算時間かかり過ぎるという虞がある。

そこで、図４に示すように、本実施形態では、自己誘導磁場の演算（上記Ｓ６）における計算領域Ｒ１を、電子流体方程式の演算（上記Ｓ５）及び移動度更新（上記Ｓ７）における計算領域Ｒ２よりも大きく設定している。つまり、上記Ｓ５，Ｓ７の計算領域Ｒ２を大きくせずに、上記Ｓ６の計算領域Ｒ１については、無限遠の境界条件を再現するべく大きい領域として計算点を生成している。

具体的には、計算領域Ｒ２は、シミュレートの対象となる所定領域であって、求めたい計算領域である真空チャンバ３内に対応している。一方で、計算領域Ｒ１については、無限遠を再現するために好ましいとして、計算領域Ｒ２に対し各座標軸において所定倍以上拡大するように設定されている。

所定倍としては、例えば計算結果の誤差が許容範囲（１０％等）に収まる値とされており、好ましいとして５倍以上とされている。図示するように、計算領域Ｒ１，Ｒ２を立方体としたとき、ここでの計算領域Ｒ１は、各座標軸の正方向及び負方向それぞれにおいて、計算領域Ｒ２よりも５Ｌ大きくなるように設定されている。
Ｌ：計算領域Ｒ２の一辺の長さ。

また、計算領域Ｒ１は、計算領域Ｒ２を含むように設定されており、計算領域Ｒ２に対し重複する重複領域を中央部分に有していると共に、計算領域Ｒ２よりも外側は重複しない非重複領域となっている。計算領域Ｒ１において外側の非重複領域では、その計算点密度（単位体積当たりの計算点の数）が重複領域の計算点密度よりも小さくされており、重複領域の格子よりも粗い格子となるように設定されている。

以上のように、本実施形態において、上記Ｓ５，Ｓ７では、シミュレートしたい領域を計算領域Ｒ２とする一方、自己誘導磁場を求める上記Ｓ６では、当該計算領域Ｒ２よりも大きい領域を計算領域Ｒ１として計算することが可能となる。よって、自己誘導磁場の計算領域Ｒ１を、境界条件として無限遠を再現するように大きくして計算精度を維持しつつ、当該計算領域Ｒ１よりも小さい計算領域Ｒ２にて電子温度及び電位の算出並びに移動度の更新を行い、その計算を簡潔化できる。従って、計算コストを低減することが可能となる。

すなわち、上記Ｓ５〜上記Ｓ７のそれぞれは、その計算成分に下式（１０）をｘ軸、ｙ軸及びｚ軸に関して有している。そのため、計算点Ｎの数に計算コストが大きく依存するところ、上記Ｓ６ではメッシュが大きい一方で上記Ｓ５，Ｓ７ではメッシュが必要最低限で小さく維持されていることから、上記Ｓ６のみで計算点Ｎの数が大きくなることとなる。従って、計算コストを低減でき、例えば計算時間を１／１００まで短縮化することが可能となる。

Ｓ：行列、
ａ，ｂ：ベクトル。

また、本実施形態の計算領域Ｒ１は、上述したように、計算領域Ｒ２を各座標軸において等距離拡大するように設定されている。これにより、例えば、計算領域Ｒ２を偏重して拡大するように計算領域Ｒ１が設定される場合に比べ、シミュレーションの精度を高めることが可能となる。

また、本実施形態では、上述したように、計算領域Ｒ１において、計算領域Ｒ２と重複しない領域の計算点の密度が、重複する領域の計算点の密度よりも小さくされている。これにより、シミュレーションの精度を維持しながら、計算コストを一層低減することができる。

以上、本発明の好適な実施形態について説明したが、本発明は上記実施形態に限られるものではなく、各請求項に記載した要旨を変更しない範囲で変形し、又は他のものに適用したものであってもよい。

例えば、本発明は、ＲＰＤ装置１のシミュレーションに適用できるだけでなく、核融合のシミュレーションやイオンエンジンのシミュレーション等にも適用することができる。特に、本発明は、プラズマＰ中を流れる電流が大きく当該電流により生じる自己誘導磁場の影響が強い場合のシミュレーションに好適に適用することができる。

また、上記計算領域Ｒ１，Ｒ２の形状は限定されるものではなく、例えば球形状であってもよい。なお、本発明は、上記プラズマシミュレーション方法を実施するプラズマシミュレータとして捉えることもでき、また、上記プラズマシミュレーションプログラム１０を記憶する記憶媒体として捉えることもできる。

１０…プラズマシミュレーションプログラム、１１…初期条件入力モジュール、１２…粒子モデル計算モジュール、１２ａ…運動方程式計算部、１２ｂ…粒子境界条件入力部、１２ｃ…衝突計算部１２ｃ、１３…流体モデル計算モジュール、１３ａ…電子流体方程式計算部、１３ｂ…自己誘導磁場計算部、１３ｃ…移動度更新部、１４…収束判定モジュール。

Claims

粒子モデル及び流体モデルを用いたプラズマシミュレーション方法であって、
コンピュータが、前記粒子モデルを用いて、イオン及び中性粒子に関する計算を行う粒子モデル計算ステップと、
コンピュータが、前記流体モデルを用いて、電子に関する計算を行う流体モデル計算ステップと、を備え、
前記流体モデル計算ステップは、
電子流体方程式を用いて電子温度及び電位を求める第１ステップと、
自己誘導磁場を求める第２ステップと、
移動度を更新する第３ステップと、を有し、
前記第２ステップにおける計算領域は、前記第１及び第３ステップにおける計算領域よりも大きい、プラズマシミュレーション方法。
前記第１及び第３ステップにおける計算領域は、シミュレートの対象となる所定領域とされている、請求項１記載のプラズマシミュレーション方法。
前記第２ステップにおける計算領域は、前記第１又は第３ステップにおける計算領域を、各座標軸において等距離拡大するように設定されてなる、請求項１又は２記載のプラズマシミュレーション方法。
前記第２ステップにおける計算領域は、前記第１又は第３ステップにおける計算領域と重複する領域と、前記第１又は第３ステップにおける計算領域と重複しない領域と、を含み、
前記重複しない領域の計算点の密度は、前記重複する領域の計算点の密度よりも小さい、請求項１〜３の何れか一項記載のプラズマシミュレーション方法。
請求項１〜４の何れか一項記載のプラズマシミュレーション方法をコンピュータに実行させるプラズマシミュレーションプログラム。