RISTニュース

No．69（2023）

全原子モデルを用いたイオン液体のシェアシニングの解明

Elucidation of Shear-thinning in Ionic Liquids Using
All-atom Models
高度情報科学技術研究機構
山田 達矢、手島 正吾
液体に与えるせん断速度を増加させた場合に、片栗粉の水溶液では粘度が増し、マヨネーズでは粘度が減り、そして、水では粘度が変わらない、などの異なる振る舞い起こすことが知られていて、それぞれ、シェアシックニング流体、シェアシニング流体、そして、ニュートン流体と呼ばれている。この様な、変形や流動に関する学問はレオロジーと呼ばれて注目されており、振る舞いの違いがどの様な物理的メカニズムに起因するのか、研究が行われている。本研究では、凝集し易いカーボンナノチューブ（CNT）を分散させることができるイオン液体について、その流動性を非平衡全原子分子動力学シミュレーションによって計算し、シェアシニングの振る舞いを起こす原因を調査した。その結果、Eyringモデルに従って、イオン分子がせん断流内で 1 つの準安定状態から別の準安定状態に移動するときに遭遇する見かけのエネルギー障壁の低下が生じることが判明し、これがシェアシニングの原因であることを突き止めた。流体の粘性の振る舞
いを分子レベルで説明できることの一例を以下で紹介する。
1 ．イントロダクション

では、陰陽両イオンの配置は、静電相互作用

イオン液体は低い蒸気圧、高い熱安定性、による秩序構造を持つ。
簡潔に言えば、
陰（陽）イオン伝導性を示す［1-3］。また、イオン

イオンの近傍には陽（陰）イオンが分布しや液体に 1wt%程度のカーボンナノチューブ

すく、その結果、各イオンは逆の電荷をもつ（CNT）を混和すると、CNTがよく分散したイオンに囲まれている［11］。Yamaguchiは、

バッキーゲル［4、5］が得られ、従来の潤滑

シェアシニングが始まるせん断速度付近にお剤では使用できない過酷な条件下での役割がいて、イオン配置の秩序構造もまた変化し始期待されている［5、6］。イオン液体は、せめることを報告した［10］。同様の結果が、ん断速度の上昇に伴って粘度が低下するとい

逆非平衡分子動力学（RNEMD）シミュレー

うシェアシニング現象を示し［7、8］、このこション及び分極力場を用いたBernardinoらとは潤滑剤としての機能に大きく関与する。

の研究によっても報告されている［13］。まイオン液体のシェアシニング現象に関した、Van-Oanhらは、イオンの回転緩和時間
̇cτrotとシェアシニングの開始せん断速度γて、平衡分子動力学（EMD）および非平衡
分子動力学（NEMD）シミュレーションを̇ c＝1/10τrot の関係にあると報告してが、γ用いた先行研究［9-13］が行われている。そいる［12］。

こでは主に、せん断速度の増大に伴ってイオ

以上で紹介した通り、イオン液体のシェアンの配置がどのように変化するかについて注

シニングとイオン配置の変化の関係性は示さ

意が向けられてきた。せん断流れのない状況

れているが、イオン液体がシェアシニングを－19－

RISTニュースNo．69（2023）

示す仕組みの物理化学的側面は未だ不明のま示すように、Eyringモデルでは液体の流動ま残されている。そこで本研究では、イオンを、各分子がエネルギー的に準安定な状態S1分子間のエネルギーの観点からイオン液体のから、他の分子との相互作用から生じる高シェアシニングを定量的に説明することを試さEaのポテンシャルエネルギー障壁Tを越えみ、そのために、Eyringモデルに注目した。

て、別の準安定状態S2に移行する遷移の集ま

マクロな粒子系や、ポリマーが示すシェア

りとして表現する。そして、液体に外部からシニングについては、Eyringモデルに基づ

