第 6 章 KAPSEL
コロイド・微粒子分散系シミュレータ
山本量一 （ヤマモト リョウイチ） 博士（工学）京都大学 工学研究科 化学工学専攻

教授

6.1

KAPSEL とは

KAPSEL（Kyoto Advanced Particle Simulator for ELectro-hydrodynamics）は，微粒子の沈降ダイナミクスやコロイド分散系のレオロジー，
荷電コロイド粒子の電気泳動など，
粒子分散系の複雑な動的現象に対して基礎方程式に基づいた高精度の直接数値計算を高効率で実現することを目的とし，独立行政法人科学技術振興機構（さきがけ，CREST）からの財政的支援のもと京都大学で開発されたソフトウェアである．シミュレーションの立場から見た微粒子分散系の特徴は，注目すべき空間・時間のスケールが原子・分子のミクロなスケールと比べて桁違いに大きいことにある．微粒子の運動だけでなく，周囲の流体やイオンによる大規模な協調運動が起こるため，微視的な時間スケール（～10-10 秒）とはか
け離れた 10-3～103 秒にも達する緩和時間を示す．このような微粒子分散系の長時間の緩和現象は，通常の分子動力学シミュレーション等で解ける最長の時間範囲（～10-6 秒）に比べて著しく長いため，ミクロなモデルに基づいた分子動力学法でコロイド分散系全体を十分な時間シミュレーションすることは不可能である．そのため微粒子・コロイド分散系
では，独自の粗視化モデルに基づいたメソスケールシミュレーション法が発展した．
コロイド・微粒子分散系のシミュレーション法としては，Brownian Dynamics（BD）
法[1]や Stokesian Dynamics（SD）法[2]など，流体運動を顕に計算しない方法が使われてきた．比較的計算が軽いことがこれらの方法の長所であるが，溶媒を単純なニュートン
流体としているため，高分子流体や電解質溶液など内部自由度のある複雑流体に粒子が分散している系への拡張が困難であった．また分散粒子の大きさが十分に小さく，粒子サイズで定義した Reynolds 数が非常に小さい状況（ Re  1 ）しか扱えなかった．1980 年代

後半に SD 法が登場して以来 10 年以上目立った変化はなかったが，2000 年を境に分散粒子のスケールに合わせた粗視化モデルを用いて溶媒流体を記述し，分散粒子の運動と連動させて顕に扱う直接数値計算（DNS）の開発が活性化し，状況が一変する．コロイド・微粒子分散系の DNS の場合，どのようなモデルで溶媒の運動を記述し，どのような方法で分散粒子と溶媒の運動を連動させるかが鍵となる．前者としては Navier-Stokes（NS）方程式や格子 Boltzmann（LB）法によって流体力学的に溶媒を扱う方法，Multi-Particle Collision（MPC）法や Dissipative Particle Dynamics（DPD）法などのように溶媒を代表粒子で置き換える方法が既に用いられている（ここで挙げた順にモデルはより粗くなる）．
溶媒流体の運動を分散粒子の運動と連動させるためには，溶媒（流体）と粒子（固体）の界面において境界条件を課す必要がある．NS 方程式で溶媒の運動を記述する場合，Finite Element 法(FEM)など粒子形状にフィットした不規則格子を用いるのが最も正当な方法であるが，計算量が膨大になってしまう．規則的な固定格子上で流体の計算を行え
ば効率はよいものの，格子の形状（立方体格子が一般的）と分散粒子（球形）の形状がフィットしないために境界条件の扱いに困難が生じる．この問題を解決するために，固体の粒子を高粘性流体としてモデル化する Fluid Particle Dynamics（FPD）法[3]や，粒子表面に近接する計算格子点上に表面からの距離に応じた補助変数を置いて界面位置を指定す

る方法[4]が提案されている．KAPSEL が実装している Smoothed Profile（SPM）では，
図 6-1 のように粒子の内側と外側を判別するために類似の補助変数を導入し，独自の方法
で分散粒子と溶媒の連動を高精度かつ高効率で実現している[5]．電解質溶液など複雑流体への拡張性も同じ補助変数を用いることで実現している[6,7]．

