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

CASES AND MATERIALS 事例/資料

キーワード・条件で検索

キーワードで検索

条件で検索

COMSOL解析事例

COMSOL紹介( 導入/検討 )

計測

サブスク限定

コロナ放電によるイオン風の解析

PDEモデルに基づくイオン風のシミュレーション
東京工業大学・竹内希先生ご提供

空気中でコロナ放電を形成すると、イオン風と呼ばれる気体の流れが発生する。本チュートリアルでは、COMSOL Multiphysics®のPDEモデルに基づいてイオン風の簡易的なシミュレーションモデルを紹介する。なお、同様の計算例が静電気学会誌の特集解説として掲載されているので、そちらも参考にしていただきたい。
イオン風のシミュレーションは、図1に示す軸対称二次元モデルで行った。計算領域は半径10mm、高さ10mmの円筒形状で、上部に設置した針電極の曲率半径はおよそ30μm,針電極先端と下部の平板電極のギャップ間距離は5mmである。平板電極を接地して針電極に負極性高電圧を印加すると、針電極先端で負極性コロナ放電が発生する。負極性コロナ放電の計算を厳密にシミュレーションする場合は、電離や付着といった電子衝突反応、分子の解離反応、生成された重粒子同士の反応、クラスターイオンの形成などを考慮する必要がある。また、負極性コロナ放電は、トリチェルパルスと呼ばれる不連続な形態となることもある。しかしながら、電子衝突反応の計算はメッシュサイズと時間刻みを極めて小さくする必要があり、メモリや計算時間が膨大となる。また、電離領域のイオン風に与える影響は小さく、ドリフト領域でのイオン分布が重要である。そこでこの簡易シミュレーションモデルでは、電子衝突反応の計算は行わずに、空気中での平均的な物性値を有する負イオンのみを荷電粒子として考慮し、針電極表面での負イオン密度を境界条件として与えることで、負極性コロナ放電のドリフト領域を模擬した。

  図1.解析モデル

本書ではモデルおよび、図2の解析結果を作成する手順を示した。
計算モデルは、まず、ポアソン方程式と負イオンの連続の式の連成計算により、電位および負イオンの分布を計算した。計算コスト削減のため、およそ5×5mm2の領域を負イオンのドリフト領域とする。ポアソン方程式は下式で表される。

     -\(∇\)⋅\(ε_rε_0∇V\)=\(ρ\)

ここで、\(ε_0\)は真空の誘電率であり、比誘電率\(ε_r\)は1とした.また、空間の電荷密度\(ρ\)は、後述する負イオンの数密度\(N_n\)と電気素量\(e_0\)を用いて、

     \(ρ\)=-\(e_0N_n\)

と表される。針電極の電位を-5kVとして与えた。下部の境界(平板電極)の電位は零とし、また、上部と側面の境界では、電束密度の法線方向成分を零とした。
負イオンの連続の式を、係数型PDEモデルによって解く。

     $\frac{∂N_n}{∂t}$+\(∇\)⋅(-\(D_n∇N_n\)-\(μ_nEN_n\))=0

ここで、\(D_n\)は空気中の負イオンの拡散定数、\(μ_n\)は負イオンの移動度、\(E\)は電界強度であり、

     \(E\)=-\(∇V\)

で表される。係数型PDEモデルにおける境界条件には、全ての境界でDirichlet条件を用いた。針電極表面での負イオン電荷密度を、対称軸からの半径方向距離rと時間tの関数としてr≤4×10–5mの範囲で

     \(N_n\)=5×1017×(1-$\frac{r^4}{(4×10^-5)^4}$)×tanh⁡(104t)

として与え、負イオンの初期数密度は宇宙線や放射線により形成された負イオンの数密度、~106m-3とした。
円筒座標系では、\(∇\)∙\(Γ\)は

     \(∇\)∙\(Γ\)=$\frac{1}{r}$$\frac{∂(rΓ_r)}{∂r}$+$\frac{∂Γ_z}{∂z}$=$\frac{∂Γ_r}{∂r}$+$\frac{∂Γ_z}{∂z}$+$\frac{Γ_r}{r}$

