第 24 回数値流体力学シンポジウム
E1-5

プラズマアクチュエータにおける推力生成の数値解析
Numerial analysis of thrust generation in plasma actuator
江波戸 翔一, 東北大工, 980-8579 仙台市青葉区荒巻字青葉 6-6-01, E-mail: ebato@rhd.mech.tohoku.ac.jp
荻野 要介, 東北大工, 980-8579 仙台市青葉区荒巻字青葉 6-6-01, E-mail: yogi@cfd.mech.tohoku.ac.jp
大西 直文, 東北大工, 980-8579 仙台市青葉区荒巻字青葉 6-6-01, E-mail: ohnishi@rhd.mech.tohoku.ac.jp

Shoichi EBATO, Department of Aerospace Engineering, Tohoku University, Sendai 980-8579, Japan
Yousuke OGINO, Department of Aerospace Engineering, Tohoku University, Sendai 980-8579, Japan
Naofumi OHNISHI, Department of Aerospace Engineering, Tohoku University, Sendai 980-8579, Japan
We have developed a PIC-MCC code to study detailed process of plasma formation and flow generation in adielectric barrier discharge (DBD) plasma actuator. Two-dimensional simulations were conducted by this codewith negative constant input voltage and positive constant one. When the negative constant input voltage is
supplied, plasma is generated on the dielectric surface because the accelerated electrons lead an electron avalanche,
and the spreading plasma distorts the initial potential distribution. The plasma is formed with sheathing theexposed electrode. In the case of the positive constant input voltage, ion cloud is generated on the exposed
electrode and spreads to the downstream followed by the plasma.

はじめに
近年, プラズマを利用した流れ場の制御や推進技術の研究が活発に行われている. 代表的な例は, 強力な磁場中にプラズマを通過させて電流を取り出す MHD 発電や,電気エネルギーによってプラズマを制御し, 推進力を得る電気推進機である. さらにここ数年, プラズマを利用した流れ場制御デバイスとしてプラズマアクチュエータが注目されている.1–4
プラズマアクチュエータが注目されている理由は, 既存の流れ場制御デバイスに比べていくつかの有利な点があるためである. 一つ目は可動部がなく構造が簡単な点である. 境界層制御に用いられている既存のデバイスは機械的な制御を必要とする. 二つ目は短時間で流れ場を制御できる点である. Singh らの研究によれば, 数マイク
ロ秒のオーダーで流れ場を制御することができる.5 三つ目は電磁的な力を比較的正確に与えることができる点である. 以上の点から, プラズマアクチュエータは実験的にも数値的にも精力的に研究されている.
現在最も研究されている誘電体バリア放電 (dielectricbarrier discharge: DBD) を用いたプラズマアクチュエー
タは, 2 枚の電極と誘電体で構成された比較的シンプルで軽量なデバイスである. 図 1 に示すように, 1 枚は気流に対して露出しており, もう 1 枚は誘電体により被覆されている. 露出電極には交流電圧を印加するが, 電子は自らの電荷により周囲の電場を弱めて加速を抑制するとともに,周囲の中性粒子と衝突を繰り返す. 衝突時のエネルギーが高いと電離衝突となり, 中性粒子から放出された電子がさらに次の電離衝突を引き起こす. この一連の過程を電子雪崩と呼ぶ. 電離で生成されたイオンが中性粒子と衝突し, 運動量を輸送することで, 流れが誘起される.
誘起された流れの有効性についてはすでに実証されている. 藤井ら6 は, プラズマアクチュエータが流体に与える力を Suzen モデルによって外力項として導入し, 流体解析を行っている. しかしながら, 何が誘起される流れに支配的な影響を及ぼしているかということや, 実際に剥離流れが付着流れに変化するメカニズムについては不明
である.7
最近では, プラズマアクチュエータを火星飛行機に適用することが検討されている. 火星飛行機とは火星探査のための手段として検討されている航空機のことで, 山岳, 渓谷の探査が可能となるため注目されている. しかし地球と火星の環境が大きく異なるため, 安定した飛行を実現することが難しい. そこで能動的な制御技術としてプラズマアクチュエータが注目されている. 高垣ら8 や安倍ら9, 10 は, 火星飛行機で使用するという目的の下, 生成される推力の圧力依存性や化学種依存性, 交流電圧の入力波形依存性, 誘電体の厚さ依存性に関してパラメータ・
スタディを行っている. しかしこれらの特性について何が支配的な因子であるのか未だ明らかになっていない.
地球大気環境下においても, より効率的な推力生成を目的として, 支配的な因子を突き止めるべく様々な研究が行われている. Kim ら11 や Enloe ら13 は入力電圧の正弦波が負の電圧である時に推力が生成されると主張している. Kim らはこの因子として負イオンの存在を指摘している. さらに Font ら12 は入力電圧の大きさと酸素含有率を変化させて, 負イオンが推力に与える影響を調べ, 負イ
オンが支配的な因子になりうる条件を示した. Singh ら5は様々な化学反応を考慮することで数値的に研究を行っている. また, よりよい入力波形や電極構成についての議論もされている. Likhanskii ら14 は正にバイアスをかけ,一定周期でバースト波を加える入力電圧がよいと報告している. Opaits ら15 は 3 つの電極を用いた配置が, 誘電体表面に堆積する荷電粒子を除去する役割を果たすため,効率的だと結論づけている.
また, プラズマアクチュエータ内でおこる放電過程についての議論も盛んである. 放電過程は推力生成に直接影響を与える荷電粒子を生成するため, 非常に重要である. 西田ら16 や Orlov ら17 は交流電圧のフェイズの違いで放電過程が違うことを指摘している. プラズマアクチュ
エータでは入力波形として正弦波がよく用いられるため,電圧が徐々に増加する場合と電圧が徐々に減少する場合

