― 32 ―SENAC Vol. 54, No. （1 2021. 1）［共同研究成果］

金蒸気－金ナノ粒子変換を伴う高エンタルピ ー水プラズマジェットの数値シミュレーション
茂田 正哉
大阪大学 接合科学研究所

金ナノ粒子の大量創製プロセスにおける熱・物質輸送過程の解明を目的とし，水を用いたプラズマジェットとその境界領域において蒸気分子からエアロゾル粒子へと集団的な変換過程にある金ナノ粒子群の対流・拡散輸送の数値シミュレーションを，東北大学サイバーサイエンスセンタ
ーのスーパーコンピュータ SX-ACE の 512 ノードを利用した並列計算によって実現した．その結果，水プラズマジェットへの周囲気体の巻込みとともに 2 層構造の渦輪列が生じ，金蒸気分子および金ナノ粒子はそれぞれ異なる対流スケールによって輸送されることが初めて明らかとなった．

1.

はじめに直径がナノメートルスケールの超微粒子（ナノ粒子）はバルクの材料とは大きく異なる物質機能を示すことが知られており，特に金ナノ粒子は生体親和性が高く，光線力学的療法，光温熱治療，X 線 CT イメージング，ドラッグデリバリー，バイオセンサーへ応用が期待されている[1]．
近年，ナノ粒子の大量創製を可能にするツールとして，高いエネルギー密度と化学的活性を有し，さらに外部電磁場によって制御可能[2, 3]な流体である高エンタルピープラズマに注目が集まっており，精力的に研究が進められている[4]．しかし高温の発光体である熱プラズマの流動場について実験によって得られる情報には限りがあるため，プロセスの効率的な制御のために必要な知見は今尚少ない．また流動場だけでなく，ナノ粒子の形成過程はマイクロ秒～ミリ秒の現象であるため，その集団的形成のメカニズムや輸送現象を直接計測することも困難である．そのため，高エンタルピープラズマによるナノ粒子創製プロセッシングは依然として現場の技術者や研究者の経験に依っているところが多く，また莫大な時間とコストを要しているのが現状である．そこで
本研究では，金ナノ粒子の大量創製プロセスにおける熱・物質輸送現象の解明を目的として，身近な物質であり，かつ熱容量が大きい水を用いたプラズマジェットとその境界領域において蒸気分子からエアロゾル粒子へと集団的な変換過程にある金ナノ粒子群の対流・拡散輸送の数値シミュレーションを行なう．

Steady parabolic
distributions withNozzle wall(300 K)+ Water(0.1 g/min) (plasma)

Au vaporyx

Max.

5754 KMax.

1119 m/s are given.
1.8

O

Plasma jet(Non-transferred)

3.6

z
Au

nanopowder
Nucleation
Condensation(Unit: mm)
Coagulation

Fig. 1 Computational domain

金蒸気―金ナノ粒子変換を伴う高エンタルピー水プラズマジェットの数値シミュレーション

2.

仮定および支配方程式

2.1

高エンタルピープラズマ－非電離気体共存系の熱流動場― 33 ―

高エンタルピープラズマ流の通常の生成条件では，圧力は大気圧と同程度で，プラズマおよび非電離気体を含む流体全域にわたって局所熱平衡が成り立ち，また光学的に薄いと仮定できる．
このとき支配方程式は以下のような質量・運動量・エネルギーに関する保存式となる．

    u   0t(1)u    uu   p     u   utr  2   uI  t3 (2)h    uh       h   p  u  p  q rad  qcon  Φt
 Cp t(3)

ここで，は密度，u は速度ベクトル，t は時間，p は圧力，は粘性係数，I は単位行列，は導電率，E は電場ベクトル，B は磁束密度ベクトル，h はエンタルピー，は熱伝導率，Cp は定圧比熱，qrad は放射損失，qcon はナノ粒子の凝縮熱，は粘性散逸である．また tr は転置を意味する．

2.2

蒸気からナノ粒子への集団的変換および輸送

