（公財）日本板硝子材料工学助成会，40（2022）

分子動力学法による環境負荷の少ない次世代燃料電池の固体電解質の研究−イオン間相互作用とダイナミカルな構造の解明−
長岡工業高等専門学校

一般教育科

松永茂樹

Research on Solid Electrolytes of Next-Generation Fuel Cells with Less Environmental
Load by Molecular Dynamics-Elucidation of Ionic Interactions and Dynamic StructureShigeki Matsunaga
Division of General Education, Nagaoka National College of Technology
次世代燃料の電解質として期待される超イオン導電体、及びその融体である溶融塩、さらに複雑な分子集合体が形成される溶液の構造や輸送現象などの物性について、分子動力学シミュレーションによる研究を継続して行ってきた。水素を燃料とする燃料電池の電解質として期待される無機固体酸の融体およびガラス中で、遮蔽された BMH 型のポテンシ
ャルを用いて分子動力学シミュレーションによる研究を行い、ガラス転移に伴って構造や輸送係数だけでなく、水素原子の結合に寄与する電子状態も変化することを見出した。また、溶融硫酸塩混合系 Li 2 SO 4 -Na 2 SO 4 系について、やはり遮蔽された BMH 型のポテン
シャルを用いて分子動力学シミュレーションを行い、イオン伝導機構の研究を行った。構造やイオン相関角の詳細な解析から、Li と Na の拡散機構に顕著な違いがあることが認められた。さらに、超イオン導電体 AgBr-CuBr 系に関して分子動力学シミュレーションを行い、Ag と Cu イオンの分布の差異に関して、単位格子面上および単位格子内部について詳細に考察した。また、電解質に関する発展的研究として、生体燃料電池の電解質に用いられる乳酸−ピルビン酸水溶液における分子動力学シミュレーションによる研究を行い、希薄濃度側での誘電率が異常に増大し極大を持つこれまで知られていない現象を見出し、その原因について議論した。
We discuss the physical properties of superionic conductors, which are expected aselectrolytes for next-generation fuels, molten salts that are their melts, and the structure andtranspor t phenomena of solutions that form complex molecular aggregates. We havecontinued my research by molecular dynamics simulation. A molecular dynamics simulation
study using the shielded BMH-type potential was conducted in a melt of an inorganic solidacid expected as an electrolyte for hydrogen-fueled fuel cells and in glass. We found that not
only the transport coefficient but also the electronic state that contributes to the bonding ofhydrogen atoms changes. In addition, for the molten sulfate mixed system Li 2 SO 4 -Na 2 SO 4
system, molecular dynamics simulation was per formed using the shielded BMH typepotential, and the ion conduction mechanism was studied. Detailed analysis of the structure− 250 −NSG Found. Mat. Sci. Eng. Rep.

