Chapter 5 Concluding Remarks
2.5 μm)
=0
εpre for a posteriori error estimator) use the first-order element.
...89 Table 4-2 Evolution of adaptive mesh refinement for computing force of interaction
between two charged spheres (radius=5) within a cylindrical pore (λ=0.416, 0003
.
=0
εpre for a posteriori error estimator) use the second-order
element...90 Table 4-3 Timing (seconds) of PAMR module for different processor numbers. ...91 Table 4-4 Evolution of force interaction between two identical charged spheres near a
like-charged flat plate (h=2.5μm, r=3μm, a=0.325μm and spheresΨ =3.0, 0 flat plateΨ =5.0). ...92 0 Table 4-5 Force of interaction with different refinement level mesh in the
simulation of two identical charged spheres near a charged flat plate
(h=2.5, r=4.5, a=0.325) for the first-order element...93 Table 4-6 Force of interaction with different refinement level mesh in the simulation
of two identical charged spheres near a charged flat plate (h=2.5, r=4.5, a=0.325) for the second-order element. ...94
List of Figures
Fig. 1-1 Electrostatic potential near a positively charged colloidal particle...95
Fig. 1-2 The illustration of colloidal particle ...96
Fig. 1-3 The Electric Double Layer distribution...97
Fig. 1-4 p-refinement ...98
Fig. 1-5 r-refinement...99
Fig. 1-6 h-refinement ...100
Fig. 2-1 Three-dimensional tetrahedral element and face...101
Fig. 2-2 Interpolation functions of ten variable-number-nodes three-dimensional tetrahedral element...102
Fig. 2-3 Traditional shared memory multiprocessor model...103
Fig. 2-4 Distributed memory multiprocessor architecture...103
Fig. 2-5 PC cluster at MuST ...104
Fig. 2-6 The Poisson-Boltzmann equation solver by subdomains...105
Fig. 2-7 The flow chart of parallel Poisson-Boltzmann equation solver ...106
Fig. 3-1 The surface mesh distribution for the potential distribution of the a sphere (176309 elements; 36396 nodes) ...107
Fig. 3-2 Comparison of potential distributions around a charged sphere (a) surface potential Ψ =1.0 (b) surface potential 0 Ψ =5.0 (18,253 nodes; 96,638 0 elements)...109
Fig. 3-3 The potential distribution for the sphere confined in a pore at separation distance r=6...110
Fig. 3-4 Potential along midplane between two spheres (r=6) ... 111
Fig. 3-5 The parallel speedup for the parallelized Poisson-Boltzmann equation solver for simulating the potential distribution around a charged sphere (radius=5) within a cylindrical pore (106,390 nodes; 612,107 elements)...112
Fig. 3-6 The distribution of potential around the charged spheres (radius a=0.325μm) near a like-charged plate (h = 2.5μm) at level-0 for solutions of (a) the first-order element (b) the second-order element...114
Fig. 3-7 The profiles of different mesh levels of the first-order element and the second-order element. the distribution of potential around the charged spheres (radius a=0.325μm) near a like-charged plate (distance h = 2.5μm) ...115 Fig. 3-8 The profiles of different mesh levels of the first-order element and the
second-order element (in local). the distribution of potential around the charged spheres (radius a=0.325μm) near a like-charged plate (distance h =
2.5μm)...116
Fig. 3-9 Potential along midplane between two spheres (r=6), compare the first-order element with the second-order element. (nodes 42,308 and elements 246,155) ...117
Fig. 4-1 Isotropic mesh refinement of tetrahedral mesh (T: Tetrahedron)...118
Fig. 4-2 Mesh refinement rules for two-dimensional triangular cell ...119
Fig. 4-3 Schematic diagram for mesh refinement rules of tetrahedron ...120
Fig. 4-4 Example of the common edge shared by the neighboring cells...121
Fig. 4-5 Schematic diagram of the proposed cell quality control ...122
Fig. 4-6 Schematic diagram of typical cell quality control...123
Fig. 4-7 A case that the proposed cell-quality-control would not affect to it...124
Fig. 4-8 Flow chart of parallel mesh refinement module...125
Fig. 4-9 Coupled PPBS-PAMR method...126
Fig. 4-10 The surface mesh distribution between two identical charged spheres (radius=5) in a cylindrical pore (λ=0.416) at (a) level-0 (initial; 924 nodes; 3,821 elements) (b) level-5 (final; 36,316 nodes; 193,368 elements) mesh refinement ...128
Fig. 4-11 Original (Level-3, 1,070,194 tetrahedral cells and 187,649 nodes) mesh distribution and next (Level-4, 2,701,178 cells and 468,944 nodes) refined mesh ...130
Fig. 4-12 Timing for different modules of PAMR at different processors...131
Fig. 4-13 Sketch of two identical charged spheres near a charged infinite flat plate 132 Fig. 4-14 The surface mesh distribution between two identical charged spheres (radius a=0.325μm) near a like-charged plate (h = 2.5μm) at the (a) level-0 (initial, 16,280 nodes; 79,535 elements) (b) level-5 (finial, 185,697 nodes; 1,014,730 elements) mesh refinement ...134
Fig. 4-15 The distribution of potential around the charged spheres (radius a=0.325μm) near a like-charged plate for (a) h = 2.5μm (b) h=9.5μm ...136
Fig. 4-16 Computed electrostatic force between charged spheres near a like-charged infinite flat plate as a function of the sphere center-to-center distance (h = 2.5μm for the second-order element. ...137
Fig. 5-1 Compare the code of tetrahedral mesh with the code of hybrid mesh of collapsed nodes ...138
Symbols
a The radius of the particle (m) E The electric field vector e The charge of a proton (C)
F The force acting on the spherical particle (Newton) h The distance between plate and particle (m)
N i FEM shape function
n i Concentration of type i ion (m ) −3
0
n i Bulk ionic concentration of type i ion (m ) −3 r The separation distance of particles (m) T Absolute temperature ( K )
z i Valence of type i ion
Greek symbols
ψ Electrical potential (V), and
T k ze
b
ψ = ψ
λ A cylindrical pore with the sphere radius to pore radius ratio
λd Debye length (m)
ρe The net charge density (C/ m3)
ε0 Permittivity of vacuum ( CmV ) κ Debye-Huckel parameter (m ) −1
κb Boltzmann constant (J /K) δ The global absolute error
δm The current mean absolute error
δpre The prescribed global relative error ΔΠ The osmotic pressure
Chapter 1 Introduction
1-1 Motivation
Understanding of the electrostatic distribution in the electrical double layer (EDL) near charged particles or charged solid wall is crucial in understanding several important physical problems, including microscopic colloidal system [Okubo and Aotani,1998; Quddus and Moussa, 2004], electrokinetics in microfluidic system [Quddus and Moussa, 2004; Shao et. al., 2006; Collins and Lee, 2004], to name a few.
Behavior of the EDL in a solution can be approximated by the nonlinear Poisson-Boltzmann equation (PBE). The Poisson-Boltzmann equation has been traditionally used to find the electrostatic potential around a macromolecule. Numerical investigation of models based on the PB equation can provide important information on effective particle interaction in colloidal systems.
The Poisson-Boltzmann equation is a second-order partial differential equation with a hyperbolic sine term that makes it exponentially nonlinear in the electrical double layer around the charged objects. Analytical solution is only possible for some problems with simple geometries and boundary conditions. However, most of the applications are three-dimensional with very complicated geometries, which necessitate the efficient and accurate numerical simulation of the Poisson-Boltzmann
three-dimensional Poisson-Boltzmann equation solver (PPBES) using the finite-element method (FEM), coupling with parallel adaptive mesh refinement (PAMR), which can be used efficiently and accurately for practical 3-D applications in the future.
1-2 Background
In this section, several important concepts about the microscopic colloidal behavior that is relevant to the Poisson-Boltzmann equation are described first. Then adaptive mesh refinement which is important to the accuracy and efficiency in solving the Poisson-Boltzmann equation is described in turn.
1-2-1 Poisson-Boltzmann Equation
1-2-1-1 Nature of Colloidal Solutions
Microscopic particles of one phase dispersed in another are generally called colloidal solution or dispersions. Three fundamental forces operate on fine particles in solution:
1. gravitational force, depending on particles density relative to the solvent;
2. viscous drag force, often described by Stoke’s Law, which depicts the force as a function of fluid viscosity, particle diameter, and particle settling
velocity;
3. the ‘natural’ kinetic energy of particles and molecules.
The ‘natural’ kinetic energy including electrostatic force, which is an important factor of interaction for particles and molecules.
Fig. 1-1 shows the typical electrostatic potential distribution around a positively charged colloidal particle. Negative counter-ions are attracted towards the particle by the electric field generated by the positively charged surface, but they are also subject to thermal motion, which tends to spread then uniformly through the surrounding medium.
1-2-1-2 Zeta Potential
Almost all particulate or macroscopic materials in contact with a liquid acquire an electronic charge on their surfaces. Zeta potential is an important and useful indicator of this charge which can be used to predict and control the stability of colloidal suspensions or emulsions.
Colloidal particles dispersed in a solution are electrically charged due to their ionic characteristics and dipolar attributes. Each particle dispersed in a solution is surrounded by oppositely charged ions called the fixed layer shown as Fig. 1-2.
Outside the fixed layer, there are varying compositions of ions of opposite polarities,
forming a cloud-like area. This area is called the diffuse double layer, and the whole area is electrically neutral in general.
When a voltage is applied to the solution in which particles are dispersed, particles are attracted to the electrode of the opposite polarity, accompanied by the fixed layer and part of the diffuse double layer, or internal side of the "sliding surface".
Zeta potential is considered to be the electric potential of this inner area including this conceptual "sliding surface". As this electric potential approaches zero, particles tend to aggregate.
1-2-1-3 Electric Double Layer
Electrokinetic phenomena are of considerable importance in many fields of science and engineering. In particular, they exert a strong influence on the flow behavior of a fluid, for example, in microchannels and capillaries.
The earliest studies of the electric double layer near a metal surface, by Helmholtz [Hunter, 1989] in the mid-nineteenth century, assumed that the electric charge on the metal surface was balanced by an equal charge which sat in the surrounding solutions a little way from the surface. Around 1910 a better model was proposed by a French scientist name Gouy and a similar treatment was developed
independently a few years later by the British scientist, Chapman. The result is known as the Gouy-Chapman model [Hunter, 1989] of the electric double layer, modified by Helmholtz model. It assumes that most solid surfaces bear electrostatic charges. When a charge surface is in contact with an electrolyte, the electrostatic charges on the solid surface will influence the distribution of nearby ions in the electrolyte solution. Ions of opposite of charges to that of the surface are attracted towards the surface, while ions of like charges are repelled from the surface; thus, an electric field is established. The charges on the solid surface and the balancing charge in the liquid are called the
“electric double layer” (EDL). The EDL potential distribution is show in Fig 1-3. The sign and magnitude of the EDL field depend on the nature of the surface and the liquid.
Immediately next to the charged solid surface, there is a layer of ions that are strongly attracted to the solid surface and are immobile. This layer is called the compact or stern layer, normally about several Angstroms thick. The charge and potential distributions in this layer are mainly determined by the geometrical restriction of ion molecule size and the short-range interactions between ions, the wall and the adjoining dipoles.
From the compact layer to the electrically neutral bulk liquid, the net charge density gradually reduces to zero. Ions in this region are affected less strongly by the electrostatic interaction and are mobile. This region is called the diffuse layer of the
EDL. As will be discussed later, an equation called the Poisson-Boltzmann equation is used to describe approximately the ion and potential distributions in the diffuse layer.
The EDL width, called Debye length and afterward noted, determines the electric interaction range between macromolecules and therefore controls the static phase behavior of these systems.
1-2-1-4 Applications of Poisson-Boltzmann Equation
The Poisson-Boltzmann equation (PBE) provides a "fast" and approximate method for calculating the effects of solvent on electrostatic interactions in the modeled system. There are at least five major applications of the PBE:
1. Calculation of the electrostatic potential at the surface of a biomolecule, which is expected to give information about the concentration of small charged solutes in the neighborhood of the molecule and whose inspection may suggest docking sites for biomolecules [Bowen and Sharif, 1998];
2. Calculation of the electrostatic potential outside the molecule, which is expected to give information on the free energy of interaction of small molecules at different positions in the surrounding of the molecule. [Kozack and Subramaniam, 1993];
3. Calculation of the free energy of a biomolecule or of different states of a
biomolecule which gives information on the stability of a biomolecule or of its different states [Sharp and Honig, 1990];
4. Calculation of the electrostatic field to derive mean forces to be added in standard molecular dynamics calculations [Gilson et. al., 1993].
5. Calculation of the electrostatic force among charged particles/walls, which is important in the study of colloidal system.
In the present thesis, we are mainly concerned about the last application of the Poisson-Boltzmann equation, to which it is not limited. The related derivation of Poisson-Boltzmann Equation will be discussed in chapter 2.
1-2-2 Adaptive Mesh Refinement (AMR)
Adaptive mesh refinement (AMR) is a computational technique for improving the efficiency and accuracy of solution of numerical simulations. In the AMR technique we start with a base coarse grid. As the solution proceeds we identify the regions requiring more resolution by some parameter characterizing the solution. We superimpose finer subgrids only on these regions. Finer and finer subgrids are added recursively until either a given maximum level of refinement is reached or the local truncation error has dropped below the desired level.
In general, mesh refinement can be categorized into three methods: (1) local
polynomial-degree-variation (p-refinement); (2) mesh movement (r-refinement); and (3) mesh enrichment (h-refinement), which will be briefly described respectively in the following.
1-2-2-1 Local Polynomial-Degree-Variation (p-refinement)
As shown in Fig. 1-4, local p-refinement scheme increases resolution of the solution by increasing the order of polynomial of the shape function in the element, which originated from finite element modeling. However, its implementation in three-dimensional finite-element modeling is rather complicated and not as popular as the mesh refinement (h-refinement) which will be introduced later.
1-2-2-2 Mesh Movement (r-refinement)
As shown in Fig. 1-5, the total mesh points remain the same in the computational domain after refinement. It is common to use a spring analogy, in which the nodes of the mesh are connected by springs whose stiffness is proportional to certain measure of solution activity over the spring. The mesh points are moved closer into the region where solution gradients are relatively large. This is often applied to the spatial adaptation of a structured mesh.
1-2-2-3 Mesh Enrichment (h-refinement)
As shown in Fig 1-6, with mesh enrichment, mesh points are added or embedded into the regions where relatively large solution gradients are detected, while the global mesh topology remains intact. It is generally regarded that mesh enrichment method has certain advantages over the first two methods. For example, the use of the
h-refinement scheme generally does not impact the coding structure of the original numerical scheme if hanging nodes are properly treated, while the use of p-refinement does, thus complicating the practical implementation. One of the most important advantages is that the mesh enrichment technique is in general many times faster and robust than the remeshing technique.
1-3 Literature Surveys
1-3-1 Numerical Simulation of Poisson-Boltzmann Equation
The Poisson-Boltzmann equation plays an important role in describing many physical phenomena, including like-charge particle interaction [Dyshlovenko, 2001;
Squires and Brenner, 2000; Tuinier, 2003], particle-membrane interaction [Bowen, 1998; Tuinier, 2003], and electrokinetic flow [Yang et. al., 2001], among others. For example, one of the applications, colloidal dispersion, requires the knowledge of electrostatic potential distribution, which can then be used to calculate other quantities,
such as particle-particle interaction force. Features of particle interaction are of great importance for the stability and properties of colloidal dispersions. Accurate numerical modeling based on the Poisson-Boltzmann equation can provide important information on effective particle interaction in colloidal systems.
Bowen and Sharif [1997] presented a two-dimensional (axisymmetric) finite-element solution combined with adaptive mesh refinement. They used the Galerkin finite-element method with nine-node quadrilateral elements to discretize the Poisson-Boltzmann equation. The Newton-Raphson iterative scheme was used to handle the nonlinear hyperbolic sine term, while the Debye-Huckel solution was used as the initial guess Later, Bowen and Sharif [1998] applied the axisymmetric Poisson-Boltzmann equation solver to simulate the interactive force between two like-charge particles within a cylindrical pore with the same sign of charge on the surface. The results showed that at some inter-particle distance the two particles exhibit attractive forces, which is similar to that observed experimentally in Larsen and Grier [1997], although the magnitude is about one order different. It was then argued that in Ref. [Bowen and Sharif, 1998] the difference might be due to the fact that the wall is planar in the experiments [Larsen and Grier, 1997], whereas in the simulation it is cylindrical. For realistic applications like this, a three-dimensional simulation is required for accurate results. However, not only are 3D simulations extremely
time-consuming in practice, but developing efficient tools for them is also quite difficult.
Dyshlovenko [2001] also presented an axisymmetric Poisson-Boltzmann equation solver using finite-element method combined with adaptive mesh refinement. Instead of using nine-node quadrilateral elements, he used six-node triangular elements with a second-order shape function. Later, Dyshlovenko [2002] simulated the same problem as Bowen and Sharif [1998], but found no like-charge attractive force at any inter-particle distance in his simulation. In addition, Das and Bhattacharjee [2005]
examined the effects of the roughness of the wall on the distribution of the electrostatic force using the two-dimensional finite-element method.
As shown in the above reviews, most existing studies use finite-element method with adaptive mesh refinement to solve the two-dimensional Poisson-Boltzmann equation. Very few studies have focused on developing a parallelized three-dimensional Poisson-Boltzmann equation solver even though the development of one would be valuable in understanding real physical phenomena. Two possible obstacles that have hindered the development of three-dimensional simulation are the time-demanding computation of a meaningful 3D simulation and the very complicated implementation of a parallel adaptive mesh refinement for a 3D unstructured mesh.
Thus, it is important that to develop a parallelized three-dimensional
Poisson-Boltzmann equation solver using the Galerkin finite-element method using the first- and the second-order shape function combined with a parallel adaptive mesh refinement (PAMR) module.
1-3-2 Mesh Refinement
1-3-2-1 Adaptive Mesh Refinement
In previous section, we have briefly described three kinds of mesh refinement scheme: local polynomial-degree-variation (p-refinement), mesh movement (r-refinement) and mesh enrichment (h-refinement). As described earlier in previous section, the first [Zienkiewicz et. al., 1983; Peano, 1976] and second [Brackbill, 1993;
Mackenzie and Robertson, 2002] methods have some disadvantages. Thus, comparatively few studies were conducted along these directions. It is generally regarded that mesh enrichment method (h-refinement) outperforms the first two methods [Rausch et. al., 1991]. One of the most important advantages is that the mesh enrichment technique is in general many times faster and robust than the re-meshing technique [Connell and Holms, 1994]. However, some disadvantages were mentioned in Ref.[Rausch et. al., 1991]. For example, appearance of the hanging nodes induces significant modification to the numerical method which couples the mesh refinement scheme. Fortunately, this can be overcome by some simple methods through the
elimination of hanging nodes, such as that proposed by Kallinderis and Vijayan [1993].
In the current study, we intend to employ similar technique to improve the accuracy of the Poisson-Boltzmann equation solver.
1-3-2-2 Parallel Adaptive Mesh Refinement (PAMR)
For treating large-scale simulation problems, not only we have to parallelize the numerical solver, but also we have to parallelize the mesh refinement scheme. There are several ways to classify PAMR algorithm, and some of them are reviewed in the following.
First, from the viewpoint of practical implementation on computers, PAMR can be generally categorized into two types: vectorized shared memory and distributed memory. In the shared memory type [Waltz, 2004], the hardware for parallel implementation is generally very expensive, and the coding is comparably easy, like using OpenMP [Chandra et. al., 2001], for instance. Memory access is direct and fast, which makes communication cost among processors minimal. Moreover, load balancing among processors is expected to be good, while memory contention may occur which possibly slows down the speed. In contrast, the hardware is more affordable for the distributed memory type than the shared memory type, while the coding is rather involved using MPI [Karypis et. al., 2003], for instance. Global data access is gained through network communication, which often becomes the bottleneck for better parallel performance if a large number of processors is used.
First, from the viewpoint of practical implementation on computers, PAMR can be generally categorized into two types: vectorized shared memory and distributed memory. In the shared memory type [Waltz, 2004], the hardware for parallel implementation is generally very expensive, and the coding is comparably easy, like using OpenMP [Chandra et. al., 2001], for instance. Memory access is direct and fast, which makes communication cost among processors minimal. Moreover, load balancing among processors is expected to be good, while memory contention may occur which possibly slows down the speed. In contrast, the hardware is more affordable for the distributed memory type than the shared memory type, while the coding is rather involved using MPI [Karypis et. al., 2003], for instance. Global data access is gained through network communication, which often becomes the bottleneck for better parallel performance if a large number of processors is used.