∑ φ (x), φ (x) = g ( x − R )(1)h[(a + ξ / 2) − x]h[(a + ξ / 2) − x] + h[ x − (a − ξ / 2)](2)

exp ( −∆ 2 / x 2 ) x ≥ 0h( x ) = 
0
x<0


(3)

φ ( x)=Rig ( x) =N

i =1 iii

図 6-1．SPM では補助関数 0 ≤ φi ≤ 1 を持いて半径 a の球状粒子に厚み ξ の境界領域を導入し，計算効率の大幅な向上を実現している． R i (i = 1, , N ) は粒子 i の座標， ∆ は流体計算の格子間隔．

本書では KAPSEL の基本原理に加え，インストールやサンプルシミュレーションの実
行方法，さらに進んだ使用方法などについて解説する．KAPSEL のソースコードやサンプ
ルファイルは，常に最新のものが公式ホームページ（KAPSEL-HP）[8]で公開されており，ライセンスに同意するものであれば誰でも自由にダウンロードして利用することが出来る．
本書とともに KAPSEL-HP の最新情報をぜひ参照して欲しい．
6.2

KAPSEL のインストールと基本操作

6.2.1

OCTA のインストール

Linux の場合，/usr/local/OCTA2013/，Windows の場合 C:¥OCTA2013¥に OCTA2013
がインストールされているものとする．以下の操作は全て，特に指示がない限りコマンド
ラインで行う．Cygwin を使用する場合，最初に以下のコマンドを実行してシンボリック
リンクを作成しておく．
> ln -s /cygdrive/c/OCTA2013/ /usr/local/.
OCTA を正常にインストールした後，以下の様にして libplatform をビルドする．
> su（Linux の場合のみ）> cd /usr/local/OCTA2013/GOURMET_2013/src> ./configure> make> make install6.2.2

KAPSEL のインストール

KAPSEL の最新のソースコード（kapsel3.*.zip）を KAPSEL-HP [8]からダウンロード
し，適当な場所で解凍する．各フォルダの内容は以下の様になっている．
解凍したフォルダ/KAPSEL のソースコード

Doxygen/

ソースコードに関するオンライン文書

Examples/

各種サンプル例題

Tools/

各種ツール

UDF/

動作確認用のファイル一式

まず，解凍したフォルダで以下の様にして，KAPSEL の実行可能ファイルをビルドする．
> unzip kapsel3.*.zip> cd kapsel3.0> make 'ENV=CYGWIN'（Cygwin の場合）

> make 'ENV=GCC'（Linux の場合）

ソースコードを解凍して出来るフォルダーUDF に簡単なサンプルがあるので，以下の様にして動作確認を行う．
> cd UDF> ../kapsel -Iinput.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
KAPSEL が正常にインストール出来ていればシミュレーションがスタートし，30 秒ほどで以下の様な出力が現れて正常終了する．
#Simulation has ended!
#Total Running Time (s):
24.77
#
(m):
0.41
#
(h):
0.01
#Average Step Time (s):
0.02
#
(m):
0.00
#
(h):
0.00
6.2.3

シミュレーション結果の確認（Gourmet を用いた可視化と解析）

上記で実行したサンプルファイルでは，流体（32×32×32 の流体格子）の中に，流体より重い粒子 5 個と流体より軽い粒子 5 個を入れた状況が設定されている．シミュレーシ
ョンの時系列データは"output.udf"に格納されており，Gourmet や Python スクリプトを用いることで容易に閲覧出来る．一例として，Gourmet を用いた粒子運動のアニメーショ
ンの作成手順を説明する．
Gourmet を起動し，File（図 6.2 の 1）→Open を選んで"output.udf"を読み込む．



下部の Python パネルに"particleshow.py"を Load する（図 6.2 の 2）．



右にある Run（図 6.2 の 3）を選択するとビューアが開き，そこで再生ボタン（図 6.2の 4）をクリックすると，アニメーションが再生される．

同様に，Gourmet から任意の Python スクリプトを実行することによって，"output.udf"に格納されたシミュレーション時系列データにアクセスし，解析に用いることが出来る．
6.2.4

Gourmet を用いないデータ解析