せん断応力が与えられている場合、この分子いた説明がなされている［14］。図 1（b）にの遷移には、
熱運動だけでなく、系のせん断応力に由来する機械的仕事Wもまた関与する。
一般的に、体積Vの流体がせん断応力σの下でひずみ量γだけせん断変形を受けるとき、この流体に与えられる機械的仕事はW＝σVγで表される。これを一つの分子に当てはめることを考える。図 1（b）に示すように、
液体中の分子層間の間隔をλ1、流速方向の分子の長さをλ2、せん断流の渦度方向の分子の長さをλ3、そして、隣接する 2 状態の間隔をλで表す。すると、 1 つの分子の占有体積はV＝λ1λ2λ3となり、分子が状態S1から障壁Tまで移動するために必要なせん断
ひずみ量はγ＝d/λ1（ここでd＝λ/2を用いた）となり、結果として、S1→Tまでの遷移の間に一つの分子に与えられる機械的仕事は
W＝σλ2λ3dで表される［14］。分子がS1の位置から見た障壁 T の見かけの大きさEa,app
は、Wによって押し下げられ、Ea,app＝Ea－
Wとなる。この結果、より大きなせん断速度
̇ では、より大きな内部応力σによって分子γに与えられる機械的仕事Wが増大し、Ea,app
がより減少し、分子の遷移が大幅に促進され、最終的に液体の粘度が減少する。以上のこと
をEyringモデルにおいて詳しく考察すると、せん断速度の増大に伴うせん断粘度の低下が

図 1 （a）単純なせん断流れ下におけるせん断
速度γ、及び、せん断応力σzxの定義。
（b）Eyringモデルにおけるせん断流れの模式図。せん断流れは、エネルギー的に準安定な状態（黒実線の円、状態S1に対応）から隣接する準安定状態（灰色破線の円、状態S2に対応）への分子の移動という、微視的な分子位置の遷移の集まり
によって表現される。

次式のように導かれる［14］。

・（式 1 ）
ここで、kB、T、及び、η0は、それぞれ、ボルツマン定数、絶対温度、及び、低せん断速度の極限におけるせん断粘度である。
－20－

RISTニュースNo．69（2023）

イオン液体のシェアシニングの仕組みがる。これらに相当する量は、せん断流れのな

Eyringモデルに基づくものであるならば、
̇ －ηの関係）は、
イオン液体の粘度曲線（γい平衡状態におけるイオン分子の平均拡散係式 1 でフィッティングできるはずである。さからも計算することができる。そこで本研究らに、そのフィッティングから、エネルギー
障壁Tの大きさ、及び、イオン分子の基礎的では、
まず、
イオン液体の非平衡せん断シミュ̇ －ηの関係式への当て
レーションを行い、γな拡散係数に関する量を求めることができはめを行い、シェアシニングに関与する障数DaveのArrheniusプロット（1/T vs lnDave）

壁の大きさE0 及び、イオン分子の基礎的な拡散係数D0 を求めた。次に、せん断流れのない平衡状態のシミュレーションを行い、イオン分子の平均拡散係数DaveのArrheniusプロットから、流動活性化エネルギーEa、及び、
イオン分子の基礎的な拡散係数Dave,0を計算した。最終的に、E0及びD0とEa及びDave,0をそれぞれ比較して、イオン液体のシェアシニ
ングの仕組みがEyringモデルに基づいて説明できることを示す。
̇ ）の両対数プロッ
なお、式 1 は、η＝f（γトにおいて、高いせん断速度領域でほぼ同じ傾きでせん断粘度が連続的に無限に減少する
ようなシェアシニングを予測する。実験的に
̇ が大きくなるとともに粘度の低下は生は、γじなくなり、最終的にプラトーまたは極小となる。この現象を考慮するため、本研究では
式１を用いた曲線当てはめの代わりに、 1 次
図2

イオン液体に適用した 1 次元Eyringモ
デルを表す模式図。（a）あるイオン（赤丸）がせん断領域（x軸と灰色点線の間）を挟んで反対側にある他のイオンに対して相対的に移動するときには、イオン間の相互作用によって形成されるポテン
シャルエネルギー障壁（青実線）を乗り越えながら移動する必要がある。（b）
実際のところ、移動するイオンによって感じられるポテンシャルエネルギー障壁
U（x）（青実線）は、機械的仕事W＝σ
lxlyx＝σSod（灰色破線）によって変化させられる（青点線）。その結果、位置エネルギーの最小値（赤点線の円）から測定されるエネルギー障壁の見かけの高さはWの増大に伴って減少する。σが十分に大きい場合、エネルギー障壁は消滅する（青一点鎖線）。

－21－