1.

Fig. 1: Schematic view of DBD plasma actuator.1

1

Copyright c 2010 by JSFM

第 24 回数値流体力学シンポジウム
E1-5の放電過程を調べているが, 前者の場合ではストリーマ
放電が見られ, 後者の場合ではグロー放電が見られることが指摘されている. ストリーマ放電が起こるメカニズ
ムについては Boeuf ら18, 19 が詳細に調べている. 彼らは露出電極付近にまずイオン雲が生成され, それを押すよ
うにしてストリーマが生成されると述べている. ところが Likhanskii20, 21 は入力電圧が低い場合には光電離過程を考慮しないとストリーマ放電が生成されず実験と一致しないことを主張しており, 放電過程の解明には議論の余地がある.
以上, 様々な特性における因子や放電過程の解明について議論の余地がある. さらにプラズマアクチュエータ
内の流れ場では熱的にも化学的にも非平衡である可能性があるが, これら研究で用いられているモデルには熱的非平衡性が考慮されていない場合が多い. そこで本研究では, プラズマアクチュエータの有効性を検討することを目標として, プラズマアクチュエータにおける流れ場
誘起を再現することのできるモデル化と数値計算コードの開発を行った. 本報では, プラズマの放電過程と, その際の運動量輸送の結果生じる力について述べる.

2.

Yj+1
yn

D

D = εE

xn

Xi

D

B

A

X i+1

(a)

(b)

Fig. 2: Weighting method of ﬁnite size particle.

る方法である. 前述のように, 超粒子は点電荷ではなく大きさを持っている. その大きさはセル程度の大きさであり, 近隣のセル内に粒子が何割存在しているか比率を計算することで, 粒子がまたがっているセルに物理量を分配する.
2 次元の場合の weighting の具体的な方法を示す. 粒子の中心位置の座標が (xn , yn ) である超粒子は, 以下の範囲に電荷密度分布を持つ.

2.1 基礎方程式
本研究では, プラズマと電磁場, および中性気体の相互作用を数値的に取り扱う. プラズマアクチュエータ内の流れ場では熱的にも化学的にも平衡が成立しない可能性があるため, 種々の非平衡過程を比較的平易に実現することができる粒子法を用いて流れ場の再現を試みた. 荷電粒子と電場との相互作用は Particle-In-Cell (PIC) 法22 を用いて評価する. その上で Monte Carlo Collision (MCC)
法23, 24 を用いて粒子間の衝突による流れ場の発展を再現する. PIC 法では各セルに配置された荷電粒子の運動を追跡する. 電磁場はセルに定義された離散量によって解き, 粒子はそれぞれ運動方程式を解く. 基礎方程式は以下に示す Poisson 方程式と粒子の運動方程式である.
dx
= qE
dt
∇ · D = −ρ

C

C

Yj

数値計算法

m

B

A

xn − ∆x/2 ≤ xn ≤ xn + ∆x/2

(4)yn − ∆y/2 ≤ yn ≤ yn + ∆y/2(5)

一方, (Xi , Yi ) に位置する格子点が持つ影響領域は図 2(a)における領域 D, (Xi+1 , Yi ) に位置する格子点が持つ影
響領域は領域 C, (Xi , Yi+1 ) に位置する格子点が持つ影響
領域は領域 A, (Xi+1 , Yi+1 ) に位置する格子点が持つ影
響領域は領域 B である. 従って面積割合より数値的には以下の電荷密度が配分される.