大規模なデータを扱う場合，以下の２つの方法で Gourmet を用いずに"output.udf"の時
系列データを解析することが出来る．


Python スクリプト（コマンドラインから実行）
インストール時のフォルダーUDF に，粒子の静的構造因子 S (k ) を計算するサンプル
スクリプト"sk.py"が同包されており，数値計算用の Python パッケージ"numpy"が使用できる環境であれば，以下の手順で使用できる（図 6-3）．
Windows のスタートメニューから，"すべてのプログラム"→"OCTA2013"→"StartGourmetTerm"を起動して，以下のコマンドを入力する．
> python sk.py> gnuplot>> plot 'sk.dat' w line
詳しくはOCTA 付属の「GOURMET PYTHON スクリプトリファレンスマニュアル」
を参照のこと．



Fortran，または C プログラム
詳しくは OCTA 付属の「プラットフォームインターフェースライブラリ libplatform
リファレンスマニュアル」を参照のこと．

図 6-2．Gourmet を用いたシミュレーション結果の可視化（アニメーションの作成）の手順．

図 6-3．粒子の静的構造因子を計算する Python スクリプトの出力例．

6.2.5

実行方法とサンプル例題について

ソースコードを解凍して出来るフォルダ Examples に種々のサンプル例題があるので，各サブフォルダの内容を列挙する．KAPSEL-HP [8]に詳しい説明があるので，詳細についてはそちらを参照して欲しい．
Examples//01

荷電コロイド粒子の電気泳動/02

球状粒子（引力なし）の重力沈降/03

球状粒子（引力なし）の拡散/04

球状粒子分散流体のレオロジー（ジグザグせん断流）/05

球状粒子（引力あり）の凝集/06

フレキシブル鎖状粒子の熱運動・せん断流下の回転運動/07

球状粒子分散流体のレオロジー（Lees-Edwards 周期境界条件）/08

非球形剛体粒子の重力沈降，せん断流下の回転運動/09

自己推進粒子の運動/10

非球形剛体粒子の重力沈降（高レイノルズ数領域）

6.3

微粒子分散系のダイナミクス

微粒子が液体(溶媒)中に分散した系は微粒子懸濁液と呼ばれ，食品・塗料・顔料・化粧品・スラリーなど我々の日常生活に数多く存在している．KAPSEL では，溶媒流体中における分散粒子の拡散，凝集[9]，沈降[10,11]など諸現象の動的シミュレーションに加え，せん断流下での微粒子分散系のレオロジー計算[12,13,14]などの直接数値シミュレーショ
ンを実行することが出来る．
6.3.1

基礎方程式

粘度η ，密度 ρ をもつニュートン流体中に分散した， N 個の半径 a の球状粒子の運動を考える．溶媒は Navier–Stokes 方程式に従う非圧縮ニュートン流体として，流速ベクトル

u(x, t ) の時間変化を追跡し，微粒子は熱揺らぎに由来するランダム力を含んだ Langevin
方程式（並進・回転）に従うものとして，並進座標 R i (t ) ，並進速度 Vi (t ) ，回転座標 Qi (t ) ，
回転角速度 Ωi (t ) の時間変化を追跡する（ i = 1, , N は粒子の番号）．この溶媒と微粒子の連成問題を固定直交格子上で高速かつ正確に解くために，SPM では図 6.1 に示すなめら
かな補助関数 φ =∑ φ （粒子領域で φ = 1 ，溶媒領域でφ = 0 ）を用いて溶媒と粒子の境
N

i =1 i

界を記述している[5,6]．溶媒，微粒子の従う基礎方程式を以下に示す．
i) 溶媒の運動方程式： Navier–Stokes 方程式
溶媒の流速ベクトル u(x, t ) ，圧力 p (x, t ) は以下の式に従う．

ρ (∂ t + u ⋅∇)u = −∇p + η∇ 2u + ρφ f p∇ ⋅u =0（Navier–Stokes 方程式）

(4)（連続の式）(5)

ここで，φ f p は分散粒子表面での境界条件を満たすための拘束力であり，SPM ではこの拘束力が固体の分散粒子の剛直性を保証している．拘束力の具体的な実装法については，既報の論文[5,6]に詳しい記述がある．
ii) 微粒子の運動方程式： Langevin 方程式

