計測エンジニアリングシステム株式会社

CASES AND MATERIALS 事例/資料

キーワード・条件で検索

キーワードで検索

条件で検索

COMSOL解析事例

COMSOL紹介( 導入/検討 )

計測

サブスク限定

ボルツマン方程式に基づくICPモデル

ボルツマン方程式とプラズマ流体モデルの連成計算手法

低温プラズマは弱電離でかつ熱的に非平衡であり、それゆえに利用価値が高く、様々な産業に幅広く応用されている。プラズマ計算用の流体モデルは、プラズマを連続体(流体)とみなし、質量保存則を表す「連続の式」、運動量保存則を表す「運動方程式」、そしてエネルギー保存則を表す「エネルギー方程式」に基づいてシミュレーションを行う。また、運動方程式を解くのは、膨大な計算時間を要するため、かわりにドリフト-拡散モデル、輸送方程式を計算する。
流体モデルによる低温プラズマシミュレーションを行うにあたっては、電子やイオンの輸送係数および電離などの衝突反応を表す衝突レート係数を事前に用意する必要である。これら係数のうち、特に電子に関する各種係数では、ボルツマン方程式 (Boltzmann equation: BE) 解析やモンテカルロシミュレーションなどの方法で様々な換算電界下での電子輸送係数や衝突レート係数を事前に求めておき、これら値をデータテーブル化して流体モデルによるプラズマシミュレーションに使用する方法が一般的である1)。しかし、プラズマ放電においては、プラズマ密度などが時間によって変わる。事前に電子輸送係数や衝突レート係数を求める手法で得たデータは、放電中のプラズマと一致していない場合がある。
そこで本例題は図1に示すGEC型ICP (Inductively Coupled Plasma: ICP) 装置2)のプラズマモデルに基づいてボルツマン方程式が流体モデルと時間に依存する連成計算を行った。
本書ではモデルの作成および、2章以後に図2と図3の解析結果を得る手順を示す。

  図1.解析モデル

計算モデルは、電子、イオンと中性粒子の輸送、空間電荷場およびコイルによる誘導放電の磁場の計算方程式を連立して計算を行うが、紙数の都合でここではこれらの方程式の説明3)を省略する。電子輸送方程式におけるソース項には電子と分子や原子の衝突によって電子密度と電子エネルギーを変化することが含まれる。反応\(j\)の反応速度係数\(k_j\)は次の式で示される。

   \(k_j\)=\(γ\)\(\int ^{\infty }_{0}\)\(εσ_j(ε)f(ε)dε\),\(γ\)=(\(2q/m\))

ここで、\(ε\)は電子エネルギー、\(σ_j(ε)\)は衝突断面積、\(f(ε)\)は電子エネルギー分布関数 (EEDF)、\(q\)は素電荷、\(m\)は電子質量である。考慮した化学種および化学反応は表1に示されている4)

 表1. 化学反応

ボルツマン方程式は、電子集団を位相空間内の連続体(流体)と見て電子分布関数の時間変化を表現したものである。ここでは、電子エネルギー分布(Electron Energy Distribution Function: EEDF)や電子輸送係数を求めることを指して、2項展開近似を用いる手法を採用し、次のような方程式となる5)

   $\frac{∂}{∂ε}$(\(Wf\)-\(D\)$\frac{∂f}{∂ε}$)=\(S\)

ここで、\(f\)はEEDF (eV-3/2)である。\(S\)はソース項であり、係数\(W\)と\(D\)の計算式は以下のように定義される。

   \(W\)=-\(γε^2σ_ε\)-3\(a\)($\frac{n_e}{N_n}$)\(A_1\)

   \(D\)=$\frac{γ}{3}$($\frac{E}{N_n}$)2($\frac{ε}{σ_m}$)+$\frac{γk_bT}{q}$\(ε^2σ_ε\)+2\(a\)($\frac{n_e}{N_n}$)(\(A_2\)+\(ε^{3/2}A_3\))

変数とその係数の扱いの詳細は例えば文献5)を参照されたい。
計算条件として、圧力は10mTorr,ICP供給パワーは159Wと245W,rf周波数は13.56MHz、温度は300Kである。245Wの供給パワーに基づくプラズマ密度を図2に示す。マクスウェル分布であるEEDFを採用した計算結果も示される。

  図2.ICPにおける電子密度

