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

KAPSEL の解説
Kyoto Advanced Particle Simulator for ELectrohydrodynamics
version 2.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 章 懸濁液のダイナミクス

7
7
8
8
8

2.1

基礎方程式

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

9
9

2.2

2.1.1 溶媒流動 : Navier–Stokes 方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.2 コロイド粒子：ニュートンの運動方程式 . . . . . . . . . . . . . . . . . . . . . . . . . .
UDF 解説 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9
9
11

2.2.1
2.2.2

入力パラメータの指定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
空間単位と時間単位 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11
13

2.2.3

分散粒子のタイプについて . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

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

15

基礎方程式

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
コロイド粒子：ニュートンの運動方程式 . . . . . . . . . . . . . . . . . . . . . . . . . .
溶媒流動 : シア流下の Navier–Stokes 方程式 . . . . . . . . . . . . . . . . . . . . . . . .

15
16
16

3.1.3 レオロジ−測定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
UDF 解説 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1720

入力パラメータの指定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20
22
22

3.1

3.1.1
3.1.2
3.2

3.2.1
3.2.2
3.2.3

空間単位と時間単位 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
分散粒子のタイプについて . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

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

4.1

基本方程式

4.1.1
4.1.2
4.2

4.3

移流拡散方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23
23

Navier–Stokes 方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24
24

4.1.3 運動方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
電気二重層の性質 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Poisson–Boltzmann 方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24
25
25

4.2.2 Debye–Hückel 近似と Debye 遮蔽長 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
電気泳動の原理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25
26

Smoluchowski の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Hückel の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Henry の式および O’Brien–White による解析 . . . . . . . . . . . . . . . . . . . . . . .

26
27
27

UDF 解説 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4.1 入力パラメータの指定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2828

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

30

4.3.1
4.3.2
4.3.3
4.4

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

4.4.2

4
第 5 章 チュートリアル

5.1

31
実行方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

5.2

計算例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31
31
32

5.2.1
5.2.2

懸濁液の沈降 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.2.3
5.2.4

定常シア流下の懸濁液の流動. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
振動シア流下の懸濁液の動的粘弾性測定 . . . . . . . . . . . . . . . . . . . . . . . . . .

32
33

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

34
35
36

5.3
5.4

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

36
37

5.5
5.6
5.7

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

38
38
38

5.2.5
5.2.6
5.2.7

懸濁液の凝集・拡散 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

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

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 のインストール

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

1.1

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

LinuxLinux 版の場合,% cd /usr/local/OCTA2007/GOURMET 2007
% tar xvfpz gourmet∗ src.tar.gz% cd src% make% make installでビルドできる. ただし標準の Makeﬁle は

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/

第 1 章 KAPSEL のインストール

8
Windows

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

% ln -s /cygdrive/c/OCTA2007/ /usr/local/.
% cd /cygdrive/c/OCTA2007/GOURMET 2007
% ./GOURMET2007 SRC.EXE (or tar xvfpz gourmet∗ src.tar.gz)% cd 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 プログラムのコンパイルは付属の Makeﬁle でおこなう.
% cd kapsel% makeとすれば Linux でも Windows(Cygwin 環境) でも GNU C コンパイラ (gcc) でコンパイルされ kapsel という
バイナリファイルができる. もし Linux 上で Intel C++コンパイラでコンパイルする場合は Makeﬁle の編集をする. この Makeﬁle では% uname -n
でホスト名を取得するがそれが, たとえば icc11.1 がインストールされた新規マシンのホスト名が hoge であれば

ifeq ($(NODE),hoge)
CC = icc
CXX = icpc
CCOPT = -O3 -xT -axT
LINKS = -lm -lplatform -lcxaguard -lstdc++
endifの部分を書き加えればよい. その後% makeとすれば icc でコンパイルされた kapsel というバイナリファイルができる.

1.4

ユーザサポート

4 http://cygwin.com/

9

第 2 章 懸濁液のダイナミクス
微粒子が液体 (溶媒) 中に分散した系は微粒子懸濁液と呼ばれ，食品・塗料・顔料・化粧品・スラリーなど我々の日常生活に数多く存在している．このような微粒子懸濁液を理解する上で，微粒子および溶媒の両者の運動を同時に把握することは非常に重要である．KAPSEL では「重力下での微粒子懸濁液の沈降現象」，「微粒子の凝集現象」，および，「微粒子の拡散現象」をシミュレートすることが可能である．

2.1

基礎方程式

粘度 η, 密度 ρ をもつニュートン流体中に分散した半径 a である N 個の球状粒子の運動を考える．溶媒であるニュートン流体は，ナヴィエ・ストークス方程式に従う粗視化された流動場 v(r) として記述し，微粒子は，熱揺らぎによるランダム力をもつニュートンの運動方程式に従う粒子描像 {X i , V i , Ωi } を採用する．ここで，X i , V i , Ωi は，それぞれ i 番目の微粒子の位置，並進速度，回転速度である．溶媒と微粒子の連成問題を，固定直交格子上で高速かつ正確に解くために，溶媒と微粒子の界面を滑らかな関数 ϕ(r) で記述する Smoothed

Proﬁle Method (SPM) を導入する [7, 8, 9]．ここで ϕ ∈ [0, 1] は，粒子領域で ϕ = 1，溶媒領域で ϕ = 0 の値をもつオーダーパラメタを表している．

2.1.1

溶媒流動 : Navier–Stokes 方程式

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

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

(2.1)

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

2.1.2

コロイド粒子：ニュートンの運動方程式

質量 Mp かつ慣性モーメント I p をもつ i 番目のコロイド粒子の位置、速度および角速度 {Ri , V i , Ω} は運
動方程式:

Ṙi

= V i,other
FH+ GVi ,
i + FiMp V̇ i=I p · Ω̇iΩ
= NH
i + Gi ,

(2.2)
(2.3)
(2.4)

によって時間発展する. I p = (2/5)Mp a2 I はコロイド球状粒子の慣性モーメントで I は単位テンソルである．
H
ここで F Hi ，N i はコロイド粒子が流体から受ける力およびトルクであり，これらは固液間での運動量収支を表す．この計算の実装法は [7, 8, 9] を参照されたい. また F otherには流動による力以外による力であり，コ
i

第2章

10

懸濁液のダイナミクス

ロイド粒子間の直接相互作用や重力，および浮力などを表す．GVi , GΩ
i は，熱揺らぎにより粒子に働くランダ
ム力およびトルクであり，白色雑音として導入されている．ずなわち, ⟨GVi ⟩ = ⟨GiΩ ⟩ = 0，かつ