高エンタルピープラズマプロセスにおけるナノ粒子形成過程の概要は次の通りである．
① プラズマの高温場において原料が蒸発する② その原料蒸気は温度低下に伴い過飽和状態となる③ 多数の臨界核が生成する（均一核生成）④ その臨界核に原料蒸気が凝縮することでナノ粒子が成長する（不均一凝縮）⑤ 同時にナノ粒子同士も衝突・合体してより大きなナノ粒子となる（粒子間凝集）
ナノ粒子は，時間スケールの異なる均一核生成や不均一凝縮を経るのみならず，２～３桁に及ぶ直径差を持つ多数の粒子同士が衝突し合体しながら集団として成長していく．これまで分子動力学に基づいた数値計算も行われてはいるものの，現在のコンピュータ性能の限界から，数十個の核の生成過程を数ナノ秒間分ほど追跡することしかできていないため[5]，ナノ粒子群全体の成長を取り扱うことは実質的に不可能である．そこでエアロゾル学に基づく理論的および数値的なアプローチが有効とされている．
本稿では，簡潔なモデルによりナノ粒子の集団的形成過程を表現するために，ナノ粒子は局所的には同じ粒径を持つ球体であるとする．また帯電の効果は無視し，粒子温度は周囲の流体の温度と等しいとする．このとき支配方程式は以下のように記述できる[3, 6-8]．

np
n 
 n 
  np / 6 1/ 6
   u   p      D p  p   J  2 2  0 n11f     K th ln T pt       (4) f  f  f f   u        D p     Jg c   0 ( nv  n s ) n 1p/ 3 f 2 / 3     K th  ln T t     
(5)



― 34 ―

SENAC Vol. 54, No. （1 2021. 1）

n 
 n 
  nv    u   v      Dv  v   Jg c   0 (nv  n s )n1p/ 3 f 2 / 3t     
(6)

ここで，n は数密度，D は拡散係数，J は均一核生成率[9]，Kth は熱泳動係数[10]，T は温度，g は１つのナノ粒子に含まれるモノマー数の平均である．添え字 p，v，c および s はそれぞれ粒子，蒸気，臨界状態および飽和状態を表している．また f は次のように定義される変数である．

f = npg(7)

また 0 は衝突頻度に関するパラメーターであり，体積 v および質量 m を用いて以下のように表される．
1/ 6

 3vv  4 0  

6k BTvvmv(8)

kB はボルツマン定数である．式(4)および式(5)の右辺第 4 項は熱泳動を表している．式(4)～(6)の右辺に含まれる粒子成長に関わる項の詳しい導出については文献[7, 11]を参照されたい．

3.

計算手法

本研究では，水プラズマと非電離気体が相互作用しながら同時に存在する熱流動場を取り扱わなければならない．両者の間には大きな温度差に起因して粘性係数・熱伝導率・定圧比熱・導電率といった物性値だけでなく密度にも大きな差が生じている[12]．その一方で，流れ場におけるマッハ数は 10-3～10-2 のオーダーにあり，工学的に有意な時間スケールでの流体運動を捉えるためには当該の系を「大きな密度変化を伴う非圧縮性流れ」として取り扱うべきである．すなわち，プラズマと周囲の低温気体との間に生じる速度，温度，密度，ナノ粒子濃度の急激な空間勾配を捉えながら，さらに時間ステップ幅を大きく取っても数値計算を安定的に進められるような計算手法が必要となる．そこで対流項をハイブリッド型 K-K スキーム[13]により，時間微分項を３次精度 Adams-Moulton-Bashforth 法により差分化し，改良型 PISO 法[14]と組み合わせることによって，上述の数値計算を実現することとする．高エンタルピープラズマ流動のシミュレーションに対するこれらの手法の有効性については文献[15]を参照されたい．なお，拡散項，圧力勾配項，熱泳動項には２次精度中心差分を用いる．

4.

計算条件

Fig. 1 に計算領域を示す．内径 1.8 mm，外径 3.6 mm の噴出孔から最高温度 5754 K，最高速度1119 m/s の放物線型の分布をもつプラズマジェットが，400 K・大気圧の非電離水蒸気雰囲気の計算領域へ定常的に噴出する．ナノ粒子の原料である金は既にプラズマ生成部で蒸気になっているものとしてプラズマジェットと共に 0.1 g/min で供給される．
これらの条件の下，３次元の計算領域を設けた．座標系の原点をプラズマジェット噴出孔の中心に取り，各軸方向にそれぞれ 0.02 mm の幅を持つスタガード格子を用いて，時間ステップ幅を0.02 ms として計算を進めた．
本計算は，東北大学サイバーサイエンスセンターの SX-ACE において MPI 並列によって 512ノードを使用して実現した．