i 番目の粒子の運動{R i (t ), Vi (t ), Qi (t ), Ωi (t )} は，粒子の質量を M p ，慣性モーメントをI p として以下の方程式に従う．

 0
 = Q ⋅  Ωz
 = V, Q
R
i
i
i
i 
i
 Ωiy


Ωiz
0
Ω

x
i

Ωiy 

Ωix  ,0 

 = N H + N other + G Ω = F H + F other + GV , I ⋅ Ω
M pV
i
i
i
i
p
i
i
i
i

(6)
(7)

H

ここで，Fi ，N iH は微粒子が流体から受ける力とトルクであり，流体・微粒子間で運動量が保存するように前述の拘束力 φ f p と関連づけられている．この実装法について興味のあ
other

る読者は既報の論文[5,6]を参照して頂きたい．Fi， N iother は流体以外からの力とトルクを表し，微粒子間の直接的な相互作用や，重力・浮力などがこれに相当する．GVi ，G iΩ は
熱揺らぎにより粒子に働くランダム力とトルクであり，平均値がゼロの白色雑音

GVi (t ) ⋅ GVj=(0) 3k BT α V δ (t )δ ij ,G iΩ (t ) ⋅ G Ωj=(0) 3k BT α Ωδ (t )δ ij ,(8)で与える．ここで  は統計平均を表す．
熱揺らぎの導入は，流体の Navier–Stokes 方程式(4)で行うのが最も正当なやり方であるが，KAPSEL では計算効率の観点から粒子の Langevin 方程式で行っている．したがって厳密には揺動散逸定理を満たさないが，それが問題になる時間スケールは粒子の拡散の時間スケールに比べて著しく短いため本質的な問題はない．ただし揺動散逸定理がそのままでは使えないため，(5)式で導入した熱揺らぎがぴったり粒子温度 T を与える保証がなく，正確な温度を決定するためには微調整が必要となる．KAPSEL では，希薄溶媒中で揺らいでいる微粒子の並進拡散係数 DV および回転拡散係数 D Ω から求まる温度を正確な粒子温度 T として定義し，そのため(5)式にランダム力の強度を微調整する補正パラメータ α V おVΩ
よび α Ω を導入している．通常は α=α=1 と設定して問題はないが，精密な温度設定が必要な場合は，以下の様にして補正パラメータを決定する．この温度決定方法の詳細については文献[15,16]を参照されたい．
ΩV
まず，希薄分散系（体積分率 ϕ  1 ）で α=α=1 として設定温度 T の元で熱平衡Ω
Vと回転拡散係数 Dsim
をそれぞ
シミュレーションを行い，微粒子の並進拡散係数 Dsim

れ求める．


Ω
V
V
Ω
， Dtheotyと比較することにより， T=
それらを同じ ϕ の理論値 Dtheoty
T=
T を満足

するように α V と α Ω を調整する．


このようにして決めた α V と α Ω を用いて，温度 T のシミュレーションを任意の ϕ で行
う．

6.3.2

単位についてKAPSEL では，長さの単位として格子幅 ∆ ，時間の単位として流体の密度 ρ と粘性率η
と格子幅で決まる粘性拡散時間τ 0 = ρ∆ 2 / η を採用している．例えば，考察したい空間スケールとして格子間隔 ∆ =1µ m ，溶媒として 水（=η 10−3 Pa ⋅=s, ρ 103 kg / m 2 ）を想定した場合，シミュレーションの単位時間はτ 0 = 10−6 s となり，現実系と対応づけられる．”input.udf”で RHO = A ，ETA = B ，DX = C と設定した場合，シミュレーションの時

2間刻み ∆t は自動的に ∆t= Tdump = ( A / B ) / kmax に設定される． kmax は流体計算の解像度を−1表す最大波数であり C に比例する．通常はこのままで問題ないが，外力やせん断速度が大きい場合には ∆t =Tdump では大きすぎてシミュレーションが異常終了することがある．
そのような場合には，調整パラメータ factor を小さくすることで，∆=t Tdump × factor と設
定する．
6.3.2