⟨GVi (t) · GVj (0)⟩Ω⟨GΩ
i (t) · Gj (0)⟩= 3kB T αV δ(t)δij ,Ω= 3kB T α δ(t)δij ,

(2.5)
(2.6)

である．ここで ⟨⟩ は統計平均をあらわしており，αV および αΩ は，粒子温度 T を決定するパラメタである．
粒子温度 T は溶媒中で揺らいでいる１粒子の並進拡散係数 DV および回転拡散係数 DΩ から決定される．
粒子温度 T の決定手順

1. ある値に αV と αΩ をそれぞれ固定して１粒子系の熱平衡計算を行う．
VΩ
2. シミュレーションから１粒子の並進拡散係数 Dsim
や回転拡散係数 Dsim
を平均二乗変位等を通してそれぞれ求める.
VΩ
3. 得られた拡散係数 Dsim
, Dsimを無限希釈系の１粒子系の拡散係数の解析解 (並進運動に対しては，D0V =kB T V /6πηa および，回転運動に対しては，D0Ω = kB T Ω /8πηa3 ) と比較することにより，それぞれの並進および回転運動に対する温度 T V , T Ω を決定することができる．
4. T Ω ̸= T V の場合，T Ω = T V (= T ) を満足するように αV または αΩ を調整する.
このアプローチによる温度決定方法の詳細については文献 [10, 11, 12] を参照されたい．

2.2. UDF 解説

2.2

UDF 解説2.2.1

入力パラメータの指定

11

constitutive eq として Navier Stokes を選ぶと Newton 流体中のコロイド粒子の運動のシミュレーショ
ンが行える．以下, 入力 UDF ファイルで指定する変数について解説する.

constitutive eq.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.kBT. . . 粒子温度.
• constitutive eq.Shear Navier Stokes.alpha v. . . 粒子の並進温度に対する補正項1 .
• constitutive eq.Shear Navier Stokes.alpha o. . . 粒子の回転温度に対する補正項.
object type.type には，球状粒子または鎖のどちらかをを選択することができる.
⃝ object type.spherical particle.Particle spec[] 以下には球状粒子のタイプの情報を入れる.
• object type.spherical particle.Particle spec[] . . . 粒子数.
• object type.spherical particle.Particle spec[].MASS RATIO. . . 粒子の流体に対する密度比.
• object type.spherical particle.Particle spec[].Surface charge RATIO. . . 表面電荷 (使用しない).
⃝ object type.chain.Chain spec[] では鎖の情報を入れる．
• object type.chain.Chain spec[].Beads number . . . 一本の鎖に属するビーズの数• object type.chain.Chain spec[].Chain number . . . 鎖の本数
• object type.chain.Chain spec[].MASS RATIO. . . ビーズの流体に対する密度比.
• object type.chain.Chain spec[].Surface charge. . . 表面電荷 (使用しない).
粒子半径, 界面厚さは以下に情報を入れる.

• A XI. . . 界面の厚さ ξ.
• A. . . 粒子半径.
gravity 以下で重力に関する情報を入れる.
• gravity.G. . . 重力加速度2 .
• gravity.direction. . . 重力を加える方向を選ぶ.
粒子間力について Lennard–Jones ポテンシャルのタイプを決める.
1 微粒子の熱拡散を導入する場合は，温度を入力する．熱平衡での粒子の拡散運動から温度は見積もられるが，界面の厚さ等の数値的な理由のために正しい温度を返さない可能性がある．この場合，もし測定された温度が 1.2 倍と低い温度が見積もられた場合，このファクターを 1.2 と設定することで既定の温度を実現することができる．回転に対しても同様である．半径 5 または 4，かつ ξ = 2 に対して
は，alpha v=alpha o=1 で既定の温度を計算できる．
2 微粒子懸濁液の沈降現象をシミュレートする場合は，gravity.G に値を入力する必要がある．

第2章

12

懸濁液のダイナミクス• 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 を選ぶ.
• switch.HYDRO int. . . 流体について Navier–Stokes 方程式を解く場合は Correct を選ぶ.
• switch.Stokes. . . Navier–Stokes 方程式を解く際, 慣性項を考慮すれば with advection を選ぶ. 慣性項を考慮しない Stokes 近似をする場合は w/o advection を選ぶ. ただし w/o advection は実装されていない.
• swithch.LJ truncate. . . 粒子間に Lennard–Jones ポテンシャルによる力について, 引力項を含む通常の形
OFF3 , 引力項を含まない斥力のみ 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 の ﬁeld ファイ
ルは./data/に data.fld というファイル名で出力される. データファイルは./data/avs/に data *.datというファイル名で出力される. *にはステップ数が入る.
3 凝集現象をシミュレートしたい場合は，swithch.LJ truncate を OFF にセットする．

2.2. UDF 解説

13• 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 として計算を再開する. 使用方法は 5.6 節で詳しく解説する.

2.2.2

空間単位と時間単位

長さの単位としてを格子幅 ∆ を採用している. 時間の単位 τ0 ついて流体の密度 ρ と粘性率 η と格子幅で決まる τ0 = ρ∆2 /η を採用している．

• 仮に入力 udf ファイルで RHO= A, ETA= B, DX= C と入力した場合, 最大波数 kmax は C を用いて与2と与えられる. 時間刻み ∆t はえられ, 時間刻みの上限は運動量の拡散時間より Tdump = (A/B)/kmax
Tdump × factor で調整する.
• 現実との対応を考察する．考察したい空間スケールとして，グリッドサイズを 1µm の長さを考える．また溶媒として水 (η = 10−3 Pa · s, ρ = 103 kg/m2 ) を対象とした場合，時間の単位は τ0 = 10−6 s となる．

2.2.3

分散粒子のタイプについて

分散粒子のタイプとして球状粒子と鎖のどちらかを選択することができる．球状粒子の場合，打ち切りレナー
ドジョーンズポテンシャルを用いることで粒子間反発のみの相互作用を導入できる．このポテンシャルは， [( )2n ( )n )]4ϵ σ− rσij+ϵ (rij < 21/n σ),
rij
(2.7)
ULJ (rij ) =0(rij > 21/n σ)であり，σ, ϵ, rij は，粒子直径，相互作用の強度，i 番目と j 番目の粒子間の距離をあらわしている．n = 6, 12, 18のいずれかを選択することができ，n が大きいほうが鋭いポテンシャルになる．これは，swithch.LJ truncate
を ON にすることで可能である．また 粒子間に引力を導入したい場合は，swithch.LJ truncate を OFF にすることで，通常のレナードジョーンズポテンシャルの計算が可能となる．
また鎖を選択した場合は，鎖内の隣接ビーズ間は FENE (ﬁnitely extensible non-linear elastic) ポテンシャル

