コロイド分散系のレオロジー・電気流体力学現象シミュレータ

KAPSEL の解説
Kyoto Advanced Particle Simulator for ELectrohydrodynamics
version 1.0
金 鋼1

名嘉山 祥也2山本 量一

京都大学大学院工学研究科化学工学専攻(独) 科学技術振興機構戦略的創造研究推進事業
21 June, 2006

1
2

現所属：分子科学研究所 計算分子科学研究系
現所属：九州大学 工学研究院 化学工学部門

2c
Copyright °2006“Hybrid Molecular Dynamics Simulation for Soft Matters” project , PRESTO-JST &Deptartment of Chemical Engineering, Kyoto University

All rights reserved.

3

目次
第 1 章 KAPSEL のインストール

1.1
1.2

OCTA のインストールと libplatform のビルド . . . . . . . . . . . . . . . . . . . . . . . . . . .
KAPSEL ソースの入手 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.3
1.4

コンパイル方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ユーザサポート . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

第 2 章 懸濁液のレオロジー

2.1

2.2

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 コロイド粒子の運動方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.2 溶媒流動 : シア下の Navier–Stokes 方程式 . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.3 レオロジ−測定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
UDF 解説 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 入力パラメータの指定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 空間単位と時間単位 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
基礎方程式

第 3 章 荷電コロイド粒子の電気泳動

3.1

3.2

3.3

3.4

11
11
11
12
12
14
14
15
17
17
18
18
18
19
19
19
20
20
21
21
22
22

3.4.2

24

空間単位と時間単位 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

実行方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
計算例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25
25

4.2.1
4.2.2

25
26

4.2.3
4.2.4
4.3
4.4

8
9

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.1 移流拡散方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.2 Navier–Stokes 方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.3 運動方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
電気二重層の性質 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Poisson–Boltzmann 方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 Debye–Hückel 近似と Debye 遮蔽長 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
電気泳動の原理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Smoluchowski の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.2 Hückel の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.3 Henry の式および O’Brien–White による解析 . . . . . . . . . . . . . . . . . . . . . . .
UDF 解説 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.1 入力パラメータの指定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
基本方程式

第 4 章 チュートリアル

4.1
4.2

7
7
8

シア下の懸濁液の流動 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1 粒子電気泳動 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
多粒子電気泳動 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

27
27

AVS/Express による可視化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Gourmet による可視化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28
28

正負に帯電した粒子混合系の電気泳動

4
4.5
4.6
4.7

gnuplot によるグラフ化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
再開 UDF の使用方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
既知の障害 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28
28
29

5

KAPSEL とは
KAPSEL(Kyoto Advanced Particle Simulator for ELectrohydrodynamics) とは, コロイド分散系のレオロジーや荷電コロイド粒子の電気泳動などの電気流体力学現象を第 1 原理 (流体力学の Navire-Stokes 方程式や静電力学の Poisson 方程式) に基づいて極力近似を導入しないで計算するために開発されたシミュレータの名称である. 2002 年 11 月から 2006 年 3 月まで, 京都大学理学研究科物理学・宇宙物理学専攻, 及び京都大学工学研究科
化学工学専攻において, 独立行政法人科学技術振興機構1 戦略的創造研究推進事業さきがけ2 「シミュレーショ
ン技術の革新と実用化基盤の構築」領域3 「ハイブリッド型分子動力学シミュレーションの開発」プロジェク
ト4 が実施された. このプロジェクトの主要な目的は, 理論的な取り扱いが困難であるコロイド溶液や懸濁液などのソフトマターに対して, スケールの異なる自由度の間の相互作用を正確に考慮しかつコンピュータ上で扱いやすい形で実現する計算機シミュレーション手法の開発であり, その成果として得られたソフトウェアがKAPSEL である.
KAPSEL のソースコードは WEB5 で完全に公開されており, ライセンス6 に同意するものであれば誰でも自由にダウンロードして利用することが出来る. 本文書では KAPSEL の基本原理に加え, インストール方法やサンプルを用いたシミュレーションの実行に関するチュートリアルを概説する. 最新の情報は WEB にあるので,この文章の記述に従ってうまくいかないときは WEB の最新情報を参照して欲しい.

1 http://www.jst.go.jp/
2 http://www.jst.go.jp/kisoken/presto/
3 http://www.simulation.jst.go.jp/
4 http://www-tph.cheme.kyoto-u.ac.jp/hymd/
5 http://www-tph.cheme.kyoto-u.ac.jp/kapsel/
6 http://octa01.t.u-tokyo.ac.jp/OCTA/license jp.html

7

第 1 章 KAPSEL のインストール
1.1

OCTA のインストールと libplatform のビルド

KAPSEL は, ソフトマテリアルに対する統合的なシミュレータとして開発された OCTA1 内部にある Gourmetとよばれるユーザインターフェースと連携し入力パラメータの管理や出力データの可視化をおこなっている.
そこで, KAPSEL を使用するためには OCTA がインストールされていることが前提となっている. OCTA の最新バージョンである OCTA2005 は OCTA ホームページ http://octa.jp/より CR-ROM の配布サービスを申し込んで入手することができる. OCTA の前バージョンである OCTA2003 は OCTA のダウンロードサイト2 よりオンラインで入手可能であるので, こちらを用いても構わない. OCTA2003 を用いる場合には, このあとに出てくるビルドの手順において, 文字列”2005” を全て”2003” と読み替える必要がある. KAPSEL は UDF
ファイルを介した変数の授受を Gourmet に付属しているライブラリ libplatform を用いておこなっており, 各
マシンで libplatform のビルドを行なわなければならない [1].

LinuxLinux 版の場合,% cd /usr/local/OCTA2005/GOURMET 2005
% tar xvfpz gourmet∗ src.tar.gz% cd src% make% make installでビルドできる. ただし標準の Makefile はCC= gcc
CXX= g++となっているので, あとで KAPSEL のビルドに Intel C++コンパイラ3 を使うようなら libplatform も Intel
C++コンパイラでビルドしなければならない. icc のバージョンによってオプションは• icc7.1 の場合% make WITHFC=no CC=icc CXX=icc• icc8.0 の場合% make WITHFC=no CC=icc CXX="icpc -U GNU "• icc8.1, icc9.0 の場合% make WITHFC=no CC=icc CXX="icpc -lcxaguard"とすればよい. make 時にたくさん警告メッセージが出るかも知れないが, 無視しても問題ないようである.

1 http://octa.jp/
2 http://octa.jp/OCTA/download.html → /pub/old/OCTA2003/
3 http://www.intel.com/cd/software/products/asmo-na/eng/linux/

8

第1章

KAPSEL のインストール

Windows
Windows 版の場合, まず UNIX を Windows 上でエミュレートするフリーソフト Cygwin4 をインストール
する. この際に gcc, gcc-mingw, make を含めることを忘れてはいけない. その後 Cygwin のシェルを起動し,以下のようにして libplatform をビルドする.

% ln -s /cygdrive/c/OCTA2005/ /usr/local/.
% cd /cygdrive/c/OCTA2005/GOURMET 2005
% ./GOURMET2005 SRC.EXE (or tar xvfpz gourmet∗ src.tar.gz)% cd src% make -f Makefile.cygmingw% make -f Makefile.cygmingw install% cd ../lib/cygwin% ln -s ../win32/libplatform.a . (Gourmet 2003 の場合は不要)

