Lattice-Boltzmann Simulations in Complex Geometries

Bachelor Thesis Lattice-Boltzmann Simulations in Complex Geometries Georg Rempfer October 1, 2010 Supervisor: Prof. Dr. Christian Holm Institute for...
Author: Oliver Hodge
0 downloads 1 Views 1MB Size
Bachelor Thesis

Lattice-Boltzmann Simulations in Complex Geometries Georg Rempfer October 1, 2010

Supervisor: Prof. Dr. Christian Holm Institute for Computational Physics , University of Stuttgart

The picture on the title page shows a snapshot from a molecular dynamics simulation of a charged slit pore containing counterions in their equilibrium distribution.

Eidesstattliche Erklärung Hiermit erkläre ich, dass ich die vorliegende Arbeit selbständig und ohne fremde Hilfe bzw. unerlaubte Hilfsmittel angefertigt, andere als die angegebenen Quellen und Hilfsmittel nicht benutzt und die den benutzten Quellen wörtlich oder inhaltlich entnommenen Stellen als solche kenntlich gemacht habe.

Stuttgart, 1. Oktober 2010

Georg Rempfer

Contents 1 Introduction

1

1.1 Objective of this thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.2 Investigated system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

2 Theoretical treatment

4

2.1 The electrokinetic equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.2 Analytical solution of the electrokinetic equations for the slit pore geometry . . . . .

5

3 Lattice-Boltzmann in ESPResSo 3.1 Brief description of the Lattice-Boltzmann-Method . . . . . . . . . . . . . . . . . . .

8 8

3.2 LB-MD coupling scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.3 Bounce-back boundaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.4 Implementation of arbitrary boundaries for the Lattice-Boltzmann fluid in ESPResSo 12 3.5 A test case - Poiseuille flow in tilted channels . . . . . . . . . . . . . . . . . . . . . . 13 4 Simulation using ESPResSo

18

4.1 Choosing a unit system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 4.2 Simulation parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.3 Detailed description of the simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 4.4 Efficiently simulating small diffusion coefficients, exploiting the special case . . . . . 22 4.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 5 Conclusion and outlook

25

5.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 5.2 Other geometries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 5.3 Electrolyte solutions with two species of ions . . . . . . . . . . . . . . . . . . . . . . 26 6 Appendices

27

6.1 Derivation of the Boltzmann distribution from the diffusion equation . . . . . . . . . . 27 6.2 The Navier-Stokes equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 6.3 Stationary state force density acting on the fluid . . . . . . . . . . . . . . . . . . . . . 29 6.4 Boundary conditions for the potential caused by the ion distribution . . . . . . . . . . 29 6.5 Detailed solution of the electrokinetic equations for the slit pore geometry . . . . . . 30 6.6 Used simulation scripts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 7 References

47

8 Acknowledgement

49

1

Introduction

1 Introduction On the microscale, biological systems consist of extensive multiparticle systems with complex interactions. Analytically determining the behaviour of such systems is an impossible task, but with the increasing calculation power of modern computer systems, this field has become one of the main areas for conducting research via computer simulations. There are all-atom simulations that essentially solve the Newtonian equations of motion by modeling the interactions between single atoms with more or less complex potentials. Unfortunately, the necessary calculation time of this approach imposes serious limitations on the possible system size and simulation times, both of which cause problems in biological systems since amino acids in biological tissue form proteins with over 100,000 atoms and the processes in such macro molecules take place on time scales ranging from femtoseconds for the movement of a single atom up to minutes or even hours for the folding of some proteins [22]. However, in most of those cases it is not necessary to account for every single particle. One approach to reduce the necessary computational effort are so-called coarse-grained simulations. In such simulations, several atoms in a macromolecule are combined into a single particle in the simulation, while preserving the original properties of the macromolecule to the extent necessary. But this approach lies beyond the scope of this thesis. In this thesis, we reduce the number of degrees of freedom by replacing the huge number of solvent molecules present in biological systems with a mesoscopic fluid model. There is a wide range of methods to model solvent properties. They range from so-called implicit solvent models [27], where hydrodynamic interactions are completely neglected (e.g. weakening occuring electrical fields to account for the dipolar nature of water) over Brownian Dynamics, where the Newtonian equations of motion are extended with frictional and stochastic terms [4], to methods that reproduce the Navier-Stokes equations on the mesoscopic scale in cases where hydrodynamics have significant influence on the investigated behaviour. There are several different algorithms that reproduce hydrodynamics on the mesoscopic scale with very different methods among which the most widely used ones are DPD, SRT and LBM. Dissipative Particle Dynamics [16] spreads momentum through the particles using additional big particles with a soft interaction potential. Stochastic Rotation Dynamics [25] mimics the effect of collisions between particles by randomly rotating parts of their velocities and the Lattice-Boltzmann-Method [12] discretises space, time and the fluid velocity to allow for an efficient calculation.