(1)
(2)
(3)

N
∑
(Xi+1 − xn )(Yi+1 − yn )
ρ
∆x∆y
n=1

(6)

ρi+1,j =N∑(xn − Xi )(Yi+1 − yn )
ρ
∆x∆y
n=1

(7)

ρi,j+1 =N∑(Xi+1 − xn )(yn − Yi )
ρ
∆x∆y
n=1

(8)

N
∑
(xn − Xi )(yn − Yi )
ρ
∆x∆y
n=1

(9)

ρi,j =

ここで m は粒子の質量, q は電荷量, ρ は電荷密度である.
本研究では静電場を考えるため, Maxwell 方程式に代わり Poisson 方程式を用いている. 運動方程式と Poisson 方
程式をつなげる物理量は電荷密度 ρ である. これらは粒子の運動を決定する位置と電荷で決まるため, 式 (1) と式(2) で, PIC 法における方程式系を閉じることができる.
以下の節で PIC 法について紹介する.

2.2 PIC 法
2.2.1 超粒子 現実のプラズマは非常に多数の粒子で構成される. 粒子 1 つ 1 つの運動をすべて追跡すること
は計算資源の観点からも非現実的である. また, 粒子の運動を追跡する場合, 任意の粒子位置において電磁場が必要となる. そこで PIC 法では多くの現実のプラズマ粒子
を代表する, 大きい電荷, 大きい質量をもつ超粒子を考え,これを用いて数値解析を行う. そしてこの超粒子に有限の大きさを持たせ, 空間格子点上に weighting25 と呼ばれる手法を用いて, 電流密度や電荷を配分する. 超粒子に作用する電磁場は, 配分した電荷をもとに Poisson 方程式を用いて更新する. そして各超粒子の運動を更新する際には, その隣接する格子点上の電磁場を超粒子位置に内挿
し, それを用いて運動方程式を解く.

ρi+1,j+1 =ここで ρ は電荷密度である.

2.3 MCC 法
気体の流れを考える時, 密度, 圧力, 温度といった状態量に着目する. これらを数値解析する際には, 流れ場を連続体で近似して連続の式やナビエ・ストークス方程式を解くのが一般的である. 気体は多数の分子の集まりであり, 流れ場のある点に多くの分子が存在して十分多くの衝突が起こっているため, 流れ場は連続体近似でき, 上記のような解析が可能となる. このとき, 微視的な分子の衝突は考慮する必要がない. しかし, 流れ場を連続体で近似しない場合, 分子の衝突を考慮して流れ場を解く必要がある. 分子運動の基礎方程式は以下に示すボルツマン方程式であり, この数値解法として MCC 法は広く知られ2.2.2 weighting weighting は, 電磁場に影響を与える超粒子の位置や速度にしたがった電流や電荷を評価す
2

Copyright c 2010 by JSFM

第 24 回数値流体力学シンポジウム
E1-5

本研究では, 一度の試行で衝突の有無と衝突形態が決定できる Nanbu 法を用いた.
粒子 i が粒子 k と衝突する場合を考える. Nanbu 法では, まず長さ 1 の線分を, 弾性衝突や電離衝突などの起こ
りうる衝突形態の数で等分する. そして等分した各小区
間左端に, 以下で示すような長さ Pe の線分をとる.

ている.

∂(nf )∂(nf )∂(nf )
+c·
+F·
=
∂t
∂c
∫ ∫∂x

[f (c0 )f (ζ 0 ) − f (c)f (ζ)]gσdΩdζ

(10)Pe (m) = Nk σm ()gik ∆t

ここで n は実空間, 速度空間における数密度, f は速度分布関数, g は相対速度の大きさ, σdΩ は微分断面積である.
MCC 法では多数のサンプル粒子を用い, これらのサン
プル粒子の挙動を追跡することで解析する. 粒子の挙動
は粒子の空間移動 (位置の変化) と粒子の衝突 (速度の変化) を微小単位時間ステップで分離して交互に計算を繰り返すことで表現する. 粒子の衝突に関して確率的な処理を行うため, 分子間の力を逐次計算する分子動力学法に比べ計算機の負荷を大幅に軽減することができる. ま
た粒子の衝突モデルを変更することによって, 単原子分
子気体だけではなく, 多原子分子気体や混合気体, また化
学反応や解離反応を含んだ流れ場へ適用可能である.
式 (10) から粒子の衝突則を導く詳細な議論は文献24 に譲る. ここでは導出される衝突則から提案された手法に関して述べる.