元上のEyringモデルをBrownian dynamics
̇－
法を用いて数値的に解くことによって、γηの関係（シェアシニング曲線）を算出し、また、フィッティングに用いた（図 2 ）。
イオン液体の構造は一般に、イオンの配置の秩序性、及び、イオンに結合した無極性基の凝集の結果生じるナノスケールの不均一性により複雑である［15］。逆に言えば、小さな非極性基のみを含むイオンで構成されるイ
オン液体の構造はかなり均一であり［10］、そのようなイオン液体を用いればEyringモ
デルに基づく解析が簡素化される。したがって本研究では、小さな非極性基を含む陽イオンと小さな球状の陰イオンのみで構成される
emi-BF4（図 3 ）を、研究とシミュレーショ

RISTニュースNo．69（2023）

Sprengerらは、RESP電荷を係数c＝0.8でスンの対象として選択した。

ケールした電荷を用いることで、

2 ．分子シミュレーションの詳細

イオン液体の密度、熱容量、
モル蒸発熱、
及び、

2.1

ゼロせん断粘度の実験測定値がGAFF力場で分子モデル作成と力場の割り当て

本 研 究 で は、 陽 イ オ ン と し てemi

よく再現されると報告している［22］。本研（C6H11N2）、陰イオンとしてBF4より成るイ

究でも事前にcを1.0から0.6までδc＝0.05刻

オン液体を扱う（図 3 ）。全原子分子動力学

みで変化させ、emi及びBF4の平均拡散係数

シミュレーションの準備として、Avogadroを計算した。その結果、Sprengerらと同様に、［16］で分子構造を作成した。分子内力場

c＝0.7～0.8において実験的な拡散係数［23］で あ る 結 合 長、 結 合 角、 結 合 二 面 角、 及が良く再現されることがわかった。
それゆえ、び、improper角ポテンシャルの定数としてシェアシニングの仕組みに対する静電相互作

GAFF力 場［17、18］、 及 び、 文 献 値［19、用の役割を定量的に評価する場合を除き、以20］ を 用 い、 フ ァ ン デ ル ワ ー ル ス（vdW）

下で行うシミュレーションでは専ら、c＝0.8及び静電相互作用のカットオフ半径としての電荷スケーリングを使用した。

0.8nmを使用した。また、
クーロンポテンシャ

イオン液体の熱平衡化されたシミュレー

ルの計算にはEwald法を用いた。原子の電ション系の準備は以下の手順で行なった。ま荷は次の順序で決定した。先ず、emiに関しず、packmol［24］を用いて、周期境界条件てb3lyp/6-31g（d,p）、BF4に 関 し てb3lyp/6-を備えた大きなシミュレーションボックス内

31+g（d,p） の 計 算 レ ベ ル でGaussian 09に陽イオンと陰イオンを2048個ずつ配置し［21］を用いて静電ポテンシャル（ESP）電た。続いて、エネルギー最小化を実施した

荷を計算し、続いてRESP電荷、即ち、emi

後、絶対温度T＝600 K、圧力P＝1000atm、

及びBF4の各 i 番目の原子の部分電荷、それ時間刻みδt＝1fsで0.1nsのNPTシミュレー

を、AmberTools内のション（粒子数N、圧力P、温度Tの全てがantechamberプログラムを用いてESP電荷一定の条件）を行なって系を圧縮した。次より計算した［17、18］。

に、T＝600K、P＝1atmに 保 ち な が ら10ns

ぞれ、

及び

イオン液体では、原子の電荷としてRESPの平衡化NPTシミュレーションを行なった。

電荷を直接用いるよりも、それらを幾らか

さらに、P＝1atmの条件のもと、100nsのシ減少させた値を用いた方がより実験結果を

ミュレーションを行い、系の温度を600Kか

正確に反映するという研究報告が多数ある。

ら250Kまで徐々に冷却し、250～350Kの温度範囲で10K刻みごとに分子構造を抽出した。
これら各温度の構造のそれぞれについて、抽出した温度T及びP＝1atmの条件下で、δ
t＝0.5fsを用いて平衡化NPTシミュレーショ
ンを10ns行ない、ボックスの平均サイズを各温度でのシミュレーション結果の最後の5nsを用いて算出した。T＝300Kの場合、ボッ

図3

本研究で計算の対象としたイオン液体分子の構造。
（左）1-ethyl-3-methylimidazolium陽イオン。
（右）tetraﬂuoroborate陰イオン。

－22－

ク ス サ イ ズ は8.25×8.25×8.25nm3と な り、
密度は1.2g/cm3となった。この値は、298K
での実験値1.279g/cm3［23］よりも約 6 %小
さいが、本研究の目的のためには十分に実験