MacOS XOCTA ホームページ内の OCTA BBS より MacOS X 版の Gourmet 2005 をダウンロードする5 . ただし
OS のバージョンによってソースファイルが異なるので注意する. MacOS X 10.4(Tiger) の場合は BBS forJapanese -> GOURMET -> GOURMET for MacOS X - >SEL[36] から, MacOS X 10.3(Panther) の場合は BBS
for Japanese -> GOURMET -> GOURMET for MacOS X -> SEL[35] からソースファイルをダウンロードできる. /Application/OCTA 以下に展開し, PATH や PF FILES など環境変数を適宜変更する. その後, Linux 版と同様に gcc を用いてビルドを行う.
% cd /Application/OCTA/GOURMET 2005/src% make% make install

1.2

KAPSEL ソースの入手

KAPSEL の最新バージョンのソースは以下の URL から入手可能である.
http://www-tph.cheme.kyoto-u.ac.jp/kapsel/

1.3

コンパイル方法

KAPSEL のソースは kapsel.tar.gz という名前の gzipped tar アーカイブファイルとして配布され,% tar xvfpz kapsel.tar.gzとすれば kapsel のディレクトリ以下にファイルが展開される.

KAPSEL プログラムのコンパイルは付属の Makefile でおこなう.
% cd kapsel% makeとすれば Linux でも Windows(Cygwin 環境) でも GNU C コンパイラ (gcc) でコンパイルされ kapsel という
バイナリファイルができる. もし Linux 上で Intel C++コンパイラでコンパイルする場合は Makefile の編集をする. この Makefile では% uname -n

4 http://cygwin.com/
5 公式には OCTA2003 は MacOS X Tiger(10.4) には対応していない.

1.4. ユーザサポート

9でホスト名を取得するがそれが, たとえば icc8.1 あるいは icc9.0 がインストールされた新規マシンのホスト
名が hoge であればifeq ($(NODE),hoge)
CC = icc
CXX = icpc
CCOPT = -tpp7 -axN -O3
LINKS = -lm -lplatform -lstdc++
endifの部分を書き加えればよい. その後% makeとすれば icc でコンパイルされた kapsel というバイナリファイルができる.
Mac OS X の場合は Makefile に以下を追加し, make を行う.
ifeq ($(NODE),hoge)
ARCH = macosx
CC = gcc
CXX = g++
CCOPT = -DMACOSX -O3 -fno-inline
LINKS = -lm -lplatform
endif

1.4

ユーザサポート

OCTA ホームページ http://octa.jp の BBS に KAPSEL の掲示板が設けてあります. KAPSEL に関するご質問, バグレポートなどがある場合はそちらに投稿してください.

11

第 2 章 懸濁液のレオロジー
懸濁液のレオロジーは応用上重要であるが，それについて理論的に一般論を構築するにはあまりに複雑である (例えばテキスト [2] や比較的最近のレビュー (日本語) [3, 4] を参照.)．粒子濃度が希薄な極限においては懸濁液のみかけの粘度 ηef f は

ηef f
η

5
= 1 + ϕ,2(2.1)

のように増加するという理論的な結果が得られている (例えば [2] を参照)．η は溶媒の粘度，ϕ はコロイド粒子の体積分率である．これは有名な Einstein の粘度式である．しかし Einstein の粘度式が成立するのは体積分率で ϕ ¿ 10% 程度である．粒子濃度が高くなると，粒子間相互作用の影響とそれによる粒子構造変化、さらにガラス化・結晶化などが起こる．そのため懸濁液のレオロジーを基本法則から理論的に予測することは途
端に難しくなる．以上のような難しさはあるものの，その重要性から実験的な研究は数多くなされ，それを整理する形でみかけの粘度についての経験式・半経験式がいくつか提案されている (Doughherty–Krieger 式 [5],文献 [6, 4] およびその引用文献参照)．
粒子濃度が高くなっていくと，コロイド体積分率以外の多様な効果が現象に関与するため，統一的な理論を作るよりも現象個別に対する解析が必要になってくる．Smoothed Profile (SP) 法は溶媒によるコロイド粒子間多体相互作用を直接計算できるシミュレーション法であり，これを用いて多様な種類の懸濁液のレオロジーを解析することが可能である．KAPSEL では，Newton 流体中の単分散系について単純せん断流下のシミュレー
ションおよびレオロジー測定が行える．

2.1

基礎方程式

Newton 流体中の N 個の球状粒子の運動を考える．

2.1.1

コロイド粒子の運動方程式

質量 Mp をもつ i 番目のコロイド粒子の位置、速度および角速度 {Ri , V i , Ω} は運動方程式:

Ṙi

=

V i,

(2.2)

Mp V̇ i

=

other
FH
,
i + Fi

(2.3)

=

NH
i ,

(2.4)

I p · Ω̇i

によって時間発展する. I p = (2/5)Mp a2 I はコロイド粒子の慣性モーメントで I は単位テンソルである．ここH
で FHi ，N i はコロイド粒子が流体から受ける力およびトルクであり，これらは固液間での運動量収支を表す．

この計算の実装法は [7, 8] を参照されたい. また F otherには流動による力以外による力であり，コロイド粒
i
子間の直接相互作用や浮力などを表す．

12

2.1.2

第 2 章 懸濁液のレオロジー

溶媒流動 : シア下の Navier–Stokes 方程式

溶媒流動は非圧縮性の流れ (∇ · v = 0) である. 溶媒およびコロイド粒子全系の速度場 v は Navier–Stokes 方
程式:

1(∂t + v · ∇)v = − ∇p + ν∇2 v + φf p + f shear ,
ρ

(2.5)

にしたがう. ρ, η = ρν, p はそれぞれ溶媒の密度, 粘性係数, 圧力である. φf p は連続体溶媒中において φ = 1 の領域にコロイド粒子が存在することによる境界条件を力として表している．φf p によってコロイド粒子表面の粘着条件が考慮される．この実装法については文献 [7, 8] を参照されたい．
単純せん断流の駆動および境界条件

外力 f shear は vx を単純せん断流のプロファイル
¶
µ
Ly
−y ,
vx = S2(2.6)に保つように制御する外力である．ここで S はシアレートである．この駆動法は一様流体に対して考案された

Schumacher & Eckhardt [9] の方法を Smoothed Profile (SP) 法に応用したものである 1 ．このとき、流れ (x−)
方向およびスパン (z−) 方向に周期境界条件とし、上下の壁には free-slip 境界条件 (∂y vx,z (y = 0, Ly ) = 0) が課されている．また上下の壁は非透過つまり ∂y vy (y = 0, LY ) = 0 である．SP 法では, 計算の毎ステップにおいてコロイド粒子および流体の移流が済んだ後に f shear がフィードバックとして印加される．以上がシミュ
レーションに必要な方程式である.

2.1.3

レオロジ−測定

懸濁液のみかけの粘度は

ηef f

=

ef f
σxy
,
S

(2.7)

ef f
によって定義される．S(シアレート) は制御パラメタであり，マクロなせん断応力 σxyは計算データから評価される．懸濁液に対する応力の計算法を導く．説明を容易にするため、まず一様流体の場合を考える．続いてコロイド分散系 (SP 法) の場合の結果を述べる.

一様流体の場合
一様流体の場合のシア流下の運動方程式を

∇·v
ρv̇

= 0,

(2.8)

= ∇ · σ + ρf

shear

,

(2.9)