扱える粒子の種類について

KAPSEL では，現在以下の種類の分散粒子を扱うことが出来る．




球状粒子


粒子間相互作用

引力 ON/OFF

2n:n の Lenard-Jones タイプポテンシャル（ n =6,12,18 ）鎖状粒子
FENE (finitely extensible non-linear elastic)ポテンシャルによるフレキシブルな直線鎖[17,18]



剛体粒子




複数の球で構成する任意形状

自己推進粒子

球状スクワマー[19]6.3.3. 入力 UDF の説明
Newton 流体中の微粒子分散系のシミュレーションは，constitutive_eq.type としてNavier_Stokes/ Shear_Navier_Stokes/ Shear_Navier_Stokes_Lees_Edward のいずれかを選択する．
constitutive_eq=Navier_Stokes
DX

長さの単位量である格子幅 ∆



RHO

溶媒の密度



ETA

溶媒の粘度



kBT

粒子温度



alpha_v

粒子の並進温度に対する補正項alpha_o

粒子の回転温度に対する補正項

constitutive_eq= Shear_Navier_Stokes（せん断流あり・ジグザグ流れ）DC

定常せん断流shear_rate

定常せん断速度AC

振動せん断流shear_rate

振動せん断流frequency

振動数

constitutive_eq= Shear_Navier_Stokes_Lees_Edwards
DC

定常せん断流shear_rate（せん断流あり・LE_PBC）
定常せん断速度分散粒子の種類は，object_type.type として spherical particle/chain/rigid のいずれかを選択する．
object_type.type= spherical particle
Particle number 粒子数



MASS_RATIO

粒子密度／流体密度janus_*

自己推進粒子で使用する（本書では省略）

object_type.type= chainBeads_number

一本の鎖に属するビーズの数

Chain_number

鎖の本数



MASS_RATIO

ビーズ密度／流体密度

object_type.type= rigid
Beads_number

一つの剛体を構成する球の数



Chain_number

剛体の数



MASS_RATIO

球密度／流体密度

Rigid_motion

free：自由運動

fix：指定した速度・角速度で運動

粒子半径や界面厚さなどは，共通パラメータとして以下で設定する．
A_XI

界面の厚さ ξ



A

粒子半径



gravity.G

重力加速度



gravity.direction 重力を加える方向

EPSILON

Lennard–Jones ポテンシャルのエネルギーの大きさ

LJ_powers

Lennard–Jones ポテンシャルのべき指数

mesh.NPX

シミュレーションセルの x 方向のサイズ： Lx = 2 NPXmesh.NPY

シミュレーションセルの y 方向のサイズ： Ly = 2 NPYmesh.NPZ

シミュレーションセルの z 方向のサイズ： Lz = 2 NPZtime_increment.type2
auto： ∆T= Tdump × factor に設定（ Tdump = ρ (η kmax ) ）．

manual：任意の値に設定
以下でシミュレーションの実行条件を入力する

switch.ROTATION

粒子の回転運動を考慮する場合は ON（デフォルト）swithch.LJ_truncate

ON：粒子間の引力なし
OFF：粒子間の引力あり
NONE：相互作用しない（幽霊粒子）

switch.INIT_distribution

粒子の初期配置を以下の様に作成する
uniform_random：ランダムに作成
random_walk：正方格子からランダムに変位

FCC：FCC 格子

BCC：BCC 格子上

user_specify：座標と速度を input.udf 内で指定．
user_specify の場合，初期位置と速度を user_specify.Particles[].R と
user_specify.Particles[].v で指定する．入力するリストの数がParticle_number で指定した数より小さければ,，Gourmet で Edit->Add an array
Element として user_specify.Particles[]を増やすか直接 UDF ファイルを編集する．剛体粒子の場合はここで粒子形状を定義する．
switch.FIX_CELL

ON：系全体のドリフト速度をゼロ（ uˆ k =0 = 0 ）とする．



output.GTS

データ出力のインターバルのステップ数



output.Num_snap

データ出力の回数，総ステップ数は GTS × Num_snap



output.AVS