RISTニュースNo．69（2023）

結果を再現している。系のサイズを各温度で̇ －ηの関係）を示す。γ̇ ≤109s-1の範囲で（γの平衡状態の値に固定した後、さらに各温は、η＝26mPa・sで粘度が一定となる領域

度において10nsの最終平衡化のためのNVTが現れ、ニュートン粘性を示した。この粘シミュレーション（粒子数N、体積V、温度度は、実験によるゼロせん断粘度17.4mPa・s

Tの全てが一定の条件）を実施した。こうして得られたイオン液体系の構造を、せん断粘［23］よりも49%大きい程度であり、実験結̇ ≤3.2果を十分に再現できたと言える。109<γ

度及び拡散係数を計算するための本番シミュ×1012s-1の領域では、emi-BF4はシェアシニ

レーションの初期状態として採用した。
なお、以上の平衡化過程では、emi-BF4中のアルキングを示し、ηは最小約0.4mPa・sまで減少
̇ ＝1.0×1013s-1では粘度の増大が観した。γ

ル鎖長が短いため、イオン液体でよく観察さ察されたが、これはレオロジー的なものではれるような不均一なナノ構造［25］は現れなく、非常に高いせん断速度におけるシミュなかった。

レーション系の温度制御に由来した数値的なものであることが判明したため、以降の解析

2.2

拡散係数とせん断粘度

ではこのデータ点を除外した。

流動活性化エネルギーの算出に用いる̇ －ηの関シミュレーションで求められたγ

emi-BF4の平均拡散係数Daveは、250～350K

係に対して、式 1 を用いたフィッティングの各温度で実施した20nsのEMDシミュレー（実際には 1 次元のEyringモデルを用いた数

ションによって生成されたNVEアンサンブ

値計算によるフィッティング）を行い、そのル（粒子数N、体積V、エネルギーEの全て

結 果 と し てD0＝2.04×10－8m2/s、 及 び、E0

が一定条件）より計算した。

＝18.7kJ/molと い う 値 が 得 ら れ た（ 図 4 ）。

せん断粘度は、SLLODアルゴリズム［26］フィッティングの精度は十分によく、それゆ

及びLees-Edwards境界条件［27］を適用したNEMD法による定常せん断流れのシミュ
レーションによって計算した。計算の実施と準 備 に はLAMMPS［28］ 及 びmoltemplate［29］を用いた。シミュレーションは、T＝300Kの条件において、時間刻みδt＝0.5fs
̇ ≤1013s-1）を用い、様々なせん断速度（108≤γで実施した。シミュレーションの長さは1.5
̇ ＝1013s-1）から10ns（γ̇ ＝108s-1）のps（γ範囲で行った。初期状態由来の影響を回避するため、シミュレーションの後半50-75%の部分で応力テンソルのzx成分の平均値σzxを
̇ の関係式を用いて最終計算し、η＝σzx /γ

図4

的なせん断粘度ηを求めた。
3 ．結果と考察
3.1

せん断粘度と拡散係数の比較

図 4 に、定常せん断流れ下のemi-BF4のシイオン液体emi-BF4のせん断粘度のせん
断速度依存性。●は、非平衡分子動力学シミュレーションの計算結果を表す。
□は、●に対して 1 次元Eyringモデルを使用したフィッティング結果を表す。
フ ィ ッ テ ィ ン グ に よ り、E0＝18.7kJ/
mol、 及 び、D0＝2.04×10-8m2/sが 得 られた。γ は、 1 次 元Eyringモ デ ル内で0.33倍にスケールされた。
・

ミュレーション結果より計算した粘度曲線－23－

RISTニュースNo．69（2023）

え、イオン液体のシェアシニングがEyring

ティングを行なった。その結果、Dave,0＝6.39

モデルに基づいて生じていることが支持され×10-6m2/s、及び、Ea＝29.4kJ/molという値た。しかしながら、シェアシニングに関与すが得られた（図 5 の●）。

るエネルギー障壁と、流体の粘度を記述する

結 果 と し て、E0＝18.7kJ/molとEa＝29.4

Eyringモデル中の流動活性化エネルギーがkJ/molの間には、無視できない不一致が存同じものであるかどうかを判断するために在することがわかった。最終的に著者らは、は、E0及びD0の値を、平衡状態のArrhenius

この不一致の原因が自由体積理論［30、31］