と書く．f shear は系にせん断をかける単位質量当たりの外力場である。これよりマクロな応力の定義は
Z1dxσ(2.10)σ ef f =Vである。これは局所的な応力のが求められればその体積積分で計算できる。
外力 ρf shear を用いて σ ef f を評価することを考える。2 階のテンソルについての恒等式，

σT= [∇ · (σx)] − x∇ · σ,(2.11)

1 単純せん断流を解く方法は，境界壁の速度／力を境界条件で与える方法や、分子動力学シミュレーションにおける Lees–Edwards 境界条件に相当する Rogallo の動格子法 [10, 11] などがあるが, ここでは格子が歪まず境界層の問題がない本方法を用いた．

2.1. 基礎方程式

13

を用いると

Z

Z
dxσ

=

h
i
T
dx [∇ · (σx)] − x∇ · σ(2.12)

dx [−x∇ · σ]

(2.13)

Z
=
Z

h ³
´i
dx x ρf shear − ρv̇
Z
Z
d
dxρv,
=
dxxρf shear −dt=(2.14)(2.15)

である．ここで表面項を落とし，また非圧縮条件 (∇ · v = 0) を用いた．定常状態を考えると，全運動量の変化はゼロなので、結局

σ ef f

1
VZ
dxxρf shear ,(2.16)∇ · σ suspension + ρt f shear ,(2.17)(1 − φ)ρf + φρp ,(2.18)=となる．

懸濁液 (Smoothed Profile 法) の場合
シア流下の懸濁液の運動量方程式を形式的にd(ρt v) =dt
ρt =と書く．ここに σ suspension は懸濁液に対する応力を表している．この σ suspension は計算できないので、応力の積分による表式 (2.10) は実行することが出来ない．
一方、シア流下の懸濁液の SP 法による速度場の方程式は

v̇

=

1
∇ · σ + φf p + f shear ,ρ(2.19)

であるが、これを式 (2.17) と比較すると

1∇ · σ suspension
ρt=

1∇ · σ + φf p ,ρ

(2.20)である．したがって懸濁液の応力の表現として

σ ef f=1
VZ
dxxρt f shear ,(2.21)を得る．この評価式には位置ベクトルが陽に入っているが、測定したいマクロなせん断応力は，座標系の原点のとり方に依らないものである．これは f shear によって全系の並進運動が駆動されないことと対応する．以下、この事実を確認する．座標原点をずらして x = x0 + x0 とすると
ZZZ
dx (x0 + x0 ) ρt f shear =dxx0 ρt f shear + x0 dxρt f shear .

(2.22)

R
外力 f shear はせん断流れを起こし，全系の並進を起こさないように制御されているので， dxρt f shear = 0 である．SP 法では式 (2.21) を用いてマクロな応力を評価している．

14

第 2 章 懸濁液のレオロジー

2.2

UDF 解説2.2.1

入力パラメータの指定

constitutive eq として Shear Navier Stokes を選ぶとシア流下における Newton 流体中のコロイド粒子の運動のシミュレーションが行える．以下, 入力 UDF ファイルで指定する変数について解説する.
constitutive eq.Shear Navier Stokes 以下で溶媒の情報を指定する.
• constitutive eq.Shear Navier Stokes.DX. . . 長さの単位量である格子幅 ∆.
• constitutive eq.Shear Navier Stokes.RHO. . . 溶媒の密度.
• constitutive eq.Shear Navier Stokes.ETA. . . 溶媒の粘度.
• constitutive eq.Shear Navier Stokes.Shear rate. . . 剪断速度 (シアレート).
particle spec[] 以下には粒子タイプの情報を入れる.
• particle spec[].Particle number. . . 粒子数.
• particle spec[].MASS RATIO. . . 粒子の流体に対する密度比.
粒子半径, 界面厚さは以下に情報を入れる.

• A XI. . . 界面の厚さ ξ. Shear Navier Stokes を選んだ場合, 少なくとも ξ ≥ 1 を選ぶ必要がある.
• A. . . 粒子半径.
gravity 以下で重力に関する情報を入れる.
• gravity.G. . . 重力加速度.
• gravity.direction. . . 重力を加える方向を選ぶ.
粒子間力について Lennard–Jones ポテンシャルのタイプを決める.

• EPSILON. . . Lennard–Jones ポテンシャルのエネルギーの単位を決める.
• LJ powers. . . Lennard–Jones ポテンシャルのべき指数を決める.
mesh 以下でシミュレーションサイズを決める (現在のバージョンは直方体のセルで周期境界条件のみ対応).
• mesh.NPX. . . x 方向のサイズは Lx = 2NPX .
• mesh.NPY. . . y 方向のサイズは Ly = 2NPY .
• mesh.NPZ. . . z 方向のサイズは Lz = 2NPZ .
time increment 以下で時間刻みを決める.
2• time increment. . . auto を選べば時間刻みの上限である Tdump = ρ/ηkmax
を決める. ここで kmax は格子幅 DX で決まる最大波数である. manual を選べば手動で値を入れる.

• time increment.auto.factor. . . 入力する値 factor を使って時間刻みは ∆t = Tdump × factor と決まる.
• time increment.manual.delta t. . . 手動で ∆t の値を入れる.
switch 以下でシミュレーションの各種条件を決める.
• switch.ROTATION. . . 粒子の回転運動の運動方程式を解く場合は ON を選ぶ.

2.2. UDF 解説

15• switch.HYDRO int. . . 流体について Navier–Stokes 方程式を解く場合は Correct を選ぶ. ただし Shear Navier Stokesの場合 Correct 以外の switch は実装されていない.
• switch.Stokes. . . Navier–Stokes 方程式を解く際, 慣性項を考慮すれば with advection を選ぶ. 慣性項を考慮しない Stokes 近似をする場合は w/o advection を選ぶ. ただし Shear Navier Stokes の場合 w/oadvection は実装されていない.
• swithch.LJ truncate. . . 粒子間に Lennard–Jones ポテンシャルによる力について, 引力項を含む通常の形
OFF, 引力項を含まない斥力のみ ON, まったく力を加えない NONE から選ぶ.
• switch.INIT distribution. . . 粒子の初期配置を, uniform random(ランダム), random walk(正方格子上からランダムにずれている), FCC(FCC 格子上), BCC(BCC 格子上), user specify(座標と速度はユー
ザによって指定される) から選ぶ. 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. . . 全速度の直流成分を 0, つまり v k=0 = 0 とする. 重力加速度を導入して沈降速度を測る場合これを ON にする.
output 以下にデータ出力の情報を入れる.
• 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. . . AVS データファイルのフォーマットを Binary か ASCII か選ぶ.
• output.UDF. . . 出力 UDF をする場合は ON を選ぶ.
出力 UDF の Particles[] 以下には各粒子の座標と速度を output.Num snap で指定した数だけレコードデータが保存される.

resume 以下で中断した計算を再計算する再開 UDF について指定する.
• resume.caclucation. . . NEW を選べば新しい計算を始める. 前回計算終了したときの計算ステップ数, 粒子
座標・速度・角速度, 溶媒速度場, イオン密度場の情報を読み込んで計算を再開させたい場合は CONTINUEを選ぶ. 再開 UDF ファイルの resume.CONTINUE.Saved Data 以下にデータはセーブされる. この場合,
前回終了時にできる再開 UDF を入力 UDF として計算を再開する. 使用方法は 4.6 節で詳しく解説する.

2.2.2

空間単位と時間単位

