静電気学会誌，43, 4（2019）164-170

J. Inst. Electrostat. Jpn.

論

文

多流体モデルへ拡張された特性曲線法による放電解析
柴田 寿一

＊, 1

（2019年3月4日受付；2019年5月27日受理）

Discharge Simulation Based on an Extended Method of Characteristics
for a Multi-fluid Model＊, 1

Hisaichi SHIBATA（Received March 4, 2019; Accepted May 27, 2019）
Electroaerodynamic thrusters can propel by utilizing ionic wind generated with corona discharges, which have multi-scale
nature both in spatial temporal directions, hence it is difficult and challenging to simulate. The method of characteristics is akind of space-marching method to simulate such the discharge and characterized with lower computational costs compared withthat of time-marching methods. So far, the method is only formulated under models excluding effects of electrons, and one of
the biggest deficits is that it is impossible to determine corona onset voltage, hence current voltage characteristics, in closed
form. To overcome the deficits, the method based on one-fluid models was extended with adopting a multi-fluid model in which
the effects of electrons, hence Townsend mechanism are considered. The consideration of Townsend mechanism in the novel
method enables self-consistent discharge simulation with holding the advantages of the space-marching method. The results
obtained with the extended method showed good agreements with the results obtained with experiments and with a one-fluid
model in which the effects of electrons are neglected.

1．序文2018 年に，マサチューセッツ工科大学（ MIT ）の研究チームがコロナ放電により生じるイオン風を応用した推進器（ Electroaerodynamic thruster: EADT ）を適用した
無人航空機の動力飛行に世界で初めて成功した ．
1）

EADT は，従来のドローンに適用されてきた電動モータ
ーと比較し，40 dB 程度も静穏，かつ可動部が全く無い
ので，将来，社会問題となり得る航空機の騒音問題を解決できる可能性や，質的に新しい航空機を構成できる可

図 1 EADT の仕組みとその多階層性
Fig.1 The mechanism of EADTs and its multiscale nature.

能性を持つ．
EADT はその曲率半径比が 30 を超えるような非対称性

電極はエミッターと呼ばれ，大きい電極はコレクターとの高い電極から構成される（図 1）．曲率半径の小さい呼ばれる．これら電極間に大気圧常温空気下で高電圧を印加すると，エミッター電極を活性電極とするコロナ放キーワード：特性曲線法，コロナ放電，電気流体力学，
シミュレーション，イオン風＊
東京大学大学院新領域創成科学研究科先端エネルギー工学専攻（〒277-8561 千葉県柏市柏の葉 5-1-5）
Department of Advanced Energy, Graduate School of
Frontier Sciences, The University of Tokyo, 5-1-5
Kashiwanoha, Kashiwa-shi, Chiba, Japan現）宇宙航空研究開発機構宇宙科学研究所学際科学
研究系（〒252-5210 神奈川県相模原市中央区由野台 3-1-1）
Department of Interdisciplinary Space Science, Institute ofSpace and Astronautical Science, Japan Aerospace
Exploration Agency, 3-1-1 Yoshinodai, Chuo-ku,
Sagamihara-shi, Kanagawa, Japan1
shibata@flab.isas.jaxa.jp

電が生じる．コロナ放電によりエミッター近傍で電離したイオンはコレクターへと電場によりドリフト運動し，空気中の中性粒子と衝突を繰り返す．この衝突により運動量がイオンから空気へと伝達（イオン風の発生）され，推進力が生み出される

．

1, 2）

EADT の実用化を加速するためには，コロナ放電とそれに誘起されるイオン風現象の数値解析が有効であると考えられる．同じくコロナ放電を応用する電気集塵機の分野では，放電を記述する方程式系として，正イオンのみを考慮した 1 流体モデルがこれまで発展してきた

．

3-5）

しかし，電子を考慮しない 1 流体モデルでは電子衝突
電離反応（ α 効果）とカソードからの 2 次電子放出（ γ効果），すなわち Townsend 機構，を考慮できない．この

