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.
However, this becomes increasingly unimportant because of the rapid increase in the
speed of network communication, and the dramatic reduction of costs in the past decade. Deciding on what type of parallel machine to implement greatly impacts the fundamental algorithm design as well as the practical implementation of the PAMR.
Second, from the viewpoint of the adaptive mesh refinement (AMR) algorithm itself, there are generally two kinds of mesh refinement: h-refinement and mesh regeneration. Mesh regeneration AMR [Wang and Harver, 1994] regenerates mesh from a mesh generator based on the distribution of error estimator obtained from the previous coarse mesh. For a large-scale 3-D computation, there exist two problems.
First, the process often takes too much time, and second, the interpolation of previous solutions onto the new mesh is rather inefficient, which may slow down the convergence for the simulation on the new mesh. However, the data structure is relatively simple as compared to the h-refinement AMR. The h-refinement AMR adds grid points onto cell edges based on the error estimator obtained from the previous coarse mesh. It is generally much faster than the mesh regeneration type, and the interpolation onto the new mesh from the previous solution is straightforward, which may speed up the convergence for the computation on the new mesh. Often, the resulting “edge-based” data structure is very complicated and memory demanding, especially for the 3-D PAMR on the memory-distributed type of machine.
Nevertheless, this type of AMR also has another advantage in that its data structure may be readily modified for mesh de-refinement, which may become important for a time-dependent simulation. Additionally, the h-refinement AMR possesses high spatial locality, which makes possible the parallel implementation on a memory-distributed machine using domain decomposition.
Third, from the viewpoint of the unstructured mesh solver itself, two types of
mesh data are required: the finite element type of connectivity (or connectivity) and cell-based connectivity (or cell neighbor-identifying information). The equation solver often needs the former, while the latter is required by some particle method for particle tracking, such as DSMC and PIC for solving the Boltzmann equation with neutral and charged species, respectively. However, the past development of the AMR scheme only produced the finite element type of connectivity. Related cell neighbor-identifying information can of course be obtained by post-processing serially (not in parallel) the node-based connectivity, while it may become time-consuming if the resulting mesh size is large. Thus, an ideal PAMR scheme may be expected to directly generate both these two sets of data.
For adaptive mesh refinement of three-dimensional tetartohedral mesh for the Poisson-Boltzmann equation, Cortis and Friesner [1997a, 1997b] proposed an automatic Finite-Element mesh generation system. Authors presented the algorithms of the construction of a tetrahedral mesh using a predetermined spatial distribution of vertices adapted to the geometry of the dielectric continuum solvent model. Besides, Baker[2000, 2001b] and his co-worker [Holst, 2000] publish the adaptive multilevel finite element treatment of the nonlinear PB equation using highly scalable parallel adaptive meshing methods which all of the numerical procedures are implemented in MANIFOLD CODE (MC) [web site, http://www.scicomp.ucsd.edu/ ~mholst/ codes/
mc/], a small self-contained parallel adaptive multilevel finite element software package.
1-4 Objectives of the Thesis
Based on previous reviews, the specific objectives of the thesis are summarized as follows:
1. To develop a parallelized 3-D Poisson-Boltzmann equation solver (PPBES) using finite element method with first- and second-order shape function on an unstructured tetrahedral mesh;
2. To develop a parallel adaptive mesh refinement (PAMR) for unstructured tetrahedral mesh;
3. To couple the PPBES and PAMR and apply it to compute a realistic 3-D problem.
1-5 Organization of the Thesis
The organization of this thesis is described as follows. Chapter 2 describes the numerical method for the parallelized Poisson-Boltzmann equation solver. Chapter 3 describes the validation and parallel performance study of the PPBES. Chapter 4 describes the proposed parallel adaptive mesh refinement, coupling scheme with
PPBES and application to a realistic 3-D problem. Finally, Chapter 5 describes the concluding remarks and recommendation of the future work. In addition, Appendix A describes the derivation of the nodal quadrature used in the present finite element
PPBES and application to a realistic 3-D problem. Finally, Chapter 5 describes the concluding remarks and recommendation of the future work. In addition, Appendix A describes the derivation of the nodal quadrature used in the present finite element