長さの単位としてを格子幅 ∆ を採用している.
時間の単位ついては以下の点に注意する.

• Navier–Stokes 方程式で ρ = η = ∆ = 1 となるような単位系を採用している.

16

第 2 章 懸濁液のレオロジー• 仮に入力 udf ファイルで RHO= A, ETA= B, DX= C と入力した場合, 最大波数 kmax は C を用いて与2
えられ, 時間刻みの上限は運動量の拡散時間より Tdump = (A/B)/kmaxと与えられる. 時間刻み ∆t はTdump × factor で調整する. ただし時間の単位自体は ρ = η = ∆ = 1 の場合と同じである.

17

第 3 章 荷電コロイド粒子の電気泳動
水のような誘電率が非常に大きい溶媒にコロイド粒子が分散しているとき, コロイド表面にある解離基からイオンが放出されて粒子表面は電荷を帯びる. 放出されたイオンは粒子表面に静電的に引き寄せられるが, 同時に熱ゆらぎによって拡散され, コロイド粒子周りに電気二重層とよばれる雲のようなイオン雰囲気を形成する. 平衡
状態のコロイド分散系の性質は Poisson–Boltzmann 方程式によって記述され, これを線型化した Debye–Hückel
近似を用いて見通しの良い理解を得ることができる. これに対して電気泳動をはじめとした界面動電現象とよばれる現象では, 粒子とイオン分布の挙動は流体力学相互作用と静電相互作用の競合によって決定される. 両者の競合の結果としてイオン分布は粒子の運動に追随できず, 電気二重層は球対称から歪み平衡状態とは異なる状況が起こる. このような複雑な状況を解析するには, 理論的には簡単な近似を導入する他なく, さらに計算機シミュレーションによっても正しく再現されたことはこれまでほとんどなかった．
本プロジェクトではこの状況を打破すべく Smoothed Profile (SP) 法という新しいシミュレーション方法の開発を行ってきた [7, 12, 8, 13]. SP 法ではコロイドは粒子描像のまま扱う.

その一方で, 溶媒およびイオンについては連続体として粗視化された密度場として与える. また, 連続体描像において通常階段関数によって表現されるコロイド粒子と溶媒の界面を, 滑らかな界面関数を用いて固定直交格子上で表わすことで数値計算上の困難を克服している. このように粒子描像と連続体描像を併用することで, 両者の自由度とそれらのあいだの相互作用を物理的に正しく記述し, マルチスケールにまたがる複雑な現象の解析を可能にするという大きな特徴を有している. また SP 法では, コロイド分散系の動的挙動を決定するうえで重要な役割を果たす流
体力学相互作用を Navier–Stokes 方程式を直接数値計算することで解決している [7, 8]. この計算の実現のために、滑らかな界面関数を前提としてコロイド粒子 (運動方程式)・イオン分布 (移流拡散方程式)・溶媒速度場(Navier–Stokes 方程式) の 3 つの自由度とそれらのカップリングをコンシステントに解くためのフォーマリズ
ムを構成した. これにより, 定量的な電気泳動シミュレーションに世界で初めて成功した [13]. 以下, KAPSEL
内部で考慮されている基本方程式を説明し, 同時に背景にある理論について簡単に触れる.

3.1

基本方程式

文献 [8, 13] にしたがいシミュレーションに必要な電気流体力学の基本方程式を概観する. 電解質溶媒中に分散している半径 a の球状コロイド粒子 N 個を考える. また溶媒の誘電率 ² は, 空間的に一様であるとする. SP
法では粒子と溶媒の界面を有限な厚さ ξ をもつ滑らかな関数 φ(r) ∈ [0, 1] によってあらわす. 固定直交格子上
で粒子領域を φ = 1, 溶媒領域を φ = 0, 界面領域を 0 < φ < 1 で表現している. この様な直交格子状で φ(r) を用いた方法としては, 他に Tanaka–Araki [14] や Kajishima ら [15] の方法がある. 界面関数を導入することによって, 非構造格子を用いた有限要素法などと比べた場合に圧倒的な計算効率の向上を実現できる. コロイド
粒子の表面は一様に帯電していると仮定し, 1 粒子あたりの帯電量は Ze であるとする. 通常の連続体描像では,
粒子表面電荷分布はデルタ関数を用いて表される. そのため有限要素法などでは, 適切な境界適合格子が用いられ計算効率の大きな枷となっている. それに対して, SP 法ではこの粒子表面電荷分布 eq(r) についても滑らかな分布を用いる. 例として界面関数 φ の 1 階微分を用いると

eq(r) =Ze|∇φ(r)|
,
4πa2

(3.1)

のように表現される. 界面関数 φ(r) が ξ → 0 でステップ関数に帰着するように構成されているのと同様に,

q(r) も ξ → 0 とすればデルタ関数に帰着するように構成されている.

18

3.1.1

第 3 章 荷電コロイド粒子の電気泳動

移流拡散方程式

価数 Zα をもつ α 種イオンの密度分布 Cα を

Cα (r, t) = (1 − φ(r, t))Cα∗ (r, t)(3.2)のように計算領域全体で定義する. イオンが存在する領域は (1 − φ) によって表現されている. Cα∗ (r, t) は計算のための補助変数であり, 計算領域全体で滑らかに定義されている. コロイドの領域 (φ = 1) における Cα∗ に物理的な意味はない. 粒子表面電荷分布を含めた全電荷分布はX
Zα Cα (r) + eq(r)ρe (r) = eである. このとき初期分布は電気的中性条件

R

(3.3)

α

ρe dr = 0 を満たすようにとられる.
∗
補助イオン密度 Cα の時間発展は移流拡散方程式:∂t Cα∗ = −∇ · Cα∗ v + Γα ∇ · (Cα∗ ∇µα ),(3.4)にしたがうとする. これは溶媒速度場 v による移流と, 化学ポテンシャル µα の勾配による拡散の 2 つの項から
R
なる. Cα∗ は移流拡散方程式にしたがうので、 drCα∗ は保存する。
イオンがコロイド粒子内部に浸透しないことは, イオン拡散流速の界面法線方向成分が 0 となる条件として考慮される。すなわち, n · ∇µα = 0 である [8, 13]. ここで n は粒子表面における (外向き) 法線方向ベクトルを表す. 界面関数を用いると n = −∇φ/|∇φ| と表すことができる. また Γα は α 種イオンの Onsager 輸送係数であり, イオンの摩擦係数および拡散係数とは fα = 1/Γα , Dα = kB T Γα のような関係にある. イオンの化学ポ
テンシャルはµα = kB T ln Cα∗ + Zα e(Ψ − E · r)

(3.5)

とする [16]. E は外部電場をあらわし, 静電ポテンシャル Ψ(r) は Poisson 方程式:²∇2 Ψ = −ρe(3.6)

を解いて得られる. この化学ポテンシャルにしたがうイオンは平衡状態において Poisson–Boltzmann 分布と
なる.

3.1.2

Navier–Stokes 方程式

溶媒流動は非圧縮性の流れ (∇ · v = 0) である. 溶媒の速度場 v は Navier–Stokes 方程式:

ρ(∂t + v · ∇)v = −∇p + η∇2 v − ρe (∇Ψ − E) + φf p

(3.7)

にしたがう. ρ, η, p はそれぞれ溶媒の密度, 粘性係数, 圧力をあらわす. ここで流体には静電気力 −ρe (∇Ψ − E)が働いていることに注意されたい. 外力項 φf p は, これは粒子の剛体性を保証するための強制力をあらわしている. すなわち φf p によって粒子表面の粘着境界条件が考慮される. 詳しくは文献 [7, 8] に定義されている.