多流体モデルへ拡張された特性曲線法による放電解析（柴田 寿一）

165（31）

ことから，放電持続条件を指定できる境界条件が別途必においては，電子の効果を無視し，正イオンのみを考慮

要となる．典型的な境界条件として，放電開始後にコロし た モ デ ル が 構 築 さ れ て き た． そ れ ら モ デ ル で は，

ナ活性電極表面の電界が一定値に維持されるとする経験

Townsend 機構を直接に再現することはできない．本節

則（ Kaptsov の仮説）を Peek の式などのコロナ放電開では，その特性曲線法を簡単にレビューし，電子の影響

始電界の経験式と組み合わせたものを挙げられる．しかを考慮した場合への拡張へ向けて準備する．

し，EADT のように電極形状が複雑な場合は，コロナ放2.1 支配方程式

電開始電界の経験式が存在しない場合がある．加え，放正イオンのドリフト方程式（1 流体モデル）は以下で電領域の非対称性が強い場合は Kaptsov の仮説による近与えられる：

似に問題が生じる可能性がある．

（1）

これら問題の解決策の一つとして，電子の効果を流体モデルに組み入れる手法を挙げられる（例えば，2 流体

ここに μp は正イオン移動度，ρp は正イオン密度，Eモデルの導入）．そのようなモデルでは，
α 効果と γ 効果，

は電場ベクトルである．正イオン密度のみを考慮したすなわち Townsend 機構，を閉じた形で再現できる．基Poisson 方程式は以下で与えられる：

本的に電極形状に依存しない，電子衝突電離係数と二次（2）電子放出係数が判明すれば，前述の仮説と経験式に依存しない，より基本原理に立ち返った解析が可能となる．
ところが，2 流体モデルを採用した EADT の解析におここに φ はスカラ電位． 0（ = 8.8×10 F/cm ）は真空-14の誘電率である．

いては，方程式系に大きく乖離する複数の時間スケール

2.2 特性曲線法の導出が含まれることが知られている

座標 t を正イオンがドリフト運動する方向に取ると，．例えば大気圧雰囲6, 7）

気下において，電子衝突電離反応はナノからマイクロ秒

以下の式を得る：

単位で生じるが，正イオンはミリ秒から秒単位で電極間（3）をドリフト運動する．このため，2 流体モデルに時間進行法を適用し EADT の問題を解く場合には，定常ソルバーを用いない限り，100 万回程度かそれ以上の時間積

ここに x は位置座標である．式
（3）
を式
（1）
へと代入し，
連鎖則による関係式

分回数が要求される．この時間積分を多次元問題におい（4）て実行することは，計算コストの面から現実的でない場合が多い．このことから，2 流体モデルの高速解法を確を利用すると，立することが望まれる．

（5）

これまで筆者らは，2 流体モデルの高速解法を摂動法によって構成してきた．先行研究では，放電開始電8-11）を得る．さらに移動度を一定と近似すれば，

圧近傍の電流電圧特性を把握するため．摂動法を基に手法を構築した．しかし，Monrolin らによって指摘

12）

（6）

さ

れたように，低次の摂動法に基づく限り，空間電荷効果を得る．もし，境界条件としてコロナ活性電極表面（こを容易に無視できない非線形領域への応用は難しい．

こではアノード電極）での正イオン密度 ρp,anode が与えら

そこで本研究では，空間進行法の一種であり，かつ非摂

れれば，式
（6）
は

動的手法である特性曲線法による定式化に注目する．特性（7）曲線法は，計算コストが時間進行法と比較し小さく，これまで電気集塵機の分野で主に適用されてきたが，筆者と解析的に解くことができる．その境界での正イオン密の知る限り，電子の効果を組み入れた多流体モデルへの適度の代表的な与え方は，コロナ活性電極表面の電界が放用例はない．本研究では，これまで正イオンのみを考慮し

電 開 始 後 も 維 持 さ れ る と 主 張 す る 経 験 則， つ ま り