ここで m は衝突形態を表し, σm () は質量中心系での衝突エネルギー  に対する衝突断面積を表している. 1 つの乱数が m 番目の小区間の線分 Pe (m) 上にあるとき, 衝突
形態 m で衝突したこととする. 乱数がどの Pe (m) 上にもない時は, 無衝突であるとする. Nanbu 法の概略図を
図 3 に示す.

Pe (1)▲

3.

入力電圧のフェイズの違いによる放電過程の変化

3.1 計算条件
推力生成の支配的な因子を調べることを目的として, プ
ラズマアクチュエータにおける放電過程を調べた. 図 5に示すように, プラズマアクチュエータを模擬した電極
間の 2 次元計算を行った. 2 mm × 2 mm の空間に, 長さ0.25 mm, 厚さ 0.5 mm の露出電極, 長さ 1.25 mm, 厚さ0.15 mm の被覆電極, 厚さ 0.5 mm の誘電体を配置した
200 × 200 個の計算格子を用いた. 被覆電極を 0 V に保ち, 露出電極には −5000 V と, 5000 V の一定電圧を印加した. 初期条件として空間全体に 107 /m3 の数密度をもつ電離した気体を与え, 温度は 300 K で圧力 15 kPa とした. 粒子法を用いているため, 電子雪崩による荷電粒子の爆発的な増大は大きな問題となる. この問題は, 粒子数

(12)

k=1

となる. ここで, n は数密度, gik は粒子 i と粒子 k の相対速度, σT は全衝突断面積である. なお全衝突断面積は相対速度に依存する. ゆえに, 非常に短い時間 ∆t の間に粒子 i が衝突する確率 νi ∆t は

Equation of motion

Nj

(13)
▲

k=1

weightingとなる. なお, ∆t の間に粒子 i と粒子 k が衝突する確率は

Δt(E,B) →MCC
v → vʼ

▼n∆tgik σtNj

Pe (3)

2.4 PIC 法と MCC 法の結合本研究では, 荷電粒子と電場との相互作用を PIC 法, 粒子間の衝突による反応や速度の更新を MCC 法を用いて実現するため, PIC 法と MCC 法を結合する必要がある.
結合には, Birdsall の方法28 を用いた. 計算サイクルの模式図を図 4 に示す. この方法では, PIC 法のサイクルの中に MCC 法による速度の更新を組み込んでいる.

Nj

∑
n
∆t
gik σTNj

Pe (2)

Fig. 3: Concept of the Nanbu method in the case of
three chemical reactions.

となる. このようにして求めたセル番号 j の粒子数を 1個ずつ増加させていくことにより, 仮に粒子 i がセル番号
j のセルの k 個目の粒子であったとすると, 「粒子 i」という情報と「セル番号 j のある k 番目の粒子」という情報が 1 対 1 に対応されることになる.
j 番目のセルに存在する粒子数を Nj とし, その中の粒子 i に注目する. この粒子が他の粒子と衝突する頻度はn ∑gik σT
Nj

Non-collision

Collision of case1

2.3.1 衝突確率 MCC 計算では衝突する粒子は同一セルに入っている粒子同士で起こると仮定され, まず粒子がどのセルに存在しているかを調べる. 1 次元的な計算領域に存在している等間隔のセルを考える. 計算領域の長さがL, セル数が N とすると, セルの一辺の長さは ∆L = L/Nとなる. よって位置座標 pi である粒子 i の存在するセル
番号は
( p )
i
j = int+1(11)∆L

νi = ngσ¯T i =

(15)

▲

n2

(14)
Maxwell equationsとなり, 式 (13) は系が非平衡状態にあっても成立する.

(ρ,J)

2.3.2 Nanbu 法 t = 0 における粒子の速度の集合 ciから, t = ∆t における粒子の速度の集合 c0i を求めるた
めの衝突判定は, Nanbu 法や衝突数法などいくつかある.

→

(E,B)

weighting(x,v) → (ρ,J)

Fig. 4: Simulation cycle of PIC-MCC method.

3

Copyright c 2010 by JSFM

第 24 回数値流体力学シンポジウム
E1-5

mm
2.0↑

1.0

exposed electrode

(a)

0.5

(b)

dielectric0buried electrode
1.5

0.25

Fig. 6: (a) Initial distribution of electric ﬁeld intensity, (b) potential and electric ﬁeld vector in the case ofnegative voltage input.

