【アプリ概要】

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

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

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

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

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

【アプリ概要】

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

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

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

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

【アプリ概要】

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

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

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

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

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

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

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

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

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

この熱源を、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変数に次式²⁾で計算される。

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

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

　　図4．温度分布

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

　ステップ１：定常　静電場（es)と伝熱（流体）(ht）
　ステップ２：定常　クリープ流 (ht）

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

参考文献

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モデルファイルと手順書をご要望のお客様は下記よりご請求ください。

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

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

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

　　　　　-\(∇\)⋅\(ε_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×10¹⁷×(1-$\frac{r^4}{(4×10^-5)^4}$)×tanh⁡(10⁴t)

として与え、負イオンの初期数密度は宇宙線や放射線により形成された負イオンの数密度、～10⁶m^-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\)

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

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

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

参考文献

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

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