た 1流体モデルで定式化されてきた特性曲線法を電子と正Kaptsov の仮説，を応用することである．しかし，そのイオンを考慮した 2 流体モデルでの定式化へと拡張する．

場合は放電開始電界を Peek の経験式，あるいは実測値

3, 4）などから決定する必要が生じる

．

3, 4）

2．特性曲線法のレビュー

2.3 計算アルゴリズム

これまで，特性曲線法（ Method of Characteristics: MoC ）Kaptsov の仮説を採用し，かつコロナ放電開始電界を経166（32）

静電気学会誌 第43巻 第 4 号（2019）

験式から求める場合，具体的な計算アルゴリズムは以下の3.1 支配方程式

ようになる．なお，本研究で提示する結果はここで解説す

局所場近似（ Local Field Approximation: LFA ）が成立

るパラメータと手法を採用して得られたものである：

すること，及び荷電粒子の拡散を十分に無視できること

1．各電極における電位を設定を要請し，以下の電子と正イオン密度の輸送方程式をス

2．アノード電極における正イオン密度を初期値に設定

カラ電位に関する Poisson 方程式と組み合わせて採用す

3．正イオン密度の修正量が許容量以下になるまで，以

る

：

6-11）

下を反復

（9）

3（a）
．電荷方程式（式
（7）
）を求解
3（b）
．Poisson 方程式（式（2））の右辺に求解した電荷

（10）

を設定

（11）

3（c）
．Poisson 方程式を全空間で解き，電場を求解
3（d）．アノード電極表面における正イオン密度の修正量を決定し，修正（12）
ここで ρe は電子密度，μe は電子移動度，α は電子衝突

以下，アルゴリズムの要点を説明する．まず 1. について，電位は放電開始電圧を超えている必要がある．次電離係数であり電場強さの関数である．
先行研究

6-11）

にならい，境界条件は以下で与える．アに，2. について，アノード電極における正イオン密度のノード電極では正イオン密度がゼロであることを要請す

初期値は現実的な値（例えば，10 cm ）を設定する．

る：

6

-3

3. の収束判定ついて，本研究では，修正量と修正される（13）

べき正イオン密度との比率が 10 以下になるときにルー-6プを打ち切ると約束する．3（a）は単純な代数計算であカソード電極で二次電子放出を考慮する：

る．3（b）
について，
座標 t と x の座標変換を行いながら，（14）

MoC で用いる格子から Poisson 方程式で用いる格子へと物理量を補間する．3（c）について，具体的には有限差分法で Poisson 方程式を空間離散化し，離散化された線

ここで γ は二次電子放出係数，は移動度比である．

形方程式を定常反復法の一種である SOR 法（緩和係数3.2 特性曲線法の導出

を 1.8 とする）で解くことができる．3（d）
について，任まず，電子密度が特性曲線上で従う式を求める．式（9）

意の特性曲線において，コロナ活性電極表面の電場がコの左辺を，式（3）及び式（12）を代入し変形する：

ロナ放電開始電界に保たれるように，任意の特性曲線に沿って（8）と更新する．ここに ω（ =1×10 ）は緩和係数，Eon は-15

コロナ放電開始電界であり実測値から決定する．Eanode は

これより，

コロナ活性電極表面での電場強さである．ここで，コロ（15）

ナ活性電極表面（アノード電極）における正イオン密度が上昇すると，アノード電極表面における電場強度は緩を得る．

和され，Eanode の値は減少する．このことから，意味のあ次に，正イオン密度が特性曲線上で従う式を求める．

るフィードバック・ループを構成するためには，ω を正

このため，電子の場合と同様に，式（10）の左辺を変形す

値に設定する必要がある．

る：

3．特性曲線法の拡張
前節では，1 流体モデルに特化した特性曲線法の導出にあたり，電子を無視した．本節では，電子の効果（ α，γ 効果）を考慮できる 2 流体モデルに対応する特性曲線法を導く．

多流体モデルへ拡張された特性曲線法による放電解析（柴田 寿一）

一方，楕円型の偏微分方程式である Poisson 方程式（式

