Confinement effects on the dynamics of a rigid particle in a nanochannel

The transport of nanoparticles in confined geometries plays a crucial role in several technological applications ranging from nanopore sensors to filtration membranes. Here we describe a Brownian approach to simulate the motion of a rigid-body nanoparticle of an arbitrary shape under confinement. A quaternion formulation is used for the nanoparticle orientation, and the corresponding overdamped Langevin equation, completed by the proper fluctuation-dissipation relation, is derived. The hydrodynamic mobility matrix is obtained via dissipative particle dynamics simulation equipped with a new method for enforcing the no-slip boundary condition for curved moving solid-liquid interfaces. As an application, we analyzed the motion of a nanoparticle in a cylindrical channel under the action of external fields. We show that both axial effective diffusion and rotational diffusion decrease with confinement.

Figure 1

Poiseuille flow simulation set-up and wall model. (a) The planar channel is made using wall particles (red spheres) distributed following Eq. (26) with parameters ${\rho }_w^{\left(0\right)}=24, a_{wf}=5$, and $\alpha =0.75$. A force parallel to the $\mathit{z}$ direction is applied to each fluid particle in a slab parallel $Oxy$ plane. Those particles are represented as blue spheres, while the remaining fluid particles are shown as gray transparent spheres. (b) Density profile of the bottom wall particles (red) and of the fluid particles (blue), as a function of the distance from the lower boundary of the simulation box, $\eta$. The horizontal dashed lines represent the quote of the hydrodynamic boundaries, and the vertical dashed line represents the bulk fluid density ${\rho }_0$. In the inset, a zoom of the density profile of the wall particles is reported. Again, the horizontal dashed line is the quote of the hydrodynamic boundary, while the dotted represents the density profile of Eq. (26). (c) Velocity profile ($\overline{g}=0.06, {\gamma }_{wf}=2\times {10}^4$) of the fluid as a function of $y$. The red areas are the regions below and above the hydrodynamic boundaries while the dashed profile is the analytic prediction (Poiseuille flow) for no-slip walls. It is apparent that, with the exception of a thin region close to the boundaries, the analytical and the numerical solutions are in very good agreement. In the inset, the velocity profile near the hydrodynamic boundary is highlighted.

Figure 2

Slip length for planar Poiseuille flow. (a) Slip length $L_s$ as a function of the dissipative force parameter ${\gamma }_{wf}$ between the wall and the fluid particles. $L_s$ is calculated from Poiseuille flow simulations via Eq. (34). (b) Mean velocity $\overline{u}$ for the planar Poiseuille flow as a function of the pressure gradient $\overline{g}$ for ${\gamma }_{wf}=2\times {10}^4$. The dashed line is the analytical expression for no-slip boundaries ($L_s=0$). The values ${\rho }_w^{\left(0\right)}=24, a_{wf}=5$, and $\alpha =0.75$ are used in all the simulations. Error bars for both graphs are of the same size of the points.

Figure 3

Wall model validation for curved surfaces. (a) Hagen-Poiseuille velocity profile for cylindrical channel of radius $R$. The flow is generated as in the planar case, i.e., an axial force $f$ is applied to every particle inside a region of the cylinder of width 2. The cylinder has length $L_c=42$. Different colors correspond to different $R$ while the black dashed line is the analytical expression. The radial position $r$ is rescaled with the radius of the cylinder $R$, and the axial velocity is rescaled with the maximum velocity $u_{\max }=\overline{g}{R}^2/\left(4\mu \right)$ for no-slip Hagen-Poiseuille flow. (b) Tangential flow velocity profile $u_{\theta }$ for the spherical particle at rest in a uniform flow. The velocity is measured in a plane normal to the streaming velocity and passing through the sphere center. The black dashed line is the analytical solution [41] while the computed velocity is in blue circles. The radial coordinate is rescaled with the particle radius $a$, and the velocity is rescaled with the velocity of the fluid far from the sphere, $u_{\infty }$. The data refer to $a=2.13$, i.e., the sphere radius is only 2 times the size of a DPD particle.

Figure 4