→
2.0
mm

Fig. 5: Simulated condition.
がある程度増大したところで, 複数の超粒子をまとめ, 再度 1 つの超粒子を作ることによって克服した. プラズマ
アクチュエータにおける電極は一方が誘電体によって被覆されており, 荷電粒子は誘電体上に堆積する. ゆえに境界条件として, 荷電粒子が誘電体に到達した時点で速度を 0 m/s とした. 露出電極では完全吸収壁を仮定し, 露
出電極に到達すると電子は電極に吸収され, イオンは中性化すると仮定する. 高エネルギーのイオンが露出電極に到達することによって電極から放出される電子を考慮するため, 2 次電子放出係数を 0.3 とした. この際放出される 2 次電子の速度は初期温度の熱速度を与える.

3.2 化学反応モデル
地球大気環境下における化学種の上位 3 種の構成は,
N2 78%, O2 21%, Ar 1% である. そこで本研究では, N2
環境下において計算を行った. N2 の衝突過程には以下の反応を考慮した.
e + N2 ←→ e + N2

(16)

e + N2 −→ 2e + N2 +(17)

e + N2 −→ e + N2 (A)N2 + N2 ←→ N2 + N2+(a) t =0.1 ns(b) t =0.2 ns(c) t =0.5 ns(d) t =0.7 ns

Fig. 7: Net charge density at early time in the case ofnegative voltage input.

(18)
+

(19)

N2 + + N2 ←→ N2 + + N2

(20)

式 (16), (20) は弾性衝突, (17) は電離衝突, (19) は電荷交換衝突, (18) は励起衝突である. 用いた衝突断面積は弾性衝突, 電荷交換衝突, A 準位への励起衝突における衝突断面積は, Phelps26 が実験で求めた値を用いた. 電離衝突
における衝突断面積は, Halas ら27 が実験で求めた値を用いた.

(a) t =1.0 ns(b) t =2.5 ns

Fig. 8: Net charge density at t =1.0 ns and t =2.5 ns
in the case of negative voltage input.

3.3 入力電圧が負の時の計算結果
初期の電場強度の空間分布, ポテンシャルを図 6, 数値計算の結果得られた, 放電開始直後の荷電粒子の正味の数密度の空間分布を図 7 に示す. t =0.1 ns では誘電体表面で荷電粒子の数密度が増大していることがわかる. これは図 6(a) より, 露出電極付近, 誘電体表面の電場強度が大きいため, 電子が場から大きなエネルギーを受け, 電
子雪崩を起こしたためである. t =0.2 ns では誘電体表面を含め, 電極間において荷電粒子の数密度が増大している. これは図 6(b) より, 電気力線上を加速する中で, 電子が大きなエネルギーを受けるためである. 図 7(c) より,
t =0.5 ns では誘電体表面の電子雪崩で生成された荷電粒

子と, 電極間を電子が加速する中で電子雪崩によって生成された荷電粒子が混在しており, 大きな数密度を持った荷電粒子が誘電体表面に生成される. 図 7(d) より, 図7(a)–(c) と比較して, 数密度にオーダーの違いがあるた
め, t =0.7 ns になると, ほとんどの電子が誘電体に到達し堆積したと考えられる.
電子が誘電体上に堆積した以後の, t =1.0 ns と t =2.5
ns における荷電粒子の正味の数密度の空間分布と電場強度の空間分布, ポテンシャルをそれぞれ図 8, 9, 10 に示

4Copyright c 2010 by JSFM

第 24 回数値流体力学シンポジウム
E1-5(a) t =1.0 ns(b) t =2.5 ns(a)Fig. 9: Spatial distribution of electric ﬁeld intensity at
t =1.0 ns and t =2.5 ns in the case of negative voltage
input.

(a) t =1.0 ns

Fig. 13: (a) Initial distribution of electric ﬁeld and (b)potential and electric ﬁeld vector in the case of positive
voltage input.

(b) t =2.5 ns(a) t =0.5 ns

Fig. 10: Spatial distribution of potential and electricﬁeld vector at t =1.0 ns and t =2.5 ns in the case ofnegative voltage input.

(a) t =1.0 ns(b) t =0.7 ns

Fig. 14: Net charge density at early time in the case of
positive voltage input.