3.1.3

運動方程式

質量 Mp をもつ i 番目コロイド粒子の位置および速度 {Ri , V i } は運動方程式:

Ṙi
Mp V̇ i

= V i,=other
FH
i + Fi

(3.8)
(3.9)

によって時間発展する. ここで F H
i は流体から受ける力をあらわし, 固液間での運動量収支を表している [8].
また F otherは Lennard–Jones ポテンシャルなど粒子間のポテンシャルによる力である. ここでは省略したがi
粒子の回転運動も同様に考慮されている [8]. 以上がシミュレーションに必要な方程式である.

3.2. 電気二重層の性質

3.2

19

電気二重層の性質

電気二重層の構造を定量的に扱うための基礎は Poisson–Boltzmann 方程式である.

3.2.1

Poisson–Boltzmann 方程式

外部電場がない (E = 0) とき, 式 (3.5) のもとでイオンの平衡分布を求める. 化学ポテンシャルが一様,µα = cst., になったとき, 平衡イオン分布として
Cα∗ (r) = C̄α expµ¶
Zα eΨ(r)−kB T(3.10)

を得る. これは静電ポテンシャル Ψ のもとでの Boltzmann 分布である. これは式 (3.6) とあわせて Poisson–

Boltzmann 方程式とよばれる.

3.2.2

Debye–Hückel 近似と Debye 遮蔽長z : z 対称電解質溶媒中にある 1 個の球状コロイド粒子を考える. Poisson–Boltzmann 方程式は
µ
¶
2zeC̄
zeΨ(r)2∇ Ψ(r) =
sinh
,
²
kB T

(3.11)となる [2, 17]. 無限遠方での境界条件は Ψ|r =∞ = 0 および C ∗ |r =∞ = C̄ である. 粒子表面における境界条
件はσe,(3.12)²
であり, ここで表面電荷密度 σe = Ze/4πa2 である. すなわち表面電荷一定の境界条件を考えている. zeΨ/kB T ¿
1 を仮定し式 (3.11) を線形化する近似は Debye–Hückel 近似と呼ばれる.
∇Ψ|surface = −∇2 Ψ(r) =2z 2 e2 C̄Ψ = κ2 Ψ,
kB T ²(3.13)となる. このとき、長さの次元をもつ定数

κ−1 = p18πλB z 2 C̄,

(3.14)

は Debye 遮蔽長とよばれる. ここで λB = e2 /4πkB T ² はビヨルン (Bjerrum) 長である. 一般の電解質の場合は,

κ−1 = q
4πλB

1
P

,

(3.15)2α Zα C̄αである. 今, 系は球対称なので動径方向 r = |r| についてのみ考えればよく,

d2 Ψ 2 dΨ+= κ2 Ψ
dr2
r dr(3.16)という式に帰着する. この一般解は Yukawa 型ポテンシャルaΨ(r) = Ψ0 exp[−κ(r − a)],r

(3.17)

である. コロイド電荷による静電力は κ−1 の程度で遮蔽される. κ−1 は, 表面電荷に異符号イオンが引き寄せられるクーロン力の効果と, イオンの局在をかき消そうとする熱拡散のつり合う距離と考えることができる. したがって Debye 遮蔽長 κ−1 は電気二重層の厚さとみなすことが出来る. 温度が高く熱エネルギー kB T が大きPい場合 κ−1 は大きくなる. またイオン強度 ( α Zα2 C̄α /2) が大きいほど, イオンによる遮蔽効果が大きくなり

20

第 3 章 荷電コロイド粒子の電気泳動

κ−1 は小さくなる. たとえば z : z 対称電解質溶媒中をとして, バルク塩濃度 C̄ を考えるとする. 媒質を 25 ◦ Cの水とするとビヨルン長は λB = 0.72nm となり, 実際の数値を入れると式 (3.14) は0.3
κ−1 = √ (nm)z C̄(3.18)となる. z = 1, C̄ = 0.1M では κ−1 = 1nm, z = 1, C̄ = 0.001M では κ−1 = 10nm となっている.
表面電位 Ψ(r = a) = Ψ0 は境界条件 dΨ/dr(r = a) = −σe/² から定まる. Ψ と表面電荷密度 σe の関係は,

σe = ²κΨ0 (1 + (κa)−1 ),(3.19)となる. すなわち表面電荷の増加に対して、表面電位が線型に増加する関係を与える. 表面電位が大きいとこ
ろでは Debye–Hückel 近似は適用できない. 実際は表面電荷の増加に対して表面電位はしだいに増加しにくくなる. このときの関係式は Poisson–Boltzmann 方程式を解くことによって得られるが, 1:1 電解質溶媒に対する近似式が Loeb–Overbeek–Wiersema [2] あるいは Ohshima–Healy–White [17, 18] によって提案されている.

3.3

電気泳動の原理

外部から電場 E をかけると, Ze に帯電したコロイド粒子は静電気力 ZeE を受けて動き出す. 粒子は静電気力によって加速するが流体から粘性抵抗力を受け, 両者がつり合ったところで定常速度 V に達し等速運動をす
る. 粘性抵抗力を半径 a の球状粒子に対する 6πηaV の Stokes 抵抗と仮定すると, 力のつりあいは,と書け, 電気泳動易動度は

ZeE = 6πηaV(3.20)

V
Ze
=
E
6πηa

(3.21)となる. しかしこれは適当ではなく, 実際の電気泳動易動度はこの値よりも小さくなる. これはコロイド周囲のイオン雰囲気にも静電力が働くため, イオン雰囲気の運動および静電力が考慮されていないことによる. コロ
イドの泳動はイオン雰囲気の運動も伴うので、泳動速度はその分遅くなる. さらに電気二重層が球対称から歪めばそれによる力も考慮しなくてはならなく, 理論的な解析は極めて複雑なものとなる. そこで簡単なモデルについて電気泳動度に対する解析がいくつかなされている [2, 17].

3.3.1

Smoluchowski の式

粒子半径 a が電気二重層の厚さ κ−1 よりも非常に大きい κa À 1 のとき Smoluchowski の式を適用することができる. この場合, 電気二重層厚さが無限に薄い極限となるので, 粒子表面の曲率を無視して, 粒子表面を平板とみなすことができる. 平板に平行な (x 方向とする) 外部電場 Ex が加わっているとする. 粒子に固定した座標系を考えて, 粒子の上の乗って液体の動きを見ると, 粒子から無限に離れたところの流体は −V の速度をもつことになる. 粘性抵抗力と静電気力のつりあいを考えると,

η∂ 2 vx X
eCα Ex = 0
+
∂y 2
α

(3.22)となる. Poisson 方程式によってイオン分布は静電ポテンシャルの 2 階微分と関係しているので,η∂ 2 vx
∂2Ψ
=
²
Ex
∂y 2
∂y 2

(3.23)となる. 無限遠方で速度勾配, ポテンシャル勾配, ポテンシャルそれぞれ 0 という境界条件のもと積分して

η[vx (y) + V ] − ²EΨ(y) = 0(3.24)

3.3. 電気泳動の原理21となる. すると粒子表面 (y = 0) で速度 0 であるので,

V
²ζ
=
E
η