となる。\(Γ\)=\(-D_n∇N_n\)-\(μ_nEN_n\)である。COMSOLのPDEインターフェースはこの式の最後の項が含まれない。次式で表す係数形式PDEの係数は表1に定義される。その以外の係数は0になる。

     \(e_a\)$\frac{∂^2 N_n}{∂t^2}$+\(d_a\)$\frac{∂N_n}{∂t}$+\(∇\)∙(-\(c∇N_n\)-\(αN_n\)+\(γ\))+\(β\)∙\(∇N_n\)+\(aN_n\)=\(f\)

 表1.PDE係数

イオン風の駆動力は、荷電粒子に働くクーロン力となる。よって、負イオン分布の計算で得られた電荷密度\(ρ\)と電界強度\(E\)の積によりクーロン力\(F\)を求め、これをNavier–Stokes方程式の外力項に代入することで、イオン風の計算が可能である。
Navier–Stokes方程式は下式で表され、変数は流速\(U\)と圧力\(P\)であり、流体力学モデルで解く。

     \(ρ_g\frac{∂U}{∂t}\)+\(ρ_g\)(\(U\)⋅\(∇\))\(U\)=-\(∇P\)+\(μ_g∇^2U\)+\(F\)

ここで、\(ρ_g\)は空気の密度、\(μ_g\)は空気の動粘性係数である。流体力学モデルにおける境界条件は、全ての境界ですべりなしとした。
計算領域における流速の初期値を零、圧力の初期値を1気圧とした。また、クーロン力を決定する電荷密度\(ρ\)と電界強度\(E\)には事前に計算した結果を用い、時間的に変化しないとして、流速および圧力の時間変化を計算した。
図2にコロナ放電によって発生した負イオン分布およびそれによるイオン風を示している。

  図2.コロナ放電における負イオン分布およびそれによるイオン風の流速

参考文献

1)竹内希, “コロナ放電により発生するイオン風のCOMSOL Multiphysicsを用いたシミュレーション[特集解説]”, 静電気学会誌, Vol. 40, No. 4, pp. 168–171, 2016.
2)竹内希, “COMSOLおよびBOLSIGによるシミュレーションの実例”, 大気圧プラズマ工学ハンドブック, 第2編第2章, エヌ・ティー・エス, pp. 273–282, 2013.

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

サブスク限定

アーク溶接におけるMHDシミュレーション

ガスタングステンアーク溶接(GTAW)、別名,タングステンイナートガス溶接(TIG)は、船舶、パイプライン、自動車など広い領域に応用されている。このプロセスの数値シミュレーションは、溶接入力パラメーターに応じて溶接部の特性を予測する大きな挑戦に臨んでいるが、それは溶接を改善し、品質と生産性を向上させる(1)。このプロセスに関係する多くの物理現象の間の強い結合があるため、数値モデリングは複雑になる。
ステンレス鋼へのGTA溶接の数値解析では、熱伝達と流れの挙動に応じて、温度場、流れ場および電磁場は、統一されたMHDの定式化でモデリングできる。アークプラズマにおける電子、イオン、中性粒子の輸送は解かず、伝熱、流体、電磁気を解く。ここでは、ガスの絶縁破壊を示すことで、放電によるガスの導電率の変化およびガスの熱伝達を計算した。図1のように、解析モデルに以下の計算条件を考慮する。

    図1. 解析モデル
  • 2D axial symmetry 2D軸対称
  • Direct current (DC) at the cathode カソードでの直流
  • Anode is grounded アノード接地
  • Electrode polarity: Cathode 電極の極性:カソード
    Direct current (DC) flows in one direction, resulting in a constant polarity
    直流が一方向に流れるため、極性が一定である。