1
UF (r) = − kc R02 ln{1 − (r/R0 )2 },2で結合される． ここで kc = 30ϵ/σ 2 , R0 = 1.5σ, r は隣接ビーズ間の距離をあらわしている

(2.8)

15

第 3 章 懸濁液のレオロジー
微粒子が液体（溶媒）中に分散している微粒子懸濁液は，化学工学や機械工学，物理学など工学・理学分野の広範囲にわたって重要であり，活発に研究が行われている．このような系は，食品・塗料・顔料・化粧品・
スラリーなどに見出すことができ，その用途および産業的な製造・加工プロセスにおいてその流動挙動の把握が重要である．しかしながら，微粒子懸濁液の動的挙動は微粒子間相互作用や熱揺らぎなどに影響され大変複雑な様相を示し，その流動挙動の把握は容易ではない．微粒子懸濁液は，異なる状況のもとで多種多様な挙動を示す．例えば，微粒子懸濁液の輸送現象においては粘度のような微粒子懸濁液の流動特性が，分散した粒子の濃度や流動の強さ，粒子界面の性質などに強く依存する．このような巨視的な流動特性と微視的なメカニズ
ムの関係を明らかにすることは，微粒子懸濁液の基礎研究における主要なテーマの一つとなっている．
上述のように懸濁液のレオロジーは応用上重要であるが，それについて理論的に一般論を構築するにはあまりに複雑である (例えばテキスト [2] や比較的最近のレビュー (日本語) [3, 4] を参照.)．粒子濃度が希薄な極限においては懸濁液のみかけの粘度 η app は

η app
η

5
= 1 + φ,2(3.1)

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

3.1

基礎方程式

粘度 η, 密度 ρ をもつニュートン流体中に分散した半径 a である N 個の球状粒子の運動を考える．溶媒であるニュートン流体は，ナヴィエ・ストークス方程式に従う粗視化された流動場 v(r) として記述し，微粒子は，熱揺らぎによるランダム力をもつニュートンの運動方程式に従う粒子描像 {X i , V i , Ωi } を採用する．ここで，X i , V i , Ωi は，それぞれ i 番目の微粒子の位置，並進速度，回転速度である．溶媒と微粒子の連成問題を，固定直交格子上で高速かつ正確に解くために，溶媒と微粒子の界面を滑らかな関数 ϕ(r) で記述する Smoothed

Proﬁle Method (SPM) を導入する [7, 8, 9]．ここで ϕ ∈ [0, 1] は，粒子領域で ϕ = 1，溶媒領域で ϕ = 0 の値をもつオーダーパラメタを表している．

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

16

3.1.1

コロイド粒子：ニュートンの運動方程式

質量 Mp かつ慣性モーメント I p をもつ i 番目のコロイド粒子の位置、速度および角速度 {Ri , V i , Ω} は運
動方程式:

Ṙi

=

V i,

(3.2)

Mp V̇ i

=

other
FH
+ GVi ,
i + Fi

(3.3)

=

Ω
NH
i + Gi ,(3.4)

I p · Ω̇i

によって時間発展する. I p = (2/5)Mp a2 I はコロイド球状粒子の慣性モーメントで I は単位テンソルである．
H
ここで F Hi ，N i はコロイド粒子が流体から受ける力およびトルクであり，これらは固液間での運動量収支を表す．この計算の実装法は [7, 8, 9] を参照されたい. また F otherには流動による力以外による力であり，コi
ロイド粒子間の直接相互作用や浮力などを表す．GVi , GΩ
i は，熱揺らぎにより粒子に働くランダム力およびトルクであり，白色雑音として導入されている．ずなわち, ⟨GVi ⟩ = ⟨GΩi ⟩ = 0，かつ

⟨GVi (t) · GVj (0)⟩Ω⟨GΩ
i (t) · Gj (0)⟩= 3kB T αV δ(t)δij ,Ω= 3kB T α δ(t)δij ,

(3.5)
(3.6)

である．ここで ⟨⟩ は統計平均をあらわしており，αV および αΩ は，粒子温度 T を決定するパラメタである．
粒子温度 T は溶媒中で揺らいでいる１粒子の並進拡散係数 DV および回転拡散係数 DΩ から決定される．
粒子温度 T の決定手順

1. ある値に αV と αΩ をそれぞれ固定して１粒子系の熱平衡計算を行う．
VΩ
2. シミュレーションから１粒子の並進拡散係数 Dsim
や回転拡散係数 Dsim
を平均二乗変位等を通してそれぞれ求める.
Ω
Vを無限希釈系の１粒子系の拡散係数の解析解 (並進運動に対しては，D0V =, Dsim3. 得られた拡散係数 Dsim

kB T V /6πηa および，回転運動に対しては，D0Ω = kB T Ω /8πηa3 ) と比較することにより，それぞれの並進および回転運動に対する温度 T V , T Ω を決定することができる．
4. T Ω ̸= T V の場合，T Ω = T V (= T ) を満足するように αV または αΩ を調整する.
このアプローチによる温度決定方法の詳細については文献 [10, 11, 12] を参照されたい．

3.1.2

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

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

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

(3.7)

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

3.1. 基礎方程式

17

単純せん断流の駆動および境界条件

外力 f shear は速度場 vx を図 3.1 に示すような下記のジグザグの速度プ

ロファイル，γ̇(−y − Ly /2), (−Ly /2 < y ≤ −Ly /4)vx (y) = γ̇y,(−Ly /4 < y ≤ Ly /4)γ̇(−y + L /2) (L /4 < y ≤ L /2)
y
y
y

(3.8)

に保つように制御する外力である．ここで γ̇ はずり速度であり，Ly は速度勾配方向の系の長さである．その流動の強さは，ずり速度 γ̇ = ∇y vx で特徴づけられる．このプロファイルの導入は，流体計算部分を高速フーリエ変換に基づいた計算流体力学 (CFD) の分野でよく使われているスペクトル法を使用しているためであり，周期境界条件下でのずり流動の流体計算が容易に可能となる 1 ．この手法の有効性は，[12] の中で詳細に議論されている．KAPSEL では, 計算の毎ステップにおいてコロイド粒子および流体の移流が済んだ後に f shear が既定の速度場になるようにフィードバックとして印加される．以上がシミュレーションに必要な方程式である.