(3.25)という Smoluchowski の式を得ることができる. ここでゼータ電位 ζ は本来はすべり面での電位として定義されるのだが粒子表面における電位 Ψ(0) と置き換えられている.

3.3.2

Hückel の式

Smoluchowski の状況とは逆に, 粒子が電気二重層の厚さに対して非常に小さい κa ¿ 1 の極限を考える. つ
まり Ze が点電荷の極限の場合は Hückel の式を適用することができる. 式 (3.21) において, 粒子表面のポテン
シャルをクーロンポテンシャル

Ze
4π²a

(3.26)

V
2 ²ζ
=
E
3 η

(3.27)

ζ=で与えるとという Hückel の式が得られる.

3.3.3

Henry の式および O’Brien–White による解析

式 (3.25) と式 (3.27) を結び付けるものとして一般の κa に対して Henry の式V²ζ= f (κa)Eη(3.28)

が知られている. f (κa) は Henry 係数といわれ

f (κa)=

1 − 5 exp(κa)E7 (κa) + 2 exp(κa)E5 (κa)·¸
2 (κa)25(κa)3
(κa)4
(κa)5
(κa)4
(κa)6
=
+
−
−
+
+
−
exp(κa)E1 (κa)
3
24
72
144
144
12
144

(3.29)(3.30)

で定義される. ここで En (κa) は n 次積分指数関数である. Smoluchowski の式は f = 1(κa → ∞), Hückel の
式は f = 2/3(κa → 0) に対応する (Fig. 3.1).
1’henry.dat’

0.95

f(κa)

0.9

0.85

0.8

0.75

0.7

0.1

1

10
κa

100

1000

Fig. 3.1: Henry 係数 f (κa).
Henry の式 (3.28) はゼータ電位に対して線形な関係式であり, ゼータ電位が低い場合にのみ適用できる. ゼー
タ電位が高くなると電気二重層の変形の効果 (緩和効果とよばれる) を考慮しなければならない. O’Brien–White は一般の κa, ζ に対して電気泳動度とゼータ電位の関係を数値解析によって提案している [19]. その後,
Ohshima–Healy–White による κa ≥ 10 に対する解析式が提案されている [20].

22

第 3 章 荷電コロイド粒子の電気泳動

3.4

UDF 解説3.4.1

入力パラメータの指定

constitutive eq として Electrolyte を選ぶと荷電コロイドの運動や, 外部電場の下での電気泳動現象をシミュレートすることができる1 . 以下, 入力 UDF ファイルで指定する変数について解説する.
constitutive eq.Electrolyte 以下には溶媒とイオン分布の情報を入れる.
• constitutive eq.Electrolyte.DX. . . 長さの単位量である格子幅 ∆.
• constitutive eq.Electrolyte.RHO. . . 溶媒の密度.
• constitutive eq.Electrolyte.ETA. . . 溶媒の粘度.
• constitutive eq.Electrolyte.kBT. . . 温度 T に対する熱エネルギー.
• constitutive eq.Electrolyte.Dielectric cst. . . 溶媒の誘電率.
• constitutive eq.Electrolyte.Init profile. . . イオン分布の初期設定について Uniform か Poisson Boltzmannを選ぶ. Poisson Boltzmann を選ぶと, 特に多粒子のときに時間がかかることがあるので注意する.
• constitutive eq.Electrolyte.Add salt. . . saltfree であれば粒子表面とは反対符号のイオン 1 種類, salt であれば正負 2 種類のイオンを考慮する.
• constitutive eq.Electrolyte.Add salt.saltfree.Valency counterion. . . 対イオンの価数.
• constitutive eq.Electrolyte.Add salt.saltfree.Onsager coeff counterion. . . 対イオンの Onsager 輸送係数.
• constitutive eq.Electrolyte.Add salt.salt.Valency positive ion. . . 正イオンの価数.
• constitutive eq.Electrolyte.Add salt.salt.Valency negative ion. . . 負イオンの価数.
• constitutive eq.Electrolyte.Add salt.salt.Onsager coeff positive ion. . . 正イオンの Onsager 輸送係数.
• constitutive eq.Electrolyte.Add salt.salt.Onsager coeff negative ion. . . 負イオンの Onsager 輸送係数.
• constitutive eq.Electrolyte.Add salt.salt.Debye length. . . Debye 遮蔽長を指定すると対応する塩濃度が設定される.
• constitutive eq.Electrolyte.Electric field. . . 外部電場の ON, OFF を選択する.
• constitutive eq.Electrolyte.Electric field.ON. . . DC(直流電場) か AC(交流電場) を選択する.
• constitutive eq.Electrolyte.Electric field.ON.DC.Ex. . . x 方向電場の強さ. 交流電場の場合も同じ.
• constitutive eq.Electrolyte.Electric field.ON.DC.Ey. . . y 方向電場の強さ. 交流電場の場合も同じ.
• constitutive eq.Electrolyte.Electric field.ON.DC.Ez. . . z 方向電場の強さ. 交流電場の場合も同じ.
• constitutive eq.Electrolyte.Electric field.ON.AC.Frequency. . . 交流電場の周波数
particle spec[] 以下には粒子タイプの情報を入れる.
• particle spec[].Particle number. . . 粒子数.
• particle spec[].MASS RATIO. . . 粒子の比重.
1 constitutive eq の他の switch として, Newton 流体中のコロイド粒子の運動は Navier Stokes を, シア流下における Newton 流
体中のコロイド粒子の運動は Shear Navier Stokes を選べばシミュレートすることができる.

3.4. UDF 解説

23• particle spec[].Surface charge. . . 粒子表面電荷の帯電量.
粒子半径, 界面厚さは以下に情報を入れる.

• A XI. . . 界面の厚さ ξ. Electrolyte を選んだ場合, 界面関数の 1 階微分 ∇φ を使用しているため, 一番最小の界面の厚さは ξ = 2 である.
• A. . . 粒子半径.
gravity 以下で重力に関する情報を入れる.
• gravity.G. . . 重力加速度.
• gravity.direction. . . 重力を加える方向を選ぶ.
粒子間力について Lennard–Jones ポテンシャルのタイプを決める.

• EPSILON. . . Lennard–Jones ポテンシャルのエネルギーの単位を決める.
• LJ powers. . . Lennard–Jones ポテンシャルのべき指数を決める.
mesh 以下でシミュレーションサイズを決める (現在のバージョンは直方体のセルで周期境界条件のみ対応).
• mesh.NPX. . . x 方向のサイズは Lx = 2NPX .
• mesh.NPY. . . y 方向のサイズは Ly = 2NPY .
• mesh.NPZ. . . z 方向のサイズは Lz = 2NPZ .
time increment 以下で時間刻みを決める.
2
2• time increment. . . auto を選べば時間刻みの上限である Tdump = min(ρ/ηkmax, 1/kB T Γα kmax
) を決め
る. ここで kmax は格子幅 DX で決まる最大波数である. manual を選べば手動で値を入れる.

• time increment.auto.factor. . . 入力する値 factor を使って時間刻みは ∆t = Tdump × factor と決まる.
• time increment.manual.delta t. . . 手動で ∆t の値を入れる.
switch 以下でシミュレーションの各種条件を決める.
• switch.ROTATION. . . 粒子の回転運動の運動方程式を解く場合は ON を選ぶ.
• switch.HYDRO int. . . 流体について Navier–Stokes 方程式を解く場合は Correct を選ぶ. ただし Electrolyteの場合 Correct 以外の switch は実装されていない.