子が誘電体表面を下流方向に広がるのと同時に, 大きな数密度を持つ荷電粒子が露出電極側に広がっていることがわかる. これは大きな数密度を持った荷電粒子が誘電体付近に存在することによってポテンシャルがゆがめられ, 露出電極側に電場強度の高い部分ができるためである. そのため電子が容易に加速され, 高密度の荷電粒子が存在するよりも露出電極側で電子雪崩を起こし, プラズ
マが広がると考えられる. 以上のことが続いていき, 露出
電極を覆うようにプラズマが生成される.
図 8(a), 図 9(a) より, 電場強度の高い部分が露出電極付近に存在したままプラズマが広がっていることがわかる. これはプラズマが誘電体表面で最初に生成され, ポテ
ンシャル分布を歪めたことに起因している. このことはグロー放電の特徴であるため, 入力電圧を負にした時はグロー放電によってプラズマが生成されることがわかる.
次に運動量輸送の結果, 中性粒子へ誘起される力の経
過した時間までの平均を図 11, 12 に示す. 図 11 は x 方向の力の時間平均, 図 12 は y 方向の力の時間平均である.
t =1.0 ns では x =1–1.25 mm 付近で右下方向の力が生成されている. これは電子が誘電体表面を移動していることに起因している. x =0.75 mm 付近では左上方向の力が生成されている. イオンは電気力線の方向に従い, 露
出電極へ向かうため, この方向の力が生成されたと考えられる.
t =2.5 ns では図 11, 12 よりイオン過多の部分で露出
電極方向の力が大きくなっていることがわかる. 一方, 被覆電極方向へ向かう力はほとんど生成されていない. これは電子と比較してイオンは質量が大きいため, イオンの加速による中性粒子への運動量の伝達が支配的になるためである. したがってプラズマが形成された後は, イオンを下流方向へ動かすように入力電圧を変化することが重要と考えられる.

(b) t =2.5 ns

Fig. 11: Spatial distribution of x-direction time averaged force at t =1.0 ns and t =2.5 ns in the case ofnegative voltage input.

(a) t =1.0 ns(b)(b) t =2.5 ns

Fig. 12: Spatial distribution of y-direction time averaged force at t =1.0 ns and t =2.5 ns in the case ofnegative voltage input.

す. 図 8(a) より, 図 7(d) と比較して, t =1.0 ns では, 電

5Copyright c 2010 by JSFM

第 24 回数値流体力学シンポジウム
E1-5(a) t =1.0 ns(a) t =1.0 ns(b) t =2.5 ns

Fig. 15: Net charge density at t =1.0 ns and t =2.5 ns
in the case of positive voltage input.

(a) t =1.0 ns

Fig. 19: Spatial distribution of y-direction time averaged force at t =1.0 ns and t =2.5 ns in the case ofnegative voltage input.

の数密度の空間分布を図 14 に示す. t =0.5 ns では露出
電極表面でイオンの数密度が増大していることがわかる.
これは図 13(b) より, 電子が場からエネルギーを受けて露出電極へ加速し電子雪崩を起こすが, 露出電極に到達す
ると電極に吸い込まれるためである. ここでこのイオンの数密度が電子の数密度よりも増大している領域を, イ
オン雲と呼ぶことにする. t =0.7 ns ではイオン雲が露出
電極から離れ, イオン雲と露出電極の間に, プラズマが生成されることがわかる.
露出電極表面にイオン雲が生成した以後の, t =1.0 nsと t =2.5 ns における荷電粒子の正味の数密度の空間分
布と電場強度の空間分布, ポテンシャルをそれぞれ図 15,
16, 17 に示す. 図 15(a) より, 図 14(b) と比較して, t =1.0
ns では, イオン雲を先頭にしてプラズマが下流方向に広がっていることがわかる. これは大きな数密度を持ったイオン雲の存在によってポテンシャルがゆがめられ, イオ
ン雲のさらに先頭に電場強度の高い部分ができるためである. 図 16(a), 図 17(a) において, このことが確認できる. そのためイオン雲の先頭にある領域で電子が露出電
極方向に電子雪崩を起こしつつ加速し, イオン雲が下流に進むと同時にプラズマが広がると考えられる. 以上のことが続いていき, 一本の線のようにプラズマが生成される. なおプラズマはこれ以上下流には成長しない. なぜなら図 17(b) より, ポテンシャル分布が歪められイオン雲の先端の電場強度が低くなるからである.
図 15(a), 図 16(a) より, 電場強度の高い部分とプラズマが同時に下流方向へ移動していることがわかる. これはイオン雲が最初に露出電極付近で生成され, ポテンシャ
ル分布を歪めたことに起因している. このことはストリー
マ放電の特徴であるため, 入力電圧を正にした時はストリーマ放電によってプラズマが生成されることがわかる.
次に運動量輸送の結果, 中性粒子へ誘起される力の経
過した時間までの平均を図 18, 19 に示す. 図 18 は x 方向の力の時間平均, 図 19 は y 方向の力の時間平均である.
t =1.0 ns では x =0.25–0.8 mm 付近, y =0.75–1.0 mm
付近で左上方向の力が生成されている. これは電子が露
出電極方向に移動していることに起因している. y =0.5–
0.75 mm 付近では右下方向の力が生成されている. イオ
ンが誘電体方向の向かうため, この方向の力が生成されたと考えられる.
t =2.5 ns では図 18, 19 より大部分で右下方向の力, つまり被覆電極へ向かう力が大きくなっていることがわかる. 一方, 露出電極方向へ向かう力はほとんど生成されていない. これはプラズマの生成が完了し, その後は電子と比較してイオンは質量が大きいため, イオンの加速による中性粒子への運動量の伝達が支配的になるからである.