-1-

1

Introduction

1.1

Objective of this thesis

1.1 Objective of this thesis The focus of this thesis lies on the simulation of processes in molecular dynamics that are governed by electrostatic and hydrodynamic interactions in volumes with boundaries. The first objective is to change the existing implementation of the Lattice-Boltzmann-Method in ESPResSo (Extensible Simulation Package for Research on Soft matter) so that it can handle systems with arbitrarily complex boundary geometries. The necessary implementations of algorithms to handle the electrostatics already exist. Additionally, a scenario involving the fluid-boundary interaction and the fluid-particle interaction is supposed to be developed. For this scenario, analytical results should be obtained with which the correctness of the implementation can be verified.

1.2 Investigated system The method of choice to treat the scenario analytically are the electrokinetic equations, which are a system of coupled partial differential equations. This is an approximation, replacing the electrostatic potential of charged particles by the potential caused by a continuous charge distribution which resembles the distribution of the particles in a timeaverage (cf. section 2.1 for a more detailed description). For us to be able to solve the electrokinetic equations analytically, the geometry of the system has to fulfill very specific conditions, which are: • Translational symmetry in two spacial coordinates, so that the electrokinetic equations consist of ordinary rather than partial differential equations, • preferrably best described in a cartesian coordinate system, because the differential operators assume the simplest possible shape, • boundaries so that there is a stationary state under the influence of homogeneous force fields, • still enough complexity to allow for nontrivial behaviour. A system which fulfills all those conditions and additionally is simple to simulate with periodic boundary conditions, is the electro-osmotic flow in an infinite slit pore. It consists of two charged, parallel, infinite planes, containing a solvent with counterions, so that the system maintains electrical neutrality, and has a constant electrical field applied along those planes (cf. figure 1). Since the walls and counterions are oppositely charged, the stationary counterion density in the proximity of the wall will be increased. The exact distribution is determined by the temperature, charge, and density of the ions and by the permittivity of the solvent.

-2-

1

Introduction

1.2

_

E

_

_

_

_

_

_ _

Investigated system

+

_ _

E

z

x

y

Figure 1: The investigated system, consisting of an infinite slit pore, bounded by charged walls, containing a fluid with counterions in solution. Additionally shown is the externally applied electrical field and the chosen cartesian coordinate system.

The ions will be pulled along the slit, in turn pulling the solvent along. A stationary state will be reached when the momentum dissipation due to friction between the solvent and the walls equals the momentum transfer to the ions caused by the electrical field. The exact shape of the profile of the fluid velocity in x-direction depends on the ion distribution but is zero at the walls due to the no-slip boundary conditions that we are going to impose and is symmetric to the plane centered in between the walls due to the symmetry of the geometry. However, it will not be a parabola since this is the result for the so-called Poiseuille flow, where there is a homogeneous force density acting on the fluid. In our case the force density is inhomogeneous due to the inhomogeneous distribution of the counterions. However, the slit pore is not the only geometry which fulfills the requirements. Another possibility would be the infinite cylinder. Section 5.1 contains the corresponding equations for that case.

-3-

2

Theoretical treatment

2 Theoretical treatment 2.1 The electrokinetic equations One way to describe the dynamics of the given system is by not looking at the trajectories of individual ions, but rather by representing them with the time average of their spacial distribution ρ. Since we are only interested in the stationary state of the system, we have no problem justifying that we can actually take such an average – we have infinite time available to average over. The equation describing the ion density ρ in this model is the diffusion equation ~j = −D∇ρ − µ∇Ψ ∇ · ~j = 0, with ~j the ion density flow, D the diffusion coefficient, µ the mobility of the ions in the solvent, and −∇Ψ a conservative force density acting on the ion distribution. The second equation is the equation of continuity ∇ · ~j = −∂t ρ in the stationary case. With the Einstein-Stokes relation D µ

= kB · T , one can show for simply connected geometries, that if a stationary state exists, the

ion density ρ resembles a Boltzmann distribution in the potential Ψ (cf. section 6.1 for the detailed derivation)   Ψ ρ = c · exp − . kB T

(1)