プロットから求められるEa及びDave,0の値とで説明できることを明らかにした。本研究で比較する必要がある。

は、せん断粘度計算のためのシミュレーショ

そこで、280～320Kの範囲の各温度で実ンは全て同一の体積で実行されたが、拡散係施されたEMDシミュレーションの結果から数計算のためのシミュレーションにおいてはイ オ ン の 拡 散 係 数 Daveを 求 め、Arrhenius

ボックスの体積が予め各温度での平衡状態のプ ロ ッ ト を 作 成 し、Dave＝Dave,0 exp（ －Ea /体積に設定された。自由体積理論によると、kBT）の式（Arrheniusの式）を用いてフィッ

せん断粘度は、流体の総体積に対する分子が移動するために利用できる自由空間の体積の割合に依存する。つまり、液体の熱膨張がシ
ミュレーションに考慮されている場合では、温度に依存せず同一の体積で行うシミュレー
ションに比べて、拡散係数の温度依存性がより大きく観察され、その結果として流動活性化エネルギーの大きさEaがより大きく評価されてしまうと考えられる。
このことを確認するため、イオン液体系の体積を温度に依存せず一定の体積（300K
かつ 1 atmで平衡化されたシミュレーショ

図5

絶対温度T＝280～320Kの範囲における拡散係数Daveの1/T 依存性。流動活性化
エネルギーEaと基本的な拡散係数Dave,0を、Dave＝Dave,0 exp（ －Ea/kBT） を 用
いて計算した。一定圧力（1atm）で予め平衡化された体積を有する系における平衡MDの結果（●）について、Dave,0＝6.39×10-6 m2/sおよびEa＝29.4kJ/molを得 た。 ま た、300Kか つ 1 atmで 予 備 平衡化されたものと同じ体積に固定した系を用いた平衡MD計算の結果では、Dave,0＝Dave,0,cstr.＝1.75×10-8m2/s、 及 び、Ea
＝Ea,cstr.＝14.6kJ/molが得られた（○）。
実験的な拡散係数（◇）は、参考文献［23］に与えられたパラメーターを使用して求
め ら れ、Dave,0＝1.46×10-6m2/s、 及 び、
Ea＝25.7kJ/molという値が得られた。
－24－

ンボックスのサイズである8.25×8.25×8.25nm3）に固定した上でEMDシミュレーショ
ンを各温度で再度実施し、こうして得られた拡散係数のArrheniusプロットを用いて、熱膨張の影響を含まない流動活性化エネル
ギーEa,cstr.を計算した（図 5 の○）。その結
果、Ea,cstr.＝14.6kJ/molと い う 値 が 得 ら れ、
E0＝18.7kJ/molにかなり近づいた。さらに、
Dave,0,cstr.＝1.75×10-8m2/sという値が得られ、この値は、D0＝2.04×10-8m2/sとほとんど同じである。
以上のことから、E0及びD0が、熱膨張の影響を取り除いたEa及びDave,0（即ち、Ea,cstr.
及びDave,0,cstr.）と十分一致していることがわかった。それ故、自由体積理論を考慮した上

RISTニュースNo．69（2023）

で、シェアシニングに関与するエネルギー障
壁と、流体の粘度を記述するEyringモデル（式 2 ）

中の流動活性化エネルギーが同じものである
ことが支持された。最終的に、温度による体こ こ で、EBは 原 子 間 の 共 有 結 合、 結 合

積膨張の影響がEa－Ea,cstr.＝14.8kJ/molの違角、 結 合 二 面 角、 及 び、improper 角 のいを生み出していると結論された。

ポ テ ン シ ャ ル エ ネ ル ギ ー の 総 和 を 表 す。
は、 静 電 相 互 作

3.2

エネルギー障壁の分解用によるエネルギーの総和を表し、εr及び引き続いて、シェアシニングを示す粘度曲ε0は、 そ れ ぞ れ、 比 誘 電 率 及 び 真 空 の 誘線から求めたエネルギー障壁の大きさE0と、電率を表す。EvdWは全ファンデルワールス

熱膨張の影響のない拡散係数から求めた流動相互作用の総和を表す。電荷スケーリング

活性化エネルギーの大きさEa,cstr.に対して、を考慮すると、静電相互作用のそれらのエネルギーがどのような相互作用に寄与はと書くことができる。

よって形成されているのかを調査するため、ここでは、原子の電荷にRESP電荷を