本書では、解析モデルおよび、2章以後に図4の解析結果を作成する手順を示す。計算用のCOMSOLインターフェースが図2に示している。プラズマ放電の熱源${Q}$を伝熱(流体)インターフェースに連成させて、支配方程式を式(1)に示す。

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

図2. COMSOLインターフェース
図2. COMSOLインターフェース

ここで、\(ρ\)は密度、\(C_p\)は比熱容量、\(T\)は温度、\(u\)は流速、\(q\)は熱流束、\(k\)は熱伝導率である。熱源\(Q\)は、マルチフィジックス機能によって以下に定義される。

  1)抵抗加熱 (オーム加熱):Q=J∙E
  2)正味放射損失 \(Q_{rad}\)
  3)エンタルピー輸送
   $\frac{∂}{∂T}$($\frac{5k_BT}{2q}$)\((∇T\)∙\(J)\)

カソード境界熱源\(Q_b\)は式(2)で表される

   \(-n\)∙(-\(k∇T\))=\(Q_b\)
   \(Q_b\)=-\(|J_{elec}|ϕ_s\)+\(|J_{ion}|V_{ion}\) 
   (2)

ここで、\(ϕ_s\)は表面仕事関数、\(V_{ion}\)はプラズマの電離電位、\(J_{elec}\)と\(J_{ion}\)はそれぞれ電子とイオンの電流密度である。
環境温度は20℃、ガスはアルゴンである。温度に依存するガスの導電率は図3に示されている。

  図3. アルゴンの導電率

放電時間は1sとします。図4はガスの導電率及び温度の分布の計算結果である。ガスの絶縁破壊を生じる時の導電率及び温度の進展は示されている。

  図4. ガスの導電率及び温度の分布

参考文献

1)A. Traidia, F. Roger, A. Chidley, J. Schroeder, T. Marlaud: “Effect of Helium-Argon Mixtures on the Heat Transfer and Fluid Flow in Gas Tungsten Arc Welding”, World Academy of Science, Engineering and Technology, 49, 2011.

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

サブスク限定

DCパルス放電プラズマシミュレーション

プラズマ窒化は古くから鉄鋼材料の表面改質用プロセスとして研究されてきた。パルスプラズマは、エッチング速度と堆積速度の向上、ダスト粒子の形成の減少、および堆積の均一性の特徴を示す。パルス放電窒化技術は、脆弱化合物層(窒化鉄化合物)の有無、窒化拡散層の深さを自在にコントロールし、理想の窒化層を形成できる。
近年、時間分解測定法による窒素とアルゴンにおけるパルスDC放電の電子密度分布の測定結果が掲載された1)。カソードに印加された電圧は-1~-2kV、パルス繰り返し周波数は0.4〜25kHz、デューティ比は10-50%である。この論文に対してDCパルス放電プラズマシミュレーションも既に報告された2)
そこで本例題は報告された解析モデル2)を参考しモデルを再構築する.モデル構造図は図1に示され、モデル領域はアノードにより囲んでカソード以上の半分の領域とする.ガスは窒素、温度は300K、圧力は10Pa、パルス繰り返し周波数は1kHz、デューティ比は25%である。アノードは接地し、カソードは電圧が-2kVで印加され、外部RC回路と繋がる。バラスト抵抗は100Ω、プロッキング容量は50nFである。
本書ではモデルおよび、図3の解析結果を作成する手順を示した。

  図1.モデル構造図1)

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

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

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

図2.COMSOLの非線形時間依存ソルバー
図2.COMSOLの非線形時間依存ソルバー

電子エネルギー分布関数はマックスウェル、カソードにイオン衝突による二次電子放出係数は0.25、二次電子平均エネルギーは2.5である。また、電極表面にイオンと励起種の付着係数は1, 窒素原子の付着係数は0.07である。
プラズマシミュレーションは、高度な非線形解析なので、ここではソルバーの追加設定する必要になっている。図2にCOMSOLの非線形時間依存ソルバーの強連成ソルバーの設定画面を示す。ヤコビアン更新は各タイムステップから「反復毎」に切り換えた。最大反復回数とトレランス因子も調整した。
図3にパルスDC放電シミュレーション結果を示している。計算時間は17時間42分でした。計算用のパソコンスペックは以下の通りである。