yzxFig. 3.1: 速度プロファイル． 矢印は，流体の速度ベクトルを示している．

3.1.3

レオロジ−測定

微粒子懸濁液の流動特性を測定するために，分散系の見かけの応力 σ app の定式化を行う．形式的に，微粒子懸濁液の運動量保存則は，以下のように記述される．

d(ρt v) = ∇ · σ dis + ρt f shear ,dtρt = (1 − ϕ)ρ + ϕρp(3.9)(3.10)

ここで σ dis は，圧力や粘性，慣性による応力への寄与を含んでいる．分散系の局所応力 s は，次のように定義される．

s = σ dis − ρt vv(3.11)

ここで，ρt vv は移流による運動量輸送を表している．この s は，直接計算することはできないが，この局所応力 s を用いてみかけの応力 σ app を求めることができる．みかけの応力は，局所応力の体積積分としてあらわされ，
1 単純せん断流を解く方法は，境界壁の速度／力を境界条件で与える方法や、分子動力学シミュレーションにおける Lees–Edwards 境界条件に相当する Rogallo の動格子法 [15, 16] などがあるが, ここでは格子が歪まない本方法を用いた．また，動的粘弾性測定の場合，
KAPSEL ではストレスの伝搬を避けるために，系全体への体積力として拘束力 f shear を導入している．一方，定常流の場合は，デルタ関数的な拘束力として与えている．

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

18

∫
1
dxs
V
∫
[
]
1
dx [∇ · (sx)]T − x∇ · s
=
V
∫
1
=
dx[−x∇ · s]
V
∫
[ (
)]
1
∂
=
dx x ρt f shear − (ρt v)
V
∂t
∫
∫
1
1
∂
shear
=
dxxρt f−dxx (ρt v),VV∂t

σ app =

(3.12)(3.13)(3.14)(3.15)(3.16)と書ける．ここで，恒等式 s = [∇ · (sx)]T − x∇ · s を用い，かつ表面項を落とした．V は系の体積を表している．

定常せん断流の場合
定常せん断流の場合のみかけの粘度の定式化を行う．定常状態では，⟨∂⟩(ρt v) = 0∂tt(3.17)

が成り立つ．ここで ⟨⟩t は時間平均を表している．したがって，時間平均されたみかけの応力 (3.16) は，∫⟩
1⟨⟨σ app ⟩t =(3.18)
dxxρt f shearVtと書ける．この評価式には位置ベクトルが陽に入っているが、測定したいマクロなせん断応力は，座標系の原点のとり方に依らないものである．これは f shear によって全系の並進運動が駆動されないことと対応する．以下、この事実を確認する．座標原点をずらして x = x′ + x0 とすると∫∫∫
dx (x′ + x0 ) ρt f shear =dxx′ ρt f shear + x0 dxρt f shear .

(3.19)∫
外力 f shear はせん断流れを起こし，全系の並進を起こさないように制御されているので， dxρt f shear = 0 である．SPM では式 (3.18) を用いて定常せん断流に対する巨視的な応力を評価している．
懸濁液のみかけの粘度は，上記のみかけの応力を用いてη app =

app
⟨σxy
⟩
γ̇

(3.20)

app
によって定義される．γ̇(ずり速度) は制御パラメタであり，マクロなせん断応力 σxyは計算データから評価される．定常せん断流における濃厚微粒子分散系に関する詳細は，文献 [13] を参照されたい．

振動せん断流の場合
振動せん断流場を系に導入することによって，懸濁液の動的粘弾性を測定することが可能である．KAPSEL
では，系に時間的に変化するせん断流を形成するために，制御パラメタとして，ずり速度 γ̇ を次のように与えている．

γ̇(t) = γ̇0 cos(ωt).

(3.21)

ここで，γ̇0 はずり速度の振幅である．また、ひずみの最大振幅 γ0 は，∫ π/2ω
γ0 =γ̇(s)ds =0γ̇0ω(3.22)

3.1. 基礎方程式

19

によって評価される．ひずみが小さい場合は線形粘弾性の領域に相当し，一方大きい場合は非線形粘弾性領域に相当する．
一般に振動外場に対するみかけの応力の応答は，
app
σxy= σ0 cos(ωt − δ)

(3.23)

appと書ける．ここで，σ0 は応力の振幅であり，δ はずり速度 γ̇ とずり応力 σxyとの間の位相差をあわらしている．この時，動的粘弾性関数である貯蔵弾性率 G′ および，損失弾性率 G′′ は，G′ (ω) =σ0 sin δ
σ0 cos δω, G′′ (ω) =
ω.
γ̇0
γ̇0

(3.24)と書ける．G′ は懸濁液の弾性的な振る舞いを特徴づけ，G′′ は粘性的な振る舞いを特徴づける．KAPSEL においては，制御パラメタはずり速度 γ̇(t) であり，出力が式 (3.16) のみかけのずり応力となっている．それぞれの変数の時間発展をモニターすることにより動的粘弾性関数をそれぞれ求めることが可能である．これらの動的粘弾性関数の導出については，レオロジーに関する教科書を参照して頂きたい．

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

20

3.2