and ionic correlation angle revealed that there was a significant difference in the diffusion
mechanism between Li and Na. Fur thermore, a molecular dynamics simulation was
per formed on the superionic conductor AgBr-CuBr system, and the dif ference in the
distribution of Ag and Cu ions was discussed in detail on the unit cell plane and inside theunit cell. In addition, as an advanced study on electrolytes, we conducted a study bymolecular dynamics simulation of lactic acid-pyr uvic acid aqueous solution used forelectrolytes of biofuel cells, and found an unknown phenomenon that the dielectric constantincreases abnormally on the dilute concentration side and has a maximum.
1. はじめに
近年の地球温暖化などの気候変動を背景に、化石燃料の燃焼による CO 2 ガスの排出をより抑制しようとする世論の動向など、国民の環境改善とエネルギー消費節約の意識は益々高まっている。我が国は産業立国のためエネルギー資源をほとんど海外からの輸入に依存しているにも拘わらず、産業や家庭で消費されるエネルギーは世界有数の消費量である。一方で国土は狭く、産業や家庭から排出される汚染は周囲に大きな影響を与える。このため、環境に影響を与えない新しいエネルギー源の開発、効率的にエネルギーを伝達・貯蔵する技術、さらに環境を改善する取り組みとそのための技術が求められている。
超イオン導電体または固体電解質と呼ばれ、固体でありながら可動陽イオンが液体と同程度の拡散係数を持つ物質がある [ 1 ]。我々は、CO 2 を排出しない次世代水素燃料電池の電解質として用いられ、水素センサーとしても期待される超プロトン伝導体の混合系（MD）と第一原理
CsHSO 4 -CsH 2 PO 4 について、融体とガラス状態において分子動力学法
計算を行い、ガラス転移に伴う電子状態の変化を認め、水素がガラス中でも結晶中と同様に SO 4 や PO 4 と相互作用しながら移動していくことを確認した [ 6 ]。
超イオン導電体中では主に一種類の可動陽イオンが伝導に寄与するが、我々はより複雑な超イオン導電体の混合系に着目し、2 種類の陽イオンによる伝導について考察した。混合系とすることで超イオン導電相への転移温度が下がり、実用化により近づくと考えられるからである。2 種類の陽イオンによる伝導が起こる超イオン導電体については、2 種類の陽イオンはほぼ同様に分布し拡散するというのが定説であった。我々は 2 種類のイオンの分布の差異、輸送係数の違い等について考察してきた。Li 2 SO 4 -Na 2 SO 4 系は高温で二種類の陽イオンが拡散する超イオン導電体で、高温で動作する二次電池、燃料電池、蓄熱器などへの応用が期待される。Li と Na では SO 4 イオンとの相互作用や拡散のメカニズムが異なることを議論した [ 11 ]。
また、貴金属ハロゲン化物には超イオン導電性を示すものがありこれまで多くの研究がなされてきた。我々は、より複雑な貴金属ハロゲン化物の混合系に着目し、超イオン導電
相や融体の構造、輸送現象やダイナミックスについて考察してきた。二種類の可動陽イオ
ン Ag と Cu を含む AgBr-CuBr 系の超イオン導電相の MD による研究を行い、単位格子内の Ag と Cu の分布の差異について詳細に検討した [ 14 , 31 ]。
さらに電解質の発展的研究として有機物の電解質溶液を取り上げた。近年、人間の体内や体表面でも使用できる新しい安全なエネルギー源として、溶液中の有機分子を燃料とする燃料電池が注目されている。有機分子として人間の体内にも存在する乳酸やグリコール酸等の生体関連物質を使うため「バイオ電池」とも呼ばれている。バイオ燃料電池の電解質− 251 −

（公財）日本板硝子材料工学助成会，40（2022）

としての有機物の水溶液の分子動力学シミュレーションを継続して行い、構造や輸送現象、[ 2 - 5 , 7 - 10 , 12 - 14 ]さらに誘電率について一連の考察を行った。
2. 無機固体酸の融体およびガラス中の水素電導のシミュレーションによる研究
無機固体酸 CsH 2 PO 4 や CsHSO 4 等は高い水素伝導度を持つ個体電解質であり、超プロ
トン伝導体とも呼ばれる。超プロトン伝導相では、水素の伝導度は液体電解質と同様の約10 - 2 S/cm 程度である。従って、無機固体酸は 400 K～500 K の温度範囲で動作する燃料電（H 2 PO 4）池の有力な候補である。無機固体酸 CsH 2 PO 4 と CsHSO 4 の混合系 Cs（HSO24）
でも転移温度 370 K 以上で高い水素伝導が認められる。この系の特徴的な性質として、高温の融体からの急冷によって高温の構造が室温でも保持され、室温において超プロトン伝導ガラスが形成される [ 16 , 17 ]。無機固体酸 CsH 2 PO 4 と CsHSO 4 についてはこれまで MDによる研究が行われてきたが、混合系 Cs（HSO（H 2 PO 4）については我々の知る限りこ2
4）
れまで MD による研究はおこなわれていない。これらの状況を踏まえて Cs（HSO2
4）（H 2 PO 4）系の MD による研究を計画した。MD には遮蔽効果とイオンの分極を考慮した
Born-Mayer-Huggins（BMH）type のポテンシャルを用いた [ 6 ]。これらは次のように表さ
れる。
（1）
（2）
（3）
方程式
（2）において f は隣接するイオンによる Coulomb 相互作用の減衰係数であり 0 . 4
方程式（1）と（2）の fij（r）と仮定する。Pij は電荷とイオンの分極 αij による双極子の相互作用、はポテンシャルの発散を防ぐための減衰項であり、以下のように表される。
（4）
ここで、r = 3 . 0 Å において f +-（r）が訳 0 . 2 となるようにパラメータ k ij= 0 . 77 とする。以上に加えて、O-S-O, S-O-S, O-P-O, P-O-P の三体相互作用（5）を導入する。
（5）
ここで θijk は O-S-O, S-O-S, O-P-O, P-O-P の結合角である [ 6 ]。MD のセルには 1875 個の（H 2 PO 4）を 含 ん で い る。Coulomb 力 の 計 算 に は原 子、 す な わ ち 125 個 の Cs（HSO
2
4）
Ewald 法を用いる。水素のイオン質量が小さいため 0 . 001 fs おきに Gear アルゴリズムを用いて数値積分を行う。温度制御には Nose の方法を用いる。系は最初 2000 K で粒子数、温度、体積（NTV）一定の条件で平衡に保たれ、次に 1K/s の割合で 1500K まで冷却される。
更に 0 . 1 K/s の割合でガラス転移温度以下の 1150 - 550 K の任意の温度まで冷却される。
急冷によって得られる以下の結果から系のガラス転移が再現されたことが推定される。
（1）
S-O、
P-O 等に対する 2 体分布関数の第二近接のピークに特徴的なショルダーが現れた。
（2）− 252 −NSG Found. Mat. Sci. Eng. Rep.