これより，
（16）

（12）
）は Deutsch の仮定

14）などの特別な仮定を置かない

限り，計算量の削減が困難である．このことから，MoC
及び Ex-MoC の計算量は，Poisson 方程式の評価回数に

を得る．
次に，

167（33）及びと定義し，式（15）と式（16）おおよそ比例すると考えられる．
具体的な Ex-MoC の計算アルゴリズムは以下の通りでに代入した方程式群，
（17）
（18）
（19）

ある：
1．各電極における電位を設定
2．アノード電極における電子密度を初期値に設定
3．電子密度の修正量が許容量以下になるまで，以下を反復を解く．アノード電極における境界条件（式
（13）
）を用

3（a ）
．電荷方程式（式
（20）
－
（22）
）を求解

い式
（17）
－
（19）を非斉次 1 階線形常微分方程式の解法に

3（b）．Poisson 方程式（式（12））の右辺に求解した電荷を設定

従って手計算で同値変形すると，
（20）
（21）
（22）

3（c）
．Poisson 方程式を全空間で解き，電場を求める
3（d）．アノード電極における電子密度の修正量の決定と修正
以下，アルゴリズムの要点を説明する．まず 1. について，電位は放電開始電圧を超えている必要がある．2．

を得る．式
（20）
－
（22）を式（17）－（19）へと代入すると，について，アノード電極における電子密度の初期値は現解の妥当性を確認できる．これに加え，カソード電極に実的な値（例えば，
10［ cm ］）を設定する．3．について，おける境界条件である式（14）を利用すると，持続放電の修正量と修正されるべき電子密度との比率が 10 以下に条件式なるときにループを打ち切ると約束する．3（a）につい
（23）

6

-3

-12

て，積分は

及び

の 2 つの量に対して実行する．3（b）について，電荷方程式（式
（20）
－
（22）
）の解

を導出できる．
式（23）は次節で述べる 2 流体モデルによる計算アルゴは電荷密度ではない点に留意して計算する．3（c）につリズムの中で修正量の決定に用いられ，1 流体モデルにいて，具体的には有限差分法や有限要素法などの空間離よる計算アルゴリズムに適用される Kaptsov の仮説に対散化が適用できる．3（d）について，修正量の決定と修応する．

正は以下の式で行う．カソード電極における二次電子放なお，式（23）を Laplace 場に適用すると，初期の放電経路を表す特性曲線が事前に分かる場合に限り，放電開出を表す境界条件が任意の特性曲線で満たされるように，各特性曲線に沿って，始電圧と二次電子放出係数を関連付けられる ．
13）

（24）

3.3 計算アルゴリズム
拡張された特性曲線法（ Ex-MoC ）により定式化された方程式群，式（20）－（22）は積分を含む代数方程式であと更新する．ここに ω（ =1×10 ）は緩和係数である．
-13

ここで，アノード電極における電子密度が上昇すると，るが，元の偏微分方程式群，式（9）－（11）と比べ容易にアノード電極表面における電場強度は緩和され，解ける．この事実は 3 次元空間で定義された偏微分方程が増大すると考えられる．このことから，意味のあるフ式の計算を 1 次元空間で定義された代数方程式に変換でィードバック・ループを構成するためには，ω を正値にきることを意味しており，特に 3 次元の実問題において設定する必要がある．

の値

計算量の削減に貢献できると考えられる．ここに，冒頭で述べたように，少なくとも EADT の解析に用いられ

4．問題設定る 2 流体モデルは，一般に硬い（ Stiff な）方程式であり，拡張された特性曲線法（ Ex-MoC ）の検証を，1 流体

1 流体モデルと比較しその求解に膨大な計算量を要求す

モデル，2 流体モデルによる数値解析結果，そして実験ることに留意されたい．

データが存在する先行研究

15）にならい，大気圧常温空168（34）

静電気学会誌 第43巻 第 4 号（2019）

気下（298 K ）における同軸円環電極間に生じる正コロ

5．結果及び議論