UDF 解説3.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.kBT. . . 粒子温度.
• constitutive eq.Shear Navier Stokes.alpha v. . . 粒子の並進温度に対する補正項2 .
• constitutive eq.Shear Navier Stokes.alpha o. . . 粒子の回転温度に対する補正項.
• constitutive eq.Shear Navier Stokes.External ﬁeld.type. . . DC (定常シア流れ）または AC (振動シア流れ) を選択する．

• constitutive eq.Shear Navier Stokes.External ﬁeld.DC.Shear rate. . . ずり速度 γ̇.
• constitutive eq.Shear Navier Stokes.External ﬁeld.AC.Shear rate. . . 振動ずり速度の振幅 γ̇0• constitutive eq.Shear Navier Stokes.External ﬁeld.AC.Frequency. . . ずり速度の周波数 ω
object type.type には，球状粒子または鎖のどちらかをを選択することができる.
⃝ object type.spherical particle.Particle spec[] 以下には球状粒子のタイプの情報を入れる.
• object type.spherical particle.Particle spec[] . . . 粒子数.
• object type.spherical particle.Particle spec[].MASS RATIO. . . 粒子の流体に対する密度比.
• object type.spherical particle.Particle spec[].Surface charge RATIO. . . 表面電荷 (使用しない).
⃝ object type.chain.Chain spec[] では鎖の情報を入れる．
• object type.chain.Chain spec[].Beads number . . . 一本の鎖に属するビーズの数• object type.chain.Chain spec[].Chain number . . . 鎖の本数
• object type.chain.Chain spec[].MASS RATIO. . . ビーズの流体に対する密度比.
• object type.chain.Chain spec[].Surface charge. . . 表面電荷 (使用しない).
粒子半径, 界面厚さは以下に情報を入れる.

• A XI. . . 界面の厚さ ξ. Shear Navier Stokes を選んだ場合, 少なくとも ξ ≥ 1 を選ぶ必要がある.
• A. . . 粒子半径.
2 熱平衡での粒子の拡散運動から温度は見積もられるが，界面の厚さ等の数値的な理由のために正しい温度を返さない可能性がある．この場合，もし測定された温度が 1.2 倍と低い温度が見積もられた場合，このファクターを 1.2 と設定することで既定の温度を実現することができる．回転に対しても同様である．半径 5 または 4，かつ ξ = 2 に対しては，alpha v=alpha o=1 で既定の温度を計算できる．

3.2. UDF 解説

21

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を決める. ここで kmax は格• time increment. . . auto を選べば時間刻みの上限である Tdump = ρ/η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 を選ぶ. ただし Shear Navier Stokesの場合 Correct 以外の switch は実装されていない.
• switch.Stokes. . . Navier–Stokes 方程式を解く際, 慣性項を考慮すれば with advection を選ぶ. 慣性項を考慮しない Stokes 近似をする場合は w/o advection を選ぶ. ただし 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 ファイルを編集する.
• switch.FIX CELL. . . 全速度の直流成分を 0, つまり v k=0 = 0 とする. 重力加速度を導入して沈降速度を測る場合これを ON にする.
output 以下にデータ出力の情報を入れる.
• output.GTS. . . データ出力のインターバルのステップ数.
• output.Num snap. . . データ出力の回数. つまり全ステップ数は GTS×Num snap で決まる.

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

22• output.AVS. . . AVS 形式のデータ出力をする場合は ON を選ぶ.

• output.AVS.ON.Out dir. . . AVS 形式のデータ出力をするディレクトリを入れる. たとえば data と指定す
れば, 事前に./data および./data/avs/のディレクトリを作成しておく必要がある. AVS の ﬁeld ファイ
ルは./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 として計算を再開する. 使用方法は 5.6 節で詳しく解説する.

3.2.2

空間単位と時間単位

長さの単位としてを格子幅 ∆ を採用している. 時間の単位 τ0 ついて流体の密度 ρ と粘性率 η と格子幅で決まる τ0 = ρ∆2 /η を採用している．

• 仮に入力 udf ファイルで RHO= A, ETA= B, DX= C と入力した場合, 最大波数 kmax は C を用いて与2と与えられる. 時間刻み ∆t はえられ, 時間刻みの上限は運動量の拡散時間より Tdump = (A/B)/kmax
Tdump × factor で調整する.
• 現実との対応を考察する．考察したい空間スケールとして，グリッドサイズを 1µm の長さを考える．また溶媒として水 (η = 10−3 Pa · s, ρ = 103 kg/m2 ) を対象とした場合，時間の単位は τ0 = 10−6 s となる．

3.2.3

分散粒子のタイプについて

分散粒子のタイプとして球状粒子と鎖のどちらかを選択することができる．球状粒子の場合，打ち切りレナー
ドジョーンズポテンシャルを用いることで粒子間反発のみの相互作用を導入できる．このポテンシャルは， [( )2n ( )n )]4ϵ σ− rσij+ϵ (rij < 21/n σ),
rij(3.25)ULJ (rij ) =0(rij > 21/n σ)であり，σ, ϵ, rij は，粒子直径，相互作用の強度，i 番目と j 番目の粒子間の距離をあらわしている．n = 6, 12, 18のいずれかを選択することができ，n が大きいほうが鋭いポテンシャルになる．これは，swithch.LJ truncate
を ON にすることで可能である．また 粒子間に引力を導入したい場合は，swithch.LJ truncate を OFF にすることで，通常のレナードジョーンズポテンシャルの計算が可能となる．
また鎖を選択した場合は，鎖内の隣接ビーズ間は FENE (ﬁnitely extensible non-linear elastic) ポテンシャル

1
UF (r) = − kc R02 ln{1 − (r/R0 )2 },2で結合される． ここで kc = 30ϵ/σ 2 , R0 = 1.5σ, r は隣接ビーズ間の距離をあらわしている(3.26)

23

第 4 章 荷電コロイド粒子の電気泳動
水のような誘電率が非常に大きい溶媒にコロイド粒子が分散しているとき, コロイド表面にある解離基からイオンが放出されて粒子表面は電荷を帯びる. 放出されたイオンは粒子表面に静電的に引き寄せられるが, 同時に熱ゆらぎによって拡散され, コロイド粒子周りに電気二重層とよばれる雲のようなイオン雰囲気を形成する. 平衡状態のコロイド分散系の性質は Poisson–Boltzmann 方程式によって記述され, これを線型化した

Debye–Hückel 近似を用いて見通しの良い理解を得ることができる. これに対して電気泳動をはじめとした界面動電現象とよばれる現象では, 粒子とイオン分布の挙動は流体力学相互作用と静電相互作用の競合によって決定される. 両者の競合の結果としてイオン分布は粒子の運動に追随できず, 電気二重層は球対称から歪み平衡状態とは異なる状況が起こる. このような複雑な状況を解析するには, 理論的には簡単な近似を導入する他なく, さらに計算機シミュレーションによっても正しく再現されたことはこれまでほとんどなかった．
本プロジェクトではこの状況を打破すべく Smoothed Proﬁle (SP) 法という新しいシミュレーション方法の開発を行ってきた [7, 17, 8, 18]. SP 法ではコロイドは粒子描像のまま扱う.

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

4.1

基本方程式

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

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

(4.1)

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

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

第4章

24

4.1.1

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

移流拡散方程式

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

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

∫

(4.3)

α

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

(4.5)

とする [21]. E は外部電場をあらわし, 静電ポテンシャル Ψ(r) は Poisson 方程式:

ϵ∇2 Ψ = −ρe(4.6)

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

4.1.2

Navier–Stokes 方程式

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

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

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

4.1.3

運動方程式

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

Ṙi
Mp V̇ i

= V i,

(4.8)

other
= FH
i + Fi(4.9)

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

4.2. 電気二重層の性質

4.2

25

電気二重層の性質

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