E0及びEa,cstr.を各相互作用成分に分割するこ

直接割り振った場合における系の静電相互作とを試みた。ここでは、E0とEa,cstr.の関係を用によるエネルギーの総和である。それゆ明確にするため、電荷スケーリング定数 c がえ、系の全体構造が電荷に鋭敏でない場合、異なるシミュレーションにおいても、同一体

EES、E0、 及 び、Ea,cstr.は、c2の 一 次 関 数 で積（8.25×8.25×8.25nm3）の条件のもとで表されると予想できる。

計算を実行した。

図 6 に、c＝0.0～0.8の範囲で計算した拡系の全ポテンシャルエネルギーの総和 Etot

散 係 数DaveのArrheniusプ ロ ッ ト を 示 し た。

は、 3 つの寄与の合計として次式で表すこと

各c値 で 計 算 し た270～330Kの 温 度 範 囲 のができる。

データ点にフィッティングを行なって、各 c値におけるEa,cstr.の値を計算した。
また、図 7 に、各 c 値において実施した定常せん断シミュレーションによる粘度曲線の結果を示した。これらの粘度曲線に 1 次元のEyringモデルによる数値計算フィッティン
グを施し、各 c 値におけるE0の値を計算した。
最 終 的 に、 図 8 に お い て、Ea,cstr.及 びE0
をc2の 関 数 と し て プ ロ ッ ト し た。Ea,cstr.及
びE0の値は、ともにc2に対して線形に増加することがわかる。c＝0では、EES＝0となるため、エネルギー障壁は内部ポテンシャ
ル 相 互 作 用 とvdW相 互 作 用 の み が 寄 与 す

図6

各電荷スケーリング係数cを用いた場合
における、emi-BF4の平均拡散係数Daveの1/T依 存 性。c＝0.75～0.8で の 計 算 上の拡散係数は、実験値［23］とほぼ定量的に一致した。

る。それらがc2と独立な変数であると想定すると、＝5.1kJ/mol、及び、＝8.7kJ/molと い う 値 が 得 ら れ る（図 8 ）。両者の差3.6 kJ/molは、c2に依存せず一定であることから、この差は静電相互作－25－

RISTニュースNo．69（2023）

用に起因するものではないと考えられる。

ている。c＝0.8の場合においては、
9.5kJ/molおよび

c＝0におけるE0とEa,cstr.の大きさの違いの～～10.3kJ/molという値

由来をより明快にするため、図 9 において、が得られた。これは、Ea,cstr.及びE0の、それ平衡状態のシミュレーション及び定常せん断

ぞれ、65%および55%に相当する。静電相互

流れのシミュレーションの間で、emiイオン作用の寄与の大きさがEa,cstr.とE0の間でほとに含まれるイミダゾリウムリングのスナップ

んど一致したことは、両者の間に密接な関係

ショットを比較した。イミダゾリウムリングがあることを反映している。即ち、シェアシの配向は、平衡状態ではランダムであるが、ニングに関与するエネルギー障壁と拡散係数

高せん断速度では流れ方向に配向する傾向がの流動活性化エネルギーの成り立ちがほとん見られる。配向の変化はせん断流れが存在す

ど同一であることを示している。そしてそのる場合のみに生じるため、先に述べた3.6kJ/

ことは、イオン液体emi-BF4のシェアシニン

molの差に関係していると考えられる。高せグの仕組みが、Eyringモデルに基づいて説ん断速度による配向変化は主に、emiイオン

明できることの証拠でもある。

が球状ではない平べったい形状を有していることに由来していると考えられる。また、

4 ．まとめemiイオンの配向の変化により、エネルギー

本研究では、イオン液体emi-BF4が示す障壁もまた変化するという考えは、先行研究

シェアシニングの仕組みを調査するために、の結果［10、32］においても主張されている。

平衡および非平衡分子動力学シミュレーショ

Ea,cstr.及びE0に対する静電相互作用の寄与

ンを実行した。先行研究においては、せん断（及び速度の上昇に伴うイオン分子の配置の変化と）は、それぞれ、c2に比例しいう、系の構造変化に焦点を当てた研究がな

図7

各 電 荷 ス ケ ー リ ン グ 係 数cを 用 い た 場
合 に お け る、emi-BF4の せ ん 断 粘 度 のせん断速度依存性。黒実線は、 1 次元
Eyringモデルを用いたフィッティングの結果を表す。シミュレーション結果（標準誤差棒付きのマーカー）に対して、どのc値 に お い て もD0＝2.04×10-8m2/sを使用した上で、E0を変化させてフィッ
ティングを行った。