金蒸気―金ナノ粒子変換を伴う高エンタルピー水プラズマジェットの数値シミュレーション― 35 ―

t = t0 + 0.64 ms

t = t0 + 0.64 ms

t = t0 + 1.28 ms

t = t0 + 1.28 ms

t = t0 + 1.92 ms

t = t0 + 1.92 ms

t = t0 + 2.56 ms

t = t0 + 2.56 ms

Fig. 2 Temperature distributions on x = 0 plane
Fig. 3 Isosurfaces of the second invariant of velocity
and isotherms of 2,000 K and 3,200 K in y < 0gradient tensor of 0.001 (normalized)
Isolines of Q on y = 0 planet = t0 + 2.56 ms

Fig. 4 Isolines of the second invariant of velocity gradient tensor on y = 0 plane

5.

計算結果

Fig. 2 に，ある時刻 t0 から 0.64 ms，1.28 ms，1.92 ms，2.56 ms 経過した温度場を示す．Fig. 3 には同時刻における渦構造を示す．渦構造は速度勾配テンソルの第 2 不変量を噴出孔出口での平均流速と直径で無次元化した等値面 Q = 0.001 で表現されている．また Fig. 4 に時刻 t0 から 2.56 ms
後の断面 y = 0 における渦構造を示す．水プラズマジェットと非電離気体の境界領域でKelvin-Helmholtz 不安定性に起因する渦生成とそれによる巻込みが生じている．温度分布の模様を流脈として捉えると境界領域における渦列の層は 1 層であるように見えるが，Fig. 3 および Fig. 4

― 36 ―

SENAC Vol. 54, No. （1 2021. 1）

t = t0 + 0.64 ms

t = t0 + 0.64 ms

t = t0 + 1.28 ms

t = t0 + 1.28 ms

t = t0 + 1.92 ms

t = t0 + 1.92 ms

t = t0 + 2.56 ms

t = t0 + 2.56 ms

Fig. 5 Number density distributions of gold vapor
molecules on x = 0 plane and isosurfaces of
3.2×1020 and 2.0×1021 in y < 0Fig. 6 Number density distributions of gold
nanoparticles on x = 0 and isosurfaces of1.0×1019 and 1.0×1020 in y < 0から高温のプラズマジェットの内部から境界領域にかけて太い渦輪列が生じ，その周囲の低温領域には細い渦輪列が生じるという，2 層構造となっていることがわかる．
Fig. 5 および Fig. 6 に金蒸気分子および金ナノ粒子それぞれの数密度分布を示す．水プラズマジ
ェット内部は高温であるため金ナノ粒子は存在せずに分子の状態をとっているが，プラズマ－非電離気体の境界領域でナノ粒子へと変換される．そのため金ナノ粒子群はプラズマ外に存在するが，上流域に多く存在するナノ粒子の数が下流域で減少している．これは拡散による希釈効果だけでなく，粒子間凝集（衝突・合体）による集団成長の表れでもある．
Fig. 5 において金蒸気の分子数密度が島状に高くなっている領域が見られるが，これは Fig. 2における周囲の低温流体の巻込みによる温度低下によって分子群が濃化しているためである．Fig.
3 および Fig. 4 に示したプラズマ内外の渦輪列の 2 層構造から，金蒸気分子は大きな渦によって，金ナノ粒子は小さな渦によって，それぞれ異なる対流輸送スケールを持つことが明らかになったといえる．

6.

まとめ

金ナノ粒子の大量創製プロセスにおける熱・物質輸送現象の解明を目的として，水を用いたプラズマジェットとその境界領域において蒸気分子からエアロゾル粒子へと集団的な変換過程にある金ナノ粒子群の対流・拡散輸送の数値シミュレーションを行なった．東北大学サイバーサイエ