1.OS: WINDOWS 10
2.CPU: Intel(R) Xeon(R) CPU E5-1660 v2 @ 3.70GHz
3.メモリ RAM 128 GB

  (c) 平均電子密度およびカソード表面のパルス電位

  図3. パルスDC放電シミュレーション結果

参考文献

1)A. Pandey, W. Sakakibara, H. Matsuoka, K. Nakamura, and H. Sugai, “Time-resolved curling-probe
measurements of electron density in high frequency pulsed DC discharges”, Jpn. J. Appl. Phys., 55, 016101 (2016).
2)L. Z. Tong, “Simulation Study of High-Frequency Pulsed DC Discharges in Nitrogen”, Proceedings of the 2017 COMSOL Conference in Boston.
https://www.comsol.jp/paper/simulation-study-of-high-frequency-pulsed-dc-discharges-in-nitrogen-51903
3)放電・プラズマ気相シミュレーション技法調査専門委員会編,「放電・プラズマ気相シミュレーション技法:佟立柱,竹内希,3.3 COMSOL Multiphysicsを用いた非熱プラズマと熱プラズマの計算」,電気学会技術報告第1488号,pp. 70-74 (2020).

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

サブスク限定

軸対称定常熱伝導解析

PDEモデルに基づく2D軸対称シミュレーション

軸対称問題は中心軸を含む任意の断面内で物理現象が同一となる問題である。この問題は二次元問題ときわめてよく似ており、円柱座標系を使う二次元的な解析手法が適用できる。COMSOL Multiphysicsは、軸対称物理インターフェースで円筒座標をサポートしている。しかしながら、COMSOL MultiphysicsのPDEモデルを利用する際、デカルト座標系から円筒座標系に手動で変更する必要がある。本チュートリアルでは、温度拘束された円筒の軸対称定常熱伝導解析を例として、COMSOL MultiphysicsのPDEモデルで座標変換する手法を紹介する。
軸対称定常熱伝導解析モデルを、図1に示す。

   図1. 解析モデル1)

寸法と物性値は表1に示される。

ここで、κは熱伝導率である。境界条件では、内壁と外壁の温度はそれぞれ10℃と0℃で、上下の境界は断熱とする。このモデルに対して、半径方向 r[mm]の温度分布の解析解を次の式により計算できる。

本書ではモデルおよび、図2の解析結果を作成する手順を示す。
定常熱伝導問題の微分方程式は次式で表される。

   \(ρC_p u\)∙\(∇T+∇\)∙\(q=Q\)
   ${q}$=-${κ∇T}$

ここで、u=0とQ=0であり、計算方程式は以下のようになる。

   ${∇}$∙${(}$-${κ∇T)=0}$

ここで、${Γ}$=-${κ∇T}$は熱流束を意味する。COMSOLのPDEインターフェースは${∇}$∙${Γ}$を

   $\frac{∂Γ_r}{∂r}$+$\frac{∂Γ_z}{∂z}$

で示す。これは${∇}$∙${Γ}$の円筒座標系の形式と一致しない。
円筒座標系では、${∇}$∙${Γ}$は

   ${∇}$∙${Γ}$=$\frac{1∂}{r}$$\frac{(rΓ_r)}{∂r}$+ $\frac{∂Γ_z}{∂z}$=$\frac{∂Γ_r}{∂r}$+$\frac{∂Γ_z}{∂z}$+ $\frac{Γ_r}{r}$