Unfortunately, this approach does not take into account the influence of the underlying solvent. To account for that, we have to introduce an additional term ρ~v into the diffusion equation, which yields the convection-diffusion equation ~j = −D∇ρ − µ∇Ψ + ρ~v , with the solvent velocity ~v . This velocity on the other hand, is determined by the Navier-Stokes equations ρfl ∂t~v + ρfl ~v ∆~v = −∇P + η ∆~v + f~fl , which can be drastically simplified in our case (cf. section 6.2 for details), yielding the Stokes equations ∆~v =

1 1 ∇P − f~fl , η η

with η the dynamic viscosity, P the pressure and f~fl the force density acting on the fluid. We also have to consider that in the case of electro-osmosis, the force density −∇Ψ acting on the ion distribution consists of two parts, the electrical field −∇Φ, caused by the ion distribution itself and

-4-

2

Theoretical treatment

2.2

Analytical solution of the electrokinetic equations for the slit pore geometry

~ ext the external field E   1 ~ −∇Ψ = ρq · −∇Φ + Eext , εr where q denotes the charge of an individual ion. The force densities f~fl and −∇Ψ are actually the same (cf. section 6.3 for the proof). Adding the Poisson equation to describe the electrostatic potential Φ caused by the ions gives us the full electrokinetic equations. This is also where the approximation happens. Instead of the charge distribution of single ions, which would essentially be several δ-peaks moving according to the Newtonian equations of motion, one uses the charge distribution coming from the distributions of ions in the time average. This approximation becomes invalid for higher ion densities [11]. For a stationary state with an isotropic, incompressible solvent of permittivity εr , the electrokinetic equations assume the following shape, ε0 being the vacuum permittivity: ~j = −D∇ρ + µρq ·



 1 ~ Eext − ∇Φ + ρ~v , εr

∇ · ~j = 0 , q · ρ, ε 0 εr   ρq 1 ~ 1 · Eext − ∇Φ , ∆~v = ∇P − η η εr

∆Φ = −

∇ · ~v = 0 . Even with the Navier-Stokes equations being linear due to drastic simplifications, this system of partial differential equations can not be solved analytically for most geometries, because of the extensive coupling between the individual equations (cf. colouring of the variables). Since these equations constitute an adequate model for processes in the field of technical biology and chemistry (i.e. substance separation, surface characterisation, or drug release), they are subject to numerical examination [18], [19].

2.2 Analytical solution of the electrokinetic equations for the slit pore geometry Making use of the properties of our slit pore geometry, we can reduce the coupling of the equations in such a way that no component of the diffusion equation is coupled to both, the Poisson equation and the Navier-Stokes equation. Elementary operations involving the no-slip boundary conditions and symmetry arguments yield for the ion density flow jx =

µρq · E + ρvx , εr

jy = −D · ∂y ρ − µρq · ∂y Φ = 0, jz = 0,

-5-

2

Theoretical treatment

2.2

Analytical solution of the electrokinetic equations for the slit pore geometry

~ ext . For the electrostatic potential we get with E the absolute of the external field intensity E ∂y2 Φ = −

q · ρ, ε0 εr

(2)

for the solvent velocity ∂y2 vx = −

qE · ρ, εr η

vy = 0,

vz = 0 ,

(3)

∂z P = 0 .

(4)

and for the pressure ∂x P = 0,

∂y P = −ρq · ∂y Φ,

This set of equations is conceptually easy to handle. One can either argue that theorem (1) applies to equation (2) or, in this special case, solve the equation for the ion density flow in ydirection directly, using separation (cf. section 6.5). The result is that the ion density ρ resembles a Boltzmann distribution in the potential qΦ   qΦ . ρ = c · exp − kB T

(5)

Combining this result with the Poisson equation (2) yields the Poisson-Boltzmann equation ∂y2 Φ = −

  qc qΦ · exp − , ε0 ε r kB T

(6)

an important ansatz, very often used to obtain electrolytic equilibrium distributions [11], [7], [26]. See section 6.1 for a more general derivation. Solving this equation (cf. section 6.5 for the used mathematical methods) and employing the obtained Φ back into ansatz (5) results in the ion distribution ε0 εr C 2 1  , ρ(y) = · 2kB T cos2 qC · y 2kB T

qC π · y 2kB T < 2 .

To fulfill the proper boundary conditions for the electrical field (cf. section 6.4 for details on those), the integration constant C has to fulfill the following transcendental equation, which has to be solved numerically:  C · tan

qd ·C 4kB T

 =−

σ , ε0 ε r

0≤C

Suggest Documents