単位温度以下ではイオンの拡散係数が急激に減少した。
（3）Voronoi 多面体解析から、多面体の最多の頂点の数は融体およびガラスともに 5 であるが、ガラス転移によって頂点の数の分布が若干大きい側にシフトする。また、面の数の分布も大きい側にシフトし、最多の面の数が融体の 12 からガラスの 13 に増加する。これらを 図 1（a）、（b）に示す。

図 1（a）

図 1（b）

次に、融体からガラス転移することに伴う電子状態の変化について考察する。これまで
一般に物質のガラス転移についての考察は主に熱力学に基づくものが主であり、電子状態の変化についての考察は我々の知る限り、我々の銀ハロゲン化物のガラス転移の研究以前（H 2 PO 4）系の融体からにこれまで見られない [ 18 ]。このような状況のもとで、Cs（HSO
24）
ガラス状態への転移における電子状態の変化について考察した。第一原理計算はこれ目での我々の研究 [ 18 ] と同様に密度汎関数法の局所密度近似（local density approximation, LDA）
[ 19 , 20 ]によって行われ、Troullier and Martins による norm 非保存方の pseudopotential [ 21 ]を用いた。周期的境界条件を用いた直方体中の平面波展開の Cut off energy は 400 eV である。Brillouin zone 中の積分は 2 x 2 x 2 Monkhorst and Pack grid を用いた [ 22 ]。これらの計算は ABINIT [ 23 ] を用いて行った。ガラス構造の 300 K における等電子密度面を 0 . 01 a.
u. について図 2 に示す。図 2 において、H-O、S-O、P-O のような陽イオンと陰イオンのペアが同一の面内に存在することが見て取れる。更に、電子の状態密度を、1150 K における融体と 300 K でのガラス状態について図 3（a）、（b）に示す。これらの図で S 3 s、P 3 s、O 2 s 状態はほとんど重なっているが、これらの電子状態が共有結合的であることを示している。融体の電子状態の図 3（a）では H 1s 状態はこれらの状態と一部が重なっているが、ガラス状態の図 3（b）では H 1 s 状態は主に O 2 p 状態と重なる。 この結果はガラス状態においても水素イオンが SO 4 および PO 4 イオンと相互作用しながら拡散する“パドルホイ
ールメカニズム”が成り立つことを示している。

− 253 −

（公財）日本板硝子材料工学助成会，40（2022）

図2

図 3（a）

図 3（b）

