A DNS simulator for the dynamics of particles in fluids
KAPSEL

Version 5.12
User’s Manual

KAPSEL Development Teamkapsel.dev@gmail.com
February 1, 2024

Contents
1

Introduction
1.1 What is KAPSEL? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 What are KAPSEL’s capabilities? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.1 Simulating particles dispersed in a Newtonian fluid . . . . . . . . . . . . . . . . . . . . . . .
1.2.2 Simulating a disperse particle system in the presence of a shear flow . . . . . . . . . . . . . .
1.2.3 Simulating charged colloidal particles dispersed in electrolytic solutions . . . . . . . . . . . .
1.2.4 Simulating particles dispersed in two-component phase-separated fluids . . . . . . . . . . . .
1.2.5 Simulating microswimmers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.6 Simulating Quincke rollers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 What is covered by this manual? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4
4
4
5
5
5
5
5
5
5

2

Installing and running KAPSEL
2.1 System-dependent procedures for installing and running KAPSEL . . . . . . . . . . . . . . . . . . .
2.1.1 The KAPSEL runtime environment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.2 Installing OCTA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.3 Installing KAPSEL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.4 Validating your KAPSEL license . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.5 Testing KAPSEL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.6 Methods for visualizing simulation data . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.7 Analyzing simulation data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Input UDF and definition UDF files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Fluid settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 Object (particle) settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.3 Common simulation settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.4 Configuring selectable features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.5 Data output settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.6 Restart settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.7 GOURMET display settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Output and restart UDF files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7
7
7
7
7
8
9
10
11
13
13
13
17
18
21
23
23
24

3

Simulating particles dispersed in Newtonian fluids
3.1 Theoretical background and basic equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.1 Basic equations for disperse particle systems . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.2 Outline of the SP method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.3 Smoothed Profile (SP) functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.4 Procedure for determining particle temperature . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Input UDF settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Fluid settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 Configuring objects (particles) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.3 Choice of length and time units . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Computational examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Sedimentation of particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.2 Particle diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.3 Particle aggregates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.4 Motion of a particle chain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.5 Motion of an arbitrarily-shaped rigid particle . . . . . . . . . . . . . . . . . . . . . . . . . .

25
25
25
27
27
29
30
30
30
30
31
31
31
32
32
33

1

CONTENTS

3.3.6

2

Sedimentation of arbitrarily-shaped rigid particles . . . . . . . . . . . . . . . . . . . . . . . .

34

4Simulating disperse particle systems in shear flows
4.1 Theoretical background and basic equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.1 Steady-state shear flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.2 Oscillatory shear flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Input UDF files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Fluid settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.2 Object (particle) settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Computational examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1 Rheology of particle suspensions under steady-state shear flow . . . . . . . . . . . . . . . . .
4.3.2 Rheology of particle suspensions under oscillatory shear flow . . . . . . . . . . . . . . . . .

37
37
39
39
40
40
40
41
41
42

5

Simulating charged particles dispersed in electrolytic solutions
5.1 Electrophoresis of charged colloidal particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Basic equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.1 The advection-diffusion equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.2 The Navier-Stokes equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.3 Equations of motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3 Properties of electric double layers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.1 The Poisson-Boltzmann equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.2 The Debye-Hückel approximation and the Debye screening length . . . . . . . . . . . . . . .
5.4 Principles of electrophoresis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4.1 Smoluchowski’s equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4.2 Hückel’s equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4.3 Henry’s equation and the O’Brien-White analysis . . . . . . . . . . . . . . . . . . . . . . . .
5.5 Input UDF files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.5.1 Fluid settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.5.2 Object (particle) settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.5.3 Choice of length and time units . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.6 Computational examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.6.1 Single-particle electrophoresis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.6.2 Many-particle electrophoresis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.6.3 Electrophoresis in a mixed particle system including both positively and negatively charged
particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.6.4 Using AVS/Express for visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.6.5 Using Gourmet for visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.6.6 Using gnuplot to plot data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45
45
45
46
46
46
47
47
47
48
48
49
49
50
50
51
51
53
53
54
55
55
55
56

Simulating particles dispersed in two-component phase-separated fluids
6.1 Adding particles to two-component phase-separated fluids . . . . . . . . . . . . . . . . . . . . . . .
6.2 Theoretical background and basic equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2.1 Basic equations for two-component phase-separated fluids . . . . . . . . . . . . . . . . . . .
6.2.2 Viscosity of suspensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3 Input UDF files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3.1 Fluid settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Settings for simulations with no imposed flow field . . . . . . . . . . . . . . . . . . . . . . .
Settings for simulations under shear flow . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3.2 Object (particle) settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3.3 Choice of length and time units . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.4 Computational examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.4.1 Motion of particles dispersed in a two-component phase-separated fluid . . . . . . . . . . . .
Running the simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Visualizing simulation data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Writing Python scripts to analyze simulation data . . . . . . . . . . . . . . . . . . . . . . . .
6.4.2 Motion of particles dispersed in a two-component phase-separated liquid under shear flow . .
6.4.3 Pickering emulsions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

57
57
57
57
59
59
59
60
63
66
66
67
67
67
68
68
70
71

6

CONTENTS

3

7

Simulating microswimmers
7.1 Theoretical background and basic equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.1.1 The squirmers model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.1.2 Basic equations for microswimmers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2 Input UDF files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2.1 Fluid settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2.2 Configuring objects (particles) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2.3 Configuring objects (planar walls) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.3 Computational examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.3.1 Motion of microswimmers with periodic boundary conditions . . . . . . . . . . . . . . . . .
7.3.2 Motion of microswimmers confined between two parallel slabs . . . . . . . . . . . . . . . . .

74
74
74
75
76
76
76
76
77
77
77

8

Simulating Quincke rollers
8.1 Theoretical background and basic equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.1.1 Quincke rollers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.1.2 Basic equations for Quincke rollers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.2 Input UDF files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.2.1 Fluid settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.2.2 Object (particle) settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.3 Computational examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.3.1 Motion of a single Quincke roller . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.3.2 Motion of multiple Quincke rollers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

79
79
79
80
81
81
81
82
82
82

A Definition of the input UDF file format
84
A.1 Fluid settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
A.2 Object settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
A.3 Common simulation settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
A.4 Configuring selectable features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
A.5 Data output settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
A.6 Definition of UDF output data classes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
A.7 Restart settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
A.8 GOURMET display settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
B Particle calculations

101

C Fluid calculations by spectral methods

102

D Fluid calculations by finite-difference methods (FDM)
104
D.1 Solving the Navier-Stokes equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
D.1.1 Explicit MAC solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
D.1.2 Implicit MAC solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
D.1.3 Solving the Navier-Stokes equations in the presence of shear flow with Lees-Edwards boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
D.2 Solving the Cahn-Hilliard equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
D.2.1 Explicit solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
D.2.2 Implicit solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
D.2.3 Solving the Cahn-Hilliard equation in the presence of shear flow with Lees-Edwards boundary
conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
E Interfacing with the Lis library107

F Compiling KAPSEL
108
F.1 The KAPSEL developer’s environment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
F.2 Installing OCTA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
F.3 Building libplatform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
F.4 Installing FFTW . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
F.5 Installing HDF5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
F.6 Installing LIS (optional) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
F.7 Building the KAPSEL executable file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

Chapter 1

Introduction
1.1

What is KAPSEL?

KAPSEL1 is a software package for simulating disperse particle systems (systems in which particles move freely

through a fluid). In simulating the motion of particles dispersed in a fluid, we must account not only for the impact
of direct interactions between particles, but also for indirect interactions—known as hydrodynamic interactions—mediated by the motion of the fluid medium. Thus, reproducing the complicated trajectories of fluid-dispersed solid
particles requires understanding how the fluid moves between them. One common computational strategy is direct numerical simulation (DNS), in which the Navier-Stokes equations are solved numerically to yield an accurate description of fluid motion on length scales smaller than the particle sizes, after which the hydrodynamic interactions between
particles at a given time step may be determined from local fluid stresses. However, DNS calculations typically entail
high computational costs, and this has spurred the development of alternative strategies that avoid numerically solving
the Navier-Stokes equations. Examples include the methods of dissipative particle dynamics (DPD) and multi-particle
collision dynamics (MPC), which model fluids as systems of particles satisfying the law of momentum conservation,and the Brownian dynamics (BD) and Stokesian dynamics (SD) approaches, in which hydrodynamic interactions between particles are determined without explicitly solving for fluid motion. In comparison to these and other simplified
approaches, DNS has the great advantage of accurately describing fluid motion based on the fundamental governing
equations; for this reason, the method is not restricted to modeling simple Newtonian fluids, but may be extended tohandle more complex fluids.
In DNS, particles are represented not by infinitesimal points but by finite-volume bodies, and momentum exchanges between particles and fluids through the solid-fluid interface are properly accounted for. The task of modelingboth the trajectories of arbitrarily-shaped solid bodies moving freely through a fluid and the flow of the fluid aroundthose bodies is a coupled multiphysics problem; the most common strategy for solving this problem via numericalcomputation is the finite-element method (FEM), in which finite bodies are represented by discretized non-uniform
meshes that are re-generated at each time step. Consequently, FEM models of disperse systems with many solid
particles involve enormous numbers of mesh nodes that must be reconfigured at each computational step, resulting
in exorbitant computational costs. To circumvent this difficulty while retaining the strengths of DNS modeling for
multi-particle disperse systems, we independently developed a smoothed-profile (SP) method, in which sharp, discontinuous interfaces between particles and fluids are replaced by finite-thickness transition regions characterized byprofile functions. The SP method, which allows simulations to achieve both high accuracy and high computational
efficiency, is implemented in KAPSEL (Kyoto Advanced Particle Simulator for ELectrohydrodynamics), a simulator
we developed for DNS modeling of disperse particle systems. This manual discusses the basic principles of KAPSEL,
explains how to install the package, and presents a number of detailed sample computations. A detailed discussion of
the SP method and its background may be found in Ref. [1].

1.2What are KAPSEL’s capabilities?

KAPSEL offers a wide range of analytical capabilities, which users may invoke and configure by preparing UDF input
1 KAPSEL is a software package developed by members of the Soft Matter Engineering Laboratory in the Department of Chemical Engineering
at Kyoto University . Early versions of KAPSEL, up to and including version 4, were distributed as open-source freeware, with users granted
the freedom to define their own specifications subject to the provisions of a license agreement. Since version 5, KAPSEL has been closed-sourcecommercial software, distributed only in the form of executable files available for purchase; viewing or editing the source code requires a separate
non-disclosure agreement (NDA). The lead developer responsible for the current version and all previous versions of the KAPSEL core is Professor
Ryoichi Yamamoto, and the copyright holders are Kyoto University and Professor Ryoichi Yamamoto.

4

CHAPTER 1. INTRODUCTION

5

files (typically given file names like input.udf). This is discussed in detail in later sections of this manual; here webriefly outline KAPSEL’s primary simulation capabilities.

1.2.1

Simulating particles dispersed in a Newtonian fluid

Particles dispersed in a Newtonian fluid may be simulated by selecting Navier Stokes as the constitutive eq.
This is useful for purposes such as predicting the trajectories of arbitrarily-shaped particles sinking due to gravity andcharacterizing the thermal-diffusion and agglomeration behavior of dispersed particles due to Brownian motion. See
Section 3 for further details. Set constitutive eq to Navier Stokes FDM to perform the analysis using a fluid
solver based on the finite-difference method (FDM) instead of the default pseudo-spectral solver.

1.2.2

Simulating a disperse particle system in the presence of a shear flow

Particles dispersed in a Newtonian fluid under shear flow may be simulated by selecting Shear Navier Stokes Lees Edwards
as the constitutive eq. This is useful for predicting the trajectories of arbitrarily-shaped particles under shear
flows or characterizing the rheological properties of disperse particle systems. See Section 4 for further details. Set
constitutive eq to Shear Navier Stokes Lees Edwards FDM to perform the analysis using an FDM-based fluid
solver.

1.2.3

Simulating charged colloidal particles dispersed in electrolytic solutions

Charged colloidal particles dispersed in an electrolytic solution may be simulated by selecting Electrolyte as theconstitutive eq. This is useful for predicting stable structures of charged colloidal systems and characterizing theelectrophoretic properties of charged colloidal particles in steady-state (DC) or oscillating (AC) electric fields. See
Section 5 for further details.

1.2.4

Simulating particles dispersed in two-component phase-separated fluids

Particles dispersed in two-component phase-separated fluids may be simulated by selecting Navier Stokes Cahn Hilliard FDM
as the constitutive eq. This is useful for purposes such as predicting fluid phase-separation structures or studying
how particle distributions are affected by changes in fluid components or in factors such as the affinity between particles and fluid phase interfaces. To simulate particles dispersed in two-component phase-separated fluids under shear
flow, set constitutive eq to Shear NS LE CH FDM. This is useful for predicting phase-separation structures in thepresence of shear flows or characterizing rheological properties of systems of particles dispersed in two-component
phase-separated fluids. See Section 6 for further details.

1.2.5

Simulating microswimmers

Microswimmer particles dispersed in Newtonian liquids may be simulated by selecting Navier Stokes as the constitutive eq.
This is useful for purposes such as predicting collective behavior resulting from the interaction among large numbersof microswimmers. See Section 7 for further details. Set constitutive eq to Navier Stokes FDM to performthe analysis using an FDM-based fluid solver. Set constitutive eq to Navier Stokes Cahn Hilliard FDM tosimulate microswimmer particles dispersed in two-component phase-separated fluids.

1.2.6

Simulating Quincke rollers

Quincke rollers—self-propelled particles produced by rolling on a slab immersed in a Newtonian fluid—may besimulated by selecting Navier Stokes as the constitutive eq. This is useful for purposes such as predicting
collective behavior resulting from the interaction among large numbers of Quincke rollers. See Section 8 for further
details. Set constitutive eq to Navier Stokes FDM to perform the analysis using an FDM-based fluid solver.

1.3

What is covered by this manual?

As a first step, Section 2 provides detailed instructions for installing and building the software and libraries needed torun KAPSEL.
The remaining six sections describe various types of simulations supported by KAPSEL. Each of these sections
begins with a thorough review of the theoretical background and basic equations relevant to the given type of simulation, then provides a concise explanation of how to create UDF input files to set simulation parameters; finally, eachCHAPTER 1. INTRODUCTION

6

section concludes with a series of detailed computational examples. The structure and format of the UDF input file isdocumented comprehensively in Appendix A.

Chapter 2

Installing and running KAPSEL
2.1

System-dependent procedures for installing and running KAPSEL

2.1.1

The KAPSEL runtime environment

The procedure for installing and running KAPSEL depends slightly on whether you are running on a Windows, Linux,
or Mac system:Linux

If you are using KAPSEL on a Linux system, you may skip to section 2.1.2 below.
Windows

To use KAPSEL on a Windows system, you must first install the Windows Subsystem for Linux (WSL)1 . OpenPowerShell or Windows Command Prompt in administrator mode by right-clicking and selecting ”Run as administrator”, enter the following command, then restart your machine. After installing WSL, skip to section 2.1.2 below.
$ wsl --installMac

To use KAPSEL on MacOS, you must first install Xcode and the command line tools.

2.1.2

Installing OCTA

For tasks such as processing input parameters and visualizing output data, KAPSEL relies on an external user-interface
module known as Gourmet. Gourmet is distributed as an internal component of the open-source package OCTA (auniversal simulator for soft materials), which must be installed before using KAPSEL.
To install OCTA, simply visit http://octa.jp/ to download and run an appropriate installer for your system.2
The instructions below assume that OCTA8.# has been installed in /usr/local/OCTA8# (for Linux or Mac) orin C:\OCTA8.# (for Windows).3

2.1.3