AVS 形式のデータ出力をする場合は ONoutput.AVS.ON.Out_dir

AVS 形式のデータ出力をするディレクトリを入れる．たとえば data と指定した場合，事前に./data および./data/avs/のディレクトリ
を作成しておく必要がある．AVS の field ファイルは./data/に data.fld というファイル名 で出力される．


データファイルは./data/avs/に data_*.dat というファイル 名で出力される．*にはステップ数が入る．



output.AVS.ON.FileType

output.UDF

AVS データファイル形式：Binary/ASCII“output.udf”に時系列データを出力する場合は ON

6.4．荷電コロイド粒子の電気泳動
水のような誘電率が非常に大きい溶媒にコロイド粒子が分散しているとき，コロイド表面にある解離基からイオンが放出されて粒子表面は電荷を帯びる．放出されたイオンは粒子表面に静電的に引き寄せられるが，同時に熱ゆらぎによって拡散し，コロイド粒子周りに電気二重層とよばれるイオン雰囲気を形成する．平衡状態の荷電コロイド分散系の性質
が Poisson–Boltzmann 方程式によって記述されるのに対して，電気泳動をはじめとした
界面動電現象では，粒子とイオン分布の挙動は流体力学と静電相互作用の兼ね合いで決まる．両者の競合の結果，コロイド周りのイオン分布が球対称から歪み分極することや，イ
オン分布が外部電場やコロイド粒子の運動に追随できないことも起こり得る．これまでこのような複雑な状況を計算機シミュレーションで正しく再現し，定量的に解析したことはほとんどなかった．
KAPSEL が実装している SPM では，コロイド分散系の動的挙動を決定するうえで重要な流体力学を Navier–Stokes 方程式に基づいて，また同時に静電相互作用を Poisson 方程式に基づいて，いずれも基礎方程式に忠実に直接数値計算を行っている．
6.4.1. 基礎原理
文献[6,7,20]にしたがいシミュレーションに必要な電気流体力学の基本方程式を解説す
る．電解質溶媒中に分散している半径 a の球状コロイド粒子 N 個を考える．また溶媒の誘電率  は，コロイド内部も含めて空間的に一様であるとする．コロイド粒子の表面は一様に帯電していると仮定し，1 粒子あたりの帯電量は Ze であるとする．通常の連続体描像では，コロイド粒子の電荷分布は表面に局在したデルタ関数で表される．そのため有限要素法などでは，適切な境界適合格子が用いられ計算効率の大きな枷となっている．それに対し，SP 法ではこの粒子表面電荷分布 eq (x) についても滑らかな分布を用いる．例えば界面補助関数 φ の 1 階微分を用いて

eq (x) =Ze | ∇φ (x) |4π a 2(9)と表現できる．界面関数 φ (x) が ξ → 0 でステップ関数に帰着するのと同様に， q (x) は

ξ → 0 でデルタ関数に帰着する．
i) イオンの流拡散方程式
価数 Zα をもつ α 種イオンの密度分布 Cα を

Cα (x, t ) = (1 − φ (x, t ))Cα∗ (x, t )(10)として計算領域全体で定義する．本来イオンが存在しないコロイド内部は (1 − φ ) ≈ 0 であり実質的に排除されている．コロイド内部における Cα∗ に物理的な意味はない．Cα∗ (x, t ) は計算効率を上げるための補助変数であり，計算領域全体（あらゆる x ）で滑らかな関数である．粒子表面電荷分布を含めた全電荷分布は

=ρe (x) e∑ Zα Cα (x) + eq (x)

(11)

α

である．補助イオン密度 Cα∗ の時間発展は移流拡散方程式:∂ t Cα∗ = −∇ ⋅ Cα∗ u + Γα ∇ ⋅ (Cα∗ ∇µα )(12)に従う．この式は，溶媒の流動 u による移流と，化学ポテンシャル µα の勾配による拡散の2 つの項からなる． Γα は α 種イオンの Onsager 輸送係数であり，イオンの摩擦係数