3. 溶融硫酸塩混合系 Li2SO4-Na2SO4系の分子動力学法によるイオン伝導機構の研究
Li 2 SO 4 -Na 2 SO 4 系は高温で二種類の陽イオンが拡散する超イオン導電体で、高温で動作する二次電池、燃料電池、蓄熱器などへの応用が期待される。超イオン導電体における陽イオンの拡散機構のモデルとしては、前章で述べたように結晶格子点で回転する正四面体の SO 4 などの陰イオンと Li や Na などの陽イオンが相互作用しながら拡散していく“パドルホイールモデル”と、格子点の SO 4 陰イオンの相関の無い振動と回転によって拡散の経路が作られるとするパーコレーションモデル”がある。超イオン導電相におけるシミュ
レーションによる先行研究としては MD [ 25 ] 及び逆モンテカルロ法（RMC）[ 24 ] による研究が行われているが、得られた構造の結果には不一致が見られる。また、融体での MD による研究はこれまで行われていない。これらを踏まえて、本研究では主に Li 2 SO 4 -Na 2 SO 4系の融体中での Li と Na の拡散機構の違いについて考察する。本研究では先行研究 [ 25 ] とは異なり、前章と同様の遮蔽された BMH 型のポテンシャルと O-S-O、S-O-S に対しては三体相関のポテンシャルを用いる [ 11 ]。MD による結果のうち、Li と Na の分布の違いはS 及び O イオンとの結合角の分布に端的に表れている。得られた結合角の分布を図 4（a）（d）に示す。本研究の MD は融体で行われたが、図 4（a）、（b）の O-S-O、S-O-S の結合角分布は SO 4 の正四面体構造がある程度残っていることを示している。MD は融体で行われているが、これらの結果は陽イオンの拡散のメカニズムがある程度固相との類推で理解できることを示唆している。図 4（c）、（d）に O-Li-O と O-Na-O, 及び S-O-Li と S-O-Na をそ− 254 −NSG Found. Mat. Sci. Eng. Rep.

れぞれ示す。図 4（c）より、Li と Na イオンは四面体位置（tetrahedral site）、及び八面体位
置（octahedral site）に分布していることが示唆されるが、Li と Na の分布には違いが見られる。Li と Na の不一致は図 4（d）において顕著に認められる。図 4（d）の S-O-Li では約60 ° 110 °にピークがあり、S-O-Na では約 100 °にピークがある。これらの結果は融体中に固相と同様の構造がある程度残っており、Na イオンは主に四面体位置に分布しパドルホ
イールメカニズムによって拡散することを示唆している。一方 Li イオンは主に八面体位置に分布し、拡散路に沿って移動する、つまりパーコレーションメカニズムによって拡散すると考えられる。

図 4（a）

図 4（b）

図 4（c）

図 4（d）

4. 超イオン導電体AgBr-CuBrにおけるAgとCuイオンの分布の差異に関する分子動力学シミュレーションによる研究
貴金属ハロゲン化物のうち AgI や CuBr は高温で超イオン導電体になるが、岩塩構造のAgBr や AgCl は融解するまで超イオン電導が現れない。そのため、これら 2 種類の異なる貴金属ハロゲン化物の混合系に特に興味が持たれる。二種類の貴金属可動陽イオンを持つ超イオン導電体の構造は主として AgI-CuI 系について考察されてきたが、Ag と Cu イ
− 255 −

（公財）日本板硝子材料工学助成会，40（2022）

オンは同じように分布し拡散するというのが定説であった [ 26 , 27 , 28 ]。 Ag − I 間と Cu − I間の二体分布関数の差異は実験から認められていたが [ 28 ]、I イオンの bcc 構造中の Ag とCu イオンの分布の差異の明確なイメージは示されていなかった。AgBr-CuBr 系について
は特に AgBr rich 側で詳細な実験的解析がなされているが [ 29 , 30 ]、特に CuBr rich 側の高温の超イオン伝導相での議論が不十分であった。これらの先行研究を踏まえて、二種類の可動陽イオン Ag と Cu を含む AgBr-CuBr 系の超イオン導電相の MD による研究を行った。
MD はイオンの分極を考慮した RVP type の二体ポテンシャルを用い、（Ag xCu 1 -x）Br,x Ag= 0 . 10 及び x Ag= 0 . 40 の組成で NTP 一定の条件で計算を行った [ 14 , 31 ]。特筆すべき結果
として、Br の bcc 格子の（1 , 0 , 0）面内における Ag と Cu 陽イオンの分布の差異が認められた。xAg= 0 . 40 の組成では Ag と Cu はどちらも主に八面体位置 6（b）site に分布し、特にCu には < 1 , 0 , 0 > 方向に帯状の分布が見られる。一方 x Ag= 0 . 10 の組成では Cu は 12（d）
site に分布するのに対して、Ag は八面体位置 6（b）site に分布しており、いずれにも< 1 , 0 , 0 > 方向に帯状の分布が見られる。x Ag= 0 . 10 の組成での格子内のイオンの分布を詳細に考察するため、（1 , 0 , 0）面を格子の一片の 1 / 16 、1 / 8 、1 / 4 、3 / 8 、1 / 2 だけセルの中
心方向にシフトさせた面でのイオン分布を求めた。Cu を図 5（a）、Ag を図 5
（b）
に示す。

図 5（a）Cu の分布