• switch.Stokes. . . Navier–Stokes 方程式を解く際, 慣性項を考慮すれば with advection を選ぶ. 慣性項を考慮しない Stokes 近似をする場合は w/o advection を選ぶ. ただし Electrolyte の場合 w/o advectionは実装されていない.

• swithch.LJ truncate. . . 粒子間に Lennard–Jones ポテンシャルによる力について, 引力項を含む通常の形
OFF, 引力項を含まない斥力のみ ON, まったく力を加えない NONE から選ぶ.
• switch.INIT distribution. . . 粒子の初期配置を, uniform random(ランダム), random walk(正方格子上からランダムにずれている), FCC(FCC 格子上), BCC(BCC 格子上), user specify(座標と速度はユー
ザによって指定される) から選ぶ. user specify を選んだ場合, 初期の粒子位置と速度をそれぞれ
user specify.Particles[].R と user specify.Particles[].v に数値を入れる. 入力するリストの数が Particle number で指定した数より小さければ, Gourmet 上で Edit->Add an array Element として
user specify.Particles[] を増やすか直接 UDF ファイルを編集する.

24

第 3 章 荷電コロイド粒子の電気泳動• switch.FIX CELL. . . 全速度の直流成分を 0, つまり v k=0 = 0 とする. 重力加速度を導入して沈降速度を測る場合これを ON にする.
output 以下にデータ出力の情報を入れる.
• 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. . . AVS データファイルのフォーマットを Binary か ASCII か選ぶ.
• output.UDF. . . 出力 UDF をする場合は ON を選ぶ.
出力 UDF の Particles[] 以下には各粒子の座標と速度を output.Num snap で指定した数だけレコードデータが保存される.

resume 以下で中断した計算を再計算する再開 UDF について指定する.
• resume.caclucation. . . NEW を選べば新しい計算を始める. 前回計算終了したときの計算ステップ数, 粒子
座標・速度・角速度, 溶媒速度場, イオン密度場の情報を読み込んで計算を再開させたい場合は CONTINUEを選ぶ. 再開 UDF ファイルの resume.CONTINUE.Saved Data 以下にデータはセーブされる. この場合,
前回終了時にできる再開 UDF を入力 UDF として計算を再開する. 使用方法は 4.6 節で詳しく解説する.

3.4.2

空間単位と時間単位

長さの単位としてを格子幅 ∆ を採用している.
時間の単位ついては以下の点に注意する.

• Navier–Stokes 方程式で ρ = η = ∆ = 1 となるような単位系を採用している.
• 仮に入力 udf ファイルで RHO= A, ETA= B, DX= C と入力した場合, 最大波数 kmax は C を用いて与2
えられ, 時間刻みの上限は運動量の拡散時間より Tdump = (A/B)/kmaxと与えられる. 時間刻み ∆t はTdump × factor で調整する. ただし時間の単位自体は ρ = η = ∆ = 1 の場合と同じである.
2
2• Electrolyte の場合、(A/B)/kmaxと (1/kB T Γα )/kmaxの小さな方を Tdump とする. この場合でも時間の単位は ρ = η = ∆ = 1 となるようにとる.

25

第 4 章 チュートリアル
4.1

実行方法

KAPSEL は,% kapsel -Iinput.udf -Ooutput.udf -Ddefine.udf -Rrestart.udfのように実行する. -I で入力 UDF, -O で出力 UDF, -D で定義 UDF, -D で再開 UDF を指定して実行する. いずれのオプションも省略することはできないことに注意する. 再開 UDF については 4.6 節で解説する.

4.2

計算例

4.2.1

シア下の懸濁液の流動

入力 UDF として shear.udf を用いればよい. 具体的には./data/と./data/avs/ が作成されていることを確認した上で以下を実行する.

% kapsel -Ishear.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
この例 (Fig. 4.1) ではメッシュサイズ 64 × 64 × 64 でシアレート S = 0.01, 粒子半径 a = 4,
粒子数 288 個で計算している. このとき粒子体積分率は ϕ = 29%.

界面厚さ ξ = 1,

Fig. 4.1: シア下の懸濁液の流動. スパン (z-) 方向中間面の速度場が矢印で表示されている．
constitutive eq として Shear Navier Stokes を選ぶと標準エラー出力に#1:dev_sigma_yx 2:shear_rate_obs 3:shear_rate ...
0.00647023 -0.00449194 0.01 ...
..
.
のようなデータが出力される．これは 1 列目がせん断応力，2 列目が実現されているシアレート，3 列目が入

26

第 4 章 チュートリアル.
力されたシアレート .. の時系列となっている．シミュレーションが終わると、最終行には上記データの過渡を除いた時間平均値が出力される．現在のところ，このデータは UDF に出力されない．

4.2.2

1 粒子電気泳動

1 粒子の電気泳動シミュレーションの計算例は, 入力 UDF として sample1.udf を用いればよい. 具体的には./data/と./data/avs/ が作成されていることを確認した上で以下を実行する.
% kapsel -Isample1.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
この例ではメッシュサイズ 64 × 64 × 64 で 1:1 電解質溶媒を考え, 半径 a = 5, 界面厚さ ξ = 2, 帯電量 Z = −100,
電場 Ex = 0.1, Debye 遮蔽長 κ−1 = 10 で計算している. Fig. 4.3 には泳動速度の時間発展を示した. 静電気力と流体抵抗力がつりあうところで定常速度を実現していることがわかる.
0.010

0.005

Z
Y

0.000

X

Fig. 4.2: +x 方向に印加された外部電場 E によって電気泳動するコロイド粒子 (青色) の様子. 濃淡は電荷密度分布を, 矢印は溶媒速度場をあらわす.

GraphSheet[]0"plot.dat" using 2:6

-0.02

-0.04

vx0

-0.06

-0.08

-0.1

-0.12

-0.14

0

2000

4000

6000

8000

time step

Fig. 4.3: コロイド粒子の泳動速度の時間発展.

10000

4.2. 計算例

4.2.3

27

多粒子電気泳動

粒子数を 32 個とした電気泳動シミュレーションは sample2.udf で計算できる. その他のパラメータは 1 粒子と同じである.
0.030

0.015

Z
Y
X

0.000

Fig. 4.4: +x 方向に印加された外部電場 E によって電気泳動するコロイド粒子 (青色) の様子.

4.2.4

正負に帯電した粒子混合系の電気泳動

正に帯電したコロイド粒子を 64 個, 負に帯電したコロイド粒子を 64 個の混合系の電気泳動は sample3.udf
で計算できる. その他のパラメータは 1 粒子と同じである.
0.050

0.000

Z
Y
X

-0.050Fig. 4.5: +x 方向に印加された外部電場 E によって電気泳動する正に帯電したコロイド粒子 (赤色) と負に帯電したコロイド粒子 (青色) の様子. 正符号コロイドは +x 方向に, 負符号コロイドは −x 方向に泳動している.

28

第 4 章 チュートリアル

4.3

AVS/Express による可視化

output.AVS を ON として AVS 形式のファイルを出力する. AVS/Express1 で avs charge.v ファイルを読み込み Read field モジュールにフィールドファイル data.fld を指定すると粒子, 溶媒速度場, 電荷分布をそれぞれ可視化することができる. 粒子は表面電荷 Ze の値をもとにカラー表示している. Fig. 4.2, Fig. 4.4, Fig. 4.5
は avs charge.v を用いて可視化している.

4.4

Gourmet による可視化