4.2.1

Poisson–Boltzmann 方程式

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

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

Boltzmann 方程式とよばれる.

4.2.2

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

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

(4.13)となる. このとき、長さの次元をもつ定数1
κ−1 = √,
8πλB z 2 C̄

(4.14)

は Debye 遮蔽長とよばれる. ここで λB = e2 /4πkB T ϵ はビヨルン (Bjerrum) 長である. 一般の電解質の場合は,1
κ−1 = √,
∑ 2
4πλB α Zα C̄α(4.15)

である. 今, 系は球対称なので動径方向 r = |r| についてのみ考えればよく,

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

(4.17)

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

第4章

26

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

κ−1 は小さくなる. たとえば z : z 対称電解質溶媒中をとして, バルク塩濃度 C̄ を考えるとする. 媒質を 25 ◦ Cの水とするとビヨルン長は λB = 0.72nm となり, 実際の数値を入れると式 (4.14) は0.3
κ−1 = √ (nm)z C̄(4.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 ),(4.19)となる. すなわち表面電荷の増加に対して、表面電位が線型に増加する関係を与える. 表面電位が大きいとこ
ろでは Debye–Hückel 近似は適用できない. 実際は表面電荷の増加に対して表面電位はしだいに増加しにくくなる. このときの関係式は Poisson–Boltzmann 方程式を解くことによって得られるが, 1:1 電解質溶媒に対する近似式が Loeb–Overbeek–Wiersema [2] あるいは Ohshima–Healy–White [22, 23] によって提案されている.

4.3

電気泳動の原理

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

ZeE = 6πηaV(4.20)

V
Ze
=
E
6πηa

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

4.3.1

Smoluchowski の式

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

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

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

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

η[vx (y) + V ] − ϵEΨ(y) = 0(4.24)

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

V
ϵζ
=
E
η

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

4.3.2

Hückel の式

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

Ze
4πϵa

(4.26)

V
2 ϵζ
=
E
3 η

(4.27)

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

4.3.3

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

式 (4.25) と式 (4.27) を結び付けるものとして一般の κa に対して Henry の式

Vϵζ= f (κa)Eη(4.28)

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

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

f (κa) =(4.29)(4.30)

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

0.95

0.9

f(κa)

0.85

0.8

0.75

0.7

0.1

1

10
κa

100

1000

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

第4章

28

4.4

UDF 解説4.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 proﬁle. . . イオン分布の初期設定について 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 coeﬀ 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 coeﬀ positive ion. . . 正イオンの Onsager 輸送係数.
• constitutive eq.Electrolyte.Add salt.salt.Onsager coeﬀ negative ion. . . 負イオンの Onsager 輸送係数.
• constitutive eq.Electrolyte.Add salt.salt.Debye length. . . Debye 遮蔽長を指定すると対応する塩濃度が設定される.
• constitutive eq.Electrolyte.Electric ﬁeld. . . 外部電場の ON, OFF を選択する.
• constitutive eq.Electrolyte.Electric ﬁeld.ON. . . DC(直流電場) か AC(交流電場) を選択する.
• constitutive eq.Electrolyte.Electric ﬁeld.ON.DC.Ex. . . x 方向電場の強さ. 交流電場の場合も同じ.
• constitutive eq.Electrolyte.Electric ﬁeld.ON.DC.Ey. . . y 方向電場の強さ. 交流電場の場合も同じ.
• constitutive eq.Electrolyte.Electric ﬁeld.ON.DC.Ez. . . z 方向電場の強さ. 交流電場の場合も同じ.
• constitutive eq.Electrolyte.Electric ﬁeld.ON.AC.Frequency. . . 交流電場の周波数
object type.type には，球状粒子または鎖を選択することができる．
⃝ object type.spherical particle.Particle spec[] 以下には球状粒子のタイプの情報を入れる.
1 constitutive eq の他の switch として, Newton 流体中のコロイド粒子の運動は Navier Stokes を, シア流下における Newton 流
体中のコロイド粒子の運動は Shear Navier Stokes を選べばシミュレートすることができる. 荷電粒子の熱揺らぎは考慮されていない

4.4. UDF 解説

29• object type.spherical particle.Particle spec[] . . . 粒子数.
• object type.spherical particle.Particle spec[].MASS RATIO. . . 粒子の流体に対する密度比（比重）.
• object type.spherical particle.Particle spec[].Surface charge RATIO. . . 粒子表面電荷の帯電量.
⃝ object type.chain.Chain spec[] では鎖のタイプの情報を入れる．
• object type.chain.Chain spec[].Beads number . . . 一本の鎖に属するビーズの数• object type.chain.Chain spec[].Chain number . . . 鎖の本数
• object type.chain.Chain spec[].MASS RATIO. . . ビーズの流体に対する密度比.
• object type.chain.Chain 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 は実装されていない.

第4章

30

荷電コロイド粒子の電気泳動• 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 ファイルを編集する.
• 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 の ﬁeld ファイ
ルは./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 として計算を再開する. 使用方法は 5.6 節で詳しく解説する.

4.4.2

空間単位と時間単位

長さの単位としてを格子幅 ∆ を採用している. 時間の単位 τ0 ついて流体の密度 ρ と粘性率 η と格子幅で決まる τ0 = ρ∆2 /η を採用している．

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

31

第 5 章 チュートリアル

5.1

実行方法

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

5.2

計算例

5.2.1

懸濁液の沈降入力 UDF として gravity.udf を用いればよい. 具体的には./avs g1/と./avs g1/avs/ が作成されていることを確認した上で以下を実行する1 .

% kapsel -Igravity.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
この例 (Fig. 5.1) ではメッシュサイズ 256 × 256 × 512 で計算している.

このとき粒子体積分率は φ = 0.05 であり，粒子数 Np = 3204, 粒子直径 10, 界面厚さ ξ = 2 である．図 5.1 には，重力下で沈降している懸濁液のスナップショットを示している．

Fig. 5.1: 重力下での懸濁液の沈降．重力は，z 方向下向きに働いている．流体の色は，z 方向の流体場の速度を示している．

1 AVS データが必要でない場合は，AVS のスイッチを OFF に切り替え，./avs g1/と./avs g1/avs/を作成する必要はない．これはその他のサンプルファイルに対しても同様である．

第5章

32

5.2.2

チュートリアル

懸濁液の凝集・拡散

入力 UDF として aggregate.udf および repulsive.udf を用いればよい. aggregate.udf では，粒子間に引力相互作用が働いた系の凝集・拡散現象を計算することができる．一方，repulsive.udf では，粒子間に反発力のみの相互作用をもつ系の拡散現象を計算することができる．