図8

－26－

NEMDシミュレーションで計算された
せん断粘度曲線から求められたシェアシ
ニングに関与するエネルギー障壁の大きさE0（●）、及び、EMDシミュレーショ
ンで計算された拡散係数から求められた流動活性化エネルギーの大きさEa,cstr.
（○）の、電荷スケーリング係数cの 2 乗に対する依存性。

RISTニュース

図9

No．69（2023）

平衡状態でのemi 陽イオンのイミダゾリウム環のスナップショット（aおよびd）、及び、定常せん断流下でのイミダゾリウム環の配向（b、c、e、及び f ）。a、b、c は、電荷が完全に取り除かれた状態（c＝0）での結果を表す。d、e、および f は、RSEP電荷がc＝0.8でスケールされた状態での結果を表す。

されてきた。一方、本研究では、温度や電荷を利用して得られたものである。筑波大学藤を変化させることによって、直接的に捉える

田淳一教授、ボノー・パトリック・アラン氏、ことの難しい障壁エネルギーや活性化エネル

及び、
牛山浩氏に多大なる支援を頂いたこと、ギーに焦点を当てて調査し、イオン液体の

ここに感謝申し上げる次第である。

シェアシニングの仕組みを明らかにすること
参考文献

ができた。
本研究の結果は、イオン液体のレオロジー［ 1 ］ K. Matsumoto et al., Syntheses,

特性の研究においては、イオン間の静電相互

structures and properties of 1-ethyl-作用に注目することが重要であることを示し

3-methylimidazoliumている。ただし、イオン分子に結合させる官

fluorocomplex anions, Dalt. Trans.

能基によって、イオン液体の種類は非常に多(2004) 144–149.

salts

of

［ 2 ］ H. Tokuda et al., Physicochemical

く存在するため、今後、様々なイオン液体に

properties and structures of room

ついて検証する余地が残されている。

temperature

本研究における以上の結果は、英文学術誌

ionic

liquids. 2.

variation of alkyl chain length inにおいて報告された［33］。

imidazolium cation, J. Phys. Chem.
B. 109 (2005) 6103–6110.

謝辞
本研究は、防衛装備庁安全保障技術研究［ 3 ］ H. Tokuda et al., Physicochemical

推進制度JPJ004596の支援を受けて実施しproperties and structures of room

た。本研究成果は、筑波大学計算科学研究セ

temperature

ンターの学際共同利用プログラム（cygnus）

Variation of anionic species, J. Phys.

－27－

ionic

liquids. 1.

RISTニュースNo．69（2023）

Chem. B. 108 (2004) 16593–16600.

［12］ N.-T. Van-Oanh et al., Viscosity of［ 4 ］ T. Fukushima et al., Molecular

the 1-ethyl-3-methylimidazolium

ordering of organic molten saltsbis (trifluoromethylsulfonyl) imide

triggered by single-walled carbon

ionic liquid from equilibrium andnanotubes, Science (80-. ). 300 (2003)

nonequilibrium molecular dynamics,

2072–2074.

Phys. Chem. Chem. Phys. 12 (2010)［ 5 ］ R. Pamies et al., Rheologicalbehavior of multiwalled carbon
nanotube-imidazolium

930–936.
［13］ K. Bernardino & M. C. C. Ribeiro,

tosylate

Relating the structure and dynamics

ionic liquid dispersions, Journal ofof ionic liquids under shear by

Rheology 61 (2017) 279–289.

means of reverse non-equilibrium［ 6 ］ M.-D. Bermúdez et al., Ionic liquids

molecular dynamics simulations,as advanced lubricant fluids,

Phys. Chem. Chem. Phys. 23 (2021)

Molecules 14 (2009) 2888–2908.

13984–13995.

［ 7 ］ G. Zarca et al., Non-Newtonian［14］ H . S p i k e s&

Z. Jie, History,

shear-thinning viscosity of carbon

originsmonoxide-selective ionic liquid

elastohydrodynamic friction, Tribol.

1-hexyl-3-methylimidazolium

Lett. 56 (2014) 1–25.

chloride

doped

Separation

with

and

prediction

of

CuCl,

［15］ Y.-L. Wang et al., Microstructural