図 5（b）Ag の分布− 256 −NSG Found. Mat. Sci. Eng. Rep.

Ag と Cu の分布の明確な差異に関する報告は我々が知る限り本研究が最初の報告である [ 14 , 31 ]。このような陽イオンの分布の差異は例えば輸送係数など他の物性に影響を与えるものと考えられる。我々は研究を継続しこれらの点を明らかにしてきた [ 32 - 36 ]。Ag とCu の分布の差異は直接的には分子動力学シミュレーションに用いた二体ポテンシャルによる二体相互作用の差異に起因する。他方、我々がこれまで行った電子状態の研究から、混成軌道 Ag 4 d, Cu 3 d, Br 4 p が形成される [ 36 ]。Cu 3 d は Ag 4 d よりも Fermi レベルに近い。換言すれば、Ag-Br は Cu-Br よりも共有結合性が強い。この結合の違いが格子中のAg と Cu の分布に影響していると考えられる。
5. 生体燃料電池の電解質に用いられる 乳酸－ピルビン酸水溶液における誘電率の異常に関する分子動力学シミュレーションによる研究
本研究の発展として、溶液中での水分子と有機分子の集合体が溶液の物性に与える影響（L-acid）を燃料とした燃料について考察する [ 2 - 5 , 7 - 10 , 12 , 13 ]。近年、人の汗にも含まれる乳酸電池が試作され、人の皮膚でも動作する携帯可能な燃料電池として注目されている。これまでの予備的研究で、生理食塩水中の乳酸とピルビン酸（P-acid）の希薄水溶液で誘電率の異常な増大を認めた [ 3 ]。誘電率は燃料電池を設計する際に重要になる物理量である。本
研究の主な目的は分子動力学法によって乳酸とピルビン酸の希薄溶液中の構造と誘電率の変化を詳細に調べるためである。誘電率εは次のように定義される [ 15 ]。
（6）
ここで k B はボルツマン定数、ε 0 は真空の誘電率、ε’はセルの外側の誘電率、T は絶対温度、V はセルの体積である。M 2 の項は次のように表される。
（7）
ここで は i 番目の分子の双極子モーメントである。周期的境界条件のもとでは ε’→∞を仮定すると ε はその結果次のように表される。
（8）
乳酸とピルビン酸の様々な組成において求めた ε の値を図 6 に示す。特徴的なピークが
約 0 . 20 mol/L の組成に認められる。このような希薄溶液側の極めて大きな ε の値は我々の研究以前に報告の無い現象である [ 4 ]。その原因としては、溶液中のイオンの集合によって大きな双極子モーメントが形成されたためと考えられる。また、その後組成の増加とともに誘電率が増大するが、これは方程式（8）に見られるように双極子モーメントを持つ[ 37 ]
イオンの組成の平方に比例する 。本研究は継続中であり、有機物水溶液の希薄溶液における大きな双極子モーメントの形成と誘電率の増大については今後さらに詳細に研究を進める予定である。

− 257 −

（公財）日本板硝子材料工学助成会，40（2022）

図6

5. 結論
無機固体酸の水素電導ガラスとその融体において、遮蔽された BMH type のポテンシャ
ルを用いて分子動力学シミュレーションを行い、急冷によってガラス転移が起きることを確認した。ガラス中で水素は SO 4 及び PO 4 イオンと相互作用しながら拡散する。また、ガラス転移に伴って電子状態の変化が起こることを第一原理計算によって確認した。 次に、燃料電池の電解質として用いられる Li 2 SO 4 -Na 2 SO 4 系の融体について Li イオンとNa イオンの分布の違いが、格子を形成する SO 4 イオンとの角度相関に顕著に現れることを確認した。これは Li イオンと Na イオンの拡散の機構が異なることを示唆している。
更に、超イオン導電体（Ag xCu 1 -x）Br における Cu rich 側の x Ag= 0 . 10 の組成で、単位格子
面（1 , 0 , 0）
上の Ag と Cu イオンの分布の顕著な差異について、単位格子の内部の分布を詳細に調べた。燃料電池の電解質としての有機物の水溶液については、希薄溶液側で顕著な誘電率の極大を認めた。これは溶液中で分子が集合し大きな双極子モーメントが生じたためと考えられる。
6. 謝辞
本研究は、2019 年度日本板硝子材料工学助成会の研究助成を受けて行ったものである。
同助成会に心より感謝いたします。
7. 参考文献
1）See for example, S. Chandra, Super Ionic Solids,（Amsterdam, North-Holland）, 19812）S. Matsunaga, IOP Conference Series: Materials Science and Engineering 369（1）012014（6 pp）, 20183）S. Matsunaga, AIP Conference Proceedings 1981 , 020115（4 pp）, 20184）S. Matsunaga, Polymer Engineering & Science 59（12）, 2474 - 2478 , 20195）田巻繁 , 松永茂樹 , 日下部征信 , 水・水溶液系の物性 , 近代科学社 Digital, 20206）S. Matsunaga, Molecular Simulation 46（17）, 1327 - 1333 , 20207）S. Tamaki, S.Matsunaga, M.Kusakabe, Electromagnetic Field Radiation in Matter, Dr.
Walter Gustavo Fano（Ed.）, IntechOpen, Chap. 5 , 69 - 101 , 20208）S. Matsunaga, Prime Archives in Molecular Sciences, Dr. Sławomir Lach（Ed）, Vide Leaf,
1 - 33 , 2020 .
− 258 −NSG Found. Mat. Sci. Eng. Rep.