fα = 1 / Γα ，あるいは拡散係数 D=
k BT Γα とそれぞれ関係がある．無次元化した摩擦係α
数を mα = 2
k BTfα 3η e 2 で定義すると，25℃の KCl 水溶液の場合は mK +  mCl+  0.184 である．イオンの化学ポテンシャルは(13)µα k BT ln Cα∗ + Zα e(Ψ − E ⋅ x)=で与えられる．ここで E は外部電場を表し，静電ポテンシャル Ψ (x) は Poisson 方程式

∇ 2 Ψ = − ρe の解である．
ii) 溶媒の Navier–Stokes 方程式
溶媒流動は非圧縮性の流れ( ∇ ⋅ u = 0 )である．溶媒の速度 u は Navier–Stokes 方程式:ρ (∂ t + u ⋅∇)u = −∇p + η∇ 2u − ρe (∇Ψ − E ) + ρφ f p(14)にしたがう．中性粒子の式(4)とは異なり，電気浸透力 − ρ e (∇Ψ − E ) が含まれることに注意する．
iii) コロイド粒子の運動方程式

i 番目のコロイド粒子の運動 {R i (t ), Vi (t ), Qi (t ), Ωi (t )} は，粒子の質量を M p ，慣性モーメントを I p として以下の方程式に従う．これらは中性粒子の場合と同じであるが，ここではコロイド粒子の熱運動（ブラウン運動）は考慮していない．

 0
 = Q ⋅  Ωz
 = V, Q
R
i
i
i
i 
i
 Ωiy


Ωiz
0
Ω

x
i

Ωiy 

Ωix  ,0 

 = N H + N other = F H + F other , I ⋅ Ω
M pV
i
i
i
p
i
i
i
6.4.2

(15)
(16)

電気二重層の性質

電気二重層の構造を定量的に扱うための基礎となる Poisson–Boltzmann 方程式について解説する．原点 r = 0 に荷電コロイド粒子が１つ存在し，外部電場がない( E = 0 )とき，式(12)のもとでイオンの平衡分布を求める．化学ポテンシャル µα が一様なとき，平衡イオン分布として Z eΨ (r ) =Cα∗ (r ) C α exp  − α(17)k BT 

を得る．これは静電ポテンシャル Ψ のもとでの Boltzmann 分布である．
正負イオンの電荷が ± z の対称電解質溶媒を考えると，Poisson–Boltzmann 方程式は
 zeΨ (r ) 
2 zeC∇ 2 Ψ (r ) =sinh 

 k BT (18)となる．無限遠方での境界条件は Ψ |r =∞ =0 および C ∗ |r =∞ =C である．コロイド粒子表面の電荷分布が一定の境界条件を用いると∇Ψ |surface = −

σe


,

σ=
e Ze / 4π a 2(19)である．いま zeΨ / k BT  1 を仮定し，式(18)を線形化した Debye–Hückel 近似を行うと∇ 2 Ψ (r )=2 z 2 e 2C
Ψ= κ 2 Ψ
k BT 

(20)

1

(21)

となる．長さの次元を持つ定数κ −1 =

8πλB z 2C

は Debye 遮蔽長と呼ばれ，電気二重層の厚みを表す．同様に，α 種のイオンの価数が Zα である一般の電解質の場合は

κ −1 =

1
4πλB ∑ α Zα2 C α(22)となる．ここで λ=e 2 / 4π k BT  は Bjerrum 長と呼ばれ，25℃の水では約 0.72nm である
B（KAPSEL では ∆ =4πλB に固定）．いま，系は球対称なので動径方向 r =| r | についてのみ考えればよいので，式(20)はd 2Ψ 2 d Ψ
+
=
κ 2Ψ
2
dr
r dr

(23)

に帰着する．この一般解は Yukawa 型ポテンシャル

Ψ (r ) =Ψ0aexp[−κ (r − a )]r

(24)であり，コロイド電荷による静電力が κ −1 程度で遮蔽されることを表している．κ −1 は温度が高いほど大きくなり，イオン強度∑α Zα C α / 2 が大きいほど小さくなる．
2

6.4.3

UDF の説明

constitutive_eq として Electrolyte を選ぶと荷電コロイドの分散・凝集や，外部電場の下での電気泳動現象をシミュレートすることができる．中性微粒子の場合（6.3.3）と共通のパラメータについては説明を省略する．
constitutive_eq=Electrolyte
kBT