Purificationand dynamical heterogeneities inTechnology 155 (2015) 96–100.
［ 8 ］ T . Yamag uc h i

and

et

al.,

ionic liquids, Chem. Rev. 120 (2020)

Shear

5798–5877.

Relaxation of Water–Ionic Liquid［16］ M.D. Hanwell et al., Avogadro: an

Mixtures, Bulletin of the Chemical

advanced semantic chemical editor,

Society of Japan 85 (2012) 701–705.

visualization, and analysis platform,［ 9 ］ T. Yamaguchi, Coupling between

J. Cheminform. 4 (2012).

the mesoscopic dynamics and shear［17］ J. Wang et al., Development andstress of a room-temperature ionic

testing of a general Amber force

liquid, Phys. Chem. Chem. Phys. 20

field, J. Comput. Chem. 25 (2004)(2018) 17809–17817.

1157–1174.

［10］ T. Yamaguchi, Shear Thinning and［18］ J. Wang et al., Automatic atom

Nonlinear Structural Deformation

type and bond type perception inof Ionic Liquids with Long Alkyl

molecular mechanical calculations,

ChainsJ. Mol. Graph. Model. 25 (2006) 247–

Studiedby

Molecular

Dynamics Simulation, J. Phys.
Chem. B. 123 (2019) 6260–6265.

260.
［19］ J. De Andrade et al., Computational［11］ R. Futamura et al., Partial breaking

study of room temperature moltenof the Coulombic ordering of ionicsalts

liquids conﬁned in carbon nanopores,

methylimidazolium cations - Force-Nat. Mater. 16 (2017) 1225–1232.

field proposal and validation, J.
－28－

composedby

1-alkyl-3-

RISTニュース

Phys. Chem. B. 106 (2002) 13344–

No．69（2023）

Cambridge, 2017.
［27］ A.W. Lees & S.F. Edwards, The

13351.
［20］ Z. Liu et al., A refined force

computer

study

of

transport

field for molecular simulation ofprocesses under extreme conditions,imidazolium-based ionic liquids, J.

J. Phys. C Solid State Phys. 5 (1972)

Phys. Chem. B. 108 (2004) 12978–

1921–1928.
［28］ S . J . P l i m p t o n ,12989.
［21］ M.J. Frisch et al., Gaussian 09,

Algorithms(2009).

Fast

for

Parallel

Short-Range

Molecular Dynamics, J. Comput.

［22］ K.G. Sprenger et al., The general

Phys. 117 (1995) 1–42.

AMBER force field (GAFF) can［29］ A.I. Jewett et al., Moltemplate: A

accurately predict thermodynamic

Tool for Coarse-Grained Modeling ofand transport properties of many

Complex Biological Matter and Softionic liquids, J. Phys. Chem. B. 119

Condensed Matter Physics, J. Mol.

(2015) 5882–5895.

Biol. 433 (2021) 166841.

［23］ A. Noda et al., Pulsed-gradient［30］ L. Atmani et al., From cellulose tospin-echo 1H and 19F NMR ionickerogen: Molecular simulation ofdiffusion coefficient, viscosity,a geological process, Chem. Sci. 8

and ionic conductivity of non-(2017) 8325–8335.

chloroaluminate room-temperature［31］ T. Shende et al., An empiricalionic liquids, J. Phys. Chem. B. 105

equation for shear viscosity of shear(2001) 4603–4610.

thickening fluids, J. Mol. Liq. 325［24］ L. Martìnez et al., Software News

(2021) 115220.

and Updates - Packmol: A package［32］ T. Yamaguchi, Stress-structure

for building initial conﬁgurations forcoupling and nonlinear rheology

molecular dynamics simulations, J.

of Lennard-Jones liquid, J. Chem.

Comput. Chem. 30 (2009) 2157–2164.

Phys. 148 (2018).

［25］ R. Hayes et al., Structure and［33］ T. Yamada, P. A. Bonnaud., S.

nanostructure in ionic liquids,

Tejima, J. Fujita, Modelling shear

Chemical reviews 115 (2015) 6357–

thinning of Imidazolium-based ionic6426.

liquids., Chemical Physics Letters［26］ B . D . T o d d&P.J. Daivis,

816 (2023).

Nonequilibrium molecular dynamics:

https://doi.org/10.1016/j. cplett. 2023.

Theory, algorithms and applications,

140387

Cambridge

University

Press,－29－