Molecular dynamics simulation of electro-osmotic flows in rough wall nanochannels

We performed equilibrium and nonequilibrium molecular dynamics simulation to study electro-osmotic flows inside charged nanochannels with different types of surface roughness. We modeled surface roughness as a sequence of two-dimensional subnanoscale grooves and ridges (step function-type roughness) along the flow direction. The amplitude, spatial period, and symmetry of surface roughness were varied. The amplitude of surface roughness was on the order of the Debye length. The walls have uniform negative charges at the interface with fluids. We included only positive ions (counterions) for simplicity of computation. For the smooth wall, we compared our molecular dynamics simulation results to the well-known Poisson-Boltzmann theory. The density profiles of water molecules showed “layering” near the wall. For the rough walls, the density profiles measured from the wall are similar to those for the smooth wall except near where the steps are located. Because of the layering of water molecules and the finite size effect of ions and the walls, the ionic distribution departs from the Boltzmann distribution. To further understand the structure of water molecules and ions, we computed the polarization density. Near the wall, its z component dominates the other components, indicating the preferred orientation (“ordering”) of water molecules. Especially, inside the groove for the rough walls, its maximum is 10% higher (stronger ordering) than for the smooth wall. The dielectric constant, computed with a Clausius-Mosotti-type equation, confirmed the ordering near the wall and the enhanced ordering inside the groove. The residence time and the diffusion coefficient, computed using the velocity autocorrelation function, showed that the diffusion of water and ions along the direction normal to the wall is significantly reduced near the wall and further decreases inside the groove. Along the flow direction, the diffusion of water and ions inside the groove is significantly lowered while it is similar to the bulk value elsewhere. We performed nonequilibrium molecular dynamics simulation to compute electro-osmotic velocities and flow rates. The velocity profiles correspond to those for overlapped electric double layers. For the rough walls, velocity inside the groove is close to zero, meaning that the channel height is effectively reduced. The flow rate was found to decrease as the period of surface roughness decreases or the amplitude of surface roughness increases. We defined the ζ potential as the electrostatic potential at the location of a slip plane. We computed the electrostatic potential with the ionic distribution and the dielectric constant both from our molecular dynamics simulation. We estimated the slip plane from the velocity profile. The ζ potential showed the same trend as the flow rate: it decreases with an increasing amplitude and a decreasing period of surface roughness.

Figure 1

(a) Model system for a smooth wall ($●$ innermost layers). It includes charged and uncharged Lennard-Jones wall atoms, sodium ions, and water molecules. The coordinate system was defined as shown. (b) Model system for a rough wall. It has two-dimensional surface roughness with a spatial period $p$ and an amplitude $A$. The flow region can be divided into the regions over the grooves and over the ridges.

Figure 2

One-dimensional profiles for a smooth wall. (a) Oxygen atom and hydrogen atom densities were normalized by their respective bulk values. (b) Water density profile. (c) Sodium ion density profiles obtained from MD simulation and PB solutions with two different boundary conditions (see text). (d) Profiles of electrostatic potential obtained from MD simulation and PB solution.

Figure 3

Density profiles for the rough wall channel with symmetric surface roughness of a period of $36.8\mathrm{\AA }$ and an amplitude of $1.8\mathrm{\AA }$. (a) Water density and (b) sodium ion density.

Figure 4

Electrostatic potential profiles in different model systems. Note that data points are more frequent than symbols. Symbols: $◯$ System 1, $*$ System 2, $∎$ System 3, $▷$ System 4, $◻$ System 5, $●$ System 6, $\times$ System 7, $▿$ System 8, $◇$ System 9, $◁$ System 10, $+$ System 11.

Figure 5

(a) Polarization density profile of water molecules in three directions and (b) dielectric constant profile, both for the smooth wall channel. We calculated the total electric dipole moment of water molecules in each bin and normalized it by the bin volume to compute the polarization density shown. We calculated the dielectric constant using Eq. (15).

Figure 6

Two-dimensional profiles for (a) the polarization density and (b) the dielectric constant of water molecules. The rough wall has symmetric surface roughness of a period of $36.8\mathrm{\AA }$ and an amplitude of $4.5\mathrm{\AA }$.

Figure 7

Velocity autocorrelation function of water molecules in the first layer near the wall for the smooth wall. We normalized it by the value at $t=0$.

Figure 8

Diffusion coefficient profiles of water molecules (a) in the $xy$ plane and (b) along the $z$ direction. The rough wall has symmetric surface roughness of a period of $36.8\mathrm{\AA }$ and an amplitude of $4.5\mathrm{\AA }$. Note that the $x$ axis in both the plots is the $z$-directional distance from the wall. The data points from the right correspond to the first, second, and third layers of water molecules and the bulk. The errors were estimated to be around 5%.

Figure 9

Velocity profiles for (a) the smooth wall and (b) the rough walls. (a) The dots are actual data points from MD and the solid line is a line obtained by smoothing these data points. The reason for this smoothing is to reduce noises of our data points (thermal noises). The dashed line was obtained by solving the Stokes' equation with the ionic distribution, denoted by “Corrected PB” in Fig. 2c. (b) All the solid lines are smoothed lines from the MD data for different types of surface roughness. Symbols are the same as in Fig. 4.

Figure 10

Velocity vector field for the rough wall with symmetric surface roughness of a period of $36.8\mathrm{\AA }$ and an amplitude of $2.7\mathrm{\AA }$.

Figure 11

Flow rate per unit width of the channel, $Q^{\prime }$, for (a) different periods of surface roughness and a fixed amplitude $\left(0.9\mathrm{\AA }\right)$ and (b) different amplitudes with a fixed period $\left(36.8\mathrm{\AA }\right)$. In (a), we simulated both the symmetric and asymmetric cases (as in Table ). The symmetry of surface roughness was defined with respect to the centerline of the channel.

Figure 12

$\zeta$ potential for (a) different periods of surface roughness and a fixed amplitude $\left(0.9\mathrm{\AA }\right)$ and (b) different amplitudes with a fixed period $\left(36.8\mathrm{\AA }\right)$. We computed the $\zeta$ potential in two ways, denoted by MD (molecular dynamics simulation) and HS (Helmholtz-Smoluchowski), which were described in more detail in the text. In (a), we simulated both the symmetric and asymmetric cases (as in Table ). The symmetry of surface roughness was defined with respect to the centerline of the channel.