Gourmet を起動して output.udf を読み込み, Gourmet の Python パネルで KAPSEL に付属の Python スク
リプト show field.py を Load してから Run するとグラフィックウインドウが開き, 粒子, 溶媒速度場, 電荷密度分布の可視化が行われる2 . そして, グラフィックウインドウの最下段にある再生ボタンをクリックすればア
ニメーションがスタートする. 粒子は出力 UDF にレコードデータとして出力される粒子座標データをもとに可視化しているので, output.UDF を ON とする必要がある. さらに溶媒速度場と電荷密度分布は AVS 形式で出力される Binary ファイルを読み込んでいるので, output.AVS を ON としてさらに output.AVS.ON.FileType
を Binary にする必要がある.

4.5

gnuplot によるグラフ化

UDF に出力されたレコードデータから粒子の位置, 速度の時系列変化を gnuplot3 を用いてグラフ化することができる [1]. Gourmet を起動して output.udf を読み込み, Gourmet の Python パネルで KAPSEL に付属のPython スクリプト plot.py を Load してから Run すると, 粒子の位置, 速度を格納しているレコードデータからグラフ描画用シートを作成することができる. データをプロットするためには “View” ボックスで Tree についているチェックを Table に変更し, 変数の一覧に GraphSheet があるので選択して GraphSheet のデータをブラウズする. そして, Gourmet の Plot パネルで Make してから Plot すると Fig. 4.3 のようなグラフを描くことができる. 線種が多すぎて見づらいときは Plot パネルのエディター中に作成された plot コマンドから不要な部分を削除すればよい. Gourmet から gnuplot を呼び出す方法は文献 [1] の第 3 章に詳しく記述されている.

4.6

再開 UDF の使用方法

入力 UDF において resume.Calculation で NEW を選び% kapsel -Iinput.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
として計算を終了すると restart.udf に計算終了時のステップ数, 粒子, 溶媒速度場, イオン分布に関する情報が出力される. これらのデータを使って計算を再開させたい場合は以下の手順でおこなえばよい.

1. restart.udf を input.udf にコピーする.
2. Gourmet で input.udf を開き, resume.Calculation で CONTINUE を選ぶ.
3. さらに output.Num step を増やしトータルで計算させたいステップ数を指定して保存する.
4. 再び% kapsel -Iinput.udf -Ooutput.udf -Ddefine.udf -Rrestart.udfとすれば input.udf からデータを読み込み計算を再開する.
1 http://www.avs.com/
2 Gourmet 2003 以前のバージョンには不具合があり, このスクリプトが正常に動作しないことが確認されている. その場合は代りに付属の Python スクリプト particleshow.py を用いて可視化する. 溶媒速度場や電荷密度分布は教示されず, 粒子のみが表示される.
3 http://www.gnuplot.info/

4.7. 既知の障害

29

5. AVS データ, output.udf ともに前回計算したデータをそのままに, 再開した計算の出力データを追加することができる.
6. restart.udf には再び計算終了時におけるデータが格納されるので, この手順を繰り返せばさらに計算を続行することも可能である.

4.7

既知の障害

1. 実行するマシンによっては, constitutevi eq に Electrolyte を選んだ場合, プログラム終了時 (再開UDF のファイル操作時) にセグメンテーションエラーやアボートして終了することがある. しかし実行結果自体には問題はないようである。
2. show field.py は Gourmet 2003 では正常に動作しないようだ.
3. Windows の Cygwin+MinGW 環境でビルドされた KAPSEL は実行時にエラーで異常終了するようだ.
MinGW を使わなければ問題ない.

31

References[1] 土井正男, 滝本淳一. 物理仮想実験室. 名古屋大学出版会, 2004.
[2] W. B. Russel, D. A. Saville, and W. R. Schowalter. Colloidal Dispersions. Cambridge University Press,
Cambridge, England, 1989.
[3] Yasufumi Otsubo. The present status and prospect of suspension rheology. Nihon Reoroji Gakkaishi,
Vol. 31, pp. 15–22, 2003.
[4] Takayoshi Matsumoto. Rheology of colloidal disperse systems. Nihon Reoroji Gakkaishi, Vol. 32, pp.
3–9, 2004.
[5] Irvin M. Krieger. Rheology of monodisperse latices. Adv. Colloid Interface Sci., Vol. 3, pp. 111–136,
1972.
[6] Chun-Hway Hsueh and Paul F. Becher. Effective viscosity of suspensions of spheres. J. Am. Ceram.
Soc., Vol. 88, pp. 1046–1049, 2005.
[7] Yasuya Nakayama and Ryoichi Yamamoto. Simulation method to resolve hydrodynamic interactions incolloidal dispersions. Phys. Rev. E, Vol. 71, pp. 036707, 2005.
[8] Y. Nakayama, K. Kim, and R. Yamamoto. Simulating (electro)hydrodynamic effects in colloidal dispersions: smoothed profile method, 2006. cond-mat/0601322.
[9] J. Schumacher and B. Eckhardt. On statistically stationary homogeneous shear turbulence. Europhys.
Lett., Vol. 52, pp. 627–632, 2000.
[10] R. S. Rogallo. Numerical experiments in homogeneous turbulence. NASA Tech. Memo., pp. 81315,
1981.
[11] Akira Onuki. A new computer method of solving dynamic equations under externally applied deformations. J. Phys. Soc. Jpn., Vol. 66, pp. 1836–1837, 1997.
[12] K. Kim and R. Yamamoto. Efficient simulations of charged colloidal dispersions: A density functional
appraoch. Macromol. Theory Simul., vol. 14, pp. 278–284, 2005.
[13] K. Kim, Y. Nakayama, and R. Yamamoto. Direct Numerical Simulations of Electrophoresis of Charged
Colloids. Phys. Rev. Lett., Vol. 96, pp. 208302, 2006.
[14] H. Tanaka and T. Araki. Simulation method of colloidal suspensions with hydrodynamic interactions:
Fluid particle dynamics. Phys. Rev. Lett., Vol. 85, pp. 1338–1341, 2000.
[15] T. Kajishima, S. Takiguchi, H. Hamasaki, and Y. Miyake. Turbulence structure of particle-laden flow
in a vertical plane channel due to vortex shedding. JSME Int. J. Ser. B, Vol. 44, pp. 526–535, 2001.
[16] J.-L. Barrat and J.-P. Hansen. Basic Concepts for Simple and Complex Liquids. Cambridge University
Press, Cambridge, 2003.

32

第 4 章 チュートリアル[17] 北原文雄, 古澤邦夫, 尾崎正孝, 大島広行. ゼータ電位 微粒子界面の物理化学. サイエンティスト社, 1995.
[18] H. Ohshima, T. W. Healy, and L. R. White. Accurate analytic expressions for the surface charge
density/surface potential relationship and double-layer potential distribution for a spherical colloidal
particle. J. Colloid Interface Sci., Vol. 90, pp. 17–26, 1982.
[19] R. W. O’Brien and L. R. White. Electrophoretic mobility of a single colloidal particle. J. Chem. Soc.
Faraday Trans. 2, Vol. 74, pp. 1607–1626, 1978.
[20] H. Ohshima, T. W. Healy, and L. R. White. Approximation analytic expressions for the electrophoretic
mobility of spherical colloidal particles and the conductivity of their dilute suspensions. J. Chem. Soc.
Faraday Trans. 2, Vol. 79, pp. 1613–1628, 1983.