Installing KAPSELFirst download the latest version of the KAPSEL package; this will be an archive file named kapsel#.#.zip.4 Then
extract the content of the archive and enter the top-level directory:$ unzip kapsel#.#.zip$ cd kapsel#.#
The top-level directory contains the following subdirectories:./bin/ KAPSEL executable files for various platforms.
1 https://learn.microsoft.com/en-us/windows/wsl
2 Answers to questions regarding OCTA/GOURMET are available upon registering as a member of the OCTA-BBS website.
3 Here the pound symbol (“#”) is to be replaced with numerals to indicate the version of your OCTA installation.
4 Here the pound symbol (“#”) is to be replaced with numerals to indicate the version of your KAPSEL installation.

7

CHAPTER 2. INSTALLING AND RUNNING KAPSEL

8./Documents/ User manuals (including this document).
./Examples/ Example calculations illustrating various KAPSEL features.
./UDF/ UDF files (define.udf/input.udf) for using all features in the latest KAPSEL release.
Next, within the ./bin/ directory, find the executable file appropriate for your runtime environment and create asymbolic link in the KAPSEL installation directory.
Linux$ ln -s ./bin/linux/* .
Windows
From Ubuntu console window (not Windows command line)$ ln -s ./bin/linux/* .
Mac (Intel or Arm CPU)$ ln -s ./bin/macOS/* .

2.1.4

Validating your KAPSEL license

Commercial KAPSEL binaries will require a valid license to run. Two type of licenses are available, node-locked andcluster licenses, both of which support offline / online environments. Node-locked licenses, issued on an individual
basis, will require a unique (anonymized) machine id before they can be activated. Likewise, cluster licenses require a
unique environment id, together with the list of allowed user groups. These identifiers can be obtained by running the
key command and checking the MACHINE FINGERPRINT code. The first time you run you should see output similar
to the following,→,→$ ./key[INFO] System name : ...
[INFO] Machine: ...
[INFO] Node name: ...
[INFO] Release: ...
[INFO] Version: ...
[INFO] Network interface (en0) : ...
[INFO] Network interface (en1) : ...
...
[INFO] Host domain : ...
[INFO] Machine UUID: ...
# [Machine Fingerprint]
8bzm9b3xae00898539e754a3adba83bb11cdcea441d06401f5f3df350fcc7705# [Environ Fingerprint]
22cca0e1fc66357d0aa8d238e7gg675a6d7ebcf323fdbbcef7e09c46c6259f17# [User Groups] user_group_1 user_group_2 ...
[ERROR] Environment variable KAPSEL_PUBLIC_KEY is missing

For node locked licenses, you will be asked to provide the MACHINE FINGERPRINT for all machines that will run
KAPSELto the license provider. For cluster licenses (Linux only), you will be asked to provide the ENVIRON FINGERPRINT,
together with the list of user groups you wish to enable on the cluster (e.g., user group 1, user group 2). To determine
the appropriate ENVIRON FINGERPRINT, be sure to run the key program on a compute node, and not on the login
node (i.e., you should submit the job as you would a KAPSEL job). You will receive a KAPSEL license key, together
with an encrypted machine/license file from the license provider. The license key uniquely identifies your license, thelicense files are linked to individual machines/clusters and allow us to support offline environments. Note that both the
key and the license are required to run KAPSEL . A license keys is like a personal password for KAPSEL, it should betreated as such, and should not be shared.
To allow KAPSEL to perform the license verification, you should set the appropriate environment variables. Thiscan be done by editting and sourcing the setvars.sh script found in the ./bin directory$ source ./bin/setvars.shCHAPTER 2. INSTALLING AND RUNNING KAPSEL

9This will set the KAPSEL ACCOUNT ID, KAPSEL PUBLIC KEY, KAPSEL LICENSE KEY, and KAPSEL LICENSE PATH
variables. Note that you should edit the script to use the KAPSEL license key you were issued as the KAPSEL LICENSE KEY
variable, and set the KAPSEL LICENSE PATH to point to your license file. The format for the license key will be of theform key/your kapsel license key. The license files are distributed as plain text files, with an encoded payloadcontained between the BEGIN/END MACHINE FILE header/footer-----BEGIN MACHINE FILE----eyJlbmMiOiJsSTc4N0QwcGZua1RvRDVOSjFpRXlaU093Q09QQ0NOdktKZHpC...
TlJWWGp6Ym5DRkF6V3lOU3NUeG9xZm9MV2FlWlhITEZnR21Ub2VBdz09Iiwi
YWxnIjoiYWVzLTI1Ni1nY20rZWQyNTUxOSJ9-----END MACHINE FILE----You should never edit the KAPSEL ACCOUNT ID or KAPSEL PUBLIC KEY variables, nor the license files, as doing so
will make it impossible to validate your license.
Alternatively you can execute the following directly in the console,→$ export KAPSEL_ACCOUNT_ID="85e8531e-915a-4795-a287-2bec0d90ae5e"$ export
KAPSEL_PUBLIC_KEY="5fae6bb532c12ef70e1ce5f69b33ecd4fea8a522658b1204a6203ee5f101f608"$ export KAPSEL_LICENSE_KEY="key/your_kapsel_licence_key"$ export KAPSEL_LICENSE_PATH="./license.lic"

With a valid/activated license, and properly defined environment variables, you can now run KAPSEL. You should see
the following message upon successful validation of your license and machine.

,→,→$ ./key# [KAPSEL LICENSE VERIFICATION]# [Machine Fingerprint]
8bzm9b3xae00898539e754a3adba83bb11cdcea441d06401f5f3df350fcc7705# [Environ Fingerprint]
22cca0e1fc66357d0aa8d238e7gg675a6d7ebcf323fdbbcef7e09c46c6259f17# [User Groups] user_group_1 user_group_2 ...
[INFO] Try as Node license[INFO] Importing license file
[INFO] Searching for file './license.lic' : YES[OK] License file successfully imported!
[INFO] Verifying...
[OK] License file successfully verified!
[INFO] Decrypting...
[OK] License file successfully decrypted!
[INFO] Parsing...
[OK] License successfully parsed!
# [INFO] Node license imported# [INFO] Attempting to connect to license server...
...
# [OK] Signature is valid!
# [OK] License key is valid! (code = "VALID")# [License ID] g9aa321f-a19f-424b-99a6-a1x9f393d69g
#
OK!

2.1.5

Testing KAPSEL

If executing KAPSEL produces the following command-line output, then KAPSEL is properly installed on your system.
$ cd UDF$ ../kapsel -Iinput.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf## OMP RUNTIME :# Number of processors: 0# Number of threads: 8# Max OMP threads: 8# Dynamic thread enabled?
: 0# Nested parallelism enabled? : 1##using input.udf as input#using output.udf as output

CHAPTER 2. INSTALLING AND RUNNING KAPSEL

10#using define.udf as definition#using restart.udf as restart# [KAPSEL LICENSE VERIFICATION]# [Machine Fingerprint]
8bzm9b3xae00898539e754a3adba83bb11cdcea441d06401f5f3df350fcc7705# [Environ Fingerprint]
22cca0e1fc66357d0aa8d238e7gg675a6d7ebcf323fdbbcef7e09c46c6259f17# [User Groups] user_group_1 user_group_2 ...
# [INFO] Node license imported# [INFO] Attempting to connect to license server...
# [INFO] Online validation OK# [OK] Signature is valid!
# [OK] License key is valid! (code = "VALID")# [License ID] g9aa321f-a19f-424b-99a6-a1x9f393d69g#...
#output.udf end.
#restart.udf end.
#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

,→

,→

The input.udf file used in the above sample describes a simulation involving 5 heavy particles and 5 light particles
sedimenting in a Newtonian fluid on a 32 × 64 × 32 CFD5 computational mesh. The simulation should complete within
1 minute; if it ran correctly, a file named output.udf will be produced.
When using the multi-core executable file, you may set environment variables before launching KAPSEL to specify
the number of CPU cores used. In the following example, we request 8 CPU cores:$ export OMP_NUM_THREADS=8

2.1.6

Methods for visualizing simulation data

The use of GOURMET with python scripts allows visualization of various types of simulation data. Below we present
some simple examples.
•

Launch GOURMET– On Windows: (OCTA install dir)\GOURMET\gourmet– On Linux: /usr/local/OCTA8#/GOURMET/gourmet– On Mac: /usr/local/OCTA8#/GOURMET/gourmet•

Open the output.udf file

File -> Open -> output.udf [1]Positions and velocities for all particles at each time step are stored in the Particle[] variables. Use the slide
bar at the bottom of the GOURMET viewer window to view variables at other time steps.
•

Create animations in GOURMET

1. Click the Load button within the Python panel at the bottom of the GOURMET viewer window. [2]2. Open the particleshow.py file and click the Run button. [3]3. A new window opens. Click the Play button in this window. [4]•Plot data in GOURMET (gnuplot)

1. Click the Load button in the Python panel at the bottom of the GOURMET viewer window. [2]5 Computational Fluid Dynamics

CHAPTER 2. INSTALLING AND RUNNING KAPSEL

11

2. Open the plot.py file and click the Run button. [3]3. From the View box at the top of the window, click Table. [4]4. Select Graph Sheet[] from the left side of the window. [5]
5. Select the Plot panel at the bottom of the window and enter the following commands into the commandbox. [6,7]
6. Click the Plot button to display the time evolution of Vx for a particle (particle #1) [8]plot 'plot.dat' using 2:7 with lines

2.1.7

Analyzing simulation data

Computational data (particle coordinates and velocities at each time step) produced by running a simulation are saved
to the file output.udf. The data in this file may be accessed by any of the following methods.
•

Using python codes to analyze data

To use python codes, import the UDFManager.py module6 to access data in the UDF files. The attachedfile sk.py is a python script that computes static structure factors S (k) from time-series data on particle positions stored in output.udf.
Before using this script, you must set up the numpy package. Then you may edit the script as appropriate to
analyze simulation data for various purposes.
6 For further details on this module, see the section named “GOURMET PYTHON Script reference manual” in the documented titled “OCTA:
A comprehensive simulator for soft materials.”CHAPTER 2. INSTALLING AND RUNNING KAPSEL

12

To use the sk.py script, first launch GOURMET:– On Windows: Start Menu > All Programs > OCTA > StartGourmetTerm– On Linux/Mac: /usr/local/OCTA8#/GOURMET/gourmetterm
Then run the following commands:
python sk.py
gnuplot>> plot 'sk.dat' w line•

Analyzing data in Jupyter Notebook

Proceed as follows to repeat the above analysis in Jupyter Notebook7$ . $PF_FILES/bin/gourmet_profile.sh$ . $PF_FILES/bin/platform_env.sh$ jupyter notebook sk.ipynb•

Running analyses in Fortran or C

To use Fortran or C, use the libplatform library8 to access data in UDF files.
7 Distributed with Anaconda and other distributions.
8 For further details on this library, see the section named “libplatform: A platform interface library reference manual” in the documented titled“OCTA: A comprehensive simulator for soft materials.”CHAPTER 2. INSTALLING AND RUNNING KAPSEL13KAPSEL writes many important data quantities (such as pressures at each time step) to stderr. Typically thesedata items will be displayed on the command line, but you may also proceed as follows to redirect stderr to a file:•

For csh or tcsh:$ ../kapsel•

For sh, bash, or Windows command prompt:$ ../kapsel

2.2-Iinput.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf >& out-Iinput.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf 2> out

Input UDF and definition UDF files

The input UDF file (conventionally named input.udf) is used to configure KAPSEL settings appropriate for a given
simulation. The structure of input.udf is specified in the definition UDF file (conventionally named define.udf),
and thus the input.udf and define.udf files must be of the same version. input.udf consists of multiple configuration sections; these are described in the next few subsections, and more detailed explanation may be found in
Appendix A.

2.2.1

Fluid settings

Each choice of constitutive eq involves its own distinct set of fluid settings. These are described in the “Input
UDF file” subsections of Sections 3-8.

2.2.2

Object (particle) settings

The property object type.type specifies the type of particle. The possible settings are spherical particle,
chain (for flexible chains), and rigid (for rigid bodies).

CHAPTER 2. INSTALLING AND RUNNING KAPSEL

14

A repulsive interaction between particles, based on the excluded volume of the particles, may be introduced by
specifying a truncated Lennard-Jones (LJ) potential. This potential takes the form h 2n  n i+ϵ (r < 21/n σ),
4ϵ σr − σr
U LJ (r) = 
(2.1)

0
(r > 21/n σ)with σ = d = 2a, ϵ, and r denoting the particle diameter, the interaction strength, and the distance between the ith
and jth particles, respectively. The parameter n may be set to 6, 12, or 18, with larger values yielding a more sharply
varying potential. To introduce only a repulsive interaction based on the excluded volume, set swithch.LJ truncate
to ON. If you wish to introduce an attractive interaction between particles, setting swithch.LJ truncate to OFF
enables computation of the usual Lennard-Jones potential with a cut-off distance rcut = 2.5σ.
When configuring the exponents of the LJ potential, the choice DLVO selects the DLVO potential, an interaction
between charged particles i and j, having charge valences of zi and z j , separated by distances that are large compared
to atomic length scales [2].
(r ≥ r∗ )VCoul (r) + VvdW (r)(2.2)
VDLVO (r) = VCoul (r) + Vexcv (r)(r < r∗ )Here,
VCoul (r) =e2
exp(κa)zi z j4πεr ε0κa + 1!21exp(−κr)r(2.3)

represents the screened Coulombic potential between charged particles,
"
!#
A
d2
d2
r2 − d2
VvdW (r) = −++2
ln
12 r2 − d2 r2
r2

(2.4)

is the van der Waals potential acting between spherical particles, at distances much larger than atomic length scales,
and
Vexcv (r) =K(r − d)2 − Vm
2

(2.5)

mimics the excluded volume potential to avoid the occurrence of overlapping particles. Vvdw and Vexcv are to besmoothly switched at an interparticle distance r = r∗ (> d ≡ 2a), thus the following two parameters are automaticallyset given the value of A!
A r∗ d211r∗2 − d2
K =
(2.6)
+
+
6 r∗ − d (r∗2 − d2 )2 r∗4
r∗6
"
!#
K ∗
A
d2
d2
r∗2 − d2Vm =(r − d)2 +
+
+
2
ln
(2.7)
2
12 r∗2 − d2 r∗2
r∗2
2

Figure 2.1 shows the schematic plots for these three potentials. The values for zi , r∗ , κa, C ≡ 4πεe r ε0 , A are set
through the following variables defined in input.udf.
•

object type.spherical particle.Particle spec[i].Surface charge = zi•

DLVO.n = r∗ /d•DLVO.kappa a = κa•

DLVO.vdw coeff = A•

DLVO.coulomb coeff = C

The choice electro osmotic flow for the exponents of the LJ potential selects an interaction potential mediatedby electro-osmotic flow, relevant for charged particles in the vicinity of a slab. This is only used for the Quincke rollers
discussed in Section 8.
Flexible-chain particles and rigid-body particles are implemented as secondary particles consisting of aggregatesof multiple spherical primary particles. Thus, initial configurations for these systems must be specified in input.udf

CHAPTER 2. INSTALLING AND RUNNING KAPSEL

15

Figure 2.1: Schematic plots for the three contributions to the DLVO potential. The green, red, and blue lines represent the screened
Coulombic VCoul (r), van der Waals (Hamarker) VvdW (r), and excluded volume Vexcv (r) potentials, respectively. We set
22
d = 1, r∗ = 1.02d, κa = 5, 4πεe r ε0 zi z j exp(κa)= 50000, A = 1 in this figure.
κa+1

by setting switch.INIT distribution=user specify. In a flexible chain, pairs of neighboring beads interact
through a finitely extensible non-linear elastic (FENE) potential of the form
1
U F (r) = − kc R20 ln{1 − (r/R0 )2 },
2

(2.8)

which is responsible for binding. Here r is the distance between neighboring beads and the parameters have values
kc = 30ϵ/σ2 , R0 = 1.5σ, If a rigid body is selected, the initial particle structure is preserved by an implicit constraint
force.
➤object type.spherical particle.Particle spec[]: Properties of spherical particles
object type.spherical particle.Particle spec[].Particle number
Number of particles.
object type.spherical particle.Particle spec[].MASS Ratio
Ratio of particle density to fluid density.
object type.spherical particle.Particle spec[].Surface charge
Surface charge (not used for simulations of two-component phase-separated fluids).
object type.spherical particle.Particle spec[].janus axis
Orientation of the Janus axis in a body-fixed coordinate system.
object type.spherical particle.Particle spec[].janus propulsion
Propulsive motion of Janus particles. Possible values: OFF (disabled), TUMBLER (particle propelled alongpropulsion axis by fixed external force), SQUIRMER (squirmer particle propelled along propulsion axis byslip boundary conditions), OBSTACLE (particle forming a fixed obstacle).
object type.spherical particle.Particle spec[].janus force.x
x component of propulsive force.
object type.spherical particle.Particle spec[].janus force.yy component of propulsive force.

CHAPTER 2. INSTALLING AND RUNNING KAPSEL

16

object type.spherical particle.Particle spec[].janus force.z
z component of propulsive force.
object type.spherical particle.Particle spec[].janus torque.x
x component of propulsive torque.
object type.spherical particle.Particle spec[].janus torque.yy component of propulsive torque.
object type.spherical particle.Particle spec[].janus torque.z
z component of propulsive torque.
object type.spherical particle.Particle spec[].janus slip vel
The parameter B1 determining the surface-slip velocity for self-propelled particles.
object type.spherical particle.Particle spec[].janus slip mode
The parameter B2 /B1 determining the type of slip motion for squirmer particles: pusher, neutral, or
puller.9
object type.spherical particle.Particle spec[].janus rotlet C1
The parameter C1 describing unipolar rotation about the propulsion axis.
object type.spherical particle.Particle spec[].janus rotlet dipole C2
The parameter C2 describing bipolar rotation about the propulsion axis.
➤object type.chain.Chain spec[]: Properties of flexible chainsobject type.chain.Chain spec[].Beads number
Number of beads associated with a single chain.
object type.chain.Chain spec[].Chain number
Number of chains.
object type.chain.Chain spec[].MASS RATIO
Ratio of bead density to fluid density.
object type.chain.Chain spec[].Surface Charge
Surface charge (not used for simulations of two-component phase-separated fluids).
object type.chain.Chain spec[].janus axis
Orientation of the Janus axis in a bead-fixed coordinate system.
➤object type.rigid.Rigid spec[]: Properties of rigid bodies
object type.rigid.Rigid spec[].Beads number
Number of beads comprising the rigid body.
object type.rigid.Rigid spec[].Rigid number
Number of rigid bodies.
9 More specifically, when this parameter is {negative, zero, positive} we have a {pusher, neutral, puller} [3, 4].

CHAPTER 2. INSTALLING AND RUNNING KAPSEL

17

object type.rigid.Rigid spec[].MASS RATIO
Ratio of bead density to fluid density. The mass of overlapping beads is not double-counted.
object type.rigid.Rigid spec[].Surface charge
Total surface charge: z (not used for simulations of two-component phase-separated fluids).
object type.rigid.Rigid spec[].Rigid motion
Selects the type of motion executed by rigid bodies: either free (free motion) or fix (motion with
fixed translational and rotational velocities). For fix, the (lab-frame) velocities are as specified by the
Rigid velocity and Rigid omega variables defined below. Note that, by default, this will fix all 6 degrees of freedom. To constrain individual translational/rotational degrees of freedom, use the switch.free rigid
options described below.
object type.rigid.Rigid spec[].Rigid velocity.x
x component of rigid body velocity (lab frame).
object type.rigid.Rigid spec[].Rigid velocity.yy component of rigid body velocity (lab frame).
object type.rigid.Rigid spec[].Rigid velocity.z
z component of rigid body velocity (lab frame).
object type.rigid.Rigid spec[].Rigid omega.x
x component of rigid body angular velocity (lab frame).
object type.rigid.Rigid spec[].Rigid omega.yy component of rigid body angular velocity (lab frame).
object type.rigid.Rigid spec[].Rigid omega.z
z component of rigid body angular velocity (lab frame).

2.2.3

Common simulation settings➤Standalone UDF properties for particle radius and interface thickness
A XIInterface thickness ξ.10
A
Particle radius.
➤gravity: Properties of gravitational forces
gravity.G
Gravitational acceleration.
gravity.G direction
Direction in which gravitational forces are applied.
➤Standalone UDF properties describing the character of interparticle forces
10 For simulations of electrolytes or two-component phase-separated fluids, we must compute the gradient of ϕ for the GL free energy F, forwhich purpose it is best to choose ξ ≥ 2.

CHAPTER 2. INSTALLING AND RUNNING KAPSEL18

EPSILON
Defines the units of energy in the Lennard-Jones potential.
LJ powers
Selects one of the various supported forms for the interparticle potential. Allowed values:
12:6
24:12
36:18
macro vdw
electro osmotic flow

Lennard-Jones potential with the given exponents
macroscopic interparticle potential
used only for Quincke rollers➤mesh: Computational mesh sizes
mesh.NPX
Base-2 logarithm of mesh size in x-direction (L x = 2NPX )mesh.NPY
Base-2 logarithm of mesh size in y-direction (Ly = 2NPY )mesh.NPZ
Base-2 logarithm of mesh size in z-direction (Lz = 2NPZ )➤time increment: Simulation timestep settings
time increment
Specifies how KAPSEL determines the simulation timestep. Allowed values:
KAPSEL will automatically set the timestep to its maxi-

automanual2
mum allowable value, given by T step = ρ/ηkmax, wherekmax is the maximum wavenumber determined by the lattice spacing ∆.
A fixed timestep will be specified by the user.

time increment.auto.factor
Specifies a multiplicative scale factor for auto-determined timesteps. The simulation timestep is given by∆t = factor × T step .
time increment.manual.delta tUser-specified timestep ∆t.

2.2.4

Configuring selectable features➤switch: Miscellaneous simulation parametersswitch.ROTATION
Set to ON to solve the equations of motion for the rotational motion of particles.
switch.LJ truncate
Selects the form of the Lennard-Jones potential acting between particles. Set to OFF for the usual formof the potential, including the attractive term. Set to ON to exclude the attractive term, retaining only therepulsive term. Set to NONE to exclude both terms, i.e. to disable the interparticle potential entirely.

CHAPTER 2. INSTALLING AND RUNNING KAPSEL

19

switch.INIT distribution
Specifies the initial particle configuration. Possible values: uniform random (random), random walk(particles randomly moved from the sites of a square lattice), FCC (FCC lattice), BCC (BCC lattice),user specify (user-specified coordinates and velocities).11
switch.INIT distribution.random walk.iteration
Number of trial iterations for the random walk setting of INIT distribution.
switch.INIT orientation
Initial particle orientation. Possible values: user specify (orientations specified by the user), random(random), space align (body-axis of particles aligned with space/lab axis).12
switch.SLIP tol
Convergence criterion for iterative calculation of fluid-flow fields when introducing slip velocities in thetangential direction at the interfaces of self-propelled particles.
switch.SLIP iter
Maximum iteration count for iterative calculation of fluid-flow fields when introducing slip velocities inthe tangential direction at the interfaces of self-propelled particles.
switch.FIX CELL.x
Set to ON to zero out the DC component of the total velocity in the x direction.
switch.FIX CELL.y
Set to ON to zero out the DC component of the total velocity in the y direction.
switch.FIX CELL.z
Set to ON to zero out the DC component of the total velocity in the z direction.
switch.pin.type
Set to YES to fix particle positions.
switch.pin.YES.pin[]
Specify the indices of any particles for which translational motion is to be prohibited.
switch.pin.YES.pin rot[]
Specify the indices of any particles for which rotational motion is to be prohibited.
switch.free rigid.type
Configure individual translational/rotation degrees of freedom (DOFs) for rigid bodies. These option aremeant to be used together with the object type.rigid.Rigid spec[].Rigid motion options, which
fix all six DOF to have specified (lab-frame) translational/rotational velocities.
To allow the six DOF to be configured independently for each rigid species, set this option to YES, and use
the vel.x|y|z and omega.x|y|z options decribed below to selectively free (YES) or fix (NO) the corresponding DOF. The fixed value of each velocity component is the one specified by the Rigid velocity
and Rigid omega variables defined via the object type.rigid.Rigid spec[] options.
For example, to specify a rigid particle with fixed vel.y and omega.z, you should specify the following free rigid options below: vel.x=YES, vel.y=NO, vel.z=YES and omega.x=YES, omega.y=YES,
omega.z=NO; i.e., the convention here is to specify which degrees of freedom are to be set free.
11 If user specify is chosen, particle positions and velocities are initialized to the values specified for user specify.Particles[].R anduser specify.Particles[].v. If the number of items in the input list is less than the value specified for Particle number, increase the lengthof user specify.Particles[], either by directly editing the UDF file or by using Edit->Add an array Element in GOURMET.
12 If user specify is selected, particle orientations will be initialized to the values specified for user specify.Particles[].q.

CHAPTER 2. INSTALLING AND RUNNING KAPSEL

20

switch.free rigid.YES.DOF[].spec id
Specify the indices of rigid bodies to configure.
switch.free rigid.YES.DOF[].vel.x
Configure degrees of freedom for x-directed translational motion of rigid bodies. Set to YES to allow
motion or NO to fix the position.
switch.free rigid.YES.DOF[].vel.y
Configure degrees of freedom for y-directed translational motion of rigid bodies. Set to YES to allow
motion or NO to fix the position.
switch.free rigid.YES.DOF[].vel.z
Configure degrees of freedom for z-directed translational motion of rigid bodies. Set to YES to allow
motion or NO to fix the position.
switch.free rigid.YES.DOF[].omega.x
Configure degrees of freedom for x-directed rotational motion of rigid bodies. Set to YES to allow motion
or NO to fix the orientation.
switch.free rigid.YES.DOF[].omega.y
Configure degrees of freedom for y-directed rotational motion of rigid bodies. Set to YES to allow motion
or NO to fix the orientation.
switch.free rigid.YES.DOF[].omega.z
Configure degrees of freedom for z-directed rotational motion of rigid bodies. Set to YES to allow motion
or NO to fix the orientation.
ns solver.OBL INT
Set to linear or spline to select the approximation function used for coordinate transformations inshear-flow simulations.
switch.wall.type
Set to FLAT to specify planar walls.
switch.wall.FLAT.axis
Set to X, Y, or Z to specify the direction normal to planar walls.
switch.wall.FLAT.DH
Thickness of planar walls, measured in units of lattice spacing.
switch.wall.FLAT.LJ Params
Specifies how the potential acting at planar walls is chosen. If set to AUTO, the potential will be the same as
the interparticle potential. If set to MANUAL, the potential is user-specified via the switch.wall.FLAT.MANUAL
property.
switch.wall.FLAT.MANUAL.truncate
Specifies the presence or absence of attractive forces at planar walls. If set to ON, the potential is purely
repulsive. If set to OFF, the potential is a sum of repulsive and attractive contributions.
switch.wall.FLAT.MANUAL.powers
Specifies the exponents in the Lennard-Jones potential at planar walls. Possible values: 12:6, 24:12,
36:18.

CHAPTER 2. INSTALLING AND RUNNING KAPSEL

21

switch.wall.FLAT.MANUAL.EPSILON
Determines the strength of the Lennard-Jones potential at planar walls.
switch.quincke.type
Set to ON when simulating Quincke rollers.
switch.quincke.ON.e dir
Set to X, Y, or Z to specify the direction of the external electric field.
switch.quincke.ON.w dir
Set to X, Y, or Z to specify the direction of the angular-velocity vector associated with rotation induced by
the Quincke effect.
switch.quincke.ON.torque amp
Magnitude of rotational torque.
switch.multipole.type
Set to ON to run simulations using the Ewald method.
switch.multipole.ON.Dipole
Set to ON to use dipolar Quincke particles.
switch.multipole.ON.Dipole.ON.magnitude
Magnitude of dipole moment.
switch.multipole.ON.Dipole.ON.type
Set to FIXED or QUINCKE to specify the type of dipole.
switch.multipole.ON.Dipole.ON.FIXED.dir
Set to X, Y, or Z to specify the direction of a fixed dipole.
switch.multipole.ON.EwaldParams.alpha
Screening parameter for Ewald method.
switch.multipole.ON.EwaldParams.delta
Convergence criterion for determination of kmax via the Ewald method.
switch.multipole.ON.EwaldParams.converge
Convergence parameter for Ewald method.
switch.multipole.ON.EwaldParams.epsilon
Dielectric permittivity at boundaries for Ewald method.

2.2.5

Data output settings➤output: Data output settings
output.GTS
Data output interval measured in number of steps.

CHAPTER 2. INSTALLING AND RUNNING KAPSEL

22

output.Num snap
Number of data outputs. The total number of timesteps is given by GTS × Num snap.
output.AVS
Set to ON to output data in AVS format.
Reference: https://ja.overleaf.com/project/5fe7389e29fad9c51e06569e
output.AVS.ON.Out dir
Specifies the directory in which AVS-format data output files are written.13
output.AVS.ON.Out nameSpecifies the filename prefix for AVS-format data output files.
output.AVS.ON.FileTypeSpecifies the format of AVS data output files. Possible values: Binary, ASCII, EXTENDED
output.ON.EXTENDED.Driver.Format
Specifies the extended output data format. At present, the only supported setting is HDF5.
output.ON.EXTENDED.Print field.Crop
Set to YES to reduce the number of output data fields.
output.ON.EXTENDED.Print field.YES.Slab x.start
Index of the first lattice point in the x direction to output.
output.ON.EXTENDED.Print field.YES.Slab x.countNumber of lattice points in the x direction to output.
output.ON.EXTENDED.Print field.YES.Slab x.stride
Interval between output lattice points in the x direction.
output.ON.EXTENDED.Print field.YES.Slab y.start
Index of the first lattice point in the y direction to output.
output.ON.EXTENDED.Print field.YES.Slab y.countNumber of lattice points in the y direction to output.
output.ON.EXTENDED.Print field.YES.Slab y.stride
Interval between output lattice points in the y direction.
output.ON.EXTENDED.Print field.YES.Slab z.start
Index of the first lattice point in the z direction to output.
output.ON.EXTENDED.Print field.YES.Slab z.countNumber of lattice points in the z direction to output.
output.ON.EXTENDED.Print field.YES.Slab z.stride
Interval between output lattice points in the z direction.
13 For example, if this property is set to data, then the directories ./data and ./data/avs must be created in advance. The AVS field file will
be written to a file named data.fld in the ./data directory. Data files will be written in the directory ./data/avs and will have filenames of theform data *.dat, where * will be replaced by the step count.

CHAPTER 2. INSTALLING AND RUNNING KAPSEL

23

output.ON.EXTENDED.Print field.Vel
Set to YES to output velocity-field values.
output.ON.EXTENDED.Print field.Phi
Set to YES to output the SP function ϕ.
output.ON.EXTENDED.Print field.Charge
Set to YES to output the charge density distribution (not used for simulations of two-component phaseseparated fluids).
output.ON.EXTENDED.Print field.Pressure
Set to YES to output the pressure field (currently not implemented).
output.ON.EXTENDED.Print field.Tau
Set to YES to output the stress tensor.
output.UDF
Set to ON for UDF output.14

2.2.6

Restart settings➤resume: Handling interrupted calculations
resume.calculation
Specifies how KAPSEL proceeds in the event of an interrupted calculation. Allowed values:
NEW
CONTINUE
CONTINUE FDM
CONTINUE FDM PHASE SEPARATION

Start new calculation.
Read data saved upon termination of previous calculation
and resume that calculation.

Which of the three CONTINUE options to use depends on the type of simulation you are running:•

Use CONTINUE for simulations using spectral methods.15•

Use CONTINUE FDM for simulations of single-component fluids via finite-difference methods.16

•

2.2.7

Use CONTINUE FDM PHASE SEPARATION for simulations of two-component fluids via finite-difference
methods.17

GOURMET display settings➤Unit Parameter: Display settings for GOURMET. Does not affect simulation results.
Unit Parameter.Name
Name of UDF file.
Unit Parameter.Comment
Comment for UDF file.
14 Coordinates and velocities for each particle will be written to the Particles[] section of the output UDF file. The number of data records
saved will be the value specified for Num snap.
15 Thisincludesthefollowing
choicesforconstitutive eq:
Navier Stokes,
Shear Navier Stokes,
Shear Navier Stokes Lees Edwards, or Electrolyte.
16 This includes the following choices for constitutive eq: Navier Stokes FDM or Shear Navier Stokes Lees Edwards FDM.
17 This includes the following choices for constitutive eq: Navier Stokes Cahn Hilliard FDM or Shear NS LE CH FDM.

CHAPTER 2. INSTALLING AND RUNNING KAPSEL

24

Unit Parameter.Temperature
Temperature unit, used as a reference when simulation parameters are reported in actual units in GOURMET.

Unit Parameter.Length
Length unit, used as a reference when simulation parameters are reported in actual units in GOURMET.
Unit Parameter.rho
Density unit, used as a reference when simulation parameters are reported in actual units in GOURMET.

2.3

Output and restart UDF files

The output UDF file (output.udf) accumulates data reported by KAPSEL simulations at fixed simulation-time intervals. The interval (number of timesteps) between data reports, and the total number of timesteps in a simulation, may
be configured via the following properties in input UDF files:
output.GTS
Data output interval measured in number of steps.
output.Num snap
Number of data outputs. The total number of timesteps is given by GTS × Num snap.
The restart UDF file (restart.udf) stores all data written at the end of a simulation.18 To use restart.udf toresume a completed calculation, follow the steps below.
1. Rename restart.udf to input2.udf.
2. Open input2.udf in Gourmet and set resume.Calculation to one of the following values as appropriate forthe simulation in question: CONTINUE, CONTINUE FDM, CONTINUE FDM PHASE SEPARATION
3. Increase the value specified for output.Num step by the number of additional timesteps you wish to compute.
Thus, the new value of output.Num step should be (number of timesteps previously computed) + (number of
additional timesteps to compute).
4. The values of various parameters and the number of particles can also be changed on restart. If particles areto be added/changed, the objects type block should be set appropriately and the additions/changes should
be applied to the Particles in the resume block. INIT distribution/INIT orientation in the switch
block is ignored on restart.
5. Execute the following command line:% kapsel -Iinput2.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
This will read data from input2.udf and resume the calculation starting at the endpoint of the previous calculation.
6. For both output.udf and AVS data files, output data from the resumed simulation may be appended to theexisting set of output data from the previous simulation.
7. When the resumed calculation is complete, all data will again be written to restart.udf, and the above procedure may be repeated to extend the calculation yet again.

18 In restart UDF files, results are saved under the item resume.CONTINUE.Saved Data. Consequently, continuation of a previous calculation
may be initiated by using the restart.udf file written by the previous calculation as the input UDF file for a new KAPSEL run.

Chapter 3

Simulating particles dispersed in Newtonian
fluidsSystems consisting of colloidal particles dispersed in liquids (solvents), commonly known as colloidal “suspensions”,
are ubiquitous in many sectors of modern daily life: foodstuffs, paints, pigments and dyes, cosmetics, slurries, andmore. Experience with microparticle suspensions shows that any hope of understanding the properties and behaviorof these systems depends crucially on having a firm grasp of the nature of the motion of both the colloidal particles
and the solvent.

3.1

Theoretical background and basic equations

KAPSELconducts simulations for several types of dispersed particle systems using a direct numerical approach known
as the smoothed profile (SP) method [5, 6], which takes into account both the solid particles and the fluid, and solves
their equations of motion simultaneously. The method is not restricted to simulations of Newtonian fluids, but may be
applied to systems in which the solvents are complex fluids with arbitrary Reynolds numbers.

3.1.1

Basic equations for disperse particle systems

The SP method allows us to solve for the motion of solid particles dispersed in a continuous medium [6]. In principle,
the SP method is applicable to solid particles of arbitrary shapes.
Also, the method is not restricted to simulating simple liquids, but can also handle complex multi-phase fluids
and other complicated systems, for which it makes liberal use of appropriate constitutive equations to enable directcomputational modeling.
To illustrate how the motion of particles and fluids ties in to the SP method, in this chapter we consider the motionof a disc-shaped rigid (solid) body moving freely in a Newtonian fluid [7]. As shown in Figure 3.1, the solid body is
described as an assembly of spherical beads; in what follows we will simply refer to this body as the particle.
The equations governing incompressible Newtonian fluids are the continuity equation and the Navier-Stokes equations:∇ · u f = 0,∂t + u f · ∇ u f = ρ−1
f ∇ · σ,t 
σ = −pI + η ∇u f + ∇u f

(3.1)
(3.2)
(3.3)

Here u f , ρ f , and η are the fluid velocity, density, and viscosity, while σ is the stress field (tensor); p is the pressure, I
is the identity tensor, and (· · · )t denotes matrix transposition. Unless stated otherwise, henceforth we assume that theparticle and fluid densities are equal, ρ f = ρ p = ρ. At solid-fluid interfaces we impose no-slip boundary conditions,
i.e., at solid surfaces we enforce u f = u p where u p is the particle velocity at the solid surface. Denoting by MI and II
the mass and inertial moment of the Ith particle (I = 1, · · · , N), the motion of a solid particle (here comprised of N
rigid spheres) is described by the Newton-Euler equations [8]:ṘI = VI ,MI V̇I=Q̇I = skew(ΩI ) · QI ,(3.4)

FIH + FIC + FIE + GVI ,

(3.5)

25

3.2

度の

己

先ほどの

度の

己

己

を

ことで円 のよりもはっきりとめた。その、以下の図に度を めた。またそのいがわかるようになると
すような円

えて

度の

に固定された座をの式を以下にす。

度の 己CHAPTER 3. SIMULATING PARTICLES DISPERSED IN NEWTONIAN FLUIDSえる26

図 7 円 に固定された座
Figure 3.1: Schematic illustration of a disc-shaped solid body moving in a three-dimensional space. The motion of the body ischaracterized by its three translational and three rotational degrees of freedom.
݀ ଶ ஶ ሺ‫ݐ‬ሻ ሺ‫ ݐ‬൅ ሻ ൅ ሺ‫ݐ‬ሻ ሺ‫ ݐ‬൅ ሻ(3.2.1)‫ כ ܥ‬ሺ߬ ‫ כ‬ሻ ൌ ൬
൰ න
‫ݐ‬
ʹ ܷ
ʹ
଴

J˙I

NIH + NIC + NIE + GΩ
I ,

=

(3.6)

ܷ
(3.2.2)
Here RI and VI are the position and velocity of the ൌcenter݀ of mass, QI and ΩI are the corresponding orientation‫כ‬matrix and angular velocity, and JI = II · ΩI is the angular momentum. The orientation matrix relates the basis
されている
y' 、z'の度である。さらに円y'、Ȧ
vectors {ei } (iここで= 1, 2,Ȧ3)
of z'aは図
fixed 中でorthogonal
coordinate
system周りの円(the laboratory
system)to the basis vectors {ẽi }Pは x' を中心にであるので、ω
とωthe
価であるとえ、それぞれのof a time-dependent
coordinate system
co-rotating y'withaccordingto ẽi = 3j=1 e己
z'は particle
j Q ji . The temporal
evolution of the rotation matrix is governed by the (skew-symmetric) angular-velocity matrix, which is determined by
P3
ௗ
kを平均したものを度の己をଶ ௎で according
化するこ(with ϵi jk the Levi–Civita
symbol)tothe angular-velocity
vectors [skew(Ω)]
i j = −己 k=1 ϵi jk Ω とした。また
 0 4−Ωz でのΩy  を す。とでそれぞれをしやすい形にした。以下に
 z


0
−Ω x  .
skew(Ω) =  Ω
(3.7)
 y

x
−Ω
Ω
0
Particles experience both forces F and torques N from their surroundings. There are three contributions to theseforces and torques: a contribution F C arising from direct interparticle interactions, a contribution F E due to gravity
or other external fields, and a contribution from the solvent. The contribution from the solvent is the sum of two
terms: one term F H representing hydrodynamic interactions between particles, and a second term F V representing
thermal fluctuations. Note that the most accurate choice of notation for the various force and torque contributions
differs depending on the system in question; we will discuss each case individually. We also introduce random forces
and torques GV and GΩ , acting on particles due to thermal fluctuations, which we model as white noise. That is, wehave ⟨GVI ⟩ = ⟨GΩ
I ⟩ = 0 and⟨GVI (t) · GVJ (0)⟩ = 3kB T αV δ(t)δI J ,(3.8)= 3kB T αΩ δ(t)δI J

(3.9)

Ω
⟨GΩ
I (t) · G J (0)⟩where ⟨⟩ denotes statistical averaging and αV and αΩ are10
adjustable parameters used to ensure that the particle temperature T accurately tracks the specified value. The particle temperature T is determined in advance by simulating asingle particle fluctuating in a solvent and adjusting T until the translational and rotational diffusion coefficients DV
and DΩ agree with the theoretical values.
We now explain how the SP method is used to solve for the motion of a system described by the above equations.
At every point x in the spatial region of the simulation—including both solvent regions and solid regions in the particle
domains—the SP approach defines a velocity field according tou(x, t) = (1 − ϕ) u f (x, t) + ϕu p (x, t),(3.10)

where
ϕu p (x, t) =NX

ϕI [VI + ΩI × rI ] ,(3.11)

I=1

expresses the rigid-body motion of solid particles in terms of the velocity field. Here rI = x − RI is the distance fromthe Ith particle’s center of mass and ϕI = ϕI (x; RI , QI ) ∈ [0, 1] is a phase-field function that distinguishes between
solid regions inside particles (ϕ ≃ 1) and fluid regions exterior to particles (ϕ ≃ 0).

CHAPTER 3. SIMULATING PARTICLES DISPERSED IN NEWTONIAN FLUIDS

27

The SP method replaces the usual sharp solid-liquid interface with a blurred finite-thickness interface region, with
the interface thickness ξ incorporated into the phase-field function ϕ; in this document we refer to this as a smoothedprofile (SP) function. Replacing discontinuous solid-fluid interfaces with finite-width transition regions yields two
key advantages for numerical calculations. First, the positions of solid-fluid interfaces may be tracked using only alimited number of lattice points. Second, interactions between the solid and fluid—such as hydrophilic or hydrophobic
behavior—may be properly specified as appropriate for the objectives of a simulation.

3.1.2

Outline of the SP method

Applying the continuity equation and the Navier-Stokes equations to the full velocity field u(r)—which includes bothparticle and fluid contributions—yields∇ · u = 0,

(3.12)(∂t + u · ∇) u = ρ ∇ · σ + ϕf p ,−1(3.13)

where ϕf p acts as a penalty force to satisfy the rigid-body constraints in the particle regions. Also, whereas equation(3.3) expresses the local stress tensor in terms of the fluid contribution to the velocity (u f ), in the SP method stress isdetermined by the total velocity (u), defined even within the solid particle domains.
The quantities ϕf p , FIH , and NIH are defined for the interval between discrete timesteps (n → n + 1). We denoteby un the velocity field at the nth timestep, i.e., at time tn = nh (with h the timestep).
We begin by neglecting the ϕf p terms in (3.13) and (3.4) and integrating those equations from tn to tn + h, yielding
the following updates to the particle position and direction:u =u +∗n

Z tn +hhids∇ · ρ−1 σ − uu ,(3.14)

dsVI ,(3.15)

ds skew(Ω) · QI .

(3.16)tn

Rn+1= RnI +I
Qn+1
= QnI +
I

Z tn +htn
Z tn +htn

At this point we must update the particle velocity field, even if the particle velocities have not changed. Assuming thatthe conservation of momentum holds for the velocity u∗ in the absence of the ϕf p term, the hydrodynamic forces andtorques acting on a particle are determined by the balance of momentum received by the fluid from the particle:"Z tn +h# Z
dsFIH =
dxρϕn+1u∗ − unp ,(3.17)
Itn"Z tn +h

#
dsNIH

=

Z

tn

h

i
dx rIn+1 × ρϕn+1u∗ − unp .
I(3.18)The result is that the particle velocities are updated to VIn+1 and Ωn+1
I . After these updates, equation (3.11) determinesn+1
the final particle velocity field ϕn+1 un+1
,
i.e.
ϕ
.
p
I
The last step is to determine the velocity field for the entire system at time tn+1 , which we do by considering thequantity ϕf p introduced to preserve particle rigidity:
un+1 = u∗ +"Z tn +h#dsϕf p ,(3.19)

tn"Z tn +h
tn

#


 h
∗
dsϕf p = ϕn+1 un+1
p − u − ∇p p ,ρ(3.20)The pressure due to the rigidity constraint is obtained from the continuity equation ∇ · un+1 = 0. Because viscousstresses are applied throughout the entire spatial region, the no-slip boundary conditions at relevant particle surfaces
are satisfied in the interior of the finite interface width ξ.

3.1.3

Smoothed Profile (SP) functions

In this section we describe the construction of SP functions for particles of arbitrary shapes. An SP function is afunction ϕ that takes the value ϕ = 1 in the interior of particles and the value ϕ = 0 in fluid regions outside of particles,CHAPTER 3. SIMULATING PARTICLES DISPERSED IN NEWTONIAN FLUIDS

28

interpolating smoothly between these limits over a finite-thickness interface region of width ξ. For the particular caseof spherical particles, several varieties of SP function with closed-form analytical expressions have been proposed [5];one example is [9, 10]"!
#
a − |rI |
1
tanh
+1 ,
(3.21)ϕI =2
ξwith a denoting the sphere radius. The interface width ξ is a freely adjustable parameter, but for numerical stability it
is generally advisable to choose ξ to be one or two times the spacing between points in the computational mesh.
The cross-sectional particle image in Figure 3.2 shows an interface region of width ξ separating the solid material
region inside the particle and the fluid region outside the particle.
∆

ξaFigure 3.2: Cross-sectional view of a particle in a fluid to illustrate the notion of an SP function. ∆ is the lattice spacing, athe particle radius, and ξ the width of the interface region. In this case the interface region occupies a volume of
approximately ∼ πa2 ξ at the particle surface. Reproduced from Phys. Rev. E 71, 036707[5], Copyright 2005, withpermission from the American Physical Society.

Non-spherical particles of complex shapes may be represented as aggregates of spherical beads. These beads may
overlap with each other, allowing arbitrary particles to be represented as assemblies of beads. Denoting by nI thenumber of spherical beads used for the Ith particle, an SP function appropriate for the overall configuration isΦI (x, t) =nI
X

ϕI,i (x, t),(3.22)

ΦI,
max (ΦI , 1)

(3.23)i=1

ϕI (x, t) =where ϕI,i is the SP function for the ith spherical bead in the bead assembly for the Ith particle.
Once an SP function has been defined for a particle—any particle, including particles of complicated shapes—thehydrodynamic impulses on the particle and the associated momentum impulse to the fluid are given by equations(3.17), (3.18), and (3.20).
Because the inertia of the fluid is determined by the SP algorithm discussed above, added mass effects are capturedby the hydrodynamic impulse in equation (3.17). To clarify the nature of the associated inertial effects, it is convenient
to decompose the integrand of expression (3.17), for the hydrodynamic impulse, in the form






∗
n+1
n
ρϕn+1
u∗ − unp = −ρϕn+1
un+1
un+1
p − u + ρϕI
p − upI
I
Z tn +h= −ρdsϕI f p − h∇p ptn

CHAPTER 3. SIMULATING PARTICLES DISPERSED IN NEWTONIAN FLUIDShi+ ρϕn+1
Vi n+1 − Vi n + Ωn+1
− ΩnI × rIn+1 ,I
I
Spatial integration subject toR

dxϕI r = 0 yields
ZZρ
H
FI = − dxρϕI f p − dxϕI ∇p p + MI
V̇I ,
ρ p,I

29

(3.24)(3.25)

with ρ p,I the mass density of the Ith particle. From this equation we see that the hydrodynamic force in (3.17) consistsof the counterreaction to the action of the rigidity constraint plus added mass effects.
The accuracy of SP calculations depends on both the timestep h and the interface width ξ. Luo et al. [11] reported
that computational accuracy does not increase monotonically as h decreases, and that optimal performance is attained
when h is of the same order of magnitude as ξ2 ρ/η (the time required for momentum to propagate a distance ξ via
viscous diffusion). This is a reflection of the fact that, in order to study the viscous Stokes layer resulting from theimpulse [equation (3.20)] of the rigidity constraints, it is necessary to select an appropriately commensurate value ofthe interface width.

3.1.4

Procedure for determining particle temperature

1. Fix values for αV and αΩ and simulate a one-particle system at thermal equilibrium.
2. From your simulation results, compute values for the translational and rotational diffusion coefficients DVsim andDΩsim by considering mean-square displacements or equivalent quantities.
3. Compare your values for DVsim and DΩsim to the analytical expressions for the diffusion coefficients in an infinitely dilute one-particle system. (These expressions are: DV0 = kB T V /6πηa for translational motion, DΩ
0 =kB T Ω /8πηa3 for rotational motion.) Based on this comparison, compute the temperatures T V and T Ω associatedwith the translational and rotational motion of the particle.
4. If T Ω , T V , adjust the values of αV and αΩ and repeat the calculation. Continue in this way until you havefound values of αV and αΩ for which T Ω = T V (= T ).
Further discussion of this temperature-determination strategy may be found in Refs.[12–14].

CHAPTER 3. SIMULATING PARTICLES DISPERSED IN NEWTONIAN FLUIDS

3.2

Input UDF settings

3.2.1

Fluid settings

30

The motion of colloidal particles in a Newtonian fluid may be simulated by setting constitutive eq to Navier Stokes.
The settings that may be configured via input UDF files for this case are described below.
➤constitutive eq.Navier Stokes: Solvent properties
constitutive eq.Shear Navier Stokes.DX
The lattice spacing ∆, which defines the length unit.
constitutive eq.Shear Navier Stokes.RHO
Density of solvent.
constitutive eq.Shear Navier Stokes.ETA
Viscosity of solvent.
constitutive eq.Shear Navier Stokes.kBT
Particle temperature.
constitutive eq.Shear Navier Stokes.alpha v
Correction term for particle temperature associated with translational motion.1
constitutive eq.Shear Navier Stokes.alpha oCorrection term for particle temperature associated with rotational motion.

3.2.2

Configuring objects (particles)

The object type.type setting specifies the type of particle. The possible values are spherical particle forspherical particles, chain for flexible chains, or rigid for rigid bodies.

3.2.3

Choice of length and time units

The length unit is given by the lattice spacing ∆. The time unit τ0 is determined by the fluid density ρ, the fluid
viscosity coefficient η, and the lattice spacing ∆ according to τ0 = ρ∆2 /η.
•••

For the Navier-Stokes equations, we choose a system of units in which ρ = η = ∆ = 1.
Assuming the input UDF file specifies RHO= A, ETA= B, and DX= C, the maximum wavenumber kmax may
be determined from C and the upper bound on the timestep is determined from the momentum diffusion time2in the form T step = (A/B)/kmax. The actual simulation timestep is related to this bound through a multiplicative
scale factor, i.e. ∆t = factor × T step .
Consider the correspondence between physical quantities in simulations and in reality. A grid spacing of 1 µm
is reasonable for the length scales we wish to consider. Assuming we use water (η = 1 × 10−3 Pa s, ρ =1 × 103 kg m−3 ) as a solvent, the time unit is then τ0 = 1 × 10−6 s.

1 When adding thermal diffusion of microparticles to simulations, you must specify a particle temperature. In principle, particle temperatures
may be estimated from the diffusive motion of particles at thermal equilibrium; however, this approach may yield incorrect temperatures due tonumerical factors such as interface thickness. In such cases, alpha v may be used to correct temperature errors. If the measured temperature
is too low by a factor of n (that is, the measured temperature is 1/n times the correct temperature), then setting alpha v to n will reproduce thecorrect temperature (and similarly with alpha o for the rotational temperature). For a particle radius of 4 or 5 and interface thickness ξ = 2, settingalpha v=alpha o=1 seems to reproduce specified temperatures.

CHAPTER 3. SIMULATING PARTICLES DISPERSED IN NEWTONIAN FLUIDS

3.3

Computational examples

3.3.1

Sedimentation of particles

31

The input UDF files for this example are gravity.udf and gravity100000.udf in the Examples/02/ folder.
These files describe simulations of particles suspended in fluid sedimenting under the influence of gravity; the numberof particles is 3204 (gravity.udf) or 100000 (gravity1000000.udf). To run these simulations, first confirm thatthe folders ./avs g1/ and ./avs g1/avs/ have been created on your system,2 then execute the following commands:$ ../../kapsel -Igravity.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf$ ../../kapsel -Igravity100000.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
For this example (Fig. 3.3), we choose a 256 × 256 × 512 computational mesh. In the gravity.udf file we specify
the following values for simulation parameters: number of particles N p = 3204, particle diameter D = 10, interface
thickness ξ = 2. With these settings, the fractional volume occupied by the particles is φ = 0.05. Figure 3.3 shows asnapshot of the 3204-particle suspension sedimenting due to gravity. of gravity.

Figure 3.3: A fluid suspension of 3204 particles sedimenting under the influence of gravity. Gravity acts in the z direction. Colors
indicate fluid field velocity in the z direction.

The file gravity100000.udf specifies the following parameter values: number of particles N p = 100000, particle
diameter D = 4, interface thickness ξ = 2. With these settings, the fractional volume occupied by the particles is
φ = 0.10. Figure 3.4 shows a snapshot of the 100000-particle suspension sedimenting due to gravity.

3.3.2

Particle diffusion

The input UDF file for this example is repulsive.udf in the Examples/03/ folder. This UDF file specifies asimulation of diffusion phenomena for a system of Brownian particles with repulsive interparticle forces but no otherinteractions.
$ ../../kapsel -Irepulsive.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
For this example (Fig. 3.5), we choose a 128 × 128 × 128 computational mesh. The fractional volume occupied by
the particles in this case is φ = 0.064. We specify the following values for the simulation parameters: number ofparticles N p = 500, particle diameter D = 8, interface thickness ξ = 2, temperature kB T = 5. The solvent has a density
ρ = 1 and a viscosity coefficient η = 1. Figure 3.5 shows a snapshot of this suspension of particles diffusing under the
influence of repulsive interparticle interactions.
2 These folders do not need to be created if you do not need AVS data; set the AVS switch to OFF in that case. This applies for all other sample
calculations as well.

CHAPTER 3. SIMULATING PARTICLES DISPERSED IN NEWTONIAN FLUIDS

32

Figure 3.4: The same as in the previous figure, but for a suspension of 100000 particles. Colors are for visualization purposes only,
all particles are identical.

3.3.3

Particle aggregates

The input UDF files for this example are aggregate.udf, aggregate LE.udf, and aggregate LE 2.udf in the
Examples/05/ folder. aggregate.udf describes a simulation of aggregation and diffusion phenomena in a systemof particles with attractive interparticle interactions.
$ ../../kapsel -Iaggregate.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
For this example (Fig. 3.6), we choose a 128 × 128 × 128 computational mesh. The fractional volume occupied by
the particles in this case is φ = 0.064, and we specify the following values for the simulation parameters: number ofparticles N p = 500, particle diameter D = 8, interface thickness ξ = 2, temperature kB T = 5. The solvent has density
ρ = 1 and viscosity coefficient η = 1. Figure 3.6 shows aggregation and diffusion phenomena for this suspension ofparticles with attractive interparticle interactions.
The alternative UDF input files aggregate LE.udf and aggregate LE 2.udf describe the same simulation, but
now in the presence of shear flow.
$ ../../kapsel -Iaggregate_LE.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf$ ../../kapsel -Iaggregate_LE_2.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
For this example (Fig. 3.7), we use the same values of computational and simulation parameters used for the simulationof Fig. 3.6: 128 × 128 × 128 computational mesh, fractional volume occupation by particles φ = 0.064, number ofparticles N p = 500, particle diameter D = 8, interface thickness ξ = 2, temperature kB T = 5. The solvent has a density
ρ = 1 and a viscosity coefficient η = 1.
The files aggregate LE.udf and aggregate LE 2.udf specify shear-flow velocities of γ̇ = 0.005 and γ̇ = 0.05,
respectively. Simulation snapshots for these two cases are shown in the left and right panels of Figure3.7.

3.3.4

Motion of a particle chain

The input UDF files for this example are s000 t010 free.udf, s001 t010 free.udf, and s001 t000 free.udf
in the Examples/06 folder.
s000 t010 free.udf describes a simulation of a single particle chain in a static fluid. s001 t010 free.udf
and s001 t000 free.udf describe similar calculations, but now in the presence of zig-zag shear flows.
$ ../../kapsel -Is000_t010_free.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf$ ../../kapsel -Is010_t010_free.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf

CHAPTER 3. SIMULATING PARTICLES DISPERSED IN NEWTONIAN FLUIDS

33

Figure 3.5: Snapshot of a suspension of particles with repulsive interparticle forces but no other interactions.

$ ../../kapsel -Is010_t000_free.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
For this example (Fig. 3.8), we choose a 32 × 64 × 16 computational mesh. The fractional volume occupied by
the particles in this case is φ = 0.005. We specify the following values for the simulation parameters: numberof particles per chain N p = 5, particle diameter D = 4, interface thickness ξ = 2. s000 t010 free.udf ands001 t010 free.udf specify a temperature of kB T = 1.0. s001 t000 free.udf specifies a temperature of kB T =0.0. s001 t010 free.udf and s001 t000 free.udf specify a shear velocity of s = 0.01.
Fig. 3.8 shows snapshots from the simulations described by these files. The left snapshot is for a particle chain in astatic fluid at temperature kB T = 1.0. The center snapshot is for a particle chain in a zig-zag shear flow at temperature
kB T = 1.0. The right snapshot is for a particle chain in a zig-zag shear flow at temperature kB T = 0.0.

3.3.5

Motion of an arbitrarily-shaped rigid particle

The input UDF files for this example are N30k3.0p52.0 vy.udf, N30k3.0p52.0 omegay.udf, N30k3.0p52.0 gravity.udf,
helixes gravity.udf, and helixes shear0005.udf in the Examples/08/ folder.
N30k3.0p52.0 vy.udf describes a simulation of a helical particle chain exhibiting translational motion at fixedvelocity.
$ ../../kapsel -IN30k3.0p52.0_vy.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
For this example we choose a 256 × 256 × 256 computational mesh. The fractional volume occupied by the particles
in this case is φ = 0.00048. We specify the following values for the simulation parameters: number of particles per
chain N p = 30, particle diameter D = 4, interface thickness ξ = 1, velocity V = −0.005
N30k3.0p52.0 omegay.udf describes a simulation of a helical particle chain exhibiting rotational motion at fixedangular velocity.
$ ../../kapsel -IN30k3.0p52.0_omegay.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
For this example we choose a 256 × 256 × 256 computational mesh. The fractional volume occupied by the particles
in this case is φ = 0.00048. We specify the following values for the simulation parameters: number of particles per
chain N p = 30, particle diameter D = 4, interface thickness ξ = 1, angular velocity Ω = 0.00041.
N30k3.0p52.0 gravity.udf describes a simulation of a helical particle chain sedimenting under the influenceof a fixed external force (gravity).

CHAPTER 3. SIMULATING PARTICLES DISPERSED IN NEWTONIAN FLUIDS

34

Figure 3.6: Snapshot illustrating aggregation and diffusion phenomena in a suspension of particles with attractive interparticle
interactions.

$ ../../kapsel -IN30k3.0p52.0_gravity.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
For this example we choose a 256 × 256 × 256 computational mesh. The fractional volume occupied by the particles
in this case is φ = 0.00048. We specify the following values for the simulation parameters: number of particles per
chain N p = 30, particle diameter D = 4, interface thickness ξ = 1, gravity g = −0.01.
helixes gravity.udf describes a simulation of three helical particle chains sedimenting due to gravity.
$ ../../kapsel -Ihelixes_gravity.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
For this example (Fig. 3.9) we choose a 256 × 256 × 256 computational mesh. The fractional volume occupied by
the particles in this case is φ = 0.0016. We specify the following values for the simulation parameters: number ofparticles per chain N p = 100, particle diameter D = 4, interface thickness ξ = 1, gravity g = −0.01. Figure 3.9 showsa snapshot of the three sedimenting helical particle chains.
helixes shear0005.udf describes a simulation of three helical particle chains in a shear flow.
$ ../../kapsel -Ihelixes_shear0005.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
For this example we choose a 256 × 256 × 256 computational mesh. The fractional volume occupied by the particles
in this case is φ = 0.0016. We specify the following values for simulation parameters: number of particles per chainN p = 100, particle diameter D = 4, interface thickness ξ = 1, shear velocity γ̇ = 0.005.

3.3.6

Sedimentation of arbitrarily-shaped rigid particles

The input UDF file for this example is input rigid free1.udf in the Examples/10/ folder. This file describes asimulation of a body falling through a fluid at Reynolds number Re = 800.
$ ../../kapsel -Iinput_rigid_free1.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
For this example (Fig. 3.10) we choose a 64 × 128 × 64 computational mesh. The density of the falling body is 1.5
times the density of the fluid, and the gravity strength is g = −50. Figure 3.10 shows a snapshot of the falling body atRe ≃ 800.

CHAPTER 3. SIMULATING PARTICLES DISPERSED IN NEWTONIAN FLUIDS

35

Figure 3.7: Snapshots from simulations like that of Fig. 3.6, but now in the presence of weak (left) or strong (right) shear flow.

Figure 3.8: Left: Snapshot of a single flexible particle chain in a static fluid at temperature kB T = 1.0. Center: Snapshot of aparticle chain in a zig-zag shear flow at temperature kB T = 1.0. Right: Snapshot of a particle chain in a zig-zag shear
flow at temperature kB T = 0.0.

CHAPTER 3. SIMULATING PARTICLES DISPERSED IN NEWTONIAN FLUIDS

36

Figure 3.9: Snapshot of three helical particle chains sedimenting due to gravity.

Figure 3.10: Snapshot of a body falling through a fluid at Reynolds number Re = 800. Light blue shading indicates fluid velocity
magnitude.

Chapter 4

Simulating disperse particle systems inshear flows
4.1

Theoretical background and basic equations

Colloidal suspensions—collections of colloidal particles dispersed in liquids (solvents)—are important and widelystudied systems in chemical engineering, mechanical engineering, physics, and many other fields of science andengineering. Applications of these systems span a wide range of products—from foodstuffs to paints, pigments,
cosmetics, slurries, and more—and a proper understanding of their fluid-mechanical behavior is important for manyindustrial production and treatment processes. However, the dynamical behavior of microparticle suspensions canbe extremely complicated due the effects of interparticle interactions, thermal fluctuations, and other factors, andaccurately modeling the fluid-mechanical behavior of these systems is a difficult challenge. Microparticle suspensions
exhibit a wide variety of behavior in various circumstances. For example, studies of transport phenomena in thesesystems reveal that their viscosities and other fluid-mechanical properties depend strongly on factors such as dispersed
particle concentrations, fluid-flow strength, and the nature of particle-fluid interfaces. The task of elucidating therelationship between such macroscopic fluid-mechanical properties and their microscopic mechanisms ranks among
the central themes of basic research on colloidal suspensions.
Despite the great practical importance of microparticle-suspension rheology, constructing a general theoreticalframework for these systems remains an extremely complex task (discussed, for example, in textbooks [15] and relatively recent (Japanese-language) reviews [16, 17]). For a solvent of viscosity η, the viscosity ηapp of a suspensionof particles in regions of low particle concentration (characterized by the particle volume fraction φ), is known on
theoretical grounds to vary linearly with φ in the form
ηapp
η

=

5
1 + φ.
2(4.1)

This is Einstein’s famous viscosity relation (see, for example, [15]). However, Einstein’s formula holds only for dilute
suspensions with particle volume fractions of φ ≪ 10% or so. At higher particle concentrations, one must accountfor the effects of interparticle interactions and the associated changes in particle structure, as well as the emergenceof phenomena such as glass transitions and crystallization; these complications pose grave challenges for the task oftheoretically predicting particle-suspension rheologies from basic physical laws. Nonetheless, despite the difficultyof this task, its practical importance has motivated many experimental studies, whose findings have been organizedin the form of various proposed empirical and quasi-empirical formulas for the viscosity of suspensions (such as theDoughherty-Krieger equation [18]; see also Refs. [17, 19] and references within).
As the particle concentration increases, phenomena begin to be affected by a variety of factors beyond the colloid
volume fraction, and the goal of establishing a unified theory must be set aside in favor of separately analyzing
individual phenomena. The SP method allows for direct computational simulation of solvent-induced multibody
interactions among colloids, enabling studies of the rheology of many types of suspensions. The key strengths ofthis direct numerical approach to rheology analysis are its ability to account for hydrodynamic interactions between
particles, thermal fluctuations of particles, and shear flows. KAPSEL allows simulation and rheology measurement(steady-state flow and dynamic viscoelasticity) for monodisperse systems in Newtonian fluids, including both simple
and oscillatory shear flows.
KAPSEL uses Lees-Edwards periodic boundary conditions to represent uniform shear flows. The advantages ofthis approach include (a) the accurate imposition of shear velocities without the use of walls or other specialized
37

CHAPTER 4. SIMULATING DISPERSE PARTICLE SYSTEMS IN SHEAR FLOWS

38

boundaries, and (b) the absence of singular points (kinks) at which the velocity gradient changes sign within the system, as in zigzag flows. With Lees-Edwards periodic boundary conditions, one solves the constitutive equations in
an oblique coordinate system in which the shear flow induces advection of lattice points with the passage of time.
The coordinate transformations used to represent shear flows with Lees-Edwards periodic boundary conditions are
depicted schematically in Figure 4.1. If there is no shear flow at time t = t0 , the collection of solid particles dispersedin the solvent is scattered over a rectangular lattice [Figure 4.1(a)]. If a shear flow develops at later times (t > t0 ),both particles and lattice points are advected by the flow; however, although the computational mesh is deformedby the shear flow, the shapes of the particles remain unchanged. More specifically, the computational mesh has two
representations—one in a fixed coordinate system [Figure 4.1(b)] and one in an oblique coordinate system distortedby the shear flow [Figure 4.1(c)]—where the shapes of the particles are distorted in the latter. The algorithm as wehave described it resembles the SLLOD algorithm used to model uniform dispersion-free flows in non-equilibrium
molecular-dynamics (MD) simulations, but implementing the algorithm in lattice-based systems requires extra caution [20–22]. Within the oblique coordinate system, which deforms as time passes due to the shear flow, we denote thevelocity by ξ = u − U . Then, using basis vectors Ê1 = e x + γ(t)ey , Ê2 = e2 , Ê3 = e3 , the appropriate contravariant
form of the Navier-Stokes equations with respect to ξ readsρ ∂tˆ + ξ̂ j ∇ˆ j ξ̂i = ∇ˆ j σ̂ ji + ϕ̂ fˆi − 2γ̇(tˆ)ξ̂2 δi,1
(4.2)
∇ˆ i ûi = ∇ˆ i ξ̂i = 0(4.3)

where hats (ˆ·) indicate tensor components in the oblique coordinate system. After solving these equations for the flow,
we compute particle-fluid interactions in the fixed (laboratory) coordinate system to solve for the particle dynamics.
Under shear flow, the spatial average of the instantaneous stress on the entire disperse particle system may be expressed

Eq.7 in [49]Figure 4.1: Schematic depiction of the coordinate transformations used to represent shear flows with Lees-Edwards periodic
boundary conditions. (a) At time t = t0 , solid particles are scattered over a rectangular lattice. The emergence of
shear flow at later times (t > t0 ) induces advection of both the particles and the solvent (lattice points); however, only
the computational mesh deforms, with particle shapes remaining unchanged. In other words, while the computationalmesh deforms in the fixed (laboratory) coordinate system (b), in the oblique coordinate system (c) particles deform in
response to the shear flow. The transformation between coordinate systems (b) and (c) may be found in equation (7) ofRef. [22]. Reproduced from J. Chem. Phys 134, 064110[22], Copyright 2011, with the permission of AIP Publishing.

in terms of the particle binding force ρϕf p in the form
Z
h
i
1
Σ=
dx σ − xρϕf p + xu · ∇(ρu) .
V

(4.4)

Subtracting from this the (spatially-averaged) contribution of the host fluid leaves just the particle contribution to thestress, s = Σ − ⟨σ⟩.

CHAPTER 4. SIMULATING DISPERSE PARTICLE SYSTEMS IN SHEAR FLOWS

4.1.1

39

Steady-state shear flow

For a steady-state shear flow, the apparent and intrinsic viscosities ηapp and [η] of the suspension are given by⟨Σ xy ⟩ ⟨σ xy ⟩ ⟨s xy ⟩⟨s xy ⟩
=
+
=η+
γ̇
γ̇
γ̇
γ̇
ηapp − η ⟨s xy ⟩
[η] =
=
.
ηφ
ηφγ̇

ηapp =

4.1.2

(4.5)
(4.6)

Oscillatory shear flow

Introducing an oscillatory shear flow allows measurement of the dynamic viscoelasticity of the suspension. In
KAPSEL, a time-varying shear flow in a system is formed by using, as a control parameter, the shear velocity γ̇,given byγ̇(t) = γ̇0 cos(ωt)(4.7)

where γ̇0 is the shear-velocity amplitude. The maximum strain amplitude γ0 may be evaluated in the form
γ0 =Z π/2ω
0γ̇(s)ds =

γ̇0
.
ω

(4.8)

Small strains correspond to the regime of linear viscoelasticity, while for larger strains we are in the nonlinear viscoelasticity regime.
For an oscillatory shear flow γ̇(t), the time-dependent stress Σ xy (t) in the suspension may, within the linear-response
regime, be expressed in the form
Σ xy (t) = σ0 cos(ωt − δ)

(4.9)

where σ0 is the stress amplitude and δ is the phase shift between the shear velocity γ̇ and the shear stress Σ xy . In
this case, the storage modulus G′ and loss modulus G′′ , which are functions of the dynamic viscoelasticity, may be
expressed in the formG′ (ω) =ωσ0 sin δωσ0 cos δ, G′′ (ω) =
.
γ̇0
γ̇0

(4.10)

G′ and G′′ respectively characterize the elastic and viscous behavior of the suspension. KAPSEL takes the shear
velocity γ̇(t) of equation (4.7) as an input and computes the shear stress Σ xy (t) of equation (4.9) as the corresponding
response (output). The moduli G′ (ω) and G′′ (ω) may be determined by conducting simulations for various values ofω and fitting the temporal evolution of variables to the above expressions. The derivations of these relations may befound in textbooks on rheology.

CHAPTER 4. SIMULATING DISPERSE PARTICLE SYSTEMS IN SHEAR FLOWS

4.2

Input UDF files

4.2.1

Fluid settings

40

The motion of colloidal particles under shear flow in a Newtonian fluid may be simulated by setting constitutive eq
to Shear Navier Stokes Lees Edwards. The settings that may be configured via input UDF files for this case aredescribed below.
➤constitutive eq.Shear Navier Stokes: Solvent properties
constitutive eq.Shear Navier Stokes.DX
The lattice spacing ∆, which defines the length unit.
constitutive eq.Shear Navier Stokes.RHO
Solvent density.
constitutive eq.Shear Navier Stokes.ETA
Solvent viscosity.
constitutive eq.Shear Navier Stokes.kBT
Particle temperature.
constitutive eq.Shear Navier Stokes.alpha v
Correction term for particle temperature associated with translational motion.1
constitutive eq.Shear Navier Stokes.alpha oCorrection term for particle temperature associated with rotational motion.
constitutive eq.Shear Navier Stokes.External field.type
Set to DC (steady-state shear flow) or AC (oscillatory shear flow).
constitutive eq.Shear Navier Stokes.External field.DC.Shear rate
Shear velocity γ̇.
constitutive eq.Shear Navier Stokes.External field.AC.Shear rate
Oscillatory shear velocity amplitude γ̇0 .
constitutive eq.Shear Navier Stokes.External field.AC.Frequency
Shear velocity frequency ω.

4.2.2

Object (particle) settings

The object type.type setting specifies the type of particle. The possible values are spherical particle forspherical particles, chain for flexible chains, or rigid for rigid bodies.

1 In principle, particle temperatures may be estimated from the diffusive motion of particles at thermal equilibrium; however, this approach may
yield incorrect temperatures due to numerical factors such as interface thickness. In such cases, alpha v may be used to correct temperature errors.
If the measured temperature is (say) too low by a factor of 1.2, then setting alpha v to 1.2 will reproduce the correct temperature (and similarly withalpha o for the rotational temperature). For a particle radius of 4 or 5 and interface thickness ξ = 2, the choice alpha v=alpha o=1 reproduces
correct temperatures.

CHAPTER 4. SIMULATING DISPERSE PARTICLE SYSTEMS IN SHEAR FLOWS

4.3

Computational examples

4.3.1

Rheology of particle suspensions under steady-state shear flow

41

The input UDF files for this example are colloid LE v31 t01.udf and colloid LE v55 t50.udf in the Examples/07/
folder.
$ ../../kapsel -Icolloid_LE_v31_t01.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf$ ../../kapsel -Icolloid_LE_v55_t50.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
For this example (Fig. 4.2), we choose a 64×64×64 computational mesh. We specify the following values for the simulation parameters: shear velocity γ̇ = 0.001, particle radius a = 4, interface thickness ξ = 2. colloid LE v31 t01.udf
specifies N p = 300 particles at temperature kB T = 0.1; for this case the fractional volume occupied by the particles
is φ = 0.307. colloid LE v55 t50.udf specifies N p = 540 particles at temperature kB T = 5.0; for this case thefractional volume occupied by the particles is φ = 0.552.

Figure 4.2: Flow snapshots for a particle suspension under steady-state shear flow using Lees-Edwards periodic boundary conditions. Left: N p = 300 particles, temperature kB T = 0.1. Right: N p = 540 particles, temperature kB T = 5.0.

Setting constitutive eq to Shear Navier Stokes Lees Edwards and External field.type to DC resultsin the following output written to standard error:#1:time 2:shear_rate 3:degree_oblique 4:shear_strain_temporal 5:lj_dev_stress_temporal 6:
0.0748173 -0.00588609 -0.000440381 ...
..
.
The significance of this output stream is as follows.
•

1:time. . . Elapsed time•

2:shear rate. . . Shear velocity•

3:degree oblique. . . Strain in oblique coordinate system2•

4:shear strain temporal. . . Applied strain (γ̇t)

2 KAPSEL represents uniform shear flows by performing calculations in an oblique coordinate system that distorts in the shear direction withthe passage of time. For computational stability, when the strain in the oblique coordinate system reaches a certain threshold, its value is remapped
into the orthogonal coordinate system. The value of 3:degree oblique is the strain with this remapping taken into account. Note that this differsfrom the applied strain (γ̇t) on the system. For further details, see Ref. [23].

CHAPTER 4. SIMULATING DISPERSE PARTICLE SYSTEMS IN SHEAR FLOWS•

5:lj dev stress temporal. . . Shear stress due to interparticle potential•

6:shear stress temporal old. . . Shear stress due to particle-fluid interactions•

7:shear stress temporal new. . . Shear stress due to particle-fluid interactions3•

8:reynolds stress. . . Reynolds shear stress•

9:fluid stress. . . Shear stress due to bulk fluid•

10:interfacial stress. . . Shear stress due to fluid-fluid interface•

11:apparent stress. . . Shear stress on suspension4•

12:viscosity. . . Viscosity of suspension5

42

Fig. 4.3 plots the temporal evolution of the stress on the suspension, while Fig. 4.4 plots the relationship between
viscosity and strain.

Figure 4.3: Temporal evolution of the stress for a suspension in steady-state shear flow. The x and y axes represent columns 1 and
11 in the output stream.

4.3.2

Rheology of particle suspensions under oscillatory shear flow

The input UDF files for this example are colloid LE v31 t01 AC.udf and colloid LE v55 t50 AC.udf in the
Examples/07 folder.
$ ../../kapsel -Icolloid_LE_v31_t01_AC.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf$ ../../kapsel -Icolloid_LE_v55_t50_AC.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
For this example (Fig. 4.5), we choose a 64×64×64 computational mesh. We specify the following values for the simulation parameters: shear velocity γ̇ = 0.001, particle radius a = 4, interface thickness ξ = 2. colloid LE v31 t01 AC.udf
specifies N p = 300 particles at temperature kB T = 0.1; for this case the fractional volume occupied by the particles
is φ = 0.307. colloid LE v55 t50 AC.udf specifies N p = 540 particles at temperature kB T = 5.0; for this case thefractional volume occupied by the particles is φ = 0.552.
Setting constitutive eq to Shear Navier Stokes Lees Edwards and External field.type to AC resultsin the following output written to standard error:
3 KAPSEL implements two methods (old and new) for computing the shear stress due to particle-fluid interactions. In most cases the new
method is the best choice.
4 This is the sum of 5:lj dev stress temporal, 7:shear stress temporal new, 8:reynolds stress, 9:fluid stress, and10:interfacial stress.
5 The result of dividing 11:apparent stress by 2:shear rate

CHAPTER 4. SIMULATING DISPERSE PARTICLE SYSTEMS IN SHEAR FLOWS

43

Figure 4.4: Viscosity-strain relationship for a suspension under steady-state shear flow. The x and y axes represent columns 4 and
12 in the output stream.

#1:time 2:shear_rate 3:degree_oblique 4:shear_strain_temporal 5:lj_dev_stress_temporal 6:
0.0748173 -0.00588609 -0.000440381 ...
The significance of the various columns in this output stream is the same as that described above for the steady-statecase.
Fig. 4.5 plots the temporal evolution of the stress on the suspension.

CHAPTER 4. SIMULATING DISPERSE PARTICLE SYSTEMS IN SHEAR FLOWS

44

Figure 4.5: Temporal evolution of the stress for a suspension in oscillatory shear flow. The x and y axes represent columns 1 and
11 in the output stream. Upper plot: particle volume fraction φ = 0.307, temperature kB T = 0.1. Lower plot: particle
volume fraction φ = 0.552, temperature kB T = 5.0.

Chapter 5

Simulating charged particles dispersed inelectrolytic solutions
5.1

Electrophoresis of charged colloidal particles

When colloidal particles are dispersed in water or another solvent with extremely high dielectric permittivity, ions are
emitted from dissociable groups on the colloid surface, causing the particle surface to acquire a net electric charge.
The emitted ions are electrostatically attracted to the particle surface, but also undergo diffusion due to thermal fluctuations, resulting in the formation of a cloud-like ionic atmosphere around the colloid known as an electric double layer.
The properties of disperse colloid systems at thermal equilibrium are described by the Poisson-Boltzmann equation;
linearizing this equation yields the Debye–Hückel approximation, from which valuable insights may be obtained. In
contrast, in systems exhibiting electrophoresis—or any of the other related effects known collectively as interfacial
electrokinetic phenomena—the behavior of particles and ion distributions is governed by a competition between hydrodynamic and electrostatic interactions. As a result of this competition, the ion distribution is unable to keep up withthe motion of the particles, causing the electric double layer to distort from a spherically symmetric distribution andyield configurations that differ from the equilibrium states. To analyze such complicated situations theoretically, there
is no choice but to introduce simple approximations. In fact, until recently, even numerical simulations were rarely, if
ever, capable of reproducing these phenomena correctly.
In applications of the SP method to simulations of electrokinetic phenomena, the colloids are treated as particles—
as in the other types of simulation—but the solvent and the ions are treated as a continuous medium and described, after
a coarse-graining procedure, by a density field. To make these calculations feasible, we have designed a computational
formalism that uses SP functions to ensure smooth interfaces and self-consistently solves a coupled set of equations
for three interacting degrees of freedom: the equations of motion for the colloids, the advection-diffusion equation for
the ion distribution, and the Navier-Stokes equations for the solvent velocity field [5, 6, 24, 25]. This approach enabled
the world’s first successful quantitative simulations of electrophoresis [25]. In the following section we describe thebasic equations considered internally by KAPSEL and briefly survey the underlying theoretical background.

5.2

Basic equations

We begin by surveying the basic equations of electrohydrodynamics needed for these simulations, following Refs. [6,25]. Consider N spherical colloidal particles of radius a dispersed in an electrolytic solvent. We assume that thedielectric permittivity ϵ of the solvent is spatially uniform, including in the interior of the colloids. In the SP method,
the interface between a particle and the solvent is represented by a smooth function ϕ(r) ∈ [0, 1] of finite width ξ.
Working on a fixed orthogonal lattice, the SP function takes the values ϕ = 1 in the interior of colloids and ϕ = 0
in solvent regions outside of colloids; the interface region is defined by the condition 0 < ϕ < 1. Other approaches
that use ϕ(r) in an orthogonal lattice include the methods of Tanaka-Araki [26] and Kajishima et al. [27]. The use ofinterface functions yields an overwhelming advantage in computational efficiency compared to finite-element methods
or other approaches based on non-structured lattices. We assume that colloid surfaces are uniformly charged, with a
total charge of Ze per particle. In typical calculations involving continuous media, particle surface charge distributions
are represented by delta functions; in finite-element and similar approaches this requires the use of an appropriate
boundary-conforming computational mesh, significantly degrading computational efficiency. In the SP method, bycontrast, the particle surface charge distribution eq(r) may itself be represented as a smooth function involving first45

CHAPTER 5. SIMULATING CHARGED PARTICLES DISPERSED IN ELECTROLYTIC SOLUTIONS

46

derivatives of the SP function ϕ:
Ze|∇ϕ(r)|
.
(5.1)
4πa2
Just as the SP function ϕ(r) is designed to approach a step function as ξ → 0, the charge density tends to a delta
function as ξ → 0.
eq(r) =

5.2.1

The advection-diffusion equation

The density distribution Cα for an ion of species α with charge zα may be defined at all points in the computational
domain by writingCα (r, t) = (1 − ϕ(r, t))Cα∗ (r, t).
(5.2)
The factor (1 − ϕ) selects only regions in which ions are present. Cα∗ (r, t) is an auxiliary variable, introduced forcomputational convenience, and defined to be smooth everywhere in the computational domain. In the interior of
colloids (ϕ = 1), the quantity Cα∗ has no physical meaning. The total charge distribution, including the colloid surface
charge, isXρe (r) = ezαCα (r) + eq(r).
(5.3)α

We constrain the initial distribution by requiring the electro-neutrality condition to be satisfied,
temporal evolution of the auxiliary ion density Cα∗ is governed by the advection-diffusion equation:Rρe dr = 0. The∂t Cα∗ = −∇ · Cα∗ v + Γα ∇ · (Cα∗ ∇µα ).

(5.4)

The two terms on the right-hand side here describe advection by the solvent velocity field v and diffusion due to thegradientof the chemical potential µα , respectively. Because Cα∗ obeys the advection-diffusion equation, the quantity
R∗
drCα is conserved.
To prevent ions from penetrating into the interior of colloids, we impose the constraint that the ion diffusion flow
velocity at particle-solvent interfaces has a vanishing component in the direction normal to the interface, i.e., thatn · ∇µα = 0 [6, 25]. Here, n is the outward-pointing normal vector at the colloid surface; in terms of the SP function
we may write n = −∇ϕ/|∇ϕ|. Also, Γα is the Onsager transport coefficient for ions of species α, related to the ion’s
friction coefficient and diffusion coefficient according to fα = 1/Γα , Dα = kB T Γα . We take the ion chemical potential
to be [28]µα = kB T ln Cα∗ + zα e(Ψ − E · r).
(5.5)
Here E represents an external electric field, and the electrostatic potential Ψ(r) solves the Poisson equation:ϵ∇2 Ψ = −ρe .

(5.6)

In an equilibrium state, ions with this chemical potential obey the Poisson-Boltzmann distribution.

5.2.2

The Navier-Stokes equations

The solvent velocity field u is described by the Navier-Stokes equations for an incompressible flow (∇ · u = 0):ρ(∂t + u · ∇)u = −∇p + η∇2 u − ρe (∇Ψ − E) + ϕ f p(5.7)

where ρ, η, and p are the density, viscosity coefficient, and pressure of the solvent. Note that the fluid experiences anelectrostatic force −ρe (∇Ψ − E). The external-force term ϕ f p represents a constraint force that preserves the rigidityof the particles. In other words, we consider adhesive boundary conditions at particle surfaces due to ϕ f p . Further
discussion of this point may be found in Refs. [5, 6].

5.2.3

Equations of motion

If the ith colloidal particle has mass M p , its time-dependent position and velocity {Ri , V i } evolve in time according to
the equations of motion:Ṙi
M p V̇ i

=

Vi,

(5.8)

=

FiH + Fother
i

(5.9)

Here FiH is the force received from the fluid, representing the balance of momentum between solid and fluid [6].
Also, Fotherrepresents forces due to interparticle potentials such as Lennard-Jones potentials. Although omitted here,isimilar considerations hold for the rotational motion of particles [6]. This completes our survey of the basic equations
governing the simulations.

CHAPTER 5. SIMULATING CHARGED PARTICLES DISPERSED IN ELECTROLYTIC SOLUTIONS

5.3

47

Properties of electric double layers

The basic tool for quantitative analysis of electric double-layer structures is the Poisson-Boltzmann equation.

5.3.1

The Poisson-Boltzmann equations

In the absence of external electric fields (E = 0), the equilibrium ion distribution may be obtained from equation (5.5).
Assuming the chemical potential is uniform (µα = const) the equilibrium ion distribution takes the form!
zα eΨ(r)∗Cα (r) = C̄α exp −.
(5.10)kB T
This is the Boltzmann distribution in the presence of an electrostatic potential Ψ. Combining this with equation (5.6)
yields the Poisson-Boltzmann equation.

5.3.2

The Debye-Hückel approximation and the Debye screening length

Consider a single spherical colloidal particle in a z : z symmetric electrolytic solvent. The Poisson-Boltzmann equation
reads [15, 29]!
2zeC̄
zeΨ(r)∇2 Ψ(r) =sinh.
(5.11)ϵkB T
At infinite distance from the system we impose the boundary conditions Ψ| r=∞ = 0 and C ∗ | r=∞ = C̄. The boundary
condition at the particle surface is such that the surface charge density remain constant and equal to σe = Ze/4πa2 ,thus requiring
σe∇Ψ|surface = − .
(5.12)ϵAssuming zeΨ/kB T ≪ 1 and linearizing equation (5.11) yields the Debye–Hückel approximation:∇2 Ψ(r) =2z2 e2C̄Ψ = κ2 Ψ,
kB T ϵwhereκ−1 = p18πλB z2C̄

(5.13)(5.14)

is a constant with dimensions of length known as the Debye screening length. Here λB = e2 /4πkB T ϵ is the Bjerrum
length. For a typical electrolyte we have1
κ−1 = q.
(5.15)
P
4πλB α z2αC̄αAs the system is spherically symmetric, we need only consider variations with respect to the radial coordinate r = |r|,yielding the simplified equation:d2 Ψ 2 dΨ+= κ2 Ψ
(5.16)
r drdr2
whose general solution takes the form of a Yukawa potential:aΨ(r) = Ψ0 exp[−κ(r − a)].
r

(5.17)The electrostatic force due to a charged colloid is screened beyond distances on the order of κ−1 . The Debye screening
length κ−1 may be thought of as the distance at which the Coulomb force attracting oppositely-charged ions to theparticle surface is balanced by the tendency of thermal diffusive motion to prevent ions from localizing. Thus κ−1
may be interpreted as the thickness of the electric double layer. At high temperatures the thermal energy kB T is largePand κ−1 is large. On the other hand, the larger the ion strength ( α z2αC̄α /2) the stronger the effect of ionic screening−1
and the shorter the screening length κ . As an example, consider the bulk salt concentration C̄ in a z : z symmetricelectrolytic solvent. Assuming the medium is water at 25 ◦ C, the Bjerrum length is λB = 0.72nm, and inserting actual
numerical values in equation (5.14) yields
0.3
κ−1 = √ (nm).
(5.18)z C̄
For z = 1, C̄ = 0.1M we have κ−1 = 1nm, while for z = 1, C̄ = 0.001M we have κ−1 = 10nm.

CHAPTER 5. SIMULATING CHARGED PARTICLES DISPERSED IN ELECTROLYTIC SOLUTIONS

48

The surface potential Ψ(r = a) = Ψ0 is determined by the boundary condition dΨ/dr(r = a) = −σe/ϵ. Ψ is related
to the surface charge density σe by
σe = ϵκΨ0 (1 + (κa)−1 ).
(5.19)
Thus the surface potential increases linearly with increasing surface charge. However, the Debye-Hückel approximation ceases to be valid in regions of high surface potential, and in practice the surface potential eventually tends toincrease by smaller and smaller increments in response to increasing surface charge. This behavior may be derived
by solving the Poisson-Boltzmann equation, but approximate formulae for 1:1 electrolytic solvents have also been
proposed by Loeb-Overbeek-Wiersema [15] and by Ohshima-Healy-White [29, 30].

5.4

Principles of electrophoresis

In the presence of an external electric field E, a colloidal particle with charge Ze feels an electrostatic force ZeE andbegins to move. The acceleration of the particle due to this electrostatic force is opposed by a viscous drag force fromthe fluid, and the balance between these two effects results in constant-velocity motion of the particle at steady-state
velocity V. Assuming the viscous drag force is the Stokes drag 6πηaV for a spherical particle of radius a, the force
balance reads
ZeE = 6πηaV(5.20)
yielding an electrophoretic mobility ofZeV=.
E 6πηa(5.21)

However, actual electrophoretic mobilities are smaller than this estimate. This is because equation (5.21) does notaccount for the effect of electrostatic forces on the ionic atmosphere surrounding the colloids, which induce motion inthat atmosphere. Colloids moving through the fluid drag their ionic atmosphere with them, reducing their velocities
correspondingly. Also, deviations from spherical symmetry in electric double layers give rise to additional forces thatmust be taken into account; this tends to make theoretical analysis prohibitively complicated, and instead a number ofsimplified models for analyzing electrophoretic mobilities have been considered [15, 29].

5.4.1

Smoluchowski’s equation

When the particle radius a is much larger than the electric double-layer thickness κ−1 , i.e. κa ≫ 1, we may apply
Smoluchowski’s equation. In this limit, the electric double layer is considered to be infinitesimally thin, so we ignore
the curvature of the particle surface and model it as a planar slab. Suppose an external electric field E x is applied
parallel to the slab (in the x direction). In a particle-fixed coordinate system, in which we ride atop a particle andobserve the motion of the surrounding fluid, the fluid velocity at infinite distance from the particle is −V; balancing
the viscous drag force and the electrostatic force on the particle yields
η∂2 v x X+
eCα E x = 0.
∂y2
α(5.22)

Poisson’s equation relates the ion distribution to the second derivative of the electrostatic potential, i.e.
η

∂2 v x∂2 Ψ= ϵ 2 Ex.
2∂y∂y

(5.23)

Integrating this expression, subject to the boundary conditions of vanishing velocity gradient, vanishing potential
gradient, and vanishing potential at spatial infinity, we find
η[v x (y) + V] − ϵEΨ(y) = 0.

(5.24)

Because the velocity at the particle surface (y = 0) is 0, we have
V
ϵζ
= .
E
η

(5.25)

This is Smoluchowski’s equation. The zeta potential ζ is conventionally defined as the potential at the slip surface, but
we have here replaced this with the potential Φ(0) at the particle surface.

CHAPTER 5. SIMULATING CHARGED PARTICLES DISPERSED IN ELECTROLYTIC SOLUTIONS

5.4.2

49

Hückel’s equation

In the opposite limit of Smoluchowski’s equation, we take the particle radius to be much smaller than the electric
double-layer thickness, i.e., κa ≪ 1. In this limit we have simply a point charge of strength Ze, and Hückel’s equation
applies. In equation (5.21), we take the potential at the particle surface to be the Coulomb potential, i.e.,Ze
4πϵa(5.26)

V
2 ϵζ
=
.
E 3 η

(5.27)

ζ=
whereupon

This is Hückel’s equation.

5.4.3

Henry’s equation and the O’Brien-White analysis

For general values of κa, a bridge between equations (5.25) and (5.27) is furnished by Henry’s equation:Vϵζ= f (κa)Eη(5.28)

where the Henry coefficient f (κa) is defined byf (κa) ==1 − 5 exp(κa)E7 (κa) + 2 exp(κa)E5 (κa)"#(κa)4 (κa)62 (κa)2 5(κa)3 (κa)4 (κa)5
+
−
−
+
+
−
exp(κa)E1 (κa).
3
24
72
144
144
12
144

(5.29)(5.30)

Here En (κa) is the nth exponential integral. Smoluchowski’s equation corresponds to f = 1(κa → ∞) while Hückel’s
equation corresponds to f = 2/3(κa → 0) (Fig. 5.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

Figure 5.1: The Henry coefficient f (κa).

Henry’s equation (5.28) depends linearly on the zeta potential and applies only for small values of this potential.
For larger zeta potentials it becomes necessary to account for the impact of deformations in the electric double-layer
potential (known as relaxation effects). O’Brien-White used numerical analysis to propose a relationship between thezeta potential and the electrophoretic mobility for arbitrary values of κa and ζ [31]. Later, an analytical formula validfor κa ≥ 10 was proposed by Ohshima-Healy-White [32].

CHAPTER 5. SIMULATING CHARGED PARTICLES DISPERSED IN ELECTROLYTIC SOLUTIONS

5.5

Input UDF files

5.5.1

Fluid settings

50

Set constitutive eq to Electrolyte to simulate charged-colloid motion or electrophoretic phenomena in thepresence of external electric fields.1 The settings that may be configured via input UDF files for this case are described
below.
➤constitutive eq.Electrolyte: Solvent and ion-distribution properties
constitutive eq.Electrolyte.DX
Lattice spacing ∆, which defines the length unit.
constitutive eq.Electrolyte.RHO
Solvent density.
constitutive eq.Electrolyte.ETA
Solvent viscosity.
constitutive eq.Electrolyte.kBT
Particle temperature.
constitutive eq.Electrolyte.alpha v
Correction term for thermal fluctuations (translational motion) of particles.
constitutive eq.Electrolyte.alpha oCorrection term for thermal fluctuations (rotational motion) of particles.
constitutive eq.Electrolyte.Dielectric cst
Dielectric permittivity of solvent.
constitutive eq.Electrolyte.Init profile
Set to Uniform or Poisson Boltzmann to select how the initial ion distribution is chosen.2
constitutive eq.Electrolyte.Add salt.type
Selects the number of ion species. Allowed values:
saltfreesalt

One ion species, with a charge opposite to that of the particle surface.
Two ion species (one positive, one negative).

constitutive eq.Electrolyte.Add salt.saltfree.Valency counterion
Counterion valency for Add salt = saltfree.
constitutive eq.Electrolyte.Add salt.saltfree.Onsager coeff counterion
Counterion Onsager transport coefficient for Add salt = saltfree.
constitutive eq.Electrolyte.Add salt.salt.Valency positive ion
Positive-ion valency for the case Add salt = salt.
1 Other possible settings for constitutive eq include Navier Stokes or Shear Navier Stokes to simulate the motion of colloidal particles
in Newtonian fluids in the absence or presence of shear flow. With these choices, thermal fluctuations of charged particles will be neglected.
2 Note that choosing Poisson Boltzmann may be rather time-consuming, particularly for large numbers of particles.

CHAPTER 5. SIMULATING CHARGED PARTICLES DISPERSED IN ELECTROLYTIC SOLUTIONS

51

constitutive eq.Electrolyte.Add salt.salt.Valency negative ion
Negative-ion valency for the case Add salt = salt.
constitutive eq.Electrolyte.Add salt.salt.Onsager coeff positive ion
Positive-ion Onsager transport coefficient for the case Add salt = salt.
constitutive eq.Electrolyte.Add salt.salt.Onsager coeff negative ion
Negative-ion Onsager transport coefficient for the case Add salt = salt.
constitutive eq.Electrolyte.Add salt.salt.Debye length
Debye screening length. If a value is specified, the salt concentration will be chosen to yield the specified
Debye length.
constitutive eq.Electrolyte.Electric field.type
Set to ON or OFF to enable or disable an external electric field.
constitutive eq.Electrolyte.Electric field.ON.type
For the case Electric field.type=ON, set to DC or AC to select a constant or oscillatory external electricfield.
constitutive eq.Electrolyte.Electric field.ON.DC.Ex
x-directed electric field strength for the case Electric field.ON.type=DC.
constitutive eq.Electrolyte.Electric field.ON.DC.Ey
y-directed electric field strength for the case Electric field.ON.type=DC.
constitutive eq.Electrolyte.Electric field.ON.DC.Ez
z-directed electric field strength for the case Electric field.ON.type=DC.
constitutive eq.Electrolyte.Electric field.ON.AC.Ex
x-directed electric field strength for the case Electric field.ON.type=AC.
constitutive eq.Electrolyte.Electric field.ON.AC.Ey
y-directed electric field strength for the case Electric field.ON.type=AC.
constitutive eq.Electrolyte.Electric field.ON.AC.Ez
z-directed electric field strength for the case Electric field.ON.type=AC.
constitutive eq.Electrolyte.Electric field.ON.AC.Frequency
Frequency of oscillatory external electric field.

5.5.2

Object (particle) settings

The object type.type setting specifies the type of particle. The possible values are spherical particle forspherical particles, or chain for flexible chains.

5.5.3

Choice of length and time units

The length unit is given by the lattice spacing ∆. The time unit τ0 is determined by the fluid density ρ, the fluid
viscosity coefficient η, and the lattice spacing ∆ according to τ0 = ρ∆2 /η.
•

For the Navier-Stokes equations, we choose a system of units in which ρ = η = ∆ = 1.

CHAPTER 5. SIMULATING CHARGED PARTICLES DISPERSED IN ELECTROLYTIC SOLUTIONS

•

•

•

52

Assuming the input UDF file specifies RHO= A, ETA= B, and DX= C, the maximum wavenumber kmax may
be determined from C and the upper bound on the timestep is determined from the momentum diffusion time2in the form T step = (A/B)/kmax. The actual simulation timestep is related to this bound through a multiplicative
scale factor, i.e., ∆t = factor × T steo .
Consider the correspondence between physical quantities in simulations and in reality. A grid spacing of 1 µm
is reasonable for the length scales we wish to consider. Assuming we use water (η = 1 × 10−3 Pa s, ρ =1 × 103 kg m−3 ) as a solvent, the time unit is then τ0 = 1 × 10−6 s.
2
2
For the case constitutive eq=Electrolyte, T step is set to the smallest value between (A/B)kmax
and (1/kB T Γα )/kmax,
and information regarding which of these two quantities is chosen is written to the console (standard error) during KAPSEL execution.

CHAPTER 5. SIMULATING CHARGED PARTICLES DISPERSED IN ELECTROLYTIC SOLUTIONS

5.6

Computational examples

5.6.1

Single-particle electrophoresis

53

The input UDF file for this example is colloid 1.udf in the folder Examples/01/. Check that the output-file folders./avs ch/ and ./avs ch/avs/ have been created on your system, then run the following command.
$ ../../kapsel -Icolloid_1.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
For this example, we choose a 64 × 64 × 64 computational mesh. We consider a 1:1 electrolytic solvent (viscosity
coefficient η = 1, density ρ = 1) and set the following values for the simulation parameters: particle radius a = 5,interface thickness ξ = 2, charge Z = −100, electric field strength E x = 0.1, Debye screening length κ−1 =10. Fig. 5.2
shows a snapshot from this simulation of single-particle electrophoresis, while Fig. 5.3 plots the time evolution of theparticle’s velocity; the velocity approaches a constant value determined by the balance of electrostatic and fluid dragforces on the particle. The plot in Fig. 5.3 is generated by the script plot.py.
0.011

0.007

0.004

Y

Z

0.000
X

Figure 5.2: Snapshot from simulation of a colloidal particle (red) exhibiting electrophoresis under the influence of a positive xdirected external electric field E. Dark/light colors indicate the charge density distribution; arrows indicate the solvent
velocity field.

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

Figure 5.3: Time evolution of electrophoretic velocity of a single colloidal particle.

CHAPTER 5. SIMULATING CHARGED PARTICLES DISPERSED IN ELECTROLYTIC SOLUTIONS

5.6.2

54

Many-particle electrophoresis

The input UDF file for this example is colloid 32.udf in the Examples/01/ folder. The simulation is similar tothat for the previous example, but now with 32 particles. Initial particle positions are chosen at random, and otherparameter values are chosen as in the previous example.
Check that the output-file folders ./avs ch/ and ./avs ch/avs/ have been created on your system, and then run
the following command.
$ ../../kapsel -Icolloid_32.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
Fig. 5.4 shows the results of this simulation of many-particle electrophoresis, while Fig. 5.5 plots the time evolutionof the average electrophoretic velocity for NP = 16, 32, 64 particles; the average velocity decreases as the number ofparticles increases.
0.04

0.03

0.01

Y

Z
X

0.00

Figure 5.4: Snapshot from a simulation of colloidal particles (green) exhibiting electrophoresis under the influence of a positive xdirected external electric field E. Dark/light colors indicate the charge density distribution; arrows indicate the solvent
velocity field.

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

Figure 5.5: Time evolution of the average electrophoretic velocity for multiparticle systems. The upper, middle, and lower plots
are for N p = 64, 32, 16 particles.

CHAPTER 5. SIMULATING CHARGED PARTICLES DISPERSED IN ELECTROLYTIC SOLUTIONS

5.6.3

55

Electrophoresis in a mixed particle system including both positively and negatively
charged particles

The input UDF file for this example is colloid p32m32.udf in the Examples/01/ folder. This file describes asimulation of electrophoresis in a mixed particle system containing 32 positively-charged colloids and 32 negativelycharged colloids. Other simulation parameters are configured as in the single-particle example above. Check thatthe output-file folders ./avs ch/ and ./avs ch/avs/ have been created on your system, and then run the following
command.
$ ../../kapsel -Icolloid_p32m32.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
Fig. 5.6 shows a snapshot from this simulation of electrophoresis in a system of 32 positively-charged and 32
negatively-charged colloids.
0.01

0.00

-0.00

Z
Y
X

-0.01

Figure 5.6: Snapshot from a simulation of a mixed particle system containing both positively-charged (red) and negatively-charged(blue) colloids exhibiting electrophoresis under the influence of a positive x-directed external electric field E. Positively
and negatively charged colloids move in the positive and negative x directions, respectively.

The input file colloid p5km5k.udf describes a similar simulation, but with the numbers of positively and negatively charged colloids increased to 5,000 each. Other parameter settings remain as in the single-particle example
above.
$ ../../kapsel -Icolloid_p5km5k.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
Fig. 5.7 shows a snapshot from this simulation of electrophoresis in a system of 5,000 positively-charged and 5,000
negatively-charged colloids.

5.6.4

Using AVS/Express for visualization

Set output.AVS to ON to enable generation of AVS output files. This will create an output file named avs charge.v,
which may be opened in the AVS/Express software package3 for visualization. In the Read field module, specify thefield information file data.fld to enable visualization of particles, solvent velocity fields, and charge distributions.
The plots in Figs. 5.2, 5.4, and 5.6 above were generated in this way using avs charge.v files.

5.6.5

Using Gourmet for visualization

Launch Gourmet, load the output file output.udf, use Gourmet’s Python panel to Load the Python script show field.py(distributed with KAPSEL), and select Run to open a new graphics window in which particles, solvent velocity fields,
and charge density distributions may be visualized.4 Click the Play button at the bottom of the screen to start an
3 http://www.avs.com/
4 Gourmet 2003 and earlier versions contain bugs that prevent this script from running properly. If you are affected by this problem, use thealternative visualization script particleshow.py. In this case the solvent velocity field and the charge density distribution will not be displayed;
only the particle visualization will be generated.

CHAPTER 5. SIMULATING CHARGED PARTICLES DISPERSED IN ELECTROLYTIC SOLUTIONS

56

Figure 5.7: Snapshot from a simulation of a mixed particle system containing both positively-charged (red) and negatively-charged(green) colloids exhibiting electrophoresis under the influence of a positive x-directed external electric field E. Positively and negatively charged colloids move in the positive and negative x directions, respectively.

animation. Particle visualizations are based on particle coordinate data written to output UDF files in the form ofdata records, so the property output.UDF must be set to ON. Also, visualizations of solvent velocity fields and charge
density distributions are produced by reading AVS-format binary output files, so output.AVS must be set to ON
and output.AVS.ON.FileType must be set to Binary. If you wish only to visualize the motion of particles, use
particleshow.py with output.UDF set to ON.

5.6.6

Using gnuplot to plot data

The gnuplot package5 may be used to plot time-series data for particle coordinates and velocities written to output
UDF files in the form of data records [33]. Launch Gourmet, load the output file output.udf, use Gourmet’s Python
panel to Load the Python script plot.py (distributed with KAPSEL), and select Run to create a new plot-generation
sheet from particle coordinates and velocities written as data records. To plot the data, in the “View” box selectTable instead of Tree; select GraphSheet from the list of variables to browse available GraphSheet data. Then, in
Gourmet’s Plot panel select Make followed by Plot to generate a plot like that in Fig. 5.3. If the plot is difficult tointerpret due to an excess of linestyles, simply use the editor in the Plot panel to edit the plot command, removing
any unneeded components. Further details on how to invoke gnuplot from Gourmet may be found in Section 3 ofRef. [33].

5 http://www.gnuplot.info/Chapter 6

Simulating particles dispersed in
two-component phase-separated fluids
6.1

Adding particles to two-component phase-separated fluids

Fluids consisting of mixtures of two immiscible components exhibit a rich variety of phase-separation structures, andtechniques for controlling these systems are extremely important for industrial applications. For example, the twocomponent systems known as emulsions—in which one component exists in the form of droplets, with sizes on theorder of microns or smaller, dispersed throughout a second fluid component—exploit special properties of low-mutualaffinity fluids in coexistence and have applications spanning a wide range of fields, from foodstuffs, to cosmetics, topetroleum recovery and beyond. However, emulsions are intrinsically thermodynamically unstable, and the stabilization and structural control of these systems pose significant practical challenges. In the early 1900s, Ramsden[34] andPickering[35] reported that emulsions could be stabilized by the addition of insoluble microparticles. This is due tothe tendency of microparticles to adhere to interfaces between two fluid components, where they play a role similar
to that of surfactants. Ever since the work of Ramsden and Pickering, an active school of research has studied various
microparticle species—including silica, carbon black, and metal complexes—as emulsion additives [36, 37]. Morerecently, the use of microparticles to stabilize fluid interfaces in two-component phase-separated mixtures with bicontinuous structures has led to the creation of new compound materials known as bijels that have been a focus of intenseinterest [38, 39]. Experience with microparticle additives for phase-separation structural control has shown that, in
addition to fluid-fluid interactions, fluid-particle interactions and the affinity between particles and fluid interfaces also
play key roles in these systems. KAPSEL is capable of taking all of these interactions into account when simulating
particle additives to two-component phase-separated fluids, enabling studies of the impact of particle additives onphase-separation structures and rheological properties.
For studies of particles in single-component fluids, the use of KAPSEL’s Smoothed Profile (SP) method is relatively straightforward: there is only one type of interface to consider, namely the particle-fluid interface. For twocomponent fluids, on the other hand, we encounter a new type of interface—the fluid-fluid interface—that must be
represented as well. Thus, for simulations of particles dispersed in two-component phase-separated fluids we introduce a new interface function ψ that distinguishes between the two fluid components that meet at fluid interfaces, andwe seek to describe how ψ evolves in time (Fig. 6.1). KAPSEL handles solid particles immersed in two-component
phase-separated fluids by combining the SP method with the well-known model H [40] approach to descriptive modeling of the dynamics of two-component phase-separated fluids. Section 6.2 describes how the SP method is extended
to handle two-component phase-separated fluids, while Section 6.3 describes various user-configurable simulation
parameters used by KAPSEL. Finally, as computational examples exploiting these new KAPSEL capabilities, Section 6.4 describes and presents results of simulations of solid particles dispersed in two-component phase-separated
fluids and Pickering emulsions.

6.2

Theoretical background and basic equations

6.2.1

Basic equations for two-component phase-separated fluids

In the model H theory [40], the equations of motion for a two-component phase-separated fluid are the Navier-Stokes(NS) equations augmented by a chemical-potential-gradient term:

57

CHAPTER 6. SIMULATING PARTICLES DISPERSED IN TWO-COMPONENT PHASE-SEPARATED FLUIDS58

Figure 6.1: Interface functions used to describe particles dispersed in a two-component phase-separated liquid. In addition to
the SP function ϕ separating particles from the surrounding fluid, we introduce an interface function ψ that separates
the two fluid components from each other. The thicknesses of the particle-fluid and fluid-fluid interfaces are denoted
respectively by ξ and ξψ . In this figure, ψ takes values of ±1 in bulk fluid regions, but this depends on the settingschosen for the double-well potential f (ψ).

∂u
1
1
ψ
ϕ
+ (u · ∇)u = − ∇p + ∇ · η(∇u + ∇ut ) − ∇µψ − ∇µϕ + ϕ f p .
∂t
ρ
ρ
ρ
ρ

(6.1)

Here the viscosity η is treated as a position-dependent physical quantity defined throughout the two-component fluid.
The superscript t indicates matrix transposition, and µψ = δF/δψ and µϕ = δF/δϕ are the locally-defined chemical
potentials with respect to ψ and ϕ, defined as functional derivatives of the Ginzburg-Landau (GL) free energy F. The
free energy F in KAPSEL is discussed in further detail below.
Next, the time evolution of the fluid-fluid interface function ψ is described by the Cahn-Hilliard (CH) equation:∂ψ= −∇ · J∂t(6.2)

where J is the mass flux, expressed here in a form that includes an advection term:
J = ψu − κ∇µψ .

(6.3)

Also, κ is the mobility of the two-component fluid. Note that the symbol κ, used in this section to denote mobility, has
no relation to the Debye length κ−1 introduced in Section 5.
We next discuss the GL free energy F. In KAPSEL we consider phase separation in two-component fluids interacting with particles, and thus F takes the form
Z


α
(6.4)
F=
dr f (ψ) + (∇ψ)2 + wξψ|∇ϕ|2 + d(ψ − ψ̄)2 ϕ + zξψ ϕ(∇ψ)2 .
2Here the function f (ψ) is a double-well potential governing the local dynamics of the fluid-fluid interface function ψ.
This potential may be chosen to be either a Landau potential:f (ψ) =

a 4 b 2
ψ − ψ
4
2

(6.5)

or a Flory-Huggins (FH) potential:"f (ψ) = kB T 0#ψ
(1 − ψ)ln ψ +ln(1 − ψ) + χψ(1 − ψ) .
NA
NB

(6.6)

CHAPTER 6. SIMULATING PARTICLES DISPERSED IN TWO-COMPONENT PHASE-SEPARATED FLUIDS59

Figure 6.2: Effect of the z parameter on the particle dispersion.

The quantities a and b in equation (6.5) are adjustable parameters used to tune the relative strengths of the twoterms in the potential. Also, in equation (6.6) the quantity kB T 0 is the reference unit of thermal energy, NA and NB are
the degrees of polymerization of the fluid components, and χ is the Flory interaction parameter.
The second term in equation (6.4) represents interfacial energy associated with fluid-fluid interfaces. The third termaψ3 −(b−2d)ψ
describes particle-fluid interactions. The fourth term is included to ensure, through the relation ψ̄ = 0 2d 0 , that ψ
takes the arbitrary fixed value ψ0 at points inside the particles.1 Finally, the fifth term describes interactions between
particles and fluid-fluid interfaces; here z is a coefficient and ξψ is the fluid-fluid interface thickness.2 Figure 6.2
illustrates how the z parameter affects particle dispersion. For z < 0, the free energy is minimized when particles
exist at fluid-fluid interfaces, so particles migrate toward those interfaces in that case. In contrast, for z > 0 particles
distance themselves from fluid-fluid interfaces. The fluid-fluid interface thickness ξψ depends on the function f (ψ);
for a Landau double-well potential it is given by [41].
r
α
ξψ =
.
(6.7)
2b
In simulating particles dispersed in two-component phase-separated fluids we successively solve the equations of
motion for the fluid and for the particles, equations (6.1) and (3.4)-(3.6), together with equation (6.2) for the time
evolution of the fluid-fluid interface function ψ, which distinguishes the two fluid components. KAPSEL solves theseequations on a semi-staggered Arakawa B lattice [42] using the marker-and-cell (MAC) method [43, 44]. Combining
this approach with the techniques of Section 4 allows for simulations of particles dispersed in two-component phaseseparated fluids in the presence of shear flows. Solution methods are discussed in detail in Appendix D.

6.2.2

Viscosity of suspensions

For a suspension of particles in a single-component fluid, the shear stress and viscosity are given by equations (4.5)
and (4.6) On the other hand, when particles are added to a two-component phase-separated fluid, there are additionalsurface-tension terms associated with fluid-fluid interfaces, whose contributions must be taken into account. The xy
component of the stress tensor due to such interfaces is given by [45]
!
Z
1
∂ψ ∂ψ
xy
σint =
dr −α
.
(6.8)
V
∂x ∂y
Consequently, the apparent viscosity of the system is obtained by adding the contribution of equation (6.8) to equation(4.6). Section 6.4.2 presents a computational example involving a viscosity calculation in the presence of steady-state
shear flow.

6.3

Input UDF files

6.3.1

Fluid settings

KAPSEL simulations of particles dispersed in two-component phase-separated fluids may be performed (a) with no

imposed flow field, or (b) in the presence of a shear flow. Set constitutive eq to Navier Stokes Cahn Hilliard FDM
1 Setting d = 0 must be avoided. To ignore this term, you should set such as d = 0.0001.
2 This term has no effect if f (ψ) is chosen to be of the Flory-Huggins form.

CHAPTER 6. SIMULATING PARTICLES DISPERSED IN TWO-COMPONENT PHASE-SEPARATED FLUIDS60

for case (a) or to Shear NS LE CH FDM for case (b).
Settings for simulations with no imposed flow field➤constitutive eq.Navier Stokes Cahn Hilliard FDM: Properties of two-component phase-separated fluids
constitutive eq.Navier Stokes Cahn Hilliard FDM.NS solver.type
Specifies the method used to solve the NS equations. The allowed values are explicit scheme or
implicit scheme.3
constitutive eq.Navier Stokes Cahn Hilliard FDM.NS solver.implicit scheme.tolerance
Convergence criterion for implicit solution of NS equations.
constitutive eq.Navier Stokes Cahn Hilliard FDM.NS solver.implicit scheme.maximum iteration
Maximum number of iterations for implicit solution of NS equations.
constitutive eq.Navier Stokes Cahn Hilliard FDM.NS solver.implicit scheme.viscosity change
Set to ON to enable simulations of two-component phase-separated fluids (A/B) with different viscosities.
constitutive eq.Navier Stokes Cahn Hilliard FDM.NS solver.implicit scheme.ON.ETA A
Viscosity of fluid component A.
constitutive eq.Navier Stokes Cahn Hilliard FDM.NS solver.implicit scheme.ON.ETA B
Viscosity of fluid component B.
constitutive eq.Navier Stokes Cahn Hilliard FDM.CH solver.type
Specifies the method used to solve the CH equation. The allowed values are explicit scheme or
implicit scheme.4
constitutive eq.Navier Stokes Cahn Hilliard FDM.CH solver.implicit scheme.tolerance
Convergence criterion for implicit solution of CH equation.
constitutive eq.Navier Stokes Cahn Hilliard FDM.CH solver.implicit scheme.maximum iteration
Maximum number of iterations for implicit solution of CH equations.
constitutive eq.Navier Stokes Cahn Hilliard FDM.DX
Lattice spacing ∆, which defines the length unit.
constitutive eq.Navier Stokes Cahn Hilliard FDM.RHO
Fluid density ρ.
3 KAPSEL implements the MAC implicit method [46] as an implicit solver for the NS equations. The use of an implicit solver increases thecomputation time, but allows simulations of two-component phase-separated fluids with different viscosities. Solution methods are discussed indetail in Appendix D.1. When implicit scheme is selected, BiCGSTAB [47] is chosen as the default iterative solver for simultaneous systemsof linear equations, and several new user-configurable simulation settings, discussed in the table above, become available. When using iterativesolvers for implicit methods, the solution is considered converged when the residual falls below the specified tolerance value. If the number of
iterations exceeds maximum iteration, KAPSEL concludes that no solution was obtained and the simulation is terminated. Iterative solvers otherthan BiCGSTAB may be used by linking with the Lis library [48]. This is discussed in more detail in Appendix E.
4 In addition to the explicit solver, methods implemented by KAPSEL for solving the CH equation include a backward-difference method for time
derivatives and a quasi-implicit solver [49] based on values extrapolated by the Adams-Bashforth method for nonlinear potential terms. Solution
methods are discussed in detail in Appendix D. For the quasi-implicit CH solver, iterative solvers other than BiCGSTAB may be used by linking
with the Lis library [48]. This is discussed in more detail in Appendix E.

CHAPTER 6. SIMULATING PARTICLES DISPERSED IN TWO-COMPONENT PHASE-SEPARATED FLUIDS61

constitutive eq.Navier Stokes Cahn Hilliard FDM.ETA
Fluid viscosity η.5
constitutive eq.Navier Stokes Cahn Hilliard FDM.kBT
Particle temperature.
constitutive eq.Navier Stokes Cahn Hilliard FDM.alpha v
Correction term for particle temperature associated with translational motion.
constitutive eq.Navier Stokes Cahn Hilliard FDM.alpha oCorrection term for particle temperature associated with rotational motion.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.type
Set to Landau or Flory Huggins to select the double-well potential for a two-component phase-separated
fluid.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.Landau.composition ratio
Fluid composition ratio for Landau potential.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.Landau.initial fluctuation
Initial concentration fluctuation for Landau potential.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.Landau.aa parameter for 4th-order term in Landau potential.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.Landau.bb parameter for 2nd-order term in Landau potential.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.Landau.dd parameter used to force the composition of fluid regions inside particles to an arbitrary value ψ0 whenusing a Landau potential.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.Landau.w
Particle-fluid interaction parameter w for Landau potential.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.Landau.z
z parameter relevant for interactions between particles and fluid-fluid interfaces when using a Landau
potential.6
constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.Landau.psi 0
Arbitrary value ψ0 for ψ function at points inside particles for Landau potential.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.Landau.alphaFluid-fluid interfacial energy parameter α for Landau potential.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.Landau.kappa
Mobility κ for Landau potential.
5 The simulation timestep is chosen based on the viscosity value specified here. Caution: When viscosity change is set to ON, the viscositiesof the two fluid components do not affect the timestep. Consequently, specifying values for ETA A or ETA B that are significantly larger than thevalue set here for ETA may degrade the accuracy of the simulations or lead to computational instabilities. This is discussed in more detail in
Section 6.3.3.
6 Particle interactions with fluid-fluid interfaces are only active when using a Landau potential.

CHAPTER 6. SIMULATING PARTICLES DISPERSED IN TWO-COMPONENT PHASE-SEPARATED FLUIDS62

constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.Flory Huggins.composition ratio
Fluid composition ratio for FH potential.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.Flory Huggins.initial fluctuation
Initial concentration fluctuation for FH potential.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.Flory Huggins.na
Degree of polymerization NA for fluid component A.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.Flory Huggins.nb
Degree of polymerization NB for fluid component B.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.Flory Huggins.chi
Flory interaction parameter χ.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.Flory Huggins.dd parameter used to force the composition of fluid regions inside particles to an arbitrary value ψ0 whenusing an FH potential.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.Flory Huggins.w
Particle-fluid interaction parameter w for FH potential.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.Flory Huggins.z
z parameter relevant for interactions between particles and fluid-fluid interfaces when using FH potential.7
constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.Flory Huggins.psi 0
Arbitrary value ψ0 for ψ function at points inside particles for FH potential.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.Flory Huggins.alphaFluid-fluid interfacial energy parameter α for FH potential.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Potential.Flory Huggins.kappa
Mobility κ for FH potential.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Wall Potential.type
Set to ON to add planar walls to simulations.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Wall Potential.ON.w
Parameter determining affinity of fluid for planar walls.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Wall Potential.ON.psi 0.magnitude
Parameter determining magnitude of affinity of fluid for planar walls.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Wall Potential.ON.psi 0.profile
Set to uniform for uniform fluid affinity on planar walls; otherwise set to user specify and configure
manually.
7 Particle interactions with fluid-fluid interfaces are only active when using a Landau potential.

CHAPTER 6. SIMULATING PARTICLES DISPERSED IN TWO-COMPONENT PHASE-SEPARATED FLUIDS63

constitutive eq.Navier Stokes Cahn Hilliard FDM.Wall Potential.ON.psi 0.user specify.PSI0[][]Manually specify a value for the affinity between fluid and planar walls.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Wall Potential.ON.DRYING.type
Set to ON to simulate drying.
constitutive eq.Navier Stokes Cahn Hilliard FDM.Wall Potential.ON.DRYING.ON.psi dry
Parameter determining rapidity of drying.
Settings for simulations under shear flow
➤Shear NS LE CH FDM: Properties of two-component phase-separated fluids
constitutive eq.Shear NS LE CH FDM.NS solver.type
Specifies the method used to solve the NS equations. The allowed values are explicit scheme or
implicit scheme.8
constitutive eq.Shear NS LE CH FDM.NS solver.implicit scheme.tolerance
Convergence criterion for implicit solution of NS equations.
constitutive eq.Shear NS LE CH FDM.NS solver.implicit scheme.maximum iteration
Maximum number of iterations for implicit solution of NS equations.
constitutive eq.Shear NS LE CH FDM.NS solver.implicit scheme.viscosity change
Set to ON to enable simulations of two-component phase-separated fluids (A/B) with different viscosities.
constitutive eq.Shear NS LE CH FDM.NS solver.implicit scheme.ON.ETA A
Viscosity of fluid component A.
constitutive eq.Shear NS LE CH FDM.NS solver.implicit scheme.ON.ETA B
Viscosity of fluid component B.
constitutive eq.Shear NS LE CH FDM.CH solver.type
Specifies the method used to solve the CH equation. The allowed values are explicit scheme or
implicit scheme.9
constitutive eq.Shear NS LE CH FDM.CH solver.implicit scheme.tolerance
Convergence criterion for implicit solution of CH equations.
constitutive eq.Shear NS LE CH FDM.CH solver.implicit scheme.maximum iteration
Maximum number of iterations for implicit solution of CH equations.
constitutive eq.Shear NS LE CH FDM.DX
Lattice spacing ∆, which defines the length unit.
constitutive eq.Shear NS LE CH FDM.RHO
Fluid density.
8 As is true for simulations with no imposed flow field, the MAC implicit method [46] is implemented as an implicit solver for the NS equations.
9 As is true for simulations with no imposed flow field, a quasi-implicit solver for the CH equation is implemented in addition to the explicit

solver. [49].

CHAPTER 6. SIMULATING PARTICLES DISPERSED IN TWO-COMPONENT PHASE-SEPARATED FLUIDS64

constitutive eq.Shear NS LE CH FDM.ETA
Fluid viscosity.10
constitutive eq.Shear NS LE CH FDM.kBT
Particle temperature.
constitutive eq.Shear NS LE CH FDM.alpha v
Correction term for particle temperature associated with translational motion.
constitutive eq.Shear NS LE CH FDM.alpha oCorrection term for particle temperature associated with rotational motion.
constitutive eq.Shear NS LE CH FDM.Potential.type
Set to Landau or Flory Huggins to select the double-well potential for a two-component phase-separated
fluid.
constitutive eq.Shear NS LE CH FDM.Potential.Landau.composition ratio
Fluid composition ratio for Landau potential.
constitutive eq.Shear NS LE CH FDM.Potential.Landau.initial fluctuation
Initial concentration fluctuation for Landau potential.
constitutive eq.Shear NS LE CH FDM.Potential.Landau.aa parameter for 4th-order term in Landau potential.
constitutive eq.Shear NS LE CH FDM.Potential.Landau.bb parameter for 2nd-order term in Landau potential.
constitutive eq.Shear NS LE CH FDM.Potential.Landau.dd parameter used to force the composition of fluid regions inside particles to an arbitrary value ψ0 whenusing a Landau potential.
constitutive eq.Shear NS LE CH FDM.Potential.Landau.w
Particle-fluid interaction parameter w for Landau potential.
constitutive eq.Shear NS LE CH FDM.Potential.Landau.z
z parameter relevant for interactions between particles and fluid-fluid interfaces when using a Landau
potential.11
constitutive eq.Shear NS LE CH FDM.Potential.Landau.psi 0
Value for ψ function at points inside particles.
constitutive eq.Shear NS LE CH FDM.Potential.Landau.alphaFluid-fluid interfacial energy parameter α for Landau potential.
constitutive eq.Shear NS LE CH FDM.Potential.Landau.kappa
Mobility κ for Landau potential.
10 The simulation timestep is chosen based on the viscosity value specified here. Caution: When viscosity change is set to ON, the viscositiesof the two fluid components do not affect the timestep. Consequently, specifying values for ETA A or ETA B that are significantly larger than the value
set here for ETA may degrade the accuracy of simulations or lead to computational instabilities. This is discussed in more detail in Section 6.3.3.
11 Particle interactions with fluid-fluid interfaces are only active when using a Landau potential.

CHAPTER 6. SIMULATING PARTICLES DISPERSED IN TWO-COMPONENT PHASE-SEPARATED FLUIDS65

constitutive eq.Shear NS LE CH FDM.Potential.Flory Huggins.composition ratio
Fluid composition ratio for FH potential.
constitutive eq.Shear NS LE CH FDM.Potential.Flory Huggins.initial fluctuation
Initial concentration fluctuation for FH potential.
constitutive eq.Shear NS LE CH FDM.Potential.Flory Huggins.na
Degree of polymerization NA for fluid component A.
constitutive eq.Shear NS LE CH FDM.Potential.Flory Huggins.nb
Degree of polymerization NB for fluid component B.
constitutive eq.Shear NS LE CH FDM.Potential.Flory Huggins.chi
Flory interaction parameter χ.
constitutive eq.Shear NS LE CH FDM.Potential.Flory Huggins.dd parameter used to force the composition of fluid regions inside particles to an arbitrary value ψ0 whenusing an FH potential.
constitutive eq.Shear NS LE CH FDM.Potential.Flory Huggins.w
Particle-fluid interaction parameter w for FH potential.
constitutive eq.Shear NS LE CH FDM.Potential.Flory Huggins.z
z parameter relevant for interactions between particles and fluid-fluid interfaces when using FH potential.12
constitutive eq.Shear NS LE CH FDM.Potential.Flory Huggins.psi 0
Value for ψ function at points inside particles for FH potential.
constitutive eq.Shear NS LE CH FDM.Potential.Flory Huggins.alphaFluid-fluid interfacial energy parameter α for FH potential.
constitutive eq.Shear NS LE CH FDM.Potential.Flory Huggins.kappa
Mobility κ for FH potential.
constitutive eq.Shear NS LE CH FDM.External field.type
Set to DC (steady-state shear flow) or AC (oscillatory shear flow).
constitutive eq.Shear NS LE CH FDM.External field.DC.Shear rate
Shear velocity γ̇.
constitutive eq.Shear NS LE CH FDM.External field.AC.Shear rate
Amplitude γ̇0 of oscillatory shear velocity.
constitutive eq.Shear NS LE CH FDM.External field.AC.Frequency
Frequency ω of oscillatory shear velocity.
12 Particle interactions with fluid-fluid interfaces are only active when using a Landau potential.

CHAPTER 6. SIMULATING PARTICLES DISPERSED IN TWO-COMPONENT PHASE-SEPARATED FLUIDS66

6.3.2

Object (particle) settings

The object type.type setting specifies the type of particle. The possible values are spherical particle forspherical particles, chain for flexible chains, or rigid for rigid bodies.

6.3.3

Choice of length and time units

The length unit is given by the lattice spacing ∆. The time unit τ0 is determined by the fluid density ρ, the fluid
viscosity coefficient η, and the lattice spacing ∆ according to τ0 = ρ∆2 /η.
Note that, for simulations of two-component phase-separated fluids with different viscosities, the viscosity η usedin computing the unit of time is the value specified for constitutive eq.*.ETA. Note carefully that the viscosities of
the individual fluid components, i.e. the values specified for constitutive eq.*.NS solver.implicit scheme.ON.ETA A
and constitutive eq.*.NS solver.implicit
scheme.ON.ETA B are not referenced when determining the simulation timestep. (Here * stands for either Navier Stokes Cahn Hillior Shear NS LE CH
FDM.)
•

•

Assuming the input UDF file specifies RHO= A, ETA= B, and DX= C, the maximum wavenumber kmax may
be determined from C and the upper bound on the timestep is determined from the momentum diffusion time2in the form T step = (A/B)/kmax. The actual simulation timestep is related to this bound through a multiplicative
scale factor, i.e. ∆t = factor × T step .
Consider the correspondence between physical quantities in simulations and in reality. A grid spacing of 1 µm
is reasonable for the length scales we wish to consider. Assuming we use water (η = 1 × 10−3 Pa s, ρ =1 × 103 kg m−3 ) as a solvent, the time unit is then τ0 = 1 × 10−6 s.

CHAPTER 6. SIMULATING PARTICLES DISPERSED IN TWO-COMPONENT PHASE-SEPARATED FLUIDS67

6.4

Computational examples

We now present examples of simulations using the methods described in the preceding subsections to study systems of particles dispersed in two-component phase-separated fluids. Section 6.4.1 describes a brief simulation thatprovides an outline of the entire KAPSEL simulation process. This simulation involves particles dispersed in twocomponent phase-separated fluids with no imposed flow field. Section 6.4.2 describes a simulation of a two-component
phase-separated fluid with dispersed particles under shear flow, in order to investigate rheological properties. Finally,
Section 6.4.3 describes a simulation of a Pickering emulsion [35].

6.4.1

Motion of particles dispersed in a two-component phase-separated fluid

Running the simulation
The UDF files for the simulations of this section are in the Examples/11a/ folder. Use the following commands toenter this folder and run the simulation.
$ cd Examples/11a$ ../../kapsel -Iinput.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
This simulation uses the Landau double-well potential of equation (6.5) to describe phase-separated two-component
fluids. Parameter values are written to the standard output when the simulation begins:################################## Landau potential is selected.
# Parameters# Composition ratio: 0# Initial fluctuation : 0.05
# a
: 1
# b
: 1
# d
: 1
# w
: -1
# z
: 0
# psi0_p : 0# alpha : 1# kappa : 1##################################For this simulation, Composition ratio is set to 0 to indicate that the two fluid components are present in equal
quantities in the system. Also, w is set to −1 to indicate that particles have greater affinity for one of the two components(component A). psi0 p is a parameter specifying the value of ψ in the particle interior.
If the simulation completes successfully, runtime information will be reported to the console:#Simulation has ended!
#Total Running Time (s):
#
(m):
#
(h):
#Average Step Time (s):
#
(m):
#
(h):

11.89
0.20
0.00
0.01
0.00
0.00

Computation times vary depending on runtime environment, but this simulation should complete in a few tens ofseconds. Upon completion of the simulation, the following output files will have been written to the Examples/11a/folder:•output.udf•

restart.udf•

fluid phase.xmf•

flux field.xmf•

particle coordinate.xmf•

particle phase.xmf•

velocity field.xmf

CHAPTER 6. SIMULATING PARTICLES DISPERSED IN TWO-COMPONENT PHASE-SEPARATED FLUIDS68•

flux *.h5•

orderparam *.h5•

particle *.h5•

particle data *.h5•

velocity *.h5

Here asterisks (*) in file names represent sequence numbers. output.udf is an output UDF file containing simulation
results that may be read by GOURMET. (Methods for viewing data in GOURMET are discussed below.) Also,
restart.udf is a restart UDF file that may be used to restart the simulation from an intermediate time. Files withextension h5 and sequence numbers in the filename store simulation data at individual timesteps in the HDF5 binary
data format. These files may be bundled for shared processing via eXtensible Data Model and Format (XDMF) control
files, which are used to manage time-series simulation data. XDMF files have extension xmf and are used by KAPSEL
to establish rules for working with groups of h5 files, as follows:
fluid phase.xmf:
flux field.xmf:
particle coordinate.xmf:
particle phase.xmf:
velocity field.xmf:

Data on the function ψ that distinguishes the two fluid components (orderparam *.h5)
Data on the mass flux J, equation (6.3) (flux *.h5)
Particle coordinate data (particle data *.h5)
Data on the function ϕ that distinguishes particles: (particle *.h5)
Data on the fluid velocity field u (velocity *.h5)

XDMF files may be displayed using visualization software such as ParaView13 or Mayavi14 , and the following
section presents an example of ParaView visualization. Data stored in .h5 files may also be viewed directly using theh5dump command-line tool, e.g.:$ h5dump orderparam_0.h5

Visualizing simulation data
In this section we describe how to use GOURMET and ParaView to visualize simulation results.
Visualizing simulation data with GOURMET
Reading an output UDF file (output.udf) into GOURMET allows browsing of simulation parameters and visualization of simulation results at various output steps [Figure (6.3)(a)]. Also, the folder Examples/11a/ contains a
Python script named gourmet particle field show.py that may be read into GOURMETto visualize the function
ψ that distinguishes between the two fluid components. More specifically, using the Python panel at the bottom of the
GOURMET screen to read and execute this Python script will open a new graphical window presenting an animatedview of simulation results [Figure (6.3)(b)]. In this figure, particles are represented by white spheres and fluid interfaces are displayed as blue constant-value surfaces, but these and other display settings may be modified directly in
the Python script.
Visualizing simulation data with ParaView
XDMF files may be read into ParaView to visualize or analyze simulation results. The file pv sample.pvsm
in the Examples/11a/ folder is a ParaView State file. Launch ParaView and select Load State Files from theFiles menu at upper left. Read in the file pv sample.pvsm to visualize the function ψ used to distinguish the two
fluid components (Figure 6.4).15 In this figure, particles are represented as white spheres; one fluid component isrepresented by red-colored (volume rendered) regions, and the other fluid component is displayed using translucent
blue-colored regions. ParaView is not only useful for data visualization, but also offers powerful analysis capabilities.
For further details, consult the ParaView documentation [50].
Writing Python scripts to analyze simulation data
Simulation data written to output files by KAPSEL may be separately read and analyzed via Python scripts. The
Examples/11a/ folder includes a Python script that computes the static structure factor S (q) [51] for a phaseseparation structure. Make sure the python packages h5py, scipy, and matplotlib are installed in your Python
13 https://www.paraview.org/
14 http://docs.enthought.com/mayavi/mayavi/
15 Choose “Search files under specified directory” for ”Load State Data File Options”, and set appropriate path for ”Data Directory”.

CHAPTER 6. SIMULATING PARTICLES DISPERSED IN TWO-COMPONENT PHASE-SEPARATED FLUIDS69(a) Listing parameter values recorded in output UDF files.

(b) Using a Python script to visualize the function ψ that distinguishes
between the two fluid components.

Figure 6.3: Visualizing simulation results in GOURMET.

Figure 6.4: Visualizing simulation results in Paraview.

CHAPTER 6. SIMULATING PARTICLES DISPERSED IN TWO-COMPONENT PHASE-SEPARATED FLUIDS70

S(q) vs q 5steps
101
100

S(q)

10 1
10 2
10 3
10 4
10 5
100

q

Figure 6.5: Plot of static structure factor S (q) vs wavenumber q, illustrating the use of python scripting for data analysis. Results
shown for the 5th simulation timestep.

environment, then execute the command below. This will produce a sequence of S (q) graphs for various timesteps(Figure 6.5).
$ python sq_analysis.py

6.4.2

Motion of particles dispersed in a two-component phase-separated liquid under shear
flow

The input file for this example is in the Examples/11/ folder. It describes a simulation of a two-component phaseseparated fluid, with dispersed particles, in the presence of a steady-state shear flow. Execute the following command
to run the simulation:$ cd Examples/11$ ../../kapsel -Iinput.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
For this example, we consider 128 particles of radius a = 4 and interface thickness ξ = 2 added to a two-component
phase-separated fluid in which a steady-state shear flow with shear velocity γ̇ = 0.01 has been imposed. We set theparticle temperature to kB T = 1 and use a 64 × 64 × 64 computational mesh. The two-component fluid exhibits phaseseparation due to the Landau double-well potential of equation (6.5). Parameter values are written to standard output
when the simulation begins:################################## Landau potential is selected.
# Parameters# Composition ratio: 0# Initial fluctuation : 0.5
# a
: 1
# b
: 1
# d
: 1
# w
: 0
# z
: 0
# psi0_p : 0# alpha : 1# kappa : 1##################################For simulations in the presence of shear flow, values for the following data items are written to standard error one line
per timestep:CHAPTER 6. SIMULATING PARTICLES DISPERSED IN TWO-COMPONENT PHASE-SEPARATED FLUIDS71
1:time2:shear rate
3:degree oblique
4:shear strain temporal
5:lj dev stress temporal
6:shear stress temporal old
7:shear stress temporal new
8:reynolds stress
9:fluid stress
10:interfacial stress
11:apparent stress
12:viscosity

Elapsed time
Shear velocity
Strain in oblique coordinate system 16
Applied strain (γ̇t)
Shear stress due to interparticle potential
Shear stress due to particle-fluid interactions
Shear stress due to particle-fluid interactions 17
Reynolds shear stress
Shear stress due to bulk fluid
Shear stress due to fluid-fluid interface
Shear stress on suspension 18
Viscosity of suspension 19

y

z

x

Figure 6.6: Snapshot of fluid interface function ψ captured on the final timestep of the simulation of Example 11 (particles in
two-component fluid under shear flow).

Figure 6.6 is a snapshot of the fluid interface function ψ captured in the final timestep (step number 20, 000) of
this simulation. In KAPSEL the shear flow is applied with the shear flow in the x direction, the shear gradient in they direction, and the z direction axis as the neutral direction (the vortex direction). Because the parameter w governing
the affinity of particles for the two-component fluid is set to 0, the particles are not preferentially located in one phaseor the other, but are instead distributed throughout the full system.
Figure 6.7 plots the time evolution of the xy component of the shear stress, with the red curve indicating the shear
stress on the suspension (11:apparent stress in the standard-error output) and the blue curve indicating the shear
stress due to the fluid-fluid interface (10:interfacial stress). The results of this simulation indicate that theshear stress due to the fluid-fluid interface contributes significantly to the overall shear stress on the suspension. Also,
whereas small fluctuations due to thermal motion of particles are visible in the curve for the overall shear stress on thesuspension, no such fluctuations are seen in the curve for the fluid-fluid-interface contribution to the shear stress. Thisis because KAPSEL introduces thermal fluctuation of particles not through the Navier-Stokes equations for the fluid,
but through the Langevin equation for the particles. A detailed discussion of particle temperatures in KAPSEL may
be found in Ref. 3.1.1.

6.4.3

Pickering emulsions

The input files for this example are in the Examples/13/ folder. They define a simulation of a Pickering emulsion,
an emulsion stabilized by the addition of microparticles as first reported by Pickering [35]. Pickering emulsions arestabilized by the adsorption of microparticles at the fluid-fluid interface. In KAPSEL the affinity of particles for
two-component fluid interfaces may be tuned via the z parameter in the 5th term of equation (6.4). This simulationconsiders 200 particles (radius a = 4 and interface thickness ξ = 2) dispersed in a two-component fluid, with theparticles having neutral affinity (w = 0) with respect to the fluid, but nonzero affinity (z = −0.6) with respect to thefluid-fluid interfaces . The two-component fluid undergoes a phase separation induced by the Landau double-wellpotential in equation (6.5), with parameter values a = b = d = κ = 1 and α = 3.

CHAPTER 6. SIMULATING PARTICLES DISPERSED IN TWO-COMPONENT PHASE-SEPARATED FLUIDS72

0.05
total stressinterfacial stressinterfacial stress

0.04

0.03

0.02

0.01

0.00

0

500

time

1000

1500

Figure 6.7: Time evolution of the shear stress (xy component) for the simulation of Example 11 (particles in two-component fluid
under shear flow). Red: shear stress on the suspension (11:apparent stress). Blue: shear stress due to the fluidfluid interface (10:interfacial stress).

The following command runs the simulation:$ cd Examples/13$ ../../kapsel -Iinput.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
In the simulation defined by the input.udf file in the Examples/13 folder, particles have nonzero affinity (z = −0.6)
for interfaces between the particle and the fluid-fluid interfaces.
Figure 6.8 shows a snapshot of the fluid interface function ψ captured in the final timestep (step 200, 000) for thissimulation. As we see in this image, the nonzero affinity of particles for fluid-fluid interfaces causes particles to be
attracted to these interfaces, where they function as surfactants.

CHAPTER 6. SIMULATING PARTICLES DISPERSED IN TWO-COMPONENT PHASE-SEPARATED FLUIDS73

Figure 6.8: Snapshot of the fluid interface function ψ captured on the final timestep of the simulation of Example 13 (Pickering
emulsion).

Chapter 7

Simulating microswimmers
7.1

Theoretical background and basic equations

Microswimmers—microscopic bodies moving autonomously through viscous fluids, with examples including spermatozoa, bacteria, and algae—are ubiquitous throughout biology and the life sciences. Microswimmers also constituteone of the most important types of active matter, soft matter exhibiting spontaneous self-propelled motion. Studies ofactive matter are extremely fruitful both for theoretical science—where active-matter research promises to make major
contributions to the development and validation of non-equilibrium physics theories—and for practical applications
such as smart-drug delivery and other cutting-edge medical innovations. However, the complex role of hydrodynamic
interactions in governing the physics of microswimmer systems remains poorly understood, with efforts to elucidate
the relevant mechanisms obstructed by difficult and longstanding challenges [52].

7.1.1

The squirmers model

KAPSEL uses a general-purpose model of microswimmers known as the squirmer model. This model, developed

by Lighthill and Blake [53, 54], is commonly used to investigate hydrodynamic interactions among spheroidal microswimmers.
KAPSEL considers spherical swimmers of radius a. As an alternative to fine-grained modeling of the motionof cilia on particle surfaces, we impose specific boundary conditions at rigid-particle surfaces, which yield similar
self-propulsion effects. We take {ẽi } to be an orthogonal system of basis vectors defining a reference system attachedto the particles and fix the direction of particle motion as ẽ3 = ẽ. Points on the particle surface are described by theirpolar and azimuthal angles θ, λ, and we denote by {θ, λ, ϱ}, the unit tangent vectors in the angular direction and theunit vector in the radial direction (so that e.g. θ = arccos ϱ · ẽ). Assuming zero velocity in the radial direction, the slip
boundary condition imposed on the squirmers reads [55]u sq =∞X∞
X
2
Bn P′n (cos θ) sin θθ +Cn P′n (cos θ) sin θλ
n(n
+
1)
n=1
n=1

(7.1)

where P′n is the derivative of the nth Legendre polynomial and Bn , Cn are the nth polar and azimuthal mode coefficients.
Neglecting all azimuthal modes (Cn = 0) and polar modes of degree 2 or more (Bn = 0 where n > 2) yields a simple
but important model capable of describing both pushers (swimmer particles for which propulsive forces arise behindthe particle) and pullers (swimmers with propulsive forces arising in front of the particle). Denoting the ratio of thefirst two polar modes as α = B2 /B1 , the surface velocity is given by


α
(7.2)
u sq (θ) = B1 sin θ + sin 2θ θ.
2The parameter B1 sets the steady-state velocity U of a single particle moving through the fluid according to U = 2/3B1 ,while B2 is the stresslet intensity. The ratio of mode coefficients α determines the type of swimmer: puller (α > 0),
pusher (α < 0), or neutral swimmer (α = 0). Figure 7.1 illustrates the motion of pushers and pullers, and the fluid
flow accompanying their propulsion. The character of a swimmer’s motion heavily influences the hydrodynamics ofthe fluid flow induced by its propulsion, which must also be taken into account when considering collective particle
behavior.

74

CHAPTER 7. SIMULATING MICROSWIMMERS

75

Figure 7.1: Schematic depiction of self-propulsion mechanisms and accompanying fluid flow for (a) pushers, (b) pullers. Describing these swimmers by Blake’s squirming model—in which microscopic propulsion mechanisms are replaced by slip
boundary conditions at particle surfaces—yields images (c) and (d). Reproduced from Soft Matter 9, 4923-4936[4],Copyright 2013, with permission from the Royal Society of Chemistry.

7.1.2

Basic equations for microswimmers

For SP-method simulations of the fluid mechanics of disperse microswimmer systems, we add a new constraint force
term ϕf sq to the continuous-medium equation (3.13) and apply slip boundary conditions at particle surfaces [4, 56]:ρ(∂t + u · ∇)u = ∇ · σ + ρϕf p + ρϕf sq .

(7.3)

To solve this equation, we add one new step to the usual procedure for SP-method calculations. After updating theadvection and viscosity terms, and computing u⋆ , the term ϕf p serves to exert rigid-body constraint conditions; at thisjuncture we compute a surface constraint force ϕf sq and use it to update a velocity field u∗∗ that satisfies the squirmer
boundary conditions. To implement this strategy, we introduce a new SP function ϕ sq , which takes nonzero values
only in the interior of particle interfaces:
ϕIsq = (1 − ϕI )|∇ϕI |.
max (|∇ϕI |)(7.4)

The new velocity field takes the formu =u +∗∗∗"Z tn +h#dsϕf sq(7.5)

tn"Z tn +h#dsϕf sq =

tn

N
X



ϕIsq VI† + Ω†I × rI + uIsq − u∗

(7.6)

I

+

N
X
I

h
ϕI (δVI + δΩi × rI ) − ∇p sq
ρThe first term on the RHS here is a term for actually fixing the slip boundary conditions, while the second term ensures
local momentum conservation (specifically, that squirmers pushing/pulling on the fluid at their boundaries experiencea counterforce). The third term ensures that the fluid velocity field ultimately remains incompressible. Because theupdated particle velocities are not yet known at this point, we solve equation iteratively, with slip boundary conditions
applied for VI† (Ω†I ). We compute hydrodynamic forces and torques as described previously (using u∗∗ instead ofu∗ ), use the results to compute particle velocities, and iterate until the particle velocities VIn=1 for the following time
step agree with the velocity VI† needed to impose slip boundary conditions. For the motion of a single swimmer, wehave validated the efficacy of this computational procedure and confirmed that simulated particle velocities and fluid
velocities are in good agreement with theoretical predictions [4].

CHAPTER 7. SIMULATING MICROSWIMMERS

7.2

Input UDF files

7.2.1

Fluid settings

76

The following choices are available for constitutive eq.
Navier Stokes:
Shear Navier Stokes:
Navier Stokes FDM:
Navier Stokes Cahn Hilliard FDM:

7.2.2

Newtonian fluid
Newtonian fluid with zigzag shear flow
Newtonian fluidTwo-component phase-separated fluid

Configuring objects (particles)

The only available option for object type.type is spherical particle.

7.2.3

Configuring objects (planar walls)

switch.wall.type may be set to NONE or FLAT.

CHAPTER 7. SIMULATING MICROSWIMMERS

7.3

Computational examples

7.3.1

Motion of microswimmers with periodic boundary conditions

77

The input UDF files for this example are squirm single a+2.udf and squirm phi0.1 a+2.udf in the Examples/09/
folder. The former of these describes a simulation with just a single particle.
$ ../../kapsel -Isquirm_single_a+2.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
In this example we use a 64×64×64 computational mesh with the the following parameter values: number of particles
N p = 1, particle diameter D = 10, and interface thickness ξ = 2, yielding a particle volume fraction of φ = 0.002.
The velocity of microswimmer motion is V = 2/3B1 = 0.01 and the character of the motion is chosen by specifying
B2 = 2.
The second file, squirm phi0.1 a+2.udf, describes a similar simulation but with multiple particles:$ ../../kapsel -Isquirm_phi0.1_a+2.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
In this case we have N p = 50 particles of diameter D = 10 and interface thickness ξ = 2, yielding a particle volume
fraction of φ = 0.1. Other parameters are set as in the single-particle simulation. Fig. 7.2 is a snapshot of the 50swimmer simulation.

Figure 7.2: Snapshot from simulation of 50 microswimmers.

7.3.2

Motion of microswimmers confined between two parallel slabs

The input UDF file for this example is squirm wall.udf in the Examples/09/ folder. This file describes a simulationin which microswimmers are confined to the region between two parallel slabs.
$ ../../kapsel -Isquirm_wall.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
For this example (Fig. 7.3) we use a 64 × 256 × 64 computational mesh. The thickness of the planar walls in thevertical direction (y axis) is 4. We specify N p = 4266 particles with particle diameter D = 4 and interface thickness
ξ = 2, yielding a particle volume fraction of φ = 0.138. The velocity of microswimmer motion is V = 2/3B1 = 0.25
and the character of the motion is chosen by specifying B2 = 0.5. We use the python script particle show wall.py
for visualization in GOURMET.

CHAPTER 7. SIMULATING MICROSWIMMERS

Figure 7.3: Motion of microswimmers confined between two parallel slabs.

78

Chapter 8

Simulating Quincke rollers
8.1

Theoretical background and basic equations

The study of active matter is not restricted to the motion of bacteria and other living organisms, but encompasses
inanimate bodies as well, with systems of self-propelled colloids furnishing an important class of examples. The studyof such systems is promising not only for deriving insight into fundamental properties, such as driving mechanisms
and collective dynamics, but also for developing practical applications involving control via external fields or othermeans; among the most important—and difficult—challenges is to understand how self-propelled colloid systems areaffected by hydrodynamic interactions.

8.1.1

Quincke rollers

Quincke rollers are spherical colloidal particles dispersed in a viscous fluid that exhibit self-propulsion initiated byrolling over the surface of a planar electrode, as illustrated in Fig. 8.1 [57]. The basic mechanism of this selfpropulsion is the Quincke effect, which induces the rotational motion of colloids in an external electrostatic field;
if the rotating colloids lie on the surface of an electrode slab, they experience additional fluid-mediated forces dueto the non-slip boundary conditions at the electrode surface. The Quincke effect is a phenomenon in which the application of a DC electric field, with a strength above a certain threshold, to a system of dielectric colloidal particles
dispersed in a conducting solvent induces spontaneous rotation of the particles [58]. The particles rotate about an axes
constrained to lie in a plane perpendicular to the applied electric field, and at rotational speeds proportional to thefield strength. Systems of multiple Quincke rollers experience not only hydrodynamic forces, but also electrostaticinteractions, excluded-volume effects, and other phenomena that affect the collective dynamics and result in highly
complex behavior [57, 59].

Figure 8.1: Schematic depiction of a Quincke roller. The application of a uniform DC electric field to a system of sphericalcolloidal particles produces torque on the particles due to the instability of the induced electric dipole moment. If
the field strength E0 exceeds a given threshold E Q , the colloids begin to rotate. The imposition of no-slip boundary
conditions at the electrode surface causes the colloids to move with a translational velocity v0 proportional to the
angular velocity ω of their rotation.

79

CHAPTER 8. SIMULATING QUINCKE ROLLERS

8.1.2

80

Basic equations for Quincke rollers

To analyze the motion of Quincke rollers using the SP method, we add three new two-body potentials to the equationsof motion (3.5) and (3.6):!
(Pi · r)(P j · r)
1 Pi · P j
UiDD
(r)
=
−
3
,
(8.1)
j
4πϵ
r3
r5
" 36  18 #σ
σ
UiLJj (r) = 4ϵ LJ
−
+ ϵ,
(8.2)
r
r
EF
UiEF
exp (−r/3σ)/r2 .
j (r) = −ϵ(8.3)

Here i, j are particle indices, ϵ is the dielectric permittivity, and U DD , U LJ , U EF describe interaction potentials due toelectric dipole moments, excluded-volume effects, and electro-osmotic flows, respectively, with ϵ LJ and ϵ EF setting themagnitudes of the latter two potentials. Assuming the external electric field points in the z direction, dipole moments
may be described by their parallel and perpendicular components Pz and P xy . For simplicity, we consider only steadystate Quincke phenomena. The directions of the colloid rotation axes are randomly distributed in a plane perpendicular
to the electric-field direction, and the constant rotation speed ω is an input parameter for which the relevant scale is
the rotational Reynolds number Rer ≡ ρ f ωσ2 /η. The xy components of a particle’s dipole moment are determined by
the cross product of the electric-field direction with the particle’s rotation axis. To ensure that colloid rotation axes do
not deviate from their plane of definition, we add restorative torques, in the form of harmonic potentials, to correct forany such deviation. Finally, for the magnitudes of the interaction potentials we choose ϵ LJ = ϵ EF .

CHAPTER 8. SIMULATING QUINCKE ROLLERS

8.2

Input UDF files

8.2.1

Fluid settings

81

The following choices are available for constitutive eq.
Navier Stokes:
Navier Stokes FDM:

8.2.2

Object (particle) settings

The only available value for object type.type is rigid.

Newtonian fluid
Newtonian fluid

CHAPTER 8. SIMULATING QUINCKE ROLLERS

8.3

Computational examples

8.3.1

Motion of a single Quincke roller

82

The input UDF file for this example is a2 N1 Rer025.udf in the folder Examples/12/.
$ cd Examples/12$ ../../kapsel -Ia2_N1_Rer025.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
In this example we use a 64 × 64 × 32 computational mesh and set the following parameter values: particle diameterσ = 4, interface thickness ξ = 2, rotational Reynolds number Rer = 0.25, gravitational acceleration g = 1. The python
script streamplot.py in the Examples/12 folder reads data from h5-format KAPSEL output files and generates
stream plots of fluid flow and simulation snapshot images.
Fig. 8.2 shows a comparison between the simulated self-propulsion velocity of a single Quincke roller and thepredictions of lubrication theory [60, 61] and fluid-flow field surrounding the Quincke roller.

Figure 8.2: Motion of a single Quincke roller. (a) Dependence of particle self-propulsion velocity v0 on separation distance δ
between the particle and planar electrode surfaces, showing both KAPSEL simulation results and the predictions oflubrication theory. The agreement between simulated and theoretical results improves with increasing particle radiusa. (b) Fluid-flow field surrounding the Quincke roller, with colors indicating the flow-velocity magnitudes.

8.3.2

Motion of multiple Quincke rollers

The input UDF file for this example is a2 N200.udf in the Examples/12/ folder.
$ cd Examples/12$ ../../kapsel -Ia2_N200.udf -Ooutput.udf -Ddefine.udf -Rrestart.udf
In this example (Fig. 8.3) we use a 128 × 128 × 32 computational mesh and consider N p = 200 particles, for a particle
volume fraction of φ = 0.00568 (the particle area fraction is 0.15). For the strength of interaction potentials we set
ϵ LJ = ϵ EF = 1.0. For dipole strengths we set P2xy /4π = 6.0, |Pz |/|P xy | = 3.0. Other parameter values are set as in the
single-particle example above. The initial particle configuration is generated by the python script init loc q.py in
the Examples/12/ folder. Snapshot images are generated by the python script snapshot.py.
Fig. 8.3 shows a simulation snapshot of the 200 Quincke rollers on the planar electrode surface.

CHAPTER 8. SIMULATING QUINCKE ROLLERS

Figure 8.3: Snapshot from a KAPSEL simulation of 200 Quincke rollers. The arrows indicate the rotation axes.

83Appendix A

Definition of the input UDF file format
This appendix documents the format of the input UDF files used to define KAPSEL simulations. The structure of
input UDF files is defined by a UDF definition file, conventionally named define.udf, and includes the following
items:•

Fluid settings•

Object settings•

Common simulation settings•

Configuring selectable features•

Data output settings•

Definition of UDF output data classes•

Restart settings•

GOURMET display settings

The following subsections will discuss each of these items in the order listed here, and additional details may be foundin the application-specific discussions of Sections 3-8 of this manual.

A.1

Fluid settings

The constitutive eq structure configures the properties of the fluid to be simulated. The type of fluid is specified
by selecting one of the 8 available options for the type selector:
type: select {'Navier_Stokes','Shear_Navier_Stokes','Shear_Navier_Stokes_Lees_Edwards','Electrolyte','Navier_Stokes_FDM','Navier_Stokes_Cahn_Hilliard_FDM','Shear_Navier_Stokes_Lees_Edwards_FDM','Shear_NS_LE_CH_FDM'}
Each of these fluid types has its own specific settings, as we now discuss.
Navier Stokes
Set constitutive eq to Navier Stokes to simulate particles dispersed in a Newtonian fluid (Section 3). The
Navier Stokes structure has the following fields describing fluid-related simulation parameters:
DX:
RHO:
ETA:
kBT:
alpha_v:
alpha_o:double [L] "lattice spacing (=1), fixed for all directions"double [rho] "mass density of solvent"double [eta] "shear viscosity of solvent"double [epsilon] "temperature"
double "correction coefficient of V"
double "correction coefficient of Omega"84

APPENDIX A. DEFINITION OF THE INPUT UDF FILE FORMAT

85

Here DX is the lattice spacing, RHO is the fluid density, ETA is the fluid viscosity, kBT is the particle temperature,alpha v is an adjustable parameter for the particle velocity, and alpha o is an adjustable parameter for the particle
angular velocity.
Shear Navier Stokes
Set constitutive eq to Shear Navier Stokes to simulate particles dispersed in a Newtonian fluid in the presence of zigzag shear flow. (This type of simulation is not discussed in this manual.) The Shear Navier Stokes
structure has the following fields describing fluid-related simulation parameters:
DX:
RHO:
ETA:
kBT:
alpha_v:
alpha_o:double [L] "lattice spacing (=1), fixed for all directions"double [rho] "mass density of solvent"double [eta] "shear viscosity of solvent"double [epsilon] "temperature"
double "correction coefficient of V"
double "correction coefficient of Omega"

Here DX is the lattice spacing, RHO is the fluid density, ETA is the fluid viscosity, kBT is the particle temperature,alpha v is an adjustable parameter for the particle velocity, and alpha o is an adjustable parameter for the particle
angular velocity.
The External field structure has the following fields describing zigzag shear-flow parameters:
External_field: {type: select {"DC","AC"}DC: {Shear_rate:
double [1/tau] "shear rate"}
AC: {Shear_rate:
double [1/tau] "shear rate"
Frequency:
double [1/tau] "alternating frequency"}}
Set the type selector to DC for steady-state shear flow. In this case, the Shear rate field sets the steady-state shear
velocity. Alternatively, set the type selector to AC for oscillatory shear flow. In this case, Shear rate sets theamplitude of the oscillatory shear velocity and Frequency sets the oscillation frequency.
Shear Navier Stokes Lees Edwards
Set constitutive eq to Shear Navier Stokes Lees Edwards to simulate particles dispersed in a Newtonian
fluid under shear flow with Lees-Edwards boundary conditions (Section 4). The Shear Navier Stokes Lees Edwards
structure has the following fields describing fluid-related simulation parameters:
DX:
RHO:
ETA:
kBT:
alpha_v:
alpha_o:double [L] "lattice spacing (=1), fixed for all directions"double [rho] "mass density of solvent"double [eta] "shear viscosity of solvent"double [epsilon] "temperature"
double "correction coefficient of V"
double "correction coefficient of Omega"

DX is the lattice spacing, RHO is the fluid density, ETA is the fluid viscosity, kBT is the particle temperature, alpha v
is an adjustable parameter for the particle velocity, and alpha o is an adjustable parameter for the particle angularvelocity.
The External field structure has the following fields describing shear-flow parameters with Lees-Edwards
boundary conditions:
External_field: {type: select {"DC","AC"}DC: {Shear_rate:
double [1/tau] "shear rate"}
AC: {Shear_rate:
double [1/tau] "shear rate"
Frequency:
double [1/tau] "alternating frequency"}}

APPENDIX A. DEFINITION OF THE INPUT UDF FILE FORMAT

86

Set the type selector to DC for steady-state shear flow. In this case, the Shear rate field sets the steady-state shear
velocity. Alternatively, set the type selector to AC for oscillatory shear flow. In this case, Shear rate sets theamplitude of the oscillatory shear velocity and Frequency sets the oscillation frequency.
Electrolyte
Set constitutive eq to Electrolyte to simulate charged colloidal particles dispersed in an electrolytic solution(Section 5). The Electrolyte structure has the following fields describing fluid-related simulation parameters:DX:
double [L] "lattice spacing (=1), fixed for all directions"
RHO:
double [rho] "mass density of solvent"
ETA:
double [eta] "shear viscosity of solvent"kBT:double [epsilon] "temperature"alpha_v:
double "correction coefficient of V"
alpha_o:double "correction coefficient of Omega"
Dielectric_cst:
double "dielectric constant"
INIT_profile: select {"Uniform","Poisson_Boltzmann"} "Initial condition for density profile of ions"
DX is the lattice spacing, RHO is the fluid density, ETA is the fluid viscosity, kBT is the particle temperature, alpha v isan adjustable parameter for the particle velocity, alpha o is an adjustable parameter for the particle angular velocity,
Dielectric cst is the dielectric permittivity of the solvent, and the INIT profile selector may be set to Uniform
or Poisson Boltzmann to select the method used to prepare the initial ion distribution.
The Add salt structure has the following fields regarding ions in the electrolytic solution:
Add_salt: {type:select {"saltfree","salt"}
saltfree: {Valency_counterion:
double "valency of counterion"
Onsager_coeff_counterion: double "Onsager coefficient of counterion"}
salt: {
Valency_positive_ion:
double "valency of positive ion"
Valency_negative_ion:
double "valency of negative ion"
Onsager_coeff_positive_ion: double "Onsager coefficient of positive ion"
Onsager_coeff_negative_ion: double "Onsager coefficient of negative ion"
Debye_length:
double "Debye screening length in the unit of DX"}}
The type selector may be set to saltfree for a simulation in which only counterions are present, or to salt fora simulation containing two species (one positive, one negative) of salt ions in addition to counterions.
For type=saltfree, Valency counterion and Onsager coeff counterion specify the valence and Onsager
transport coefficient for counterions.
For type=salt, the options Valency positive ion, Valency negative ion, Onsager coeff positive ion,and Onsager coeff negative ion specify valences and Osager coefficients for positive and negative ions, andDebye length specifies the Debye screening length.
The Electric field structure offers the following options to describe external electric fields:
Electric_field: {type: select {"ON","OFF"}ON: {type: select {"DC","AC"}
DC: {Ex: double
Ey: double
Ez: double}
AC: {Ex: double
Ey: double
Ez: double
Frequency:
double}}}APPENDIX A. DEFINITION OF THE INPUT UDF FILE FORMAT

87

The type selector may be set to ON or OFF for simulations with or without the application of an external electric field.
In the former case, the choice ON.type=DC selects a constant electric field with Cartesian vector components specified
by the Ex, Ey, and Ez options; the choice ON.type=AC selects an oscillating electric field with frequency set by the
Frequency option and amplitudes given by the Ex, Ey, and Ez options.
Navier Stokes FDM
Set constitutive eq to Navier Stokes FDM to simulate particles dispersed in a Newtonian fluid using finitedifference methods (Appendix D). In this case, the NS solver field of the Navier Stokes FDM structure configures
the method used to solve the Navier-Stokes equations:NS_solver: {type: select {'explicit_scheme','implicit_scheme'} "explicit_scheme: explicit MAC scheme, ON: implicit MAC scheme"
implicit_scheme: {tolerance:
double "stopping criteria"
maximum_iteration: int "number of maximum iteration"}}
Set the type selector to explicit scheme or implicit scheme to specify an explicit or implicit solver. For the
implicit case, use tolerance and maximum iteration to specify the convergence criterion and maximum numberof iterations.
The Navier Stokes FDM structure also includes the following fields describing fluid-related simulation parameters:
DX:
RHO:
ETA:
kBT:
alpha_v:
alpha_o:double [L] "lattice spacing (=1), fixed for all directions"double [rho] "mass density of solvent"double [eta] "shear viscosity of solvent"double [epsilon] "temperature"
double "correction coefficient of V"
double "correction coefficient of Omega"

Here DX is the lattice spacing, RHO is the fluid density, ETA is the fluid viscosity, kBT is the particle temperature,alpha v is an adjustable parameter for the particle velocity, and alpha o is an adjustable parameter for the particle
angular velocity.
Navier Stokes Cahn Hilliard FDM
Set constitutive eq to Navier Stokes Cahn Hilliard FDM to simulate particles dispersed in a two-component
phase-separated fluid (Section 6). In this case, the NS solver field of the Navier Stokes Cahn Hilliard FDM
structure configures the method used to solve the Navier-Stokes equations:NS_solver: {type: select {'explicit_scheme','implicit_scheme'} "explicit_scheme: explicit MAC scheme, ON: implicit MAC scheme"
implicit_scheme: {tolerance:
double "stopping criteria"
maximum_iteration: int "number of maximum iteration"
viscosity_change: select {'ON','OFF'}ON: {ETA_A:
double [eta] "shear viscosity of solvent A"
ETA_B:
double [eta] "shear viscosity of solvent B"}}}
Set the type selector to explicit scheme or implicit scheme to specify an explicit or implicit solver. For the
implicit case, use tolerance and maximum iteration to specify the convergence criterion and maximum numberof iterations. The viscosity change selector allows the two fluid components to be assigned different viscosities. If
viscosity change is set to ON, then the settings ETA A and ETA B specify separate viscosities for the A and B fluid
components. If viscosity change is set to OFF, then the single setting ETA discussed in detail later specifies thecommon viscosity of both fluid components.
The CH solver structure describes the method used by KAPSEL to solve the Cahn-Hilliard equations:

APPENDIX A. DEFINITION OF THE INPUT UDF FILE FORMAT

88

CH_solver: {type: select {'explicit_scheme','implicit_scheme'} "explicit_scheme: explicit Euler scheme, ON: implicit BDFAB scheme"
implicit_scheme: {tolerance:
double "stopping criteria"
maximum_iteration: int "number of maximum iteration"}}
Set the type selector to explicit scheme or implicit scheme to specify an explicit or implicit solver. For the
implicit case, use tolerance and maximum iteration to specify the convergence criterion and maximum numberof iterations.
The Navier Stokes Cahn Hilliard FDM structure also includes the following fields describing fluid-related
simulation parameters:
DX:
RHO:
ETA:
kBT:
alpha_v:
alpha_o:double [L] "lattice spacing (=1), fixed for all directions"double [rho] "mass density of solvent"double [eta] "shear viscosity of solvent"double [epsilon] "temperature"
double "correction coefficient of V"
double "correction coefficient of Omega"

DX is the lattice spacing, RHO is the fluid density, ETA is the fluid viscosity, kBT is the particle temperature, alpha v
is an adjustable parameter for the particle velocity, and alpha o is an adjustable parameter for the particle angularvelocity.
The Potential structure describes the phase-separation potential. Its type selector may be set to Landau or
Flory Huggins:
type: select {'Landau','Flory_Huggins'}
Selecting Landau chooses a Landau double-well potential, in which case the following settings may be configured
to specify parameter values:
composition_ratio:
initial_fluctuation:
a:
b:
d:
w:
z:
psi_0:
alpha:
kappa:

double "composition ratio of A and B fluids"
double "initial fluctuation of concentration"
double "GL parameter (third order term)"
double "GL parameter (first order term)"
double "penalty factor in fictitious particle domain"
double "penalty factor on particle surface domain"
double "penalty factor on fluid interface"double "psi_0 for particle"
double "surface parameter on fluid-fluid surface"double "mobility parameter"

Here composition ratio sets the composition ratio of the two-component phase-separated fluid, initial fluctuation
sets the initial magnitude of concentration fluctuations, a is the coefficient of the cubic term, b is the coefficient of thelinear term, d is a parameter governing the magnitude of the term responsible for ensuring that the fluid composition
inside the particle domains approaches psi 0, w is a parameter governing the affinity of particle interfaces, z is aparameter governing the affinity between particle interfaces and fluid-fluid interfaces, psi 0 is the value of ψ in theinterior of particles, alpha is a fluid interface parameter, and kappa is a mobility parameter.
On the other hand, selecting Flory Huggins for type chooses a Flory-Huggins double-well potential, in whichcase the following settings may be configured to specify parameter values:
composition_ratio:
double "composition ratio of A and B fluids"
initial_fluctuation:
double "initial fluctuation of concentration"na:double "(reduced) number of polymerization of A component"nb:double "(reduced) number of polymerization of B component"chi:double "Flory's interaction parameter (chi parameter)"d:double "penalty factor in fictitious particle domain"w:
double "penalty factor on particle surface domain"z:
double "penalty factor on fluid interface"
psi_0:
double "psi_0 for particle"alpha:
double "surface parameter on fluid-fluid surface"kappa:
double "mobility parameter"

APPENDIX A. DEFINITION OF THE INPUT UDF FILE FORMAT

89

Here composition ratio sets the composition ratio of the two-component phase-separated fluid, initial fluctuation
sets the initial magnitude of concentration fluctuations, na and nb are the degrees of polymerization of the A and B
fluid components, chi is the Flory interaction parameter, d is a parameter governing the magnitude of the term responsible for ensuring that the fluid composition inside the particle domains approaches psi 0, w is a parameter governing
the affinity between fluid and particle interfaces , z is a parameter governing the affinity between particle interfaces
and fluid-fluid interfaces, psi 0 is the value of ψ in the interior of particles, alpha is a fluid interface parameter, andkappa is a mobility parameter.
The Wall Potential structure describes planar walls:
type: select {'ON', 'OFF'}ON:{w :double "penalty factor on wall surface domain"
psi_0: {magnitude:
double "psi_0 magnitude"
profile : select {'uniform', 'user_specify'}
user_specify:{
PSI0[][]:{
value : double}}}
DRYING:{
type: select {'ON', 'OFF'}ON:{
psi_dry : double "psi_dry"}}}
Setting the type selector to ON enables planar walls. In this case, w is a parameter governing the affinity between
planar walls and fluids, and psi 0 controls the value of ψ in the interior of planar walls. psi 0 contains the field’smagnitude, specifying the magnitude of the affinity between planar walls and fluids, and profile, which may be set
to uniform (uniform affinity over surfaces of planar walls) or user specify (user-specified variation of affinity overplanar walls). If profile is set to user specify, use PSI0[][] to configure values by hand. Within the DRYING
structure, set the type selector to ON to specify a drying simulation. In this case, psi dry is a parameter governing
the speed of drying.
Shear Navier Stokes Lees Edwards FDM
Set constitutive eq to Shear Navier Stokes Lees Edwards FDM for finite-difference simulations of particles dispersed in a Newtonian fluid under shear flow with Lees-Edwards boundary conditions (Appendix D). In thiscase, the NS solver field of the Shear Navier Stokes Lees Edwards FDM structure configures the method used tosolve the Navier-Stokes equations:NS_solver: {type: select {'explicit_scheme','implicit_scheme'} "explicit_scheme: explicit MAC scheme, ON: implicit MAC scheme"
implicit_scheme: {tolerance:
double "stopping criteria"
maximum_iteration: int "number of maximum iteration"}}
Set the type selector to explicit scheme or implicit scheme to specify an explicit or implicit solver. For the
implicit case, use tolerance and maximum iteration to specify the convergence criterion and maximum numberof iterations.
The Shear Navier Stokes Lees Edwards FDM structure also includes the following fields describing fluidrelated simulation parameters:
DX:
RHO:
ETA:
kBT:
alpha_v:
alpha_o:double [L] "lattice spacing (=1), fixed for all directions"double [rho] "mass density of solvent"double [eta] "shear viscosity of solvent"double [epsilon] "temperature"
double "correction coefficient of V"
double "correction coefficient of Omega"

APPENDIX A. DEFINITION OF THE INPUT UDF FILE FORMAT

90

Here DX is the lattice spacing, RHO is the fluid density, ETA is the fluid viscosity, kBT is the particle temperature,alpha v is an adjustable parameter for the particle velocity, and alpha o is an adjustable parameter for the particle
angular velocity.
The External field structure contains parameters affecting shear flows under Lees-Edwards boundary conditions:
External_field: {type: select {"DC","AC"}DC: {Shear_rate:
double [1/tau] "shear rate"}
AC: {Shear_rate:
double [1/tau] "shear rate"
Frequency:
double [1/tau] "alternating frequency"}}
Set the type selector to DC for steady-state shear flow. In this case, the Shear rate field sets the steady-state shear
velocity. Alternatively, set the type selector to AC for oscillatory shear flow. In this case, Shear rate sets theamplitude of the oscillatory shear velocity and Frequency sets the oscillation frequency.
Shear NS LE CH FDM
Set constitutive eq to Shear NS LE CH FDM to simulate particles dispersed in two-component phase-separated
fluids under shear flow with Lees-Edwards boundary conditions (Section 6). In this case, the NS solver field of theShear NS LE CH FDM structure configures the method used to solve the Navier-Stokes equations:NS_solver: {type: select {'explicit_scheme','implicit_scheme'} "explicit_scheme: explicit MAC scheme, ON: implicit MAC scheme"
implicit_scheme: {tolerance:
double "stopping criteria"
maximum_iteration: int "number of maximum iteration"
viscosity_change: select {'ON','OFF'}ON: {ETA_A:
double [eta] "shear viscosity of solvent A"
ETA_B:
double [eta] "shear viscosity of solvent B"}}}
Set the type selector to explicit scheme or implicit scheme to specify an explicit or implicit solver. For the
implicit case, use tolerance and maximum iteration to specify the convergence criterion and maximum numberof iterations. The viscosity change selector allows the two fluid components to be assigned different viscosities. If
viscosity change is set to ON, then the settings ETA A and ETA B specify separate viscosities for the A and B fluid
components. If viscosity change is set to OFF, then the single setting ETA discussed in detail later specifies thecommon viscosity of both fluid components.
The CH solver structure describes the method used by KAPSEL to solve the Cahn-Hilliard equations:
CH_solver: {type: select {'explicit_scheme','implicit_scheme'} "explicit_scheme: explicit Euler scheme, ON: implicit BDFAB scheme"
implicit_scheme: {tolerance:
double "stopping criteria"
maximum_iteration: int "number of maximum iteration"}}
Set the type selector to explicit scheme or implicit scheme to specify an explicit or implicit solver. For the
implicit case, use tolerance and maximum iteration to specify the convergence criterion and maximum numberof iterations.
The Shear NS LE CH FDM structure also includes the following fields describing fluid-related simulation parameters:

APPENDIX A. DEFINITION OF THE INPUT UDF FILE FORMAT

DX:
RHO:
ETA:
kBT:
alpha_v:
alpha_o:91

double [L] "lattice spacing (=1), fixed for all directions"double [rho] "mass density of solvent"double [eta] "shear viscosity of solvent"double [epsilon] "temperature"
double "correction coefficient of V"
double "correction coefficient of Omega"

Here DX is the lattice spacing, RHO is the fluid density, ETA is the fluid viscosity, kBT is the particle temperature,alpha v is an adjustable parameter for the particle velocity, and alpha o is an adjustable parameter for the particle
angular velocity.
The Potential structure describes the phase-separation potential. Its type selector may be set to Landau or
Flory Huggins:
type: select {'Landau','Flory_Huggins'}
Selecting Landau chooses a Landau double-well potential, in which case the following settings may be configured
to specify parameter values:
composition_ratio:
double "composition ratio of A and B fluids"
initial_fluctuation:
double "initial fluctuation of concentration"a:double "GL parameter (third order term)"b:double "GL parameter (first order term)"d:double "penalty factor in fictitious particle domain"w:
double "penalty factor on particle surface domain"z:
double "penalty factor on fluid interface"
psi_0: double "psi_0 for particle"alpha: double "surface parameter on fluid-fluid surface"kappa: double "mobility parameter"
composition ratio sets the composition ratio of the two-component phase-separated fluid, initial fluctuation
sets the initial magnitude of concentration fluctuations, a is the coefficient of the cubic term, b is the coefficient of thelinear term, d is a parameter governing the magnitude of the term responsible for ensuring that the fluid composition
inside the particle domains approaches psi 0, w is a parameter governing the affinity of particle interfaces, z is aparameter governing the affinity between particle interfaces and fluid-fluid interfaces, psi 0 is the value of ψ in theinterior of particles, alpha is a fluid interface parameter, and kappa is a mobility parameter.
On the other hand, selecting Flory Huggins for type chooses a Flory-Huggins double-well potential, in whichcase the following settings may be configured to specify parameter values:
composition_ratio:
double "composition ratio of A and B fluids"
initial_fluctuation:
double "initial fluctuation of concentration"na:double "(reduced) number of polymerization of A component"nb:double "(reduced) number of polymerization of B component"chi:double "Flory's interaction parameter (chi parameter)"d:double "penalty factor in fictitious particle domain"w:
double "penalty factor on particle surface domain"z:
double "penalty factor on fluid interface"
psi_0: double "psi_0 for particle"alpha: double "surface parameter on fluid-fluid surface"kappa: double "mobility parameter"
Here composition ratio sets the composition ratio of the two-component phase-separated fluid, initial fluctuation
sets the initial magnitude of concentration fluctuations, na and nb are the degrees of polymerization of the A and B
fluid components, chi is the Flory interaction parameter, d is a parameter governing the magnitude of the term responsible for ensuring that the fluid composition inside the particle domains approaches psi 0, w is a parameter governing
the affinity between fluid and particle interfaces , z is a parameter governing the affinity between particle interfaces
and fluid-fluid interfaces, psi 0 is the value of ψ in the interior of particles, alpha is a fluid interface parameter, andkappa is a mobility parameter.
The External field structure has the following fields describing shear-flow parameters under Lees-Edwards
boundary conditions:
External_field: {type: select {"DC","AC"}DC: {Shear_rate:
double [1/tau] "shear rate"}
AC: {APPENDIX A. DEFINITION OF THE INPUT UDF FILE FORMAT

Shear_rate:
Frequency:

92

double [1/tau] "shear rate"
double [1/tau] "alternating frequency"}}
Set the type selector to DC for steady-state shear flow. In this case, the Shear rate field sets the steady-state shear
velocity. Alternatively, set the type selector to AC for oscillatory shear flow. In this case, Shear rate sets theamplitude of the oscillatory shear velocity and Frequency sets the oscillation frequency.

A.2

Object settings

The object type structure describes the properties of objects dispersed in fluids. Its type selector defines the typeof object:
type: select {'spherical_particle','chain','rigid'}

spherical particle
Set object type.type to spherical particle to simulate spherical particles dispersed in fluids. In this casethe following parameters may be configured:
Particle_spec[]:{Particle_number: int "number of colloidal particles"
MASS_RATIO: double "mass density ratio colloid/solvent"Surface_charge: double "surface charge of colloid"janus_axis: select {'NONE', 'X', 'Y', 'Z'} "janus axis in body_fixed frame"janus_propulsion: select{'OFF', 'TUMBLER', 'SQUIRMER', 'OBSTACLE'}
janus_force: Vector3d "self-propulsion force"
janus_torque: Vector3d "self-propulsion torque"janus_slip_vel: float "Slip velocity coeff B1"
janus_slip_mode: float "Blake squirmer mode B2/B1"
janus_rotlet_C1: float "rotlet coefficient C1"
janus_rotlet_dipole_C2: float "rotlet dipole C2"}
Particle number sets the number of particles. MASS RATIO specifies the ratio of particle density to fluid density. Surface charge sets the particle surface charge; this is only applicable if constitutive eq.type is set toElectrolyte.
Properties of the form janus * specify properties of Janus particles. The janus axis selector may be set to NONE,
X, Y, or Z to set the orientation of the Janus axis in a particle-fixed coordinate system. The propulsive motion of Janus
particles is specified by the janus propulsion selector, which may be set to OFF (disabled), TUMBLER (particle propelled along propulsion axis by fixed external force), SQUIRMER (squirmer particle propelled along propulsion axis byslip boundary conditions), or OBSTACLE (particle forming a fixed obstacle). The Cartesian components of the propulsive force and torque on Janus particles may be specified by setting values for janus force.x, janus force.y,janus force.z, janus torque.x, janus torque.y, and janus torque.z. janus slip vel sets the parameterB1 determining the surface-slip velocity for Janus particles, while janus slip mode sets the parameter B2 /B1 determining the type of slip motion for Janus-particles; finally, janus rotlet C1 and janus rotlet dipole C2 are theJanus-particle parameters C1 and C2 describing rotlet dipole exerted on particles/fluids.
chain
Set object type.type to chain to simulate flexible particle chains dispersed in fluids. In this case the following
parameter settings are available:
chain:{
Chain_spec[]:{Beads_number: int "number of beads in a chain"
Chain_number: int "number of chains"
MASS_RATIO: double "mass density ratio chain/solvent"Surface_charge: double "surface charge of colloid"janus_axis: select {'NONE', 'X', 'Y', 'Z'} "janus axis in body_fixed frame"}}

APPENDIX A. DEFINITION OF THE INPUT UDF FILE FORMAT

93

Beads number is the number of beads associated with a single chain, while Chain number is the number of chains.
MASS RATIO is the ratio of bead density to fluid density. Surface charge sets the particle surface charge; this is only
applicable if constitutive eq.type is set to Electrolyte. The janus axis selector may be set to NONE, X, Y, or
Z to specify the orientation of the Janus axis in a bead-fixed coordinate system.
rigid
Set object type.type to rigid to simulate rigid bodies comprised of spherical beads dispersed in fluids. In thiscase the following parameter settings are available:
rigid:{
Rigid_spec[]:{Beads_number: int "number of beads in a rigid"
Rigid_number: int "number of rigids"
MASS_RATIO: double "mass density ratio rigid/solvent"Surface_charge: double "surface charge of particle"
Rigid_motion: select {'fix','free'}Rigid_velocity: Vector3d "speed of translation ### fix only ###"
Rigid_omega: Vector3d "angular velocity ### fix only ###"}}
Here Beads number is the number of beads associated with a rigid body, while Rigid number is the number of rigid
bodies. MASS RATIO is the ratio of bead density to fluid density. Surface charge sets the particle surface charge;
this is only applicable if constitutive eq.type is set to Electrolyte. The motion of rigid bodies is described by
the Rigid motion selector, which may be set to free (for free motion) or fix (for motion with fixed translationaland rotational velocities). For the choice Rigid motion=fix, all six the Cartesian components of the translationaland rotational velocities of rigid bodies (in the lab frame) are specified via Rigid velocity.x, Rigid velocity.y,
Rigid velocity.z, Rigid omega.x, Rigid omega.y, and Rigid omega.z. To free/fix individual degrees of freedom use this option together with the switch.free rigid option detailed below ‘Configuring selectable features’.

A.3

Common simulation settings

Sections A.1 and A.2 described configurable settings for fluids and objects. In this section we describe common
simulation settings applicable in general.
Objects in KAPSEL simulations are comprised of particles. The following parameters affect common propertiesof these primary particles:
A_XI: double "interface thickness in the unit of DX"
A: double "colloid radius in the unit of DX"
A XI and A respectively specify the particle interface thickness and the particle radius in units of the lattice spacingDX.
The gravity structure describes gravitational forces on the system simulated by KAPSEL :gravity: {G: double [L*tauˆ-2] "gravitational acceleration constant"
G_direction: select {'-X','-Y','-Z'} "direction of gravitational acceleration"}
Here G is the gravitational acceleration and the G direction selector may be set to -X, -Y, or -Z to specify the
direction of gravitational forces.
The following common parameters affect interparticle interaction potentials:
EPSILON: double [epsilon] "Lennard-Jones depth"
LJ_powers: select {'12:6','24:12','36:18','macro_vdw','electro_osmotic_flow'} "set of,→
power exponents of LJ potential"
EPSILON specifies the unit of energy in the Lennard-Jones potential, while LJ powers may be set to one of thefollowing 5 values:•

12:6, 24:12 or 36:18 to select a Lennard-Jones potential with the given exponents,•

macro vdw to select a macroscopic interparticle potential, orAPPENDIX A. DEFINITION OF THE INPUT UDF FILE FORMAT•

94electro osmotic flow, a choice used only for Quincke rollers.

The mesh structure describes the computational mesh used for fluid calculations:mesh: {
NPX: int "number of mesh in x-direction = 2ˆNPX"
NPY: int "number of mesh in y-direction = 2ˆNPY"
NPZ: int "number of mesh in z-direction = 2ˆNPZ"}
Here NPX, NPY, and NPZ are the base-2 logarithms of the numbers of mesh points in the x, y, and z directions; that is,
the numbers of mesh points are L x = 2NPX , Ly = 2NPY , and Lz = 2NPZ .
Timestamp settings
The time increment structure describes simulation timesteps:
time_increment: {type: select {"auto","manual"}
auto: {
factor: double "delta_t = factor * h(determined by system parameters)"}
manual: {delta_t: double [tau]
}
}
2
where
Set the type selector to auto to set the timestamp to its maximum allowed value, given by T step = ρ/ηkmax
kmax is the maximum wavenumber determined by the lattice spacing ∆. In this case, factor specifies a multiplicative
scale factor, so the actual simulation timestep is given by ∆t = T step × factor. Alternatively, set the type selector to
manual and set delta t to your preferred timestep value.

A.4

Configuring selectable features

This section describes features defined in the switch structure.
Set the ROTATION selector to ON or OFF to include or exclude particle rotational motion from consideration in
KAPSEL simulations.
ROTATION: select {'ON','OFF'} "OFF: not solve rotation, ON: solve rotation"
Use the LJ truncate selector to specify the form of the Lennard-Jones potential acting between particles. Set to
OFF for the usual form of the potential, including the attractive term. Set to ON to exclude the attractive term, retaining
only the repulsive term. Set to NONE to exclude both terms, i.e., to disable the interparticle potential entirely.
LJ_truncate: select {'ON','OFF','NONE'} "OFF:normal LJ, ON:WCA, NONE: no-interaction,→
at all"The INIT distribution structure specifies the initial configuration of particles.
INIT_distribution: {type: select {'uniform_random','random_walk','FCC','BCC','user_specify'}"uniform_random:distributed uniformly in box, random_walk:distributed uniformly inbox, FCC:distributed on FCC lattice, BCC:distributed on BCC lattice,,→,→user_specify:configuration and velocity specified by user"
random_walk: {iteration: int}
user_specify: {Particles[]: Particle}}

APPENDIX A. DEFINITION OF THE INPUT UDF FILE FORMAT

95

Allowed values for the type selector are uniform random (random), random walk (particles randomly moved fromthe sites of a square lattice), FCC (FCC lattice), BCC (BCC lattice), or user specify (user-specified coordinatesand velocities). For random walk, iteration specifies the number of random-walk iterations. For user specify,
Particles[] sets the particle positions and velocities.
Similarly, the INIT orientation selector describes the initial orientation of particles:
INIT_orientation: select {'user_specify', 'random', 'space_align'}
The allowed values are user specify (coordinates and velocities specified by user), random (random), or space align(particles aligned with the lab-frace axes). For the case user specify, initial particle orientations are specified via
Particles[].
SLIP tol and SLIP iter set the convergence criterion and maximum iteration count for iterative calculation of
fluid-flow fields when introducing slip velocities in the tangential direction at Janus particle interfaces:
SLIP_tol: float "Tolerance for iterative slip convergence"
SLIP_iter: int "Maximum number of iterations for iterative slip convergence"
The FIX CELL structure configures properties of the simulation cells:
FIX_CELL: {x: select{'ON','OFF'}"OFF:w/o DC current, ON:with DC current"
y: select{'ON','OFF'}"OFF:w/o DC current, ON:with DC current"
z: select{'ON','OFF'}"OFF:w/o DC current, ON:with DC current"}
Set the x selector to ON to zero out the DC component of the total velocity in the x direction, or OFF to omit thisadjustment (and similarly for the y and z selectors).
The pin structure configures particle degrees of freedom:
pin: {
type: select{"NO","YES"}
YES:{
pin[]: int
pin_rot[]: int}}
Set the type selector to YES to fix (“pin”) particle positions. In this case, the pin[] / pin rot[] fields lists the indicesof particles whose translational/rotational motion is to be prohibited.
The free rigid structure configures individual degrees of freedom (DOFs) for the translational/rotational motionof rigid bodies. This should be used together with the object type.rigid.Rigid spec option, which will fix all
six components of the linear and angular velocities. free rigid can then be used to selectively free any of the DOFs.
free_rigid:{
type: select{'NO', 'YES'} "Free rigid degrees of freedom"
YES:{
DOF[]:{
spec_id: int"Rigid body species id"vel:{
x:select{'NO', 'YES'},y:select{'NO', 'YES'},z:select{'NO', 'YES'}} "Free velocity components"
omega:{
x:select{'NO', 'YES'},y:select{'NO', 'YES'},z:select{'NO', 'YES'}} "Free omega components"}}}
Set the type selector to YES to allow DOFs to be configured separately for each rigid body in the DOF[] field. In
this case, spec id specifies the index of a rigid body, and the x, y, and z selectors of the vel and omega structures
may be set to YES (free) or NO (fixed) to enable or disable translational and rotational motion in the various Cartesian
directions.

APPENDIX A. DEFINITION OF THE INPUT UDF FILE FORMAT

96

The ns solver structure configures the method used to solve the Navier-Stokes equations (for shear-flow simulations under Lees-Edwards boundary conditions):
ns_solver:{OBL_INT: select {'linear', 'spline'} "interpolation scheme for Oblique/Rectangular,→transform"}
The OBL INT selector may be set to linear or spline to select the approximation method used for the coordinatetransformations between the cartesian and oblique grids in the shear-flow simulations.
The wall structure specifies properties of planar walls:
wall:{
type: select{'NONE', 'FLAT'}
FLAT: {
axis: select{'X', 'Y', 'Z'}"perpendicular axis to flat parallel walls"
DH: int "wall thickness in number of grid points"
LJ_Params: select{"AUTO", "MANUAL"}MANUAL:{
truncate: select{'ON', 'OFF', 'NONE'} "Truncate OFF: attractive force, ON: no,→
attractive force"
powers: select{'12:6','24:12','36:18'} "type of LJ potential"
EPSILON: double "LJ parameter, default (Basically, a large value.)"}}}
Set the type selector to FLAT to specify planar walls. In this case, set the axis selector to X, Y, or Z to specify thenormal direction to the planar walls. The DH field specifies the thickness of the planar walls, measured in units of thelattice spacing. The LJ Params selector describes the potential acting at planar walls; the choice LJ Params=AUTO
sets the potential to be the same as the interparticle potential, while the choice LJ Params=MANUAL allows a different
potential to be specified. For the MANUAL case, the truncate selector specifies the presence or absence of attractive
forces at planar walls. If set to ON, the potential is purely repulsive. If set to OFF, the potential is a sum of repulsive
and attractive contributions. The powers selector may be set to 12:6, 24:12, or 36:18 to specify the exponents in the
Lennard-Jones potential at planar walls, and EPSILON sets the energy unit for the Lennard-Jones potential at planar
walls.
The quincke structure configures properties of Quincke rollers:quincke:{type: select{'ON','OFF'}
ON: {e_dir: select{'X', 'Y', 'Z'}"constraint axis(the direction of external electric,→field E)"w_dir: select{'X', 'Y', 'Z'}"the direction of constant angular velocity vector(body,→frame)"
torque_amp: double "the amplitude of constraint torque"}}
Set the type selector to ON to simulate Quincke rollers. In this case, set the e dir selector to X, Y, or Z to specify the
direction of the external electric field, set the w dir selector to X, Y, or Z to specify the direction of the angular-velocity
vector associated with rotation induced by the Quincke effect, and set torque amp to the torque magnitude.
The multipole structure configures Ewald-method simulations:multipole:{
type: select{'ON', 'OFF'}ON:{
Dipole:{type : select{'ON', 'OFF'}ON:{
magnitude: double "the dipole strength"
type: select{'FIXED', 'QUINCKE'} "type of dipole"
FIXED:{
dir: select{'X', 'Y', 'Z'} "direction of dipolar axis"}
QUINCKE:{
type : select{'with_mirror_image','no_mirror_image'} "Mirror image component,→of electrode surface"
Pz_factor : double "Pz strength += Pz_factor*magnitude"

APPENDIX A. DEFINITION OF THE INPUT UDF FILE FORMAT

97

}
}
}
EwaldParams:{alpha: double "Ewald screening parameter"delta: double "Tolerance parameter, used to determine k_max"
converge: double "Convergence parameter (fraction of k vectors to consider)"epsilon : double "Permittivity at boundary ( if negative, set to tinfoil)"}}}
Set the type selector to ON to run simulations using the Ewald method. The Dipole structure describes properties
of electric dipoles; set its type selector to ON to enable dipoles. In this case, magnitude sets the magnitude of thedipole moment, while the type selector may be set to FIXED for fixed dipoles or QUINCKE for Quincke rollers. For
type=FIXED, dir may be set to X, Y, or Z to specify the direction of the fixed dipoles. For type=QUINCKE, the typeselector may be set to with mirror image or no mirror image to enable or disable mirror-image components ofthe electrode surface, while Pz factor is a parameter specifying the magnitude of the z-component of the dipole
moment. Finally, the EwaldParams structure specifies Ewald-method parameters: alpha is the screening parameter,delta is a convergence criterion for determination of kmax , converge is a convergence parameter, and epsilon is thedielectric permittivity of the surrounding medium (negative values are used to specify tin-foil boundary conditions,
ϵ = ∞).

A.5

Data output settings

The output structure contains fields that may be used to configure KAPSEL’s data output settings.
output: {GTS: int "interval between snapshots"
Num_snap: int "number of snapshots"
AVS: select {"ON","OFF"}ON:{
Out_dir: string "directory name"Out_name: string "prefix name for data file"
FileType: select {"BINARY","ASCII","EXTENDED"} "output data type"
EXTENDED:{
Driver:{Format: select {"HDF5"}}
Print_field:{
Crop:select{"YES", "NO"}"Crop Field Data to Hyperslab"
YES:{
Slab_x: SlabSelection
Slab_y: SlabSelection
Slab_z: SlabSelection}
Vel: select{"YES", "NO"}"Print velocity field"
Phi: select{"YES", "NO"}"Print phi field"
Charge: select{"YES", "NO"}"Print charge fields (surface & solute charge &,→potential)"
Pressure: select{"YES", "NO"} "Print pressure field"
Tau: select{"YES", "NO"}"Print stress tensor"}}}
UDF: select {"ON","OFF"}}GTS sets the data output interval measured in number of steps. Num snap sets the number of data outputs. Thus, thetotal number of timesteps in the simulation is given by GTS × Num snap. Set the AVS selector to ON to enable data
output in AVS format. For AVS=ON, the Out dir and Out name options set the directory in which AVS-format data
output files are written and the filename prefix for those files. The format of AVS data output files is determined by
the FileType selector, which may be set to Binary, ASCII, or EXTENDED. The choice EXTENDED specifies HDF5 as
the extended data format in KAPSEL. The Print field structure specifies properties of the output data fields. Set
the Crop selector to YES to reduce the number of output data fields. Use Slab x, Slab y, and Slab z to specify
the reduction in the x, y, and z directions. Use start, stride, and count to set the first lattice-point index, theinterval (number of lattice points) between lattice-point outputs, and the number of lattice points to output. Set theAPPENDIX A. DEFINITION OF THE INPUT UDF FILE FORMAT

98

Vel selector to YES to output values of velocity fields. Set the Phi selector to YES to output values of the SP function
ϕ. Set the Charge selector to YES to output the charge density distribution (only applicable if the type selector forconstitutive eq is set to Electrolyte). Set the Pressure selector to YES to output values of pressure fields.
(Currently not implemented). Set the Tau selector to YES to output values of the stress tensor. Finally, set the UDF
selector to ON to enable UDF output.

A.6

Definition of UDF output data classes

The following UDF output data classes are defined in define.udf. As these are not configurable settings, we omit
detailed descriptions.
\begin{def}class outParticle:{
R:Vector3d [L],R_raw:Vector3d [L],v:Vector3d [L/tau],
q:Quaternion,omega:Vector3d,f_hydro:Vector3d [mass*L*tauˆ-2],
torque_hydro:Vector3d,
f_r:Vector3d [mass*L*tauˆ-2],
torque_r:Vector3d,
f_slip:Vector3d [mass*L*tauˆ-2],
torque_slip:Vector3d}
E: float [epsilon] "total kinetic energy of the system"
t: float "total time"
Particles[]: outParticle
RigidParticles[]: outParticle
PSI[][][]: {psi: float}\end{def}\begin{def}class sParticle:{
R:Vector3d [L],R_raw:Vector3d [L],v:Vector3d [L/tau],v_old:Vector3d [L/tau],f_hydro:Vector3d [mass*L*tauˆ-2],
f_hydro_previous:Vector3d [mass*L*tauˆ-2],
f_hydro1:Vector3d [mass*L*tauˆ-2],
f_slip:Vector3d [mass*L*tauˆ-2],
f_slip_previous:Vector3d [mass*L*tauˆ-2],
fr:Vector3d [mass*L*tauˆ-2],
fr_previous:Vector3d [mass*L*tauˆ-2],
omega:Vector3d,
omega_old:Vector3d,
torque_hydro:Vector3d,
torque_hydro_previous:Vector3d,
torque_hydro1:Vector3d,
torque_slip:Vector3d,
torque_slip_previous:Vector3d,
torque_r:Vector3d,
torque_r_previous:Vector3d,
q:Quaternion,q_old:Quaternion}
class Matrix3d:{
xx: float,
xy: float,
xz: float,yx: float,
yy: float,
yz: float,zx: float,zy: float,zz: float}
class CTime:{
ts:inttime:float [tau]APPENDIX A. DEFINITION OF THE INPUT UDF FILE FORMAT

99}\end{def}

A.7

Restart settings

The resume structure contains settings that control restart settings for resuming calculations. The Calculation
selector describes how to proceed from a previous calculation:
Calculation: select {'NEW','CONTINUE','CONTINUE_FDM','CONTINUE_FDM_PHASE_SEPARATION'} "flg in order to specify resumed simulation or not"The choice NEW starts a new calculation. The remaining three choices (CONTINUE, CONTINUE FDM, and CONTINUE FDM PHASE SEPARATI
specifies that data saved upon termination of previous calculation should be read in and the calculation resumed. Whichof the three CONTINUE options to use depends on the type of simulation you are running:•

Use CONTINUE for simulations using spectral methods, i.e., constitutive eq set to Navier Stokes, Shear Navier Stokes,
Shear Navier Stokes Lees Edwards, or Electrolyte.

•

Use CONTINUE FDM for simulations of single-component fluids via finite-difference methods, i.e., constitutive eq
set to Navier Stokes FDM or Shear Navier Stokes Lees Edwards FDM.

•

Use CONTINUE FDM PHASE SEPARATION for simulations of two-component fluids via finite-difference methods,i.e. constitutive eq set to Navier Stokes Cahn Hilliard FDM or Shear NS LE CH FDM.

In each case, results needed to resume calculations are stored under Saved Data:
CONTINUE:{
Saved_Data:{
jikan: CTime
Particles[] : sParticle
GR_body[]: Vector3d
GR_masses[] : float
GR_moments_body[]: Matrix3d
Zeta[][][]:{
zeta0: float
zeta1: float}
uk_dc: Vector3d
Concentration[][][][]: {ck:float}
oblique: {degree_oblique: float}}}
CONTINUE_FDM: {Saved_Data: {jikan: CTime
Particles[]: sParticle
GR_body[]: Vector3d
GR_masses[]: float
GR_moments_body[]: Matrix3d
U[][][]: {u0: float
u1: float
u2: float}
U_OLD[][][]: {u_old_0: floatu_old_1: floatu_old_2: float}
oblique: {degree_oblique: float}}

APPENDIX A. DEFINITION OF THE INPUT UDF FILE FORMAT

100}
CONTINUE_FDM_PHASE_SEPARATION: {Saved_Data: {jikan: CTime
Particles[]: sParticle
GR_body[]: Vector3d
GR_masses[]: float
GR_moments_body[]: Matrix3d
U[][][]: {u0: float
u1: float
u2: float}
U_OLD[][][]: {u_old_0: floatu_old_1: floatu_old_2: float}
PSI[][][]: {psi: float}
PSI_OLD[][][]: {psi_old: float}
STRESS_OLD[][][]: {stress_old_0: float
stress_old_1: float
stress_old_2: float}
oblique: {degree_oblique: float}}}
As these settings are not configured manually, we omit detailed descriptions.

A.8

GOURMET display settings

The Unit Parameter structure contains fields used to configure settings relevant for GOURMET displays. Thesesettings do not affect simulations or simulation results.
Unit_Parameter:{Name: string"Name"
Comment:string"Comment"
Temperature:double [K] "Temperature"
Length:double [nm]"Grid spacing"
rho:double [g/cmˆ3]"Fluid density"} "Parameters for unit conversion"
Name sets the name of the UDF file. Comment adds a comment to the UDF file. Temperature specifies a temperatureunit used as a reference when simulation parameters are reported in actual units in GOURMET. Length and rho
similarly specify length and density units.

Appendix B

Particle calculations
Particle configurations (positions, orientations, velocities, and angular velocities) are updated using a 2-step AdamsBashforth scheme, with the exception of the initial step for which Euler’s method is used. Denoting the dynamic
variables of interest by YI and putting YIn = YI (tn ), we haveY n+1 = Y n +h n
3Ẏ − Ẏ n−1
2

(B.1)

In the body-fixed (primary-axis) reference coordinate system (in which the moment-of-intertia tensor Ie is diagonal),
we solve Euler’s equation in the form [8] 1 ωe˙   e
e2 ω
e3 (e
τ1 + ω
I 22 − e
I 33 )  ˙ 2   2
e3 ω
e1 (e
τ +ω
I 33 − e
I 11 ) 
ωe  =  e 3   3
e
e1 ω
e2 (e
τ +ω
I 11 − e
I 22 ).
e˙
ω

(B.2)

For better precision in numerical calculations, particle orientations are updated using a rotational quaternion q, normalized at each timestep, instead of the rotation matrix R; the dynamical equation for the orientation then reads [62].
q̇ =

1
1
e
w◦q= q◦w
2
2

(B.3)

e = (0, ωe ) and ω = (0, ω) are quaternions corresponding
Here ◦ denotes quaternion multiplication and the quantities w
to the angular velocity in the laboratory system (ω) and in the body-fixed reference coordinate system (eω = Rt · ω).

101Appendix C

Fluid calculations by spectral methods
As noted above in our discussion of simulation methods, KAPSEL uses the fractional-step method. To begin, thismethod determines u∗ by time-evolving a modified version of the Navier-Stokes equations in which the advection
and viscous-stress terms are present but the ϕf p term is omitted. To simplify computations, KAPSEL does not solve
directly for u, but rather for ω = ∇ × u. In a fixed Cartesian coordinate system with periodic boundary conditions,the vorticity equations satisfying the incompressibility condition (∇ · u = 0) read∂t ω = −∇ × ∇ · (uu) + ρ−1 ∇ × ∇ · σ,b = k × [k · Fk (uu)] − ρ−1 k × k · σ
b
∂t ω

(C.1)

where Fk ( f ) = b
f (k) denotes the Fourier transform and k is a wave-vector. The fluid velocity field is obtained fromthe definition of vorticity and the incompressibility condition:
b=i
u

b
k×ω
|k|2

(C.2)

The advantage of solving an equation for the vorticity ω is that there is no need to solve the pressure Poisson equation
to ensure incompressibility. Instead, incompressibility is satisfied via a Fourier-space projection of the form"#kkb −→ I − 2 · u
b.
u
(C.3)
|k|
To reduce the amount of memory needed to store program code, KAPSEL further exploits the fact that the threevorticity components are not linearly independent, whereupon it suffices to consider only two components. Denotingb the Fourier transform of the curl of an arbitrary vector field A, for all wavevectors (withby b
a = Fk (∇ × A) = ik × A
k , 0) we define an invertible map († : C 3 → C 2 ) of the form


(b
a2 ,b
a3 ) k1 , 0







b
a1 = −(k3ba3 + k2ba2 )/k1






(b
a
,b
a
)
k
=
0,
k
,
0
†
 3 1
1
2
b
a = (ba1 ,ba2 ,ba3 )† : 
(C.4)


b
b
a
=
−k
a3 /k2

2
3






(b
a1 ,b
a2 ) k1 = k2 = 0, k3 , 0





b
a3 = 0Because ba(k = 0)—which corresponds to a volume integral over the entire vector field—cannot be determined frombb† ; this reduces bya† , we compute it separately. Applying this mapping to equations (C.1) yields a formula for ζb = ω
1/3 the computational cost of updating the velocity field.
The form of the equations that are integrated depends on the stress tensor σ that is used. For a Newtonian fluidwith ∇ · σ = η∇2 u, the equation for the vorticity ω readsb = −νk2 ω
b + k × [k · Fk (uu)] .
∂t ω

(C.5)

∂t a = La − N(t, a(t))(C.6)

We have expressed this in the form

102

APPENDIX C. FLUID CALCULATIONS BY SPECTRAL METHODS

103

where L is a time-independent linear operator and N is a nonlinear operator. The general solution of this equation isa(tn + h) = e a(tn ) − e
Lh

Z h

Lh

dτe−Lτ N(tn + τ, a(tn + τ))

(C.7)

0

and the linear portion may be solved exactly. The integral of the nonlinear portion may be approximated by various
methods [63, 64]. In particular, the first-order approximation N(tn + τ, a(tn + τ)) = N(tn , a(tn )) yieldsh
ia(tn + h) = eLh a(tn ) − L−1 1 − e−hL N(tn , a(tn ))
(C.8)



= a(tn ) + e−hL − 1 a(tn ) + L−1 N(tn , a(tn ))(C.9)
which reduces to Euler’s method in the limit |L| → 0.

Appendix D

Fluid calculations by finite-difference
methods (FDM)
For simulations of particles dispersed in two-component phase-separated fluids, discussed in Section 6, KAPSEL
uses finite-difference methods for fluid calculations. This appendix discusses the implementations of these methods.

D.1

Solving the Navier-Stokes equations

KAPSEL solves the Navier-Stokes (NS) equations (6.1) using explicit and implicit marker-and-cell (MAC) solvers [43,46]. Variables are arranged on the computational mesh in a semi-staggered Arakawa-B lattice configuration, with velocity variables associated with lattice points and pressure variables associated with the centers of lattice cells [42].
The explicit and implicit MAC solvers are described in the following subsections.

D.1.1

Explicit MAC solverFirst, from equation (6.1), omitting the term describing the volume force exerted by particles on fluids, we discretize
the pressure term implicitly and the remaining terms explicitly:ũn+1 − un1ψnϕn
+ (un · ∇)un + ∇pn+1 − ν∇2 un + ∇µnψ + ∇µnϕ = 0.
∆t
ρ
ρ
ρ

(D.1)

Here quantities with superscript n+1 are unknowns, while quantities with superscript n are known at step n. The tildeon the unknown flow velocity indicates that the velocity field obtained here is a predicted velocity field that must
be corrected at a later stage. Assuming the continuity equation holds for the unknown flow velocity, we obtain thefollowing Poisson equation for the pressure:"#∆t 2 n+1
ψnϕn∇ p = ∇ · un − ∆t∇ · (un · ∇)un − ν∇2 un + ∇µnψ + ∇µnϕ .
(D.2)
ρ
ρ
ρ
The first term on the right-hand side here is expected to vanish under ordinary circumstances, but we retain it for
self-adjustment of cycle errors. Solving this Poisson equation yields the unknown pressure pn+1 , which we insert into(D.1) to obtain the unknown flow velocity ũn+1 . The value thus obtained for ũn+1 fails to distinguish between sites
inside particles and in fluid regions, and for this reason it cannot be used in its present form as the value for the next
timestep. Instead, we modify the velocity field by applying a constraint force in particle regions to ensure that thevelocity field in particle regions properly accounts for the motion of particles:
un+1 = ũn+1 + ϕ f p ∆t.

(D.3)

The constraint force applied here is divergence-free, so the resulting velocity field satisfies the continuity equation.
This yields the velocity field for the next step, un+1 .

D.1.2

Implicit MAC solver

By separating velocities from pressures and applying an implicit solver, we achieve greater computational stability
than what is possible with an explicit solver alone, as we now explain [46].
104APPENDIX D. FLUID CALCULATIONS BY FINITE-DIFFERENCE METHODS (FDM)

105First, we discretize the NS equation in the form
1
1
1
ψ∗
ũn+1 − unη∗  n+ 1
ϕn+1 ∗∇u 2 + ∇uTn+ 2 + ∇µ p si∗ ++ u∗ · ∇un+ 2 + ∇pn+1 − ∇ ·∇µϕ = 0.
∆t
ρ
ρ
ρ
ρ

(D.4)

Here quantities with superscript n+1/2 are unknowns, at timestep n + 1/2, and are determined by a Crank-Nicolson
scheme, while quantities with superscript ∗ are known quantities, extrapolated using an Adams-Bashforth (AB)method. For example, flow velocities are given by
1 n(u + ũn+1 ),21u∗ = (3un − un−1 ).
2
1

un+ 2 =

(D.5)
(D.6)

Also, here the viscosity η is a position-dependent variable. If the continuity equation holds for the unknown flow
velocity un+1 in equation (D.4), we obtain the following Poisson equation for the pressure:"# ψ∗∆t 2 n+1η∗  ∗
ϕn+1 ∗∇u + ∇uT∗ + ∇µ p si∗ +∇ p = ∇ · un − ∆t∇ · u∗ · ∇u∗ − ∇ ·
∇µϕ .
(D.7)
ρ
ρ
ρ
ρ
Note that all unknown terms on the right-hand side here are replaced by known terms obtained via the AB method, so
this Poisson equation may be solved independently. The pressure thus obtained is inserted into equation (D.4) to yielda simultaneous system of linear equations that may be solved for the unknown flow velocity. By default, KAPSEL
solves this system using BiCGSTAB [47] with no preprocessing; however, as discussed in Appendix E, the parallel
iterative solver library Lis may be used as an alternative to allow various preprocessing steps and various iterativesolvers.

D.1.3

Solving the Navier-Stokes equations in the presence of shear flow with Lees-Edwards
boundary conditions

The combination of explicit and implicit MAC solvers may also be used for computations involving shear flow with
Lees-Edwards boundary conditions imposed via coordinate transformations based on tensor analysis [23], as we now
discuss. This method is described in detail in Section 4. Consider a shear flow with shear velocity γ̇(t) imposedin the e x direction, where e x belongs to a set of basis vectors for an orthogonal coordinate system. Then the shear
velocity at coordinate y along the ey axis—the direction of the shear-flow gradient—is given by U = γ̇(t)ye x . Next,
for a shear flow subject to these conditions we consider transforming to an oblique coordinate system (covariant basis)that varies in time as the system evolves. The basis vectors for this oblique coordinate system are described in the
text preceding equation (4.2). For KAPSEL simulations in the presence of shear flows with Lees-Edwards boundary
conditions, working in the oblique coordinate system we subtract the shear-flow contribution from the flow velocity toyield ξ̂ = u − U, then solve an NS equation for the contravariant components ξ̂i :∂t ξ̂i + ξ̂ j ∂ˆ j ξ̂i = −ρ−1Gi j ∂ˆ j p̂ + νG jk ∂ˆ j ∂ˆ k ξ̂i − 2γ̇ξ̂2 δi1 ,∂ˆ j ξ̂ j = 0.

(D.8)
(D.9)

where p̂ is the pressure defined in the oblique coordinate system. In equations (D.8) and (D.9) we have adopted
the Einstein summation convention; also, ∂t denotes differentiation with respect to time, while ∂ˆ i denotes spatial
derivatives in the oblique coordinate system (covariant basis), i.e. ∂ˆ i = ∂/∂ x̂i , and Gi j are the contravariant componentsof the metric tensor. Then for the shear flow we have [23]1 + (γ̇t)2 −γ̇t 0

1 0
Gi j =  −γ̇t(D.10)

0
0 1
KAPSEL solves for equations (D.8) and (D.9) via the same procedure used for computations in orthogonal coordinate
systems, discretizing time and applying explicit and implicit MAC solvers in the oblique coordinate system; we omit
further details.

D.2

Solving the Cahn-Hilliard equation

We next describe KAPSEL’s strategy for solving the Cahn-Hilliard (CH) equation given by equation (6.2); as in thecases discussed above, KAPSEL implements both explicit and quasi-implicit solvers for this purpose.

APPENDIX D. FLUID CALCULATIONS BY FINITE-DIFFERENCE METHODS (FDM)

D.2.1

106

Explicit solver

KAPSEL implements Euler’s method as an explicit solver for the CH equation:hiψn+1 − ψn+ ∇ · (ψn un ) = κ∇2 f ′ (ψn ) − α∇2 ψn + wξ|∇ϕn+1 |2 + 2dϕn+1 (ψn − ψ̄) − 2zξψ ∇ϕn+1 · ∇ψn .
∆t(D.11)

Here we have used the conservative form of the convection term in the CH equation to ensure good conservation
behavior for ψ. Note that, in this equation, the quantity ϕn+1 , denoting the value of ϕ on the next timestep, is a known
quantity because the weak-coupling solver evolves ϕ before evolving ψ.

D.2.2

Implicit solver

Numerical approaches to the CH equation involve as many as 4 layers of partial-derivative operators, and are thus
prone to computational instability in general. KAPSEL addresses this difficulty by implementing the following quasiimplicit solution strategy [49]:h3ψn+1 − 4ψn + ψn−1
+ ∇ · (ψn+1 un+1 ) = κ∇2 2 f ′ (ψn ) − f ′ (ψn−1 ) − α∇2 ψn+1 + wξ|∇ϕn+1 |2 + 2dϕn+1 (ψn+1 − ψ̄)2∆ti−2zξψ ∇ϕn+1 · ∇ψn+1 .
(D.12)
This formulation uses a backward-difference scheme for time-derivative terms, and also performs an implicit solveby using values extrapolated via the AB method to linearize the nonlinear potential term on the right-hand side.
By default, KAPSEL uses BiCGSTAB with no preprocessing to solve the resulting simultaneous system of linear
equations; however, as discussed in Appendix E, the parallel iterative solver library Lis may be used as an alternative
to allow various preprocessing steps and various iterative solvers.

D.2.3

Solving the Cahn-Hilliard equation in the presence of shear flow with Lees-Edwards
boundary conditions

The CH equation in the oblique coordinate system readsno∂t ψ̂ + ∂ˆ i (ûi ψ̂) = κGi j ∂ˆ i ∂ˆ j f ′ (ψ̂) − αGkl ∂ˆ k ∂ˆ l ψ̂ + wξGkl ∂ˆ k ϕ̂∂ˆ l ϕ̂ + 2dϕ̂(ψ̂ − ψ̄) − 2zξψGkl [∂ˆ k ϕ̂∂ˆ l ψ̂ + ϕ̂∂ˆ k ∂ˆ l ψ̂](D.13)

where ϕ̂ and ψ̂ are the SP functions, defined in the oblique coordinate system, that respectively identify particle-fluid
interfaces and fluid-fluid interfaces . As in the cases discussed above, both explicit and quasi-implicit solvers are
implemented for the oblique coordinate system.

Appendix E

Interfacing with the Lis library
In KAPSEL the default method for solving simultaneous systems of linear equations in finite-difference method(FDM) calculations is BiCGSTAB [47]. BiCGSTAB is a stable and high-speed algorithm that is widely used for fluid
calculations and is the best choice for almost all KAPSEL calculations, but the possibility that this method may—inextremely rare cases—fail to achieve successful computations cannot be entirely excluded. To allow BiCGSTAB to
be replaced by more robust methods (such as GMRES [65]) in such situations, KAPSEL implements the ability tointerface with the Lis [48] library for general-purpose iterative solvers. When building the KAPSEL executable file,setting the LIS option in the Makefile to ON will automatically select Lis as the default choice of implicit solver
for the Navier-Stokes (NS) or Cahn-Hilliard (CH) equations; this enables a variety of combinations of preprocessing
techniques and iterative solvers. (Note that KAPSEL must be recompiled to use Lis.)
Preprocessing steps and iterative solvers are configured within the program. The Init lis function in the source
file fdm matrix solver.cxx refers to lis solver set option in two places, corresponding to implicit-solver
settings for the NS and CH equations respectively. Lis settings are briefly discussed in comments at relevant places in
the code; for further details, consult the freely available Lis user manual [66].

107Appendix F

Compiling KAPSELKAPSEL version 5 is distributed in the form of binary executable files, so there is no need for users to compile thecode from source. For the benefit of developers interested in the compilation process, this appendix briefly discusses
the procedure for compiling KAPSEL. The following description assumes that you are working as the root user. When
working with a user account that has administrative rights, use the sudo command as appropriate. The most recent
information may be found at the KAPSEL website: https://kapsel-dns.com.

F.1

The KAPSEL developer’s environment

The procedure for compiling and running KAPSEL depends slightly on whether you are running on a Windows(+Cygwin),
Linux, or Mac system:Linux
If you are using KAPSEL on a Linux system, you may skip to section F.2 below.
Windows
To use KAPSEL on a Windows system, you must first install the Cygwin1 package. In the Select Package section,
set View “Category” and be sure that your Cygwin installation includes the following Packages.
•

all Packages in the Devel category•

all Packages related to fftw3 in the Libs category•

all Packages related to hdf5 in the Libs category•

python3 Package in the Python category

After installing Cygwin, skip to section F.2 below.
Mac

To use KAPSEL on MacOS, you must first install Xcode and the command line tools. By default, the gcc
command invokes the clang compiler; to extract maximal performance from your hardware it is best to install gcc-12
explicitly using Homebrew.

F.2

Installing OCTA

For tasks such as processing input parameters and visualizing output data, KAPSEL relies on an external user-interface
module known as Gourmet. Gourmet is distributed as an internal component of the open-source package OCTA (auniversal simulator for soft materials), which must be installed before using KAPSEL.
To install OCTA, simply visit http://octa.jp/ to download and run an appropriate installer for your system.2
The instructions below assume that OCTA8.# has been installed in /usr/local/OCTA8# (for Linux or Mac) orin C:\OCTA8.# (for Windows).3
1 http://cygwin.com/
2 Answers to questions regarding OCTA/GOURMET are available upon registering as a member of the OCTA-BBS website.
3 Here the pound symbol (“#”) is to be replaced with numerals to indicate the version of your OCTA installation.

108

APPENDIX F. COMPILING KAPSEL

F.3

109

Building libplatform

libplatform is an I/O library used by OCTA to access data stored in UDF files. Depending on your environment,
building libplatform may require installing additional software tools such as OpenJDK4 or Python35 .
Windows
Log on to Windows with administrator privileges, open a Cygwin terminal window, and execute the following
commands.
$ ln -s /cygdrive/c/OCTA8.# /usr/local/OCTA8#$ cd /usr/local/OCTA8#/GOURMET/src$ make distclean$ ./configure --with-python$ make$ make installLinux•

If using gcc:6$ cd /usr/local/OCTA8#/GOURMET/src$ make distclean$ ./configure --with-python$ make$ make install$ mv ../lib/linux64/libplatform.a ../lib/linux64/libplatform_gcc.a•

If using icc:7$ cd /usr/local/OCTA8#/GOURMET/src$ make distclean$ ./configure CC=icc CXX=icpc --with-python$ make$ make install$ mv ../lib/linux64/libplatform.a ../lib/linux64/libplatform_icc.a

Mac

The build process varies depending on the compiler you are using.
•If using clang:8$ cd /usr/local/OCTA8#/GOURMET/src$ make distclean$ ./configure --with-python$ make$ make install$ mv ../lib/macosx/libplatform.a ../lib/macosx/libplatform_clang.a•

If using gcc:9$ cd /usr/local/OCTA8#/GOURMET/src$ make distclean$ ./configure CC=gcc-12 CXX=g++-12 --with-python$ make$ make install$ mv ../lib/macosx/libplatform.a ../lib/macosx/libplatform_gcc.a

4 AdoptOpenJDK: https://adoptopenjdk.net / For arm64 (Apple silicon): https://www.azul.com/downloads/
5 Anaconda: https://www.anaconda.com/
6 gcc may be installed by running brew install gcc-12 or a similar command.
7 The Intel C compiler, available as part of Intel oneAPI Toolkits
8 clang is the default C compiler on macOS
9 gcc may be installed by running brew install gcc-12 or a similar command.

APPENDIX F. COMPILING KAPSEL

110

After building libplatform, there is a possibility that GOURMET will not run correctly for one reason or another.
In this case, execute one of the the following commands:$ cd /usr/local/OCTA8#/GOURMET$ ./Make-All.sh$ cd /usr/local/OCTA8#/GOURMET/src$ ./build-gourmet
For further details regarding libplatform, consult the OCTA manual.

F.4

Installing FFTW

Windows

The source code must be downloaded and installed manually.
$ ./configure --prefix=/opt/fftw/latest.gcc CFLAGS="-O3" FFLAGS="-O3"--enable-openmp --enable-threads --enable-shared --disable-fortran$ make$ make install

,→

Linux

The source code must be downloaded and installed manually.
•

If using gcc:,→•

$ ./configure --prefix=/opt/fftw/latest.gcc CFLAGS="-O3" FFLAGS="-O3"--enable-openmp --enable-threads --enable-shared --disable-fortran$ make$ make install

If using icc:,→

$ ./configure --prefix=/opt/fftw/latest.gcc CC=icc CXX=icpc CFLAGS="-O3"FFLAGS="-O3" --enable-openmp --enable-threads --enable-shared --disable-fortran$ make$ make installMac•

If using clang:brew install fftw•

If using gcc:

,→

F.5

$ ./configure --prefix=/opt/fftw/latest.gcc CC=gcc-12 CXX=g++-12 CFLAGS="-O3"FFLAGS="-O3" --enable-openmp --enable-threads --enable-shared --disable-fortran$ make$ make install

Installing HDF5

Windows

The source code must be downloaded and installed manually.

APPENDIX F. COMPILING KAPSEL$ ./configure --prefix=/opt/hdf5/latest.gcc CFLAGS="-O3" FFLAGS="-O3"--enable-fortran --enable-cxx$ make$ make install

,→

Linux

The source code must be downloaded and installed manually.
•

If using gcc:,→•

$ ./configure --prefix=/opt/hdf5/latest.gcc CFLAGS="-O3" FFLAGS="-O3"--enable-fortran --enable-cxx$ make$ make install

If using icc:,→

$ ./configure --prefix=/opt/hdf5/latest.icc CC=icc CXX=icpc CFLAGS="-O3"FFLAGS="-O3" --enable-fortran --enable-cxx$ make$ make installMac•

If using clang:$ brew install hdf5•

If using gcc:,→

F.6

$ ./configure --prefix=/opt/hdf5/latest.gcc CC=gcc-12 CXX=g++-12--enable-fortran --enable-cxx$ make$ make install

Installing LIS (optional)

Windows

The source code must be downloaded and installed manually.
$ $./configure --prefix=/opt/lis/latest.gcc$ make$ make installLinuxThe source code must be downloaded and installed manually.
•

If using gcc:
$ $./configure --prefix=/opt/lis/latest.gcc$ make$ make install•

If using icc:
$ ./configure --prefix=/opt/lis/latest.icc CC=icc CXX=icpc$ make$ make install

111

APPENDIX F. COMPILING KAPSEL

112

Mac

The source code must be downloaded and installed manually.
•

If using clang :$ ./configure --prefix=/opt/lis/latest.clang CC=clang CXX=clang++$ make$ make install•

If using gcc:
$ ./configure --prefix=/opt/lis/latest.gcc CC=gcc-12 CXX=g++-12$ make$ make install

F.7

Building the KAPSEL executable file

Download the latest version of the KAPSEL source-code archive file kapsel#.#.zip10 and uncompress the archive:$ unzip kapsel#.#.zip$ cd kapsel#.#/src
To ensure that the library file libplatform.a built in Section F.3 is properly linked, modify the GOURMET HOME PATH
variable in your Makefile. Depending on your environment, it may also be necessary to modify the -I (search pathfor include files) and/or -L (search path for library files) flags. Once you have correctly configured all paths, delete
temporary working files.
$ make clean
Then select and execute one of the following make commands to build the KAPSEL binary executable file appropriate
for your computing environment:
Windows
Use Cygwin:$ make ENV=CYGWIN$ make ENV=CYGWIN FFT=FFTW$ make ENV=CYGWIN_OMP FFT=FFTW
Linux•

If using gcc:
$ make ENV=GCC$ make ENV=GCC FFT=FFTW$ make ENV=GCC_OMP FFT=FFTW•

If using icc:
$ make ENV=ICC$ make ENV=ICC FFT=IMKL$ make ENV=ICC_OMP FFT=IMKLMac•

If using clang:

10 Here the pound symbol (“#”) is to be replaced with numerals to indicate the version of your KAPSEL installation.

APPENDIX F. COMPILING KAPSEL$ make ENV=CLANG$ make ENV=CLANG FFT=FFTW$ make ENV=CLANG_OMP FFT=FFTW•

If using gcc:
$ make ENV=GCC_MAC$ make ENV=GCC_MAC FFT=FFTW$ make ENV=GCC_MAC_OMP FFT=FFTW

All platforms

Copy the executable kapsel to bin, and create a symbolic link in the KAPSELinstall directory.
$ cd ..
$ cp ./src/kapsel ./bin/kapsel_self_made$ ln -s ./bin/kapsel_self_made ./kapsel
Set appropriate environment variables, and run KAPSEL.

,→$ export DYLD_LIBRARY_PATH = "/opt/fftw/latest.gcc/lib: /opt/hdf5/latest.gcc/lib:/opt/lis/latest.gcc/lib: DYLD_LIBRARY_PATH"$ ./kapsel
Usage:> ./kapsel -I[input UDF] -O[output UDF] -D[define UDF] -R[restart UDF]

113

Acknowledgements
This work was supported by JST PRESTO, JST-CREST, Grant-in-Aid for Scientific Research, New Energy and Industrial Technology Development Organisation (NEDO), and by the National Institute of Advanced Industrial Scienceand Technology (AIST), Kyoto University and Kyushu University.

114References[1]

R. Yamamoto, J. J. Molina, and Y. Nakayama, “Smoothed profile method for direct numerical simulations ofhydrodynamically interacting particles”, Soft Matter 17, 4226 (2021).

[2]

J. N. Israelachvili, Intermolecular and Surface Forces, 3rd Edition (Academic Press, 2011).

[3]

J. R. Blake, “A spherical envelope approach to ciliary propulsion”, Journal of Fluid Mechanics 46, 199 (1971).

[4]

J. J. Molina, Y. Nakayama, and R. Yamamoto, “Hydrodynamic interactions of self-propelled swimmers”, Soft
Matter 9, 4923 (2013).

[5]

Y. Nakayama and R. Yamamoto, “Simulation method to resolve hydrodynamic interactions in colloidal dispersions”, Physical Review E 71, 036707 (2005).

[6]

Y. Nakayama, K. Kim, and R. Yamamoto, “Simulating (electro)hydrodynamic effects in colloidal dispersions:
Smoothed profile method”, The European Physical Journal E 26, 361 (2008).

[7]

J. J. Molina and R. Yamamoto, “Direct numerical simulations of rigid body dispersions. I. Mobility/friction
tensors of assemblies of spheres”, The Journal of Chemical Physics 139, 234105 (2013).

[8]

J. V. José and E. J. Saletan, Classical Dynamics: A Contemporary Approach (Cambridge University Press, Aug.
1998).

[9]

H. Tanaka and T. Araki, “Simulation Method of Colloidal Suspensions with Hydrodynamic Interactions: Fluid
Particle Dynamics”, Physical Review Letters 85, 1338 (2000).

[10]

R. Yamamoto, “Simulating Particlçe Dispersions in Nematic Liquid-Crystal Solvents”, Physical Review Letters
87, 075502 (2001).

[11]

X. Luo, M. R. Maxey, and G. E. Karniadakis, “Smoothed profile method for particulate flows: Error analysis
and simulations”, Journal of Computational Physics 228, 1750 (2009).

[12]

T. Iwashita, Y. Nakayama, and R. Yamamoto, “A numerical model for brownian particles fluctuating in incompressible fluids”, Journal of the Physical Society of Japan 77, 10.1143/JPSJ.77.074007 (2008).

[13]

T. Iwashita, Y. Nakayama, and R. Yamamoto, “Velocity autocorrelation function of fluctuating particles in incompressible fluids: toward direct numerical simulation of particle dispersions”, Progress of Theoretical Physics
Supplement 178, 10.1143/PTPS.178.86 (2008).

[14]

T. Iwashita and R. Yamamoto, “Short-time motion of brownian particles in a shear flow”, Physical review. E,
Statistical, nonlinear, and soft matter physics 79, 031401 (2009).

[15]

D. A. S. W. B. Russel and W. R. Schowalter, Colloidal dispersions (Cambridge University Press, Cambridge,
England, 1989).

[16]

Y. Otsubo, “The present status and prospect of suspension rheology”, Nihon Reoroji Gakkaishi 31, 15 (2003).

[17]

T. Matsumoto, “Rheology of colloidal disperse systems”, Nihon Reoroji Gakkaishi 32, 3 (2004).

[18]

I. M. Krieger, “Rheology of monodisperse latices”, Advances in Colloid and Interface Science 3, 111 (1972).

[19]

C.-H. Hsueh and P. Becher, “Effective viscosity of suspensions of spheres”, Journal of the American CeramicSociety 88, 1046 (2005).

[20]

A. Onuki, “Phase transitions of fluids in shear flow”, Journal of Physics: Condensed Matter 9, 6119 (1997).

[21]

S. Toh, K. Ohkitani, and M. Yamada, “Enstrophy and momentum fluxes in two-dimensional shear flow turbulence”, Physica D: Nonlinear Phenomena 51, 569 (1991).

[22]

H. Kobayashi and R. Yamamoto, “Implementation of Lees–Edwards periodic boundary conditions for direct
numerical simulations of particle dispersions under shear flow”, The Journal of Chemical Physics 134, 064110(2011).
115

REFERENCES

116

[23]

H. Kobayashi and R. Yamamoto, “Implementation of lees–edwards periodic boundary conditions for direct
numerical simulations of particle dispersions under shear flow”, The Journal of Chemical Physics 134, 064110(2011).

[24]

K. Kim and R. Yamamoto, “Efficient simulations of charged colloidal dispersions: a density functional approach”, Macromolecular Theory and Simulations 14, 278 (2005).

[25]

K. Kim, Y. Nakayama, and R. Yamamoto, “Direct numerical simulations of electrophoresis of charged colloids”,
Phys. Rev. Lett. 96, 208302 (2006).

[26]

H. Tanaka and T. Araki, “Simulation method of colloidal suspensions with hydrodynamic interactions: fluid
particle dynamics”, Phys. Rev. Lett. 85, 1338 (2000).

[27]

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 International Journal Series B Fluids and Thermal Engineering
44, 526 (2001).

[28]

J.-L. Barrat and J.-P. Hansen, Basic concepts for simple and complex liquids (Mar. 2003).

[29] A Book in Japanese (1995).
[30]

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”, Journal of Colloid and Interface Science 90, 17 (1982).

[31]

R. O’Brien and L. White, “Electrophoretic mobility of a spherical colloidal particle”, Journal of The Chemical
Society, Faraday Transactions 2 74, 1607 (1978).

[32]

H. Ohshima, T. W. Healy, and L. R. White, “Approximate analytic expressions for the electrophoretic mobilityof spherical colloidal particles and the conductivity of their dilute suspensions”, J. Chem. Soc., Faraday Trans.
2 79, 1613 (1983).

[33] A Book in Japanese (2004).
[34]

W. Ramsden, “Separation of solids in the surface-layers of solutions and ‘suspensions’(observations on surfacemembranes, bubbles, emulsions, and mechanical coagulation).—preliminary account”, Proceedings of the royalSociety of London 72, 156 (1904).

[35]

S. U. Pickering, “CXCVI.-Emulsions”, Journal of the Chemical Society, Transactions 91, 2001 (1907).

[36]

F. Fenouillot, P. Cassagnau, and J.-C. Majesté, “Uneven distribution of nanoparticles in immiscible fluids: morphology development in polymer blends”, Polymer 50, 1333 (2009).

[37]

Y. Yang, Z. Fang, X. Chen, W. Zhang, Y. Xie, Y. Chen, Z. Liu, and W. Yuan, “An overview of pickering
emulsions: solid-particle materials, classification, morphology, and applications”, Frontiers in pharmacology 8,287 (2017).

[38]

M. E. Cates and P. S. Clegg, “Bijels: a new class of soft materials”, Soft Matter 4, 2132 (2008).

[39]

M. N. Lee and A. Mohraz, “Bicontinuous macroporous materials from bijel templates”, Advanced Materials
22, 4836 (2010).

[40]

P. C. Hohenberg and B. I. Halperin, “Theory of dynamic critical phenomena”, Reviews of Modern Physics 49,
435 (1977).

[41]

M. Doi, Soft Matter Physics (Oxford University Press, 2013).

[42]

A. Arakawa and V. R. Lamb, “Computational design of the basic dynamical processes of the ucla generalcirculation model”, General circulation models of the atmosphere 17, 173 (1977).

[43]

F. H. Harlow and J. E. Welch, “Numerical calculation of time-dependent viscous incompressible flow of fluidwith free surface”, The physics of fluids 8, 2182 (1965).

[44]

J. H. Ferziger and M. Peric, Computational methods for fluid dynamics (Springer Science & Business Media,
2012).

[45]

D. Jacqmin, “Calculation of two-phase navier–stokes flows using phase-field modeling”, Journal of Computational Physics 155, 96 (1999).

[46]

A. Maruoka, J. Matsumoto, and M. Kawahara, “A fractional step finite element method for incompressiblenavier–stokes equations using quadrilateral scaled bubble function”, Journal of Structural Engineering A 44,
383 (1998).

[47]

H. A. Van der Vorst, “Bi-cgstab: a fast and smoothly converging variant of bi-cg for the solution of nonsymmetric linear systems”, SIAM Journal on scientific and Statistical Computing 13, 631 (1992).

REFERENCES

117[48] Lis: Library of Iterative Solvers for Linear Systems, https : / / www . ssisc . org / lis / index . ja . html,
Accessed: 2019-01-28.
[49]

Y. He, Y. Liu, and T. Tang, “On large time-stepping methods for the Cahn–Hilliard equation”, Applied Numerical Mathematics 57, 616 (2007).

[50] ParaView Documentation, https://docs.paraview.org/en/latest/, Accessed: 2021-09-30.
[51]

C. Domb, Phase transitions and critical phenomena (Elsevier, 2000).

[52]

M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, “Hydrodynamics of soft active matter”, Reviews of Modern Physics 85, 1143 (2013).

[53]

M. J. Lighthill, “Hydromechanics of Aquatic Animal Propulsion”, Annual Review of Fluid Mechanics 1, 413(1969).

[54]

J. R. Blake, “A spherical envelope approach to ciliary propulsion”, Journal of Fluid Mechanics 46, 199 (1971).

[55]

O. S. Pak and E. Lauga, “Generalized squirming motion of a sphere”, Journal of Engineering Mathematics 88,
1 (2014).

[56]

N. Oyama, J. J. Molina, and R. Yamamoto, “Simulations of Model Microswimmers with Fully Resolved Hydrodynamics”, Journal of the Physical Society of Japan 86, 101008 (2017).

[57]

A. Bricard, J.-B. Caussin, N. Desreumaux, O. Dauchot, and D. Bartolo, “Emergence of macroscopic directedmotion in populations of motile colloids”, Nature 503, 95 (2013).

[58]

G. Quincke, “Ueber rotationen im constanten electrischen felde”, Annalen der Physik 295, 417 (1896).

[59]

A. Mauleon-Amieva, M. Mosayebi, J. E. Hallett, F. Turci, T. B. Liverpool, J. S. Van Duijneveldt, and C. P.
Royall, “Competing active and passive interactions drive amoebalike crystallites and ordered bands in active
colloids”, Physical Review E 102, 032609 (2020).

[60]

A. J. Goldman, R. G. Cox, and H. Brenner, “Slow viscous motion of a sphere parallel to a plane wall—i motion
through a quiescent fluid”, Chemical engineering science 22, 637 (1967).

[61]

B. Cichocki and R. Jones, “Image representation of a spherical particle near a hard wall”, Physica A: Statistical
Mechanics and its Applications 258, 273 (1998).

[62]

M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Oxford Science Publications, 1987).

[63]

S. M. Cox and P. C. Matthews, “Exponential time differencing for stiff systems”, Journal of Computational
Physics 176, 430 (2002).

[64]

M. Hochbruck and A. Ostermann, “Exponential integrators”, Acta Numerica 19, 209 (2010).

[65]

Y. Saad and M. H. Schultz, “Gmres: a generalized minimal residual algorithm for solving nonsymmetric linear
systems”, SIAM Journal on scientific and statistical computing 7, 856 (1986).

[66] Lis User Guide, https://www.ssisc.org/lis/lis-ug-en.pdf, Accessed: 2019-01-28.