Particle-PLUS コラム - プラズマモデリング（連載） - 第７回
ウェーブフロントでは「粒子法」でプラズマをモデル化したシミュレーションソフトウェア『Particle-PLUS』や 「流体法」のプラズマシミュレーションソフトウェア『VizGlow』を取り扱っており， これらを適宜使い分けながら皆様方へ長年ソリューションを提供しております． 本コラムを読んでプラズマシミュレーションに少しでも興味をお持ちになりましたら 資料請求や 無料セミナーなど承りますので， いつでもお気軽にお問い合わせください．
プラズマの流体方程式の応力および熱流束を決定するための熱力学的な近似モデルを3種類紹介しました． 今回は，異なる視点で「ドリフト拡散近似」（あるいは「ドリフト拡散モデル」）と呼ばれる定式化の方法を紹介します． ドリフト拡散近似はプラズマを流体法でシミュレーションする場合に最もよく使われる近似ですので，ぜひとも本コラムが理解の一助となりますと幸いです．
2.4. ドリフト拡散近似Drift-diffusion approximation
2.4.1. ドリフト拡散方程式
荷電粒子の平均速度は電場に対して即座に応答して，その慣性が無視できると考えるドリフト拡散近似と呼ばれるモデル化手法があります．
ドリフト拡散近似では運動量方程式が後述のドリフト拡散方程式に置き換わります．
ドリフト拡散近似では，等方性応力に加えて外部磁場が無視できる状況を考えます*1．
さらに，運動量交換項 $\vec{K}_{s}$ について，運動量移行衝突周波数 ${\nu_{\mathrm{m}}}_{ss^{\prime}}=1/{\tau_{\mathrm{m}}}_{ss^{\prime}}$
（${\tau_{\mathrm{m}}}_{ss^{\prime}}$ は運動量移行衝突の平均自由時間）を導入して，
\begin{align}
\vec{K}_{s} = \sum_{s^{\prime}\neq s} \frac{\tilde{m}_{ss^{\prime}}n_{s}\left(\vec{u}_{s^{\prime}}-\vec{u}_{s}\right)}{{\tau_{\mathrm{m}}}_{ss^{\prime}}}
= \sum_{s^{\prime}\neq s} \tilde{m}_{ss^{\prime}}n_{s}\left(\vec{u}_{s^{\prime}}-\vec{u}_{s}\right){\nu_{\mathrm{m}}}_{ss^{\prime}}
\label{eq_definition_K}
\end{align}
と書き直します．
ここで，
\begin{align}
\tilde{m}_{ss^{\prime}}\equiv \frac{m_{s}m_{s^{\prime}}}{m_{s}+m_{s^{\prime}}}
\label{eq_reduced_mass}
\end{align}
は換算質量です．
なお，この表式(\ref{eq_definition_K})については後に第3.2節で議論しますので，ここでは一種のansatzだと考えてください．
さて，以上の仮定と連続方程式
\begin{align*}
\frac{\partial n_{s}}{\partial t}+\nabla\cdot\vec{\varGamma}_{s}=S_{s}
\end{align*}
を用いて運動量方程式
\begin{align*}
\frac{\partial}{\partial t}(\rho_{s}\vec{u}_{s})
+ \nabla \cdot \left\{ \rho_{s}\left(\vec{u}_{s}\otimes\vec{u}_{s}\right) \right\}
= q_{s} n_{s} \left(\vec{E}+\vec{u}_{s}\times\vec{B}\right)
- \nabla p_{s}
+ \vec{K}_{s}
\end{align*}
を変形すると，
\begin{align*}
m_{s}n_{s}\left(\frac{\partial \vec{u}_{s}}{\partial t}\right) + m_{s}S_{s}\vec{u}_{s}
+ m_{s}n_{s}\vec{u}_{s}\cdot\left(\nabla\otimes\vec{u}_{s}\right)
= q_{s}n_{s}\vec{E} - \nabla\left(n_{s}k_{\mathrm{B}} T_{s}\right)
+ \sum_{s^{\prime}\neq s} \tilde{m}_{ss^{\prime}}n_{s}\left(\vec{u}_{s^{\prime}}-\vec{u}_{s}\right){\nu_{\mathrm{m}}}_{ss^{\prime}}
\end{align*}
となります．
左辺第1項は平均運動量の時間変化を，左辺第2項は粒子の生成・消滅に応じた運動量変化を表していますが，
これらの変化の時間スケール*2が運動量移行（右辺第3項）の時間スケール ${\forall\tau_{\mathrm{m}}}_{ss^{\prime}}$ と比較して十分遅い場合
（即ち，熱化の時間スケールからするとほとんど定常流と見做せて，さらに電離や結合を伴う衝突よりも熱化のみが起きるような衝突が支配的である場合）には，
左辺第1,2項は無視できます．
また，熱速度 $\vec{w}_{s}(\simeq \sqrt{k_{\mathrm{B}} T_{s}/m_{s}})$ と比較して流速 $\vec{u}_{s}$ が十分小さければ，
右辺第2項と比べて左辺第3項は無視できるようになります．
以上の仮定が成り立つ場合，運動量方程式は以下のドリフト拡散方程式に置き換えられます．
ここで，$\mu_{s}$ と $D_{s}$ はそれぞれ移動度と拡散係数と呼ばれる量で， まとめて輸送パラメータ（輸送係数やスウォームパラメータ (swarm parameter) とも）と呼ぶこともあります。 移動度と拡散係数は， \begin{align} \mu_{s} \equiv \left(\sum_{s^{\prime}\neq s}\tilde{m}_{ss^{\prime}}{\nu_{\mathrm{m}}}_{ss^{\prime}}\right)^{-1}q_{s}, \qquad D_{s} \equiv \left(\sum_{s^{\prime}\neq s}\tilde{m}_{ss^{\prime}}{\nu_{\mathrm{m}}}_{ss^{\prime}}\right)^{-1}k_{\mathrm{B}} T_{s} \end{align} のように運動量移行衝突周波数 ${\nu_{\mathrm{m}}}_{ss^{\prime}}$ によって決定することができます． これらの輸送パラメータについては，拡散項と共に次回以降のコラムで説明する予定です．
以上より，ドリフト拡散近似において解くべき方程式が以下であることが分かります．
ドリフト拡散近似では，熱速度 $\vec{w}_{s}(\simeq \sqrt{ k_{\mathrm{B}} T_{s}/m_{s}})$ と比較して流速 $\vec{u}_{s}$ が十分小さいと仮定していました． 即ち，$\varepsilon_{s}\ll k_{\mathrm{B}} T_{s}$ となりますので，$\varepsilon_{s}$ の項を無視しています． また，熱流束 $\vec{Q}_{s}$ は座標に依存しないスカラー熱伝導率 $\kappa_{s}$ を用いて \begin{align} \vec{Q}_{s} = -\kappa_{s}\nabla T_{s} = -\frac{5}{2} k_{\mathrm{B}} D_{s} n_{s} \nabla T_{s}\qquad (\kappa_{s}=\frac{5}{2} k_{\mathrm{B}} D_{s}n_{s}). \label{eq_relation_kappa_and_D} \end{align} のように書けるとしました*3．
ドリフト拡散モデルの本質的な変数は $n_{s}$，$\vec{\varGamma}_{s}$，$T_{s}$ であり，電場および散乱項を含む輸送パラメータを除けば方程式系が閉じています．
*1 磁場があると，後述の移動度や拡散係数が非等方的になることもあり問題が複雑になるため，近似のメリットが薄れてしまいます．
*2 平均運動量変化の時間スケールとしては，例えば印加電圧の周期を想定します．粒子増減の時間スケールは，例えば電離衝突の平均自由時間です．
*3 熱伝導率 $\kappa$ と拡散係数 $D$ の関係式 (\ref{eq_relation_kappa_and_D}) について，拡散係数の定義 \begin{align} D=\frac{\kappa}{\rho c_{\mathrm{p}}} \label{eq_definition_D} \end{align} から確認しておきましょう． 粒子種を表す下付き添え字の $s$ は省略してあります． ここで，$c_{\mathrm{p}}$ は定圧比熱（定圧モル比熱ではないことに注意）です． さらに，定常状態における熱伝導率 $\kappa$ はFourierの法則 $\vec{Q}=-\kappa\nabla T$ で定義されるスカラー定数としておきます． さて，自由度 $d$ の理想気体の定圧モル比熱は $\tilde{c}_{\mathrm{p}}=\{(d+2)/2\}N_{\mathrm{A}} k_{\mathrm{B}}$ （$N_{\mathrm{A}}$ はAvogadro数）です． （この関係式の導出は熱力学の典型的な演習問題ですので，本コラムでは導出は省略します．） また，モル質量 $\tilde{m}=mN_{\mathrm{A}}$ を使うと，定圧比熱 $c_{\mathrm{p}}$ と定圧モル比熱 $\tilde{c}_{\mathrm{p}}$ の関係は $c_{p}=\tilde{c}_{\mathrm{p}}/\tilde{m}$ と書けます． したがって，拡散係数の定義 (\ref{eq_definition_D}) から逆算すると， \begin{align} \kappa = \rho c_{\mathrm{p}} D = \frac{d+2}{2} k_{\mathrm{B}} D n \end{align} であることが分かります． なお，第2.1.1項で議論したように， ここでは並進自由度を考えれば十分なので $d=3$ とします．
2.4.2. ドリフト拡散方程式の意味と有効範囲
ドリフト拡散方程式 (\ref{eq_drift_diffusion}) の右辺第1項は電場によるドリフト速度を，
右辺第2項は圧力（あるいは密度と温度）勾配による拡散速度を，
右辺第3項は背景流れによる速度修正を表しています．
ドリフト拡散近似を用いるための時間スケールに関する仮定
（(平均運動量の時間変化及び粒子増減による運動量変化の時間スケール)$\,\ll{\tau_{\mathrm{m}}}_{ss^{\prime}}$）
は，経験的にガス圧力がおよそ$10$–$100\ \mathrm{mTorr}(=1.33$–$13.3\ \mathrm{Pa})$以上
の場合に成り立ちます．（もちろん，印可電圧の周波数やガス種によっても異なるため，あくまでも目安です．）
熱速度が流速よりも十分大きくて対流が無視できるという過程に関しては，プロセスプラズマであればおよそ成り立ちます．
また，第2.1節でも述べたように，等方性を仮定するためには電場が十分小さくなければなりません．
さらに，背景ガス密度よりも荷電粒子密度が非常に小さい弱電離プラズマであることや，
荷電粒子と背景ガスとの間の平均自由行程が系の特徴長さよりも十分短いような衝突性プラズマであることも必要です．
したがって，ドリフト拡散近似はプラズマシミュレーションにおいて実用上重要な近似ではありますが，
高圧であってもマイクロ放電のような現象では近似が不適になるため，適用範囲には気を付けなければならないことが分かります．
第７回まとめ
第７回ではドリフト拡散近似を紹介しました． 前回までの熱力学的な近似も含めて，プラズマの流体方程式系中に現れる余分な未知数である応力と熱流束を決定するためのモデル化をいくつか紹介したことになります． 残りの余分な未知数は電場と衝突パラメータですが，次回からは衝突パラメータについて議論していきたいと思います．
ウェーブフロントでは「粒子法」でプラズマをモデル化したシミュレーションソフトウェア『Particle-PLUS』や 「流体法」のプラズマシミュレーションソフトウェア『VizGlow』を取り扱っており， これらを適宜使い分けながら皆様方へ長年ソリューションを提供しております． 本コラムを読んでプラズマシミュレーションに少しでも興味をお持ちになりましたら 資料請求や 無料セミナーなど承りますので， いつでもお気軽にお問い合わせください．
 |
上野 崚一郎 | 博士（理学） 1991年 広島県生まれ 2019年 ウェーブフロント入社 2019年 広島大学大学院理学研究科 博士後期課程修了 学生時代は数値シミュレーションを使った素粒子論（格子ゲージ理論）の研究に従事． 入社後は，専門職（エンジニア）としてプラズマ解析ソフトウェアの開発をはじめとして， 解析コンサルティング業務や国内外のユーザー向けの技術サポート・トレーニングなどを担当． 最後までお読みいただきありがとうございます． お気づきの点や扱ってほしい話題がございましたらお気軽にお問い合わせください． |