ナ放電で実行する．ここでは簡単化のために対称性の高

本節では，Ex-MoC（ここでは 2 流体モデル）の結果をい電極形状を選択するが，Ex-MoC は任意形状電極間に示し，MoC（従来の 1流体モデル）の結果と比較検討する．

生じるコロナ放電を扱えるように定式化されていること

5.1 2流体モデルの結果に留意されたい．

印加電圧を 40 kV に設定し，Ex-MoC でのシミュレー

4.1 計算パラメータの設定

ションを実行した．図 2 にシミュレーションで得られたにならい，輸送係数を以下のように設定

電荷密度分布を示す．電子密度は，コロナ活性電極（アする．電子移動度を 5.00×10 m V s ，正イオン移動度

ノード）近傍では正イオン密度に比べて優勢であるが，を 1.86×10 m V s ，電子衝突電離係数を電場強さの関カソードに近づくにつれて急速に減衰する．一方，正イ

先行研究

15）

-2

-4

2

2

-1 -1

-1 -1

数として，

と設定す

オン密度は，アノードではゼロであるが，その後極大値を取った後，カソードに近づくにつれて減衰する．カソる．ここに電場強さ |E| の単位は V/m である．2 流体モードでの正イオン密度は電子密度の約 10 万倍である．

デルでパラメータとなる二次電子放出係数（ γ ）は，1

これは移動度比（ μ ）と二次電子放出係数（ γ ）の積と流体モデルと可能な限りフェアに予測性能を比較するた

同じオーダーであり，境界条件（式（14））と一貫した結め，実験から得られた放電開始電圧を再現できるように果が得られている．

逆算し，4.0×10 と設定する．また，同軸円環電極の形-3状は内径 0.07 cm，外径 10.35 cm である．本研究では，実験で適用された 3 次元構造の電極形状を 1 次元軸対称のそれに帰着させて解く．

4.2 計算格子の設定
アノード電極及びカソード電極近傍で格子が密になる
不等間隔格子を採用する． 1 流体モデルでは，シースを解像する必要がない．しかし，2 流体モデルでは，精度よく計算を実行するためには，シースを解像する必要がある．ここにデバイ長さは，電子温度を 1 eV，コロナ
活性電極表面近傍での電子密度を 10 cm と見積もると，
7-3220 μm のオーダーである．このことから，22 μm 程度の格子分解能が必要と考えられる．本研究では，この事
実，そして格子収束性を考慮しながら節点数を検討し，
1 流体モデル及び 2 流体モデルの両者において，収束性が十分に確認できた 2,000 節点と設定し，最小格子幅は

図 2 Ex-MoC により得られた電荷密度分布
Fig.2 Distribution of charge densities obtained by Ex-MoC.

5.2 1流体モデル，2流体モデル，実験結果の比較

共に安全側で約 1 μm と設定した．一般に，シースを解

印加電圧を 40 kV に設定した場合の空間電荷分布及び像する必要がない観点から，1 流体モデルで用いる格子

電場分布を比較検討する．図 3 に 1 流体モデル（ MoC ）は 2 流体モデルのそれより粗くとれると考えられる．し及び 2 流体モデル（ Ex-MoC ）を用いて計算した正イオ

かし，本研究における問題設定においては，コロナ活性

ン分布の結果を示す．カソード電極近傍では，物理量分電極近傍の急峻な電場分布を精度よく解像する観点か

布が定性的に一致しているため，結果は類似したものにら，このシースを解像するために必要な格子と同程度のなると考えられる．しかし，電子を陽に考慮したので，幅の格子が双方のモデルで必要になるので，注意が必要

アノード電極近傍で両者の分布は大きく乖離している．

である．

図 4 に 1 流体モデル及び 2 流体モデルを用いて計算され座標系は動径方向に r を取る．ここに座標 t と座標 rとの関係は

い一致を見ている．1 流体モデルと 2 流体モデルにおけ（25）であり，計算過程で利用する．