Scheme of the simulated system and calculation of the resistance matrix. (a) A spherical particle of radius $a$ is located in a cylindrical channel of radius $R$. The particle reference frame has origin in the particle center and it is constituted by the radial axis $\vec{\rho }$, the transversal axis $\vec{\theta }$, and the axial direction $\vec{z}$. The distance between the particle center and the cylinder center is $r$. (b) Snapshot of the DPD simulation used to calculate the resistance matrix of the particle inside the cylinder. To determine the resistance matrix for a given particle position, six simulations have been carried out. In each simulation one of the six velocity components, three translational in the $\left(\vec{\rho },\vec{\theta },\vec{z}\right)$ directions and three rotational, is imposed while the other five are zero. For each run, the six components of the generalized forces, three forces and the three torques, acting on the particle are measured, leading to the 36 values of the resistance matrix Eq. (37). In panel (b), the procedure to calculate ${\gamma }_{zz}=-F_z/V_z$, with $V_z$ the imposed translational velocity and $F_z$ the average measured force along the $z$ axis, is sketched.

Figure 5

Diagonal components of the resistance matrix, Eq. (37), as a function of the normalized position $\chi =r/(R-a)$ for different confinement ratios $\lambda =R/a$. The resistance values are normalized with respect to the bulk resistance, ${\gamma }_t=6\pi \mu a$ for the translational components and ${\gamma }_r=8\pi \mu a^3$ for the rotational components.

Figure 6

Results of the Langevin simulations. (a) Mean-squared displacement $C_{zz}$, Eq. (39), for different confinement $\lambda$. The axial position is normalized with the radius of the particle $a$, and the reference time is the time necessary for the particle to move in the bulk of a distance equal to its radius ${\tau }_{\mathrm{ref}}=\frac{a^2}{2D_t}$. The dashed line corresponds to the bulk diffusion slope. (b) Effective diffusion coefficient as a function of $\lambda$, as computed from the Langevin simulations (red squares). The dashed line is the analytical result calculated with the centerline approximation [43], and the dotted line represents the curve obtained by Dechadilok and Deen fitting computational results [44]. (c) Axial rotational displacement covariance $C_{\zeta \zeta }$ for $\lambda =3$ (red points) and for $\lambda =10$ (yellow points). The horizontal dotted line represents the value reached for large $\tau$ in the case of free diffusion, Eq. (44) [18]. The inset is a zoom of the curve for low $\tau$, where $C_{\zeta \zeta }$ is almost linear. The rotational diffusion coefficient was computed considering the slope of $C_{\zeta \zeta }$ for small $\tau$. The dashed line in the inset corresponds to the slope in the bulk diffusion case. (d) Axial rotational diffusion coefficient $D_{\zeta \zeta }$ as a function of $\lambda$.

Figure 7

Distribution of particle position when a constant force orthogonal to the cylinder axis is applied. The position is projected into a section of the cylinder to obtain a 2D map. (a) Sketch of the applied force and of the corresponding 2D distribution for $\lambda =10$ and $f=0.125k_{B}T/a$. The probability distribution $P(x,y)$ is normalized with respect to the uniform distribution $P_0=1/\left(\pi R^2\right)$. (b) Particle equilibrium distribution along the $y$ axis. The red and blue histograms correspond to $\lambda =2.87$ and $\lambda =9.87$, respectively. For both cases, $f=0.125k_{B}T/a$. The dashed lines represent the predictions from Eq. (48).

Figure 8

Results of the Langevin simulations with a constant force orthogonal to the cylinder axis applied on the particle. (a) Effective diffusion coefficient $D_{\mathrm{eff}}$ as a function of the external force $f$ for different confinement parameters $\lambda$. For each value of $\lambda$, the effective $D_{\mathrm{eff}}$ is normalized with respect to $D_{\mathrm{eff}}\left(0\right)$, i.e., the effective diffusion in absence of external forces, but in the presence of confinement. (b) Axial rotational diffusion coefficient, Eq. (44), as a function of $f$ for different $\lambda$. Again, the values are normalized with respect to the diffusion in absence of external force.

Figure 9

Results of the Langevin simulations with a constant electric field $\vec{E}=E_0\widehat{z}$ parallel to the cylinder axis applied on a particle with an electric dipole. (a) Sketch of the system. A particle with an electric dipole $\vec{p}$ moves in a cylindrical channel in which a constant electric field $\vec{E}$ is applied. The relevant angle for the potential energy of the particle, i.e., Eq. (49), is the angle $\alpha$ between $\vec{p}$ and $\vec{E}$. (b) Probability distribution of $\alpha$ for three different values of the electric field intensity $E_0$. The numerical results, squares and triangles, are compared with the analytical formula, Eq. (54), dashed lines.

Figure 10

Sketch of the relative orientation between particle dipole $\vec{p}$ and electric field $\vec{E}$. Two angles, $\alpha$ and $\theta$, specify the orientation of the dipole direction $\widehat{p}=\vec{p}/p_0$ while the angle $\psi$ specifies the rotation of the particle around $\widehat{p}$. The electric field is directed along $\widehat{z}$.