(b) t =2.5 ns

Fig. 16: Spatial distribution of electric ﬁeld intensity at
t =1.0 ns and t =2.5 ns in the case of negative voltage
input.

(a) t =1.0 ns(b) t =2.5 ns

Fig. 17: Spatial distribution of potential and electricﬁeld vector at t =1.0 ns and t =2.5 ns in the case ofnegative voltage input.

(a) t =1.0 ns(b) t =2.5 ns(b) t =2.5 ns

Fig. 18: Spatial distribution of x-direction time averaged force at t =1.0 ns and t =2.5 ns in the case ofnegative voltage input.

まとめ
本研究では, プラズマアクチュエータにおける流れ場の誘起を再現することができるモデル化と数値計算コードの開発を行い, プラズマの生成過程と運動量輸送の結4.
3.4 入力電圧が正の時の計算結果
初期の電場強度の空間分布, ポテンシャルを図 13, 数
値計算の結果得られた, 放電開始直後の荷電粒子の正味6

Copyright c 2010 by JSFM

第 24 回数値流体力学シンポジウム
E1-5

果生じる力を計算した. 入力電圧が負の場合は, プラズマは誘電体表面で最初に生成された. 誘電体表面に堆積した荷電粒子が電位分布を歪めることでプラズマが大きく広がっていき, 最終的に誘電体を覆うように形成された.
そして, 電場強度の高い部分は露出電極付近に保たれ, グ
ロー放電によってプラズマが生成されることがわかった.
一方, 入力電圧が正の場合は, イオン雲が露出電極表面で最初に生成される. そしてイオン雲を先頭にしてプラズマが広がった. 電場強度の高い部分がイオンと共に移動するため, ストリーマ放電によってプラズマが生成されることがわかった. なお, 運動量輸送にはイオンの加速が重要となるので, プラズマ形成後はイオンの動きを考慮
した入力電圧をかけることが必要である.

14. A. V. Likhanskii, M. N. Shneider, S. O. Macheret
and R. B. Miles, “Modeling of dielectric barrier discharge plasma actuator in air”, Journal of Applied
Physics, 103, 2008, 053305.
15. D. F. Opaits, M. N. Shneider, R. B. Miles, A. V.
Likhanskii, and S. O. Macheret, “Surface charge
in dielectric barrier discharge plasma actuators”,
Physics of Plasmas, 15, 2008, 073505.
16. H. Nishida, and T. Abe, “Numerical Analysis
for Plasma Dynamics in SDBD Plasma Actuator”,
AIAA Paper 2010-4634, 2010.
17. D. M. Orlov, G. I. Font, and D. Edelstein, “Characterization of Discharge Modes of Plasma Actuators”,
AIAA Journal, 46, 2008, pp. 3142–3148.

参考文献
1. G. I. Font, C. L. Enloe, T. E. McLaughlin, and D.
Orlov, “Plasma Discharge Characteristics and Experimentally Determined Boundary Conditions for
a Plasma Actuator”, AIAA Paper 2007-188, 2007.

18. J. P. Bouef, Y. Lagmich, Th. Callegari and L. C.
Pitchford, “EHD Force in Dielectric Barrier Discharges Parametric Study and Inﬂuence of NegativeIons”, AIAA Paper 2007-183, 2007.

2. S. Roy, “Flow actuation using radio frequency
in partially ionized collisional plasmas”, Applied
Physics Letters, 86, 2005, 101502.