9）S. Matsunaga, Journal of Physics: Conference Series 1730 , 012045（8 pp）, 2021
10）S.Matsunaga, AIP Conference Proceedings, 2380（1）, 200001（6 pp）, 2021
11）S.Matsunaga, Physics and Chemistry of Liquids 60（2）, 190 - 202 , 2021
12）S.Matsunaga, 29th European Biomass Conference and Exhibition Proceedings 5CV.2.20
1392 - 1395 , 2021
13）S.Matsunaga, New Ideas Concerning Science and Technology, Vol. 9 , B P International,
125 - 132 , 2021
14）S.Matsunaga, Newest Updates in Physical Science Research Vol. 11 , B P International,
1 - 9 , 2021
15）J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids, 2 nd Ed., London: Academic
Press, 1986 .
16）S. M. Haile, D.A. Boysen, C.R.I. Chisholm, R.B. Merle, Nature 410 , 910 , 2001 .
17）C. R. I. Chisholm, S.M. Haile, Acta Crystallogr. B 55 , 937 , 1999 .
18）S. Matsunaga, Ionics, 24（5）, 1371 - 1376 , 201819）P. Hohenberg and W. Kohn, Physical Review, 136 , B 864 , 196420）W. Khon and L. J. Sham, Physical Review, 140 , A 1133 , 196521）N. Troullier amd J. L. Martins, Physical Review B, 43 , 1993 , 1991
22）H.J. Monkhorst, J.D. Pack, Phys Rev B 13 , 5188 , 1976
23）ABINIT code is a common project of the University Catholique de Louvain, Corning
Incorporation, and other contributions（URL http://www.abinit.org ）; X. Gonze et al.,Zeitschrift für Kristallographie, 220 , 558 , 200524）L. Karlsson, R. L. McGreevy, Solid State Ionics. 76 , 301 - 308 , 199525）M. Ferrario, M. L. Klein, I. R. McDonald, Mol. Phys. 86 , 923 - 938 , 1995
26）A. K. Ivanov-Shitz, B. Yu. Mazniker, E. S. Povolotskaya. Crystallogr. Rep. 46 , 292 - 295 ,200127）J. Bośko, J. Rybicki, Solid State Ionics. 157 , 227 , 2003 .
28）J. Kimura, T. Ida, M. Mizuno, K. Endo, M. Suhara, and K. Kihara, J. Mol. Struct. 522 , 61 ,200029）M. Saito, S. Tamaki, Solid State Ionics. 60 , 237 , 199330）J. Teltow, Z. Phys. Chem. 18 , 134 , 195831）S. Matsunaga, J. Phys. :Conf. Ser. 144 , 012011 , 200932）S. Matsunaga, S. Tamaki, J. Phys.: Condens. Matter. 20 , 114116 , 200833）S. Matsunaga, Molecular Simulation, 38（5）, 384 , 201234）S. Matsunaga, Proc. of 2 nd Int. Conf. Smar t Materials and Nanotechnology inEngineering（SMNE 2012）, 20 , 309 , 201235）S. Matsunaga, Trans. Mater. Res. Soc. Jpn. 38（2）, 177 , 201336）S. Matsunaga, Ionics, 21 , 161 , 201537）P. Wang and A. Anderko, Fluid Phase Equilibria, 186 , 103 - 122 , 2001

− 259 −