金蒸気―金ナノ粒子変換を伴う高エンタルピー水プラズマジェットの数値シミュレーション― 37 ―

ンスセンターのスーパーコンピュータ SX-ACE の 512 ノードを利用した並列計算によって数値シミュレーションが達成された．その結果，水プラズマジェットへの周囲気体の巻込みとともに 2層構造の渦輪列が生じ，金蒸気分子および金ナノ粒子はそれぞれ異なる対流スケールによって輸送されることが初めて明らかとなった．

謝辞
本研究は，東北大学サイバーサイエンスセンターのスーパーコンピュータを利用することで実現することができました．計算コードの高速化およびデータ整理にあたっては同センター関係各位よりご協力をいただきました．シミュレーションに要する費用の一部は科学研究費補助金（基盤研究(B)：19H01887）によって賄われました．水プラズマの物性値モデルおよび条件設定にあたっては，渡辺隆行先生（九州大学 教授）より多大なる協力をいただきましたので，ここに感謝の意を表します．

参考文献[1] Elahi, N., Kamali, M., Baghersad, M.H., “Recent biomedical applications of gold nanoparticles: Areview,” Talanta, 184 (2018), pp. 537-556.
[2] Sato, T., Shigeta, M., Kato, D., Nishiyama, H., “Mixing and magnetic effects on a nonequilibriumargon plasma jet,” Int. J. Thermal Sci., 40 (2001), pp. 273-278.
[3] Shigeta, M., “Numerical Study of Axial Magnetic Effects on a Turbulent Thermal Plasma Jet forNanopowder Production Using 3D Time-Dependent Simulation,” J. Flow Control, Measure. Visual., 6(2018), pp. 107-123.
[4] Shigeta, M. and Murphy, A.B., “Thermal plasmas for nanofabrication,” J. Phys. D: Appl. Phys., 44(2011), pp. 174025-(16 pages).
[5] Lümmen, N. and Kraska, T., “Homogeneous nucleation and growth in iron-platinum vapour
investigated by molecular dynamics simulation,” Euro. Phys. J. D, 41 (2007), pp. 247-260.
[6] Shigeta, M., “Simple nonequilibrium model of collective growth and transport of metal nanomist in athermal plasma process,” Theor. Appl. Mech. Jpn., 63 (2015), pp. 147-154.
[7] Shigeta, M., “Modeling and Simulation of a Turbulent-like Thermal Plasma Jet for Nanopowder
Production,” IEEJ Trans. Electrical Electronic Eng., 14, (2019), pp. 16-28.
[8] Shigeta, M., “Simulating Turbulent Thermal Plasma Flows for Nanopowder Fabrication,” Plasma
Chem. Plasma Process., 40 (2020), pp. 775-794.
[9] Girshick, S.L., Chiu, C.P. and McMurry, P.H., “Time-dependent aerosol models and homogeneousnucleation rates,” Aerosol Sci. Tech., 13 (1990), 465-477.
[10] Talbot, L., Cheng, R.K., Schefer, R.W. and Willis, D.R., “Thermophoresis of particles in a heatedboundary layer,” J. Fluid Mech., 101 (1980), pp. 737-758.
[11] Nemchinsky, V.A. and Shigeta, M., “Simple equations to describe aerosol growth,” Modelling Simul.
Mater. Sci. Eng., 20 (2012), pp. 045017-(11 pages).
[12] Aubreton, J., Elchinger, M. F., Vinson, J. M., “Transport Coefficients in Water Plasma: Part I:Equilibrium Plasma,” Plasma Chem. Plasma Process., 29 (2009), pp. 149-171.

― 38 ―

SENAC Vol. 54, No. （1 2021. 1）[13] Komurasaki, S., “A Hydrothermal Convective Flow at Extremely High Temperature,” 7th Int. Conf.
Comp. Fluid Dynamics, (2012), ICCFD7-3001.
[14] Oliveira, P.J. and Issa, R.I., “An improved PISO algorithm for the computation of buoyancy-drivenflows,” Numer. Heat Transfer B 40 (2001), pp. 473-493.
[15] Shigeta, M, “Turbulence modelling of thermal plasmas flows,” J. Phys. D: Appl. Phys., 49 (2016), pp.

493001-(18 pages).