た電場を示す．両者モデルにおける電場の計算結果はよる放電現象のモデル化はシース領域の考慮の有無を含め大きく異なるが，類似する電場分布が得られることは注目に値する．電流密度の比較結果を図 5 に示す．1 流体
モデル，2 流体モデルの相対誤差は 30 kV から 50 kV の

多流体モデルへ拡張された特性曲線法による放電解析（柴田 寿一）

169（35）

5.3 イオン風現象の予測精度
イオン風を駆動する体積力（クーロン力）を以下の式10）で計算し，その空間分布を図 6 に示す．
F=（ρp - ρe ）
・E

（26）

コロナ活性電極の半径の 2 倍程度までの領域においては，1 流体モデルで得られる体積力が 2 流体モデルのそれを上回る．一方，計算空間の大部分を占めるそれ以外の領域においては，両者の結果は定量的に一致している．
このことから，本研究で扱った問題に限れば，体積力予測に 1 流体モデルを用いても問題がないと考えられる．
図 3 MoC 及び Ex-MoC から計算された正イオン密度の分布
Fig.3 Distribution of positive ion density obtained with MoC and
Ex-MoC.

図 6 体積力分布のモデルによる比較
Fig.6 Comparison of body force distribution between the different
models.

図 4 MoC 及び Ex-MoC から計算された電場の分布
Fig.4 Distribution of electric field obtained with MoC and
Ex-MoC.

5.4 計算コスト
図 7 に，
収束までに要した SOR 反復回数の総計を示す．
多くの場合で，1 流体モデルと比較して 2 流体モデルでは収束までに多くの反復回数を要するが，アノードにおける電子密度の初期値が収束値に近い場合は，直ちに収束することが分かる．
陽解法による時間進行法で本問題を解いた場合，1 流体モデルにおいては，正イオンがギャップ d ~ 0.1 m の電極間をドリフト速度 vp~10 m/s で移動する時間尺度 τdp3

はおおよそ τdp ≈ d/vp ≈ 0.1/10 =100 μs である．空間格子
3幅 Δx = 1 μm，CFL 数を 0.1 と設定すれば，時間刻み幅の最大値は Δt = CFL・Δx /vp ≈ 100 ps と概算できる．このため，少なくとも 100 万回の時間積分が必要と考えられる．一方，2 流体モデルの場合は電子ドリフト現象を含
図 5 実験，MoC，及び Ex-MoC で得られた電流電圧特性
Fig.5 Current-Voltage Characteristics obtained with theexperiments, MoC, and Ex-MoC.

むので，さらに制約が厳しくなる．本研究では，電子の移動度は正イオンのそれの約 300 倍であるので，CFL 条
件から少なくとも 3 億回の時間進行が必要と考えられる．また，Poisson 方程式ソルバーの呼び出し回数は，1

区間では 2％以内である．加え，実験データともよい一流体モデルの場合と比較して約 300 倍となり，その呼び致を見ている．

出し回数の下限は 3 億回と概算できる．

170（36）

静電気学会誌 第43巻 第 4 号（2019）
ら，より複雑な電極でのベンチマークを行う．
・本研究では，Ex-MoC を正イオンと電子を考慮した 2
流体モデルで定式化したが，応用先を広げるため，負イオンや準安定種などを加えた定式化を試みる．
謝辞
本研究を遂行する機会を与えてくださった，東京大学の鈴木宏二郎教授に感謝いたします．
参考文献

図 7 全 SOR 反復回数のモデルによる比較
Fig.7 Comparison of the total number of SOR iteration between
the different models.