系の温度 T （イオンの分布に関係する）



Dielectric_cst

溶媒の誘電率

Init_profile

イオン分布の初期設定 Uniform／Poisson_Boltzmann

Add_salt

saltfree：対イオンのみ
salt：対イオン＋正負 2 種類の塩イオン



Add_salt.saltfree.Valency_counterion

対イオンの価数

Add_salt.saltfree.Onsager_coeff_counterion

対イオンの Onsager 輸送係数



Add_salt.salt.Valency_positive_ion

正イオンの価数

Add_salt.salt.Valency_negative_ion

負イオンの価数

Add_salt.salt.Onsager_coeff_positive_ion

正イオンの Onsager 輸送係数



Add_salt.salt.Onsager_coeff_negative_ion

負イオンの Onsager 輸送係数



Add_salt.salt.Debye_length

Debye 遮蔽長Electric_field

外部電場 ON／OFF

Electric_field.ON

DC：直流電場 ／AC：交流電場

Electric_field.ON.DC.Ex

x 方向電場の強さElectric_field.ON.DC.Eyy 方向電場の強さElectric_field.ON.DC.Ez

z 方向電場の強さElectric_field.ON.AC.Frequency

交流電場の周波数

object_type.type には，球状粒子／鎖／剛体を選択することができる．
object_type.type= spherical particle/chain/rigidSurface_charge

time_increment.type

粒子表面電荷の帯電量
auto： ∆T= Tdump × factor に設定（ Tdump は，粘性拡散時間22
ρ (η kmax) とイオン拡散時間1 (k BT Γα kmax
) の小さい方）

manual：任意の値に設定

参考文献[1] D. L. Ermak and J. A. McCammon, J. Chem. Phys. 69, 1352 (1978).
[2] J. F. Brady and G. Bossis, Ann. Rev. Fluid Mech. 20, 111 (1988).
[3] H. Tanaka and T. Araki, Phys. Rev. Lett. 85, 1338 (2000).

[4] T. Kajishima, S. Takiguchi, H. Hamasaki, and Y. Miyake, JSME Int. J., Ser. B 44,
526 (2001).
[5] Y. Nakayama and R. Yamamoto, Phys. Rev. E 71, 036707 (2005).
[6] Y. Nakayama, K. Kim, and R. Yamamoto, Eur. Phys. J. E. 26, 361 (2008).

[7] K. Kim, Y. Nakayama, and R. Yamamoto, Phys. Rev. Lett. 96, 208302 (2006).
[8] KAPSEL-HP: http://www-tph.cheme.kyoto-u.ac.jp/kapsel/[9] Y. Matsuoka, T. Fukasawa, K. Higashitani, and R. Yamamoto, Phys. Rev. E 86,
051403 (2012).

[10] A. Hamid and R. Yamamoto, J. Phys. Soc. Jap. 82, 024004 (2013).
[11] A. Hamid and R. Yamamoto, Phys. Rev. E 87, 022310 (2013).
[12] T. Iwashita and R. Yamamoto, Phys. Rev. E 80, 061402 (2009).
[13] T. Iwashita, T. Kumagai, and R. Yamamoto, Eur. Phys. J. E, 32, 357 (2010).
[14] H. Kobayashi and R. Yamamoto, J. Chem. Phys. 134, 064110 (2011).
[15] T. Iwashita, Y. Nakayama, and R. Yamamoto, J. Phys. Soc. Jap. 77, 074007 (2008).
[16] T. Iwashita and R. Yamamoto, Phys. Rev. E 79, 031401 (2009).
[17] H. Kobayashi and R. Yamamoto, Phys. Rev. E 81, 041807 (2010).
[18] H. Kobayashi and R. Yamamoto, Phys. Rev. E 84, 051404 (2011).
[19] J. J. Molina, Y. Nakayama, and R. Yamamoto, Soft Matter 9, 4923 (2013).
[20] 北原文雄，古澤邦夫，尾崎正孝，大島広行，「ゼータ電位：微粒子界面の物理化学」（サイエンティスト社，1995）．