ボルツマン方程式と流体モデルの連成計算によって得たプラズマ分布はマクスウェル分布であるEEDFに基づいた計算結果に比べて、放電中心に集中し、電子密度の最大値は3倍以上高くあった。ボルツマン方程式との連成計算を採用して得た電子密度はMillerら2)が提示した実測値と比較し、図3に示された。EEDFをマクスウェル分布からボルツマン分布に切り替えることで、電子密度分布の計算値は実測値と大体一致することが分かった。

 図3.実測値と比較する電子密度分布(z=1.5cm)

参考文献

1)菅原広剛,志村尚彦,シミュレーションの諸技法:電子エネルギー分布・速度分布計算法1:ボルツマン方程式解析,放電・プラズマ気相シミュレーション技法調査専門委員会編,電気学会技術報告第1488号,pp. 4-7, 2020.
2)Miller, P.A., Hebner, G.A., Greenberg, K.E., Pochan, P.D., Aragon, B.P., An Inductively Coupled Plasma Source for the Gaseous Electronics Conference RF Reference Cell, J. Res. Natl. Inst. Stand. Technol., Vol.100, No.4, pp. 427-439, 1995.
3)放電・プラズマ気相シミュレーション技法調査専門委員会編,「放電・プラズマ気相シミュレーション技法:佟立柱,竹内希,3.3 COMSOL Multiphysics®を用いた非熱プラズマと熱プラズマの計算」,電気学会技術報告第1488号,pp. 70-74, 2020.
4)Lymberopoulos, D.P., Economou, D.J., Fluid Simulations of Glow Discharge& Effect of Metastable Atoms in Argon, J. Appl. Phys., Vol.73, No.8, pp.3668-3679, 1993.
5)Hagelaar, G.J.M., Pitchford, L.C.,Solving the Boltzmann Equation to Obtain Electron Transport Coefficients and Rate Coefficients for Fluid Models, Plasma Sources Sci. Technol., Vol.14, No.4, pp.722–733, 2005.

*該当のCOMSOLモデルファイルと手順書をご要望のお客様は下記よりご請求ください。

一般

粒子トレーシングモデリング

【アプリ概要】

静電場存在下の荷電粒子の挙動をシミュレートしています。この間電極は移動しています。

粒子間相互作用のない質点モデルでは粒子同士は干渉しますが、剛体球モデルでは干渉しません。これはレナード・ジョーンズ型ポテンシャルの斥力項のみを考慮することで実装できます。剛体球モデルは砂粒や金属粒子等が対象となります。

一方ソフト粒子モデルでは双極子間または誘起双極子間の引力、すなわちロンドン力を含む完全なレナード・ジョーンズ型ポテンシャルにより実装します。凝集系では斥力項だけでなく引力項を含める必要があります。

このアプリはCOMSOL Multiphysicsセミナーツアーコース例題11に準拠しています。

*CAEアプリをお試しいただける特別URLを用意しております。画面右上の「お問い合わせ」からご要望ください。

一般

電気化学インピーダンス分光法

【アプリ概要】

このアプリの目的はEIS、ナイキスト、ボード線図を理解することです。このアプリではバルク濃度、拡散係数、交換電流密度、二重層容量および最大周波数と最小周波数を変更できます。

電気化学インピーダンス分光法(EIS)は電気分析の一般的な手法であり、電気化学系の調和応答を解析します。作用電極の電位に小さな正弦波の変化が適用され、結果として生じる電流が周波数領域で解析されます

インピーダンスの実数成分と虚数成分は、二重層静電容量による表面特性と同様に、セルの動的および質量輸送特性に関する情報を提供します。

*CAEアプリをお試しいただける特別URLを用意しております。画面右上の「お問い合わせ」からご要望ください。

一般

粒あん羊羹をつくろう!

【アプリ概要】

指定数の小豆(回転楕円体)を位置および方向をランダムに指定寸法の羊羹(円筒または直方体)内に配置させます。
小豆同士の干渉チェックは含まれていません。
Javaメソッドだけで作成されています。
これはCOMSOLカンファレンス2020 東京ライブセミナー用モデル最終更新版です。

使用モジュール:COMSOL Multiphysicsモジュール・デザインモジュール

*CAEアプリをお試しいただける特別URLを用意しております。画面右上の「お問い合わせ」からご要望ください。

サブスク限定

マイクロ流体モデル