19. Y. Lagmich, Th. Callegari, Th. Unfer, L. C. Pitchford and J. P. Boeuf, “Electrohydrodynamic force
and scaling laws in surface dielectric barrier discharges”, Applied Physics Letters, 90, 2007, 051502.

3. T. C. Corke, “Overview of Plasma Flow Control:
Concepts, Optimization, and Applications”, AIAAPaper 2005-563, 2005.

20. A. V. Likhanskii, V. V. Semak, M. N. Shneider, D.
F. Opaits, R. B. Miles and S. O. Macheret, “The roleof the photoionization in the numerical modeling of
the DBD plasma actuator”, AIAA Paper 2009-841,
2009.

4. G. I. Font, “Boundary Layer Control with Atmospheric Plasma Discharge”, AIAA Paper 2004-3574,
2004.

21. A. V. Likhanskii M. N. Shneider, S. O. Macheret
and R. B. Miles, “Modeling of dielectric barrierdischarge plasma actuators driven by repetitive
nanosecond pulses”, Physics of Plasmas, 14, 2007,
073501.

5. K. P. Singh and S. Roy, “Modeling plasma actuatorswith air chemistry for eﬀective ﬂow control”, Journalof Applied Physics, 101, 2007, 123308.
6. 藤井 考蔵, 浅田 健吾, “DBD プラズマアクチュエー
タによる翼剥離制御メカニズム–バースト波解析が語るもの–”, ながれ, 29, 2010, pp. 259–270.

22. C. K. Birdsall and A. B. Langdon, Plasma Physics
via Computer Simulation, Institute of Physics Publishing, London, 2004.

7. 藤井 考蔵, 大山 聖, “低騒音・低エネルギーを実現するプラズマ利用の流体制御–輸送機・流体機器の設計
概念革新を目指して–”, 第 71 回応用物理学会学術講
演会 講演予稿集, 2010.

23. G. A. Bird, Molecular Gas Dynamics and Direct
Simulation of Gas Flows, Oxford Science Publications, New York, 2003.

8. 高垣 雅彦, “プラズマアクチュエータの作動原理の解
明と火星大気への適用,” 東北大学大学院修士学位論
文, 2009.

24. 南部 健一, ボルツマン方程式の確率解法, 数値流体力学–基礎と応用, 東京大学出版会, 東京, 1992.
25. H. Matsumoto and Y. Omura, Computer Space
Plasma Physics: Simulation Techiques and Software,Radio Atmospheric Science Center, Kyoto, 1993, pp.
21–65.

9. T. Abe, Y. Takizawa, S. Sato and N. Kimura, “Experimental Study for Momentum Transfer in a Dielectric Barrier Discharge Plasma Actuator”, AIAA
Journal, 46, 2008, pp. 2248–2256.

26. A. V. Phelps, “Cross Sections and Swarm Coeﬃcients for Nitrogen Ions and Neutrals in N2 and Argon Ions and Neutrals in Ar for Energies from 0.1
eV to 10 keV”, Journal of Physical and Chemical
Reference Data, 20, 1991, pp. 557–573.

10. T. Abe, Y. Takizawa, S. Sato and N. Kimura,“A Parametric Experimental Study for Momentum
Transfer by Plasma Actuator”, AIAA Paper 2007187, 2007.
11. W. Kim, H. Do, M. G. Mungal and M. A. Cappelli,“On the role of oxygen in dielectric barrier discharge
actuation of aerodynamic ﬂows”, Applied Physics
Letters, 91, 2007, 181501.

27. S. T. Halas and B. Adamczyk, “Cross sections forthe production of N2+ , N+ , and N2++ from nitrogen by electrons in the energy range 16–600eV”, International Journal of Mass Spectrometry and Ion
Physics, 10, 1972, pp. 157–160.

12. G. I. Font, C. L. Enloe, J. Y. Newcomb, A. L.
Teaque, A. R. Vasso and T. E. McLaughlin, “Effects of oxygen content on the behavior of the dielectric barrier discharge aerodynamic plasma actuator”, AIAA Paper 2010-545, 2010.

28. C. K. Birdsall, “Particle-in-Cell Charged-Particle
Simulations, Plus Monte Carlo Collisions with Neutral Atoms, PIC–MCC”, IEEE Trans. on Plasma
Science, 19, 1991, pp. 65–85.

13. C. L. Enloe, M. G. Mcharg, G. I. Font and T. E.
Mclaughlin, “Plasma-induced force and self-induced
drag in the dielectric barrier discharge aerodynamic
plasma actuator”, AIAA Paper 2009-1622, 2009.
7

Copyright c 2010 by JSFM