6．結論
正イオンのみが考慮される従来の特性曲線法を拡張し，正イオンに加えて電子の影響を考慮できる新しい特性曲線法を定式化した．同軸円環電極内部に生じるコロ
ナ放電を対象に，正イオンのみを考慮した従来の特性曲線法（ MoC ），および電子も考慮できる特性曲線法（ ExMoC ）での予測結果を比較検討した．これにより，以下の知見を得た．
・Ex-MoC では最終的に 2 本の常微分方程式が得られ，これらは MoC と同じく解析的に解ける．
・Ex-MoC を用いた数値解析により，少なくとも本研究で扱った条件下では，電流予測に際し，実験結果をよく説明できることが分かった．
・カソード電極近傍での正イオン密度の空間分布について，MoC と Ex-MoC による解析結果の間で定性的な差異は見受けられない．しかし，コロナ活性電極近傍におけるその分布は大きく異なる．
・Ex-MoC と MoC の間で，電場分布は大局的によい一致を見た．
・Poisson 方程式ソルバーに含まれる SOR 反復ルーチンの呼び出し回数を基準として考えた場合，特性曲線法
かつ 2 流体モデル（ Ex-MoC ）は特性曲線法かつ 1 流体モデル（ MoC ）と比較してより高い計算コストを要する．一方，時間進行法かつ 2 流体モデルの計算コストは，時間進行法かつ 1 流体モデルのそれと比較し，電子移動度を正イオン移動度で除算した程度，つまり少なくとも約 300 倍の計算コストを要する．
今後は以下に挙げる研究を遂行する計画である．
・少なくとも定式化の観点からは，電極形状を制限していないので，Ex-MoC は MoC と同じく任意形状電極
間の放電解析に応用できると考えられる．このことか

1） H. Xu, Y. He, K. L. Strobel, C. K. Gilmore, S. P. Kelley, C. C.
Hennick, and S. R. H. Barrett: Flight of an aeroplane withsolid-state propulsion. Nature, 563（2018）532
2） F. Plouraboué: Flying with ionic wind. Nature, 563（2018）476
3） K. Adamiak: Numerical models in simulating wire-plateelectrostatic precipitators: a review. J. Electrostat., 71（2013）,
673
4） K. Adamiak and P. Atten: Simulation of corona discharge inpoint–plane configuration. J. Electrostat., 61（2004）85
5） 伊藤航平，森 悠真，瑞慶覧章朝，川田吉弘，田岡智浩，
柴田憲司：線対平板型電気集塵機におけるイオン風の解析と測定．静電気学会誌，43（2019）25
6） U. Ebert, C. Montijn, T. M. Briels, W. Hundsdorfer, B.
Meulenbroek, A. Rocco, and van E. M. Veldhuizen: The
multiscale nature of streamers. Plasma Sources Sci. Technol.,15（2006）S118
7） S. T. Surzhikov, and J. S. Shang: Two-component plasma
model for two-dimensional glow discharge in magnetic field.
J. Comput. Phys., 199（2004）437
8） H. Shibata: Electric breakdown detection of gases based on
global linear stability analysis. University of Tokyo Ph.D.
thesis（2016）9） H. Shibata, Y. Ohmichi, Y. Watanabe, and K. Suzuki: Global
stability analysis method to numerically predict precursor of
breakdown voltage. Plasma Sources Sci. Technol. 24（2015）
05501410）H. Shibata, Y. Watanabe, and K. Suzuki: Performance
prediction of electrohydrodynamic thrusters by the
perturbation method. Phys. Plasmas, 23（2016）053512
11） H. Shibata and R. Takaki: A novel method to predict current
voltage characteristics of positive corona discharges based ona perturbation technique. I. Local analysis. AIP Adv., 7（2017）
11502612）N. Monrolin, O. Praud, and F. Plouraboué: Revisiting thepositive DC corona discharge theory: Beyond Peek's andTownsend's law. Phys. Plasmas, 25（2018）063503
13）L. K. Warne, and R. E. Jorgenson: Ionization CoefficientApproach to Modeling Breakdown in Nonuniform Geometries.
SANDIA Report（2003）SAND 2003-4078
14）A. Bouziane, K. Hidaka, M. C. Taplamacioglu, and R. T.
Waters: Assessment of corona models based on the Deutsch
approximation. J. Phys. D: Appl. Phys., 27（1994）320
15）Y. Zheng, Bo. Zhang, and J. He: Current-voltage characteristics
of dc corona discharges in air between coaxial cylinders. Phys.
Plasmas 22（2015）023501