% kapsel -Iaggregate.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf% kapsel -Irepulsive.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
この例 (Fig. 5.2) ではメッシュサイズ 128 × 128 × 128 で計算している.

このとき粒子体積分率は φ = 0.063

であり，粒子数 Np = 500, 粒子直径 8, 界面厚さ ξ = 2，温度 kB T = 5 である．溶媒は，密度 ρ = 1，粘性係数

η = 1 である．図 5.2(左) には，粒子間に引力相互作用をもった懸濁液の凝集・拡散の様子を示している．また，図 5.2(右) には，反発相互作用をもった懸濁液の拡散のスナップショットを示している．

Fig. 5.2: 左図は，粒子間相互作用に引力相互作用をもつ懸濁液の凝集・拡散の様子のスナップショットを示している．右図は，反発相互作用のみの粒子間相互作用をもつ懸濁液のスナップショットを示している．

5.2.3

定常シア流下の懸濁液の流動

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

% kapsel -Ishear dc.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
この例 (Fig. 5.3) ではメッシュサイズ 128 × 128 × 128 でずり速度 γ̇ = 0.01, 粒子半径 a = 4,
粒子数 500 個で計算している.

界面厚さ ξ = 2,

このとき粒子体積分率は φ = 0.064%.

constitutive eq として Shear Navier Stokes を選び，Externa field として DC を選択すると標準エラー
出力に
2 AVS データが必要でない場合は，AVS のスイッチを OFF に切り替え，./shear dc/と./shear dc/avs/を作成する必要はない．これはその他のサンプルファイルに対しても同様である．

5.2. 計算例

33
0.32

0.11

-0.11

Y
X Z
-0.32

Fig. 5.3: シア流下の懸濁液の流動. スパン (z-) 方向中間面の流れ (x-) 方向の速度場が色で表示されている．また，速度勾配 (y-) 方向下面の速度が矢印で表示されている．
#1:time 2:shear_rate_temporal 3:shear_strain_temporal 4:shear_stress_temporal 5:viscosity
0.0748173 -0.00588609 -0.000440381

...

..
.
のようなデータが出力される．この時系列出力の意味は以下のとおり。
1 列目 経過時間
2 列目 実現されているずり速度 γ̇ sim (t).
3 列目 ひずみ (γ̇ × t)
app
4 列目 計算される懸濁液のみかけのずり応力 σxy(t)app
5 列目 懸濁液の粘度 η app = σxy(t)/γ̇ sim (t)現在のところ，このデータは UDF に出力されない．Fig. 5.4 は，いくつかの体積分率に対する粘度の時間発展を示している．分散粒子の濃度が高いほど，定常シア流粘度が高くなっていることがわかる．

5.2.4

振動シア流下の懸濁液の動的粘弾性測定

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

% kapsel -Ishear ac.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
constitutive eq として Shear Navier Stokes を選び，Externa field として AC を選択すると標準エラー
出力に

第5章

34

チュートリアル10
Volume fraction=0.383
Volume fraction=0.191
Volume fraction=0.064

987Viscosity

6
5
4
3
2
1
0
0

100

200

300

400

500

600

700

time

Fig. 5.4: 定常シア流下の懸濁液の粘度の時間発展. システムサイズ Lx = Ly = Lz = 128. 体積分率 φ は上か
1 列目と 5 列目のプロットを表示

ら 0.383(N = 3000), 0.191(N = 1500), 0.064(N = 500) のデータである.
している．

#1:time 2:shear_rate_temporal 3:shear_strain_temporal 4:shear_stress_temporal 5:
0.0748173 -0.00988654 -0.000739684 -68.4334 68.5048 0.0714706...
..
.
のようなデータが出力される．この時系列出力の意味は以下のとおり。
1 列目 経過時間
2 列目 実現されているずり速度 γ̇ sim (t).
3 列目 ひずみ (γ̇ × t)
4 列目 外力 f shear から計算される応力成分
5 列目 懸濁液の加速からくる応力成分
6 列目 懸濁液のみかけのずり応力
現在のところ，このデータは UDF に出力されない．Fig. 5.5 ではずり速度 (1 列目) とみかけのずり応力 (6 列目) の時間発展を示している．メッシュサイズ 128 × 128 × 128 でずり速度の振幅 γ̇ = 0.01,

振動数 ω = 0.05,

粒子半径 a = 4, 界面厚さ ξ = 2, 温度 kB T = 0.01, 粒子数 500 個で計算している.

このとき粒子体積分率は
φ = 0.064. これらのデータからずり速度と応力の位相差 δ と振幅 σ0 , γ̇0 を見積もり，動的粘弾性関数 G′ および G′′ を計算することができる.

5.2.5

1 粒子電気泳動

1 粒子の電気泳動シミュレーションの計算例は, 入力 UDF として colloid 1.udf を用いればよい. 具体的には./avs ch/と./avs ch/avs/ が作成されていることを確認した上で以下を実行する.
% kapsel -Icolloid 1.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
この例ではメッシュサイズ 64 × 64 × 64 で 1:1 電解質溶媒 (粘性係数 η = 1, 密度 ρ = 1) を考え, 半径 a = 5, 界
面厚さ ξ = 2, 帯電量 Z = −100, 電場 Ex = 0.1, Debye 遮蔽長 κ−1 = 10 で計算している. Fig. 5.6 には 1 粒子

5.2. 計算例

350.04
Shear stress
Shear rate
0.03

Shear rate and Shear stress

0.02

0.01

0

-0.01

-0.02

-0.03

-0.04
0

50

100

150

200

250

300

time

Fig. 5.5: 振動シア流下の懸濁液のずり速度と応力の時間発展. システムサイズ Lx = Ly = Lz = 128.
電気泳動のスナップショットを示している．Fig. 5.7 には泳動速度の時間発展を示した. 静電気力と流体抵抗力
がつりあうところで定常速度を実現していることがわかる. Fig. 5.7 のプロットには，plot.py を用いている.
0.011

0.007

0.004

Y

Z

0.000
X

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

5.2.6

多粒子電気泳動

粒子数を 32 個とした電気泳動シミュレーションは colloid 32.udf で計算できる. 初期配置は，ランダムであり，その他のパラメータは 1 粒子と同じである. 具体的には./avs ch/と./avs ch/avs/ が作成されていることを確認した上で以下を実行する.