となる。COMSOLのPDEインターフェースはこの式の最後の項が含まれない。
COMSOLの物理インターフェースはデカルト座標系または曲線座標系に対応する微分演算子(共変微分. 物理法則は座標系に関わらず成り立つ)を利用するため、デカルト座標系から曲線座標系に変更する必要がなくなるが、PDE方程式ベースのモデリングでは、偏微分が使われるため、この座標系の変更が必要になる。本チュートリアルでは、COMSOLの係数形式PDEおよび一般形式PDEインターフェースを利用する。支配方程式は次式で表される。

ここで、${∇}$=[$\frac{∂}{∂r}$,$\frac{∂}{∂z}$]、変数とその係数の扱いの詳細は例えば文献2)を参照されたい。${∇}$∙${Γ}$の発散式によって、係数形式PDEの係数は以下に定義できる。その以外の係数は0になる。

一般形式PDEの係数は${Γ}$に、-${κ×T_r}$と-${κ×T_z}$、${f}$に${κ×T_r/r}$で、以下のように設定される。

  図2. 温度分布

COMSOLの伝熱(固体)、係数形式PDEおよび一般形式PDEのインターフェースによって計算した結果を図2に示した。温度分布の計算結果は解析解とよく一致したことは示された。

参考文献

1)田中正隆,松本敏郎,中村正行,計算力学とCAEシリーズ: 境界要素法,培風館,1991.
2)COMSOL Multiphysics 5.6 – http://www.comsol.com/products and COMSOL Multiphysics Reference Manual.

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

サブスク限定

ノズル噴流のシミュレーション

CFD解析のベンチマーク分析

噴流現象は、ノズル、スリットなどの小孔から速度を持った流体が空間中に噴出する現象であり、加熱・冷却、洗浄、消防、加工、噴霧、燃焼、など極めて多様な分野で使用されており工業的にも極めて重要な事象である。
ノズルによる噴流はノズル出口近傍の速度分布、拡がり、巻き込み特性、などの流動特性がノズル形状によって大きく影響される。そこで本例題は図1に示す四分円ノズルの軸対称円形ノズルを用い、圧縮性を考慮しなくてよい範囲で十分発達した乱流の場合(Re=1.5×1041)について、空気噴流を計算し、測定値と比較して噴流の流動特性を検討する。CRはノズル絞り面積比である。

  図1. ノズル形状

解析モデルは図2に示す。軸対称二次元モデルであり、直管部の長さL=500mm,\(L/d_o\)=50、ノズル出口からモデルの流出口までの距離は400mm,流出口幅は80.75mmである。温度は20℃,圧力は1atmである。
本書ではモデルおよび、図3の解析結果を作成する手順を示した。

  図2.解析モデル

計算には、${k}$-${ε}$乱流モデルを用いて圧縮流と非圧縮流を考慮する。レイノルズ数${Re}$は次の式で表される1)
 \(Re=u_m d_o/v\)   (1)
ここで、\(u_m\)はノズル出口平均速度、\(d_o\)はノズル出口直径、${v}$は動粘性係数である。本モデルは、流入口の流速が式(1)に示したノズル出口のレイノルズ数${Re}$により決まる。COMSOLスタディの補助スイープを用いて、Re=1.5×104に満たす流入口の流速\(u_m\)=2.565m/sを得た。
${k}$-${ε}$乱流モデルの支配方程式は次の式で示される。

   (2)

紙数の都合でここではこれらの方程式の説明を省略する。
COMSOLスタディの補助スイープによって、計算したノズル出口平均速度は以下のテーブルに示され、理論値とフィッティングしたことが明らかになった。

図3(a)にノズル噴流の軸対称二次元モデルの計算結果の3D表示である。最大の流速は24.8 m/sになった。測定値1)と比較した結果は図3(b)に示され、大体一致することが分かった。

  図3. ノズルの噴流速度および測定値との比較

参考文献

1)鬼頭みずき,社河内敏彦,軸対称噴流の流動特性に対するノズル形状の影響 奈良工業高等専門学校研究紀要 , No.45, pp.19-24 (2009).

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