電気流体力学(EHD)解析
東京理科大学・元祐昌廣先生ご提供

マイクロ流体チップは、掌サイズの基板上に指紋サイズの極めて微細な流路を集積した器具で、少量の細胞・試薬での実験や分析時間の短縮が可能であり、化学・生物の実験効率化が期待されている。しかし、微小デバイス内部の流動は、高い比界面積があり、壁面との相互作用や電気的な力などの影響が強く、現象は複雑になる。ここでは、交流電場の印加によるマイクロ誘起流れ現象の解析を着目する。
図1にはタンパク質分析チップおよび解析モデル1)を示す。計算は静電場、温度場および流動場を連成している。各モデルは図2に示されている。ここでは、AC信号のRMS値を利用し、同じ熱効果を生み出すのであると考えられる。

  図1.タンパク質分析チップおよび解析モデル
  図2.電場および熱流体モデル

本書では、静電場、伝熱およびクリープ流(ほふく流)を連成する電場誘起流れ現象を示すモデルおよび、2章以後に図4と図5の解析結果を得る手順を示す。
計算モデルは、以下に示した計算手順
交流電場=>温度上昇(ジュール加熱)=>流れ物性勾配=>電磁力が作用=>誘起流れの発生
を考慮する。
AC信号のRMS値を利用する静電場の支配方程式は次式で表される。

    \(∇\)∙\({\bf D}\)=ρ  \({\bf D}\)=ϵ\({\bf E}\)
    \({\bf E}\)=-\(∇\)\(V\)

ここで、は比誘導率である。ジュール熱源は

    \(Q\)=$\frac{1}{2}$σ|\({\bf E}\)|2

この熱源を、COMSOLの伝熱 (流体)インターフェースに導入する。

    \(ρC_p\frac{∂T}{∂t}\)+\(∇\)∙\({\bf q}\)=\(Q\)+\(q_0\)
    \({\bf q}\)=-\(k∇T\)

ここで、\(ρ\)は密度、\(C_p\)は比熱容量、\(T\)は温度、 \(k\)は熱伝導率、\(Q\)は熱源、\(q_0\)は面外熱流束である。流体の物性は図3に示した温度に依存するデータであると考えられる。
環境温度は20℃である。静電場と伝熱を連成計算して、流れ計算を行う。クリープ流れインターフェース(spf)を利用する。支配方程式は次式で表される。

    0=\(∇\)∙[-\(p\)\({\bf I}\)+\(u\)+(\(∇\)\({\bf u}\)+(\(∇\)\({\bf u}\))\(T\))+\({\bf F}\)
    \(ρ∇\)∙\({\bf u}\)=0

ここで、\({\bf u}\)は速度、\(p\)は圧力、\(μ\)は粘性係数、\({\bf F}\)は電磁力である。COMSOLには\({\bf F}\)が体積力で設定される。COMSOL変数に次式2)で計算される。

    \({\bf F}\)=-$\frac{1}{2}$[($\frac{∇σ}{σ}$-$\frac{∇ϵ}{ϵ}$)∙\({\bf E}\)$\frac{ϵ\bf E}{1+\left( ωτ\right) ^{2}}$]+$\frac{1}{2}$ |\({\bf E}\)|2\(∇ϵ\)

ここで、\(τ\)=\(ϵ/σ\)は緩和時間である。

  図3.液体試料の熱伝導率と密度

図4.温度分布
  図4.温度分布

液体試料の熱伝導率と密度などの物性を補間関数として設定して、以下の定常計算を行う。

 ステップ1:定常 静電場(es)と伝熱(流体)(ht)
 ステップ2:定常 クリープ流 (ht)

図4に全体の温度分布を示している.誘起流れ場は、図5に示される。図5(a)は温度場および速度場である.図5(b)と図5(c)は電極ギャップ中心線上の流速分布および下流側電極上方の流速分布である。

  図5.電場誘起流れ場の計算結果

参考文献

1)元祐昌廣、数値シミュレーションと実験を駆使した最先端熱流体システム開発、COMSOL Conference 2018 Tokyo.
2)A. Ramos, H. Morgan, N. G. Green and A. Castellanos, “Ac electrokinetics: a review of forces in microelectrode structures”, J. Phys. D: Appl. Phys. 31 (1998) 2338–2353.

*該当のCOMSOLモデルファイルと手順書をご要望のお客様は下記よりご請求ください。