% kapsel -Icolloid 32.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
Fig. 5.8 には多粒子系での電気泳動のスナップショットを示している．Fig. 5.9 には，平均泳動速度の時間発展を示している．この図では，粒子数 NP = 16, 32, 64 の場合のデータがプロットされている．粒子数が増加するにつれて粒子の平均泳動速度が減少することがわかる．

第5章

36

チュートリアル0"plot.dat" u 2:6

-0.01

-0.02

Vx

-0.03

-0.04

-0.05

-0.06

-0.07
0

1000

2000

3000

4000

5000

time step

Fig. 5.7: １個のコロイド粒子の泳動速度の時間発展.
0.04

0.03

0.01

Y

Z
X

0.00

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

5.2.7

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

正に帯電したコロイド粒子を 32 個, 負に帯電したコロイド粒子を 32 個の混合系の電気泳動は colloid p32m32.udf
で計算できる. その他のパラメータは 1 粒子と同じである. 具体的には./avs ch/と./avs ch/avs/ が作成されていることを確認した上で以下を実行する.

% kapsel -Icolloid p32m32.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
Fig. 5.10 には，正負に帯電したコロイド粒子の電気泳動のスナップショットを表している．

5.3

AVS/Express による可視化

output.AVS を ON として AVS 形式のファイルを出力する. AVS/Express3 で avs charge.v ファイルを読み込み Read field モジュールにフィールドファイル data.fld を指定すると粒子, 溶媒速度場, 電荷分布をそれ
3 http://www.avs.com/

5.4. Gourmet による可視化

37
0
Np=16
Np=32
Np=64

-0.005-0.01-0.015Vx-0.02-0.025-0.03-0.035-0.04-0.045
-0.05
0

5000

10000

15000
time

20000

25000

30000

Fig. 5.9: 多粒子系での平均泳動速度の時間発展. 上から粒子数 Np = 64, 32, 16 のデータである.
0.01

0.00

-0.00

Z
Y
X

-0.01

Fig. 5.10: +x 方向に印加された外部電場 E によって電気泳動する正に帯電したコロイド粒子 (赤色) と負に帯電したコロイド粒子 (青色) の様子. 正符号コロイドは +x 方向に, 負符号コロイドは −x 方向に泳動している.
ぞれ可視化することができる. Fig. 5.6, Fig. 5.8, Fig. 5.10 は avs charge.v を用いて可視化している.

5.4

Gourmet による可視化

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

4 Gourmet 2003 以前のバージョンには不具合があり, このスクリプトが正常に動作しないことが確認されている. その場合は代りに付属の Python スクリプト particleshow.py を用いて可視化する. 溶媒速度場や電荷密度分布は表示されず, 粒子のみが表示される.

第5章

38

5.5

チュートリアルgnuplot によるグラフ化

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

5.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 からデータを読み込み計算を再開する.

5. AVS データ, output.udf ともに前回計算したデータをそのままに, 再開した計算の出力データを追加することができる.

6. restart.udf には再び計算終了時におけるデータが格納されるので, この手順を繰り返せばさらに計算を続行することも可能である.

5.7

既知の障害

1. 実行するマシンによっては, constitutevi eq に Electrolyte を選んだ場合, プログラム終了時 (再開UDF のファイル操作時) にセグメンテーションエラーやアボートして終了することがある. しかし実行結果自体には問題はないようである。

2. show field.py は Gourmet 2003 では正常に動作しないようだ.
3. Windows の Cygwin+MinGW 環境でビルドされた KAPSEL は実行時にエラーで異常終了するようだ.
MinGW を使わなければ問題ない.

5 http://www.gnuplot.info/

39

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. Eﬀective 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 eﬀects in colloidal dispersions: smoothed proﬁle method. Eur. Phys. J. E., Vol. 26, pp. 361, (2008).
[9] X. Lio, M. R. Maxey, and G. E. Karniadakis. Smoothed proﬁle method for particulate ﬂows: Error
analysis and simulations. J. Comput. Phys., Vol. 228, pp. 1750, (2009).
[10] T. Iwashita, Y. Nakayama, and R. Yamamoto. A numerical model for Brownian particles ﬂuctuating
in incompressible ﬂuids. J. Phys. Soc. Jpn., Vol. 77, pp. 074007, (2008).
[11] T. Iwashita, Y. Nakayama, and R. Yamamoto. Velocity autocorrelation function of ﬂuctuating particles
in incompressible ﬂuids. -Toward direct numerical simulation of particle dispersions- Prog. Theor. Phys.
Suppl, Vol. 178, pp. 86, (2009).
[12] Takuya Iwashita and Ryoichi Yamamoto. Short-time motion of Brownian particles in a shear ﬂow Phys.
Rev. E, Vol. 79, pp. 031401, (2009).
[13] Takuya Iwashita and Ryoichi Yamamoto. Direct numerical simulations for non-Newtonian rheology ofconcentrated particle dispersions Phys. Rev. E, 2009. cond-mat/0905.0130.
[14] J. Schumacher and B. Eckhardt. On statistically stationary homogeneous shear turbulence. Europhys.
Lett., Vol. 52, pp. 627–632, 2000.
[15] R. S. Rogallo. Numerical experiments in homogeneous turbulence. NASA Tech. Memo., pp. 81315,
1981.
[16] Akira Onuki. A new computer method of solving dynamic equations under externally applied deformations. J. Phys. Soc. Jpn., Vol. 66, pp. 1836–1837, 1997.

第5章

40

チュートリアル[17] K. Kim and R. Yamamoto. Eﬃcient simulations of charged colloidal dispersions: A density functional
appraoch. Macromol. Theory Simul., vol. 14, pp. 278–284, 2005.
[18] K. Kim, Y. Nakayama, and R. Yamamoto. Direct Numerical Simulations of Electrophoresis of Charged
Colloids. Phys. Rev. Lett., Vol. 96, pp. 208302, 2006.
[19] 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.
[20] T. Kajishima, S. Takiguchi, H. Hamasaki, and Y. Miyake. Turbulence structure of particle-laden ﬂow
in a vertical plane channel due to vortex shedding. JSME Int. J. Ser. B, Vol. 44, pp. 526–535, 2001.
[21] J.-L. Barrat and J.-P. Hansen. Basic Concepts for Simple and Complex Liquids. Cambridge University
Press, Cambridge, 2003.
[22] 北原文雄, 古澤邦夫, 尾崎正孝, 大島広行. ゼータ電位 微粒子界面の物理化学. サイエンティスト社, 1995.
[23] 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.
[24] 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.
[25] 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.