Flow of colloidal solids and fluids through constrictions: dynamical density functional theory versus simulation

We consider point-like Brownian particles in two spatial dimensions which interact via a pairwise potential

$$\begin{eqnarray}&&u(r)=\frac{{{u}_{0}}}{{{r}^{3}}},\end{eqnarray} \tag{ 1 }$$

where r is the distance between two particles and the amplitude u₀  >  0 sets the interaction strength. A real-world analogue of this system is given by superparamagnetic particles that are confined in a 2d plane with an uniform external magnetic field ${{\mathbf{B}}_{\text{ext}}}$ applied perpendicular to the plane. The external magnetic field ${{\mathbf{B}}_{\text{ext}}}$ induces a dipole–dipole interaction between the colloidal particles, which can be tuned by changing its strength. In bulk, the only relevant length scale present in this system is the typical interparticle distance, which is given by

$$\begin{eqnarray}&&l=\rho _{0}^{-1/2},\end{eqnarray} \tag{ 2 }$$

with

$$\begin{eqnarray}&&{{\rho}_{0}}=N/{{A}_{0}}\end{eqnarray} \tag{ 3 }$$

the number density of the system, N the number of particles, and A₀ the accessible area, which will be defined later. Due to the inverse power law scaling of equation (1), a change in density of the system is equivalent to a change in the interaction strength u₀. It is therefore convenient to rewrite equation (1) as

$$\begin{eqnarray}&&\frac{u(r)}{{{k}_{\text{B}}}T}=\frac{\Gamma}{{{(r/l)}^{3}}}\end{eqnarray} \tag{ 4 }$$

where $\Gamma={{u}_{0}}\rho _{0}^{3/2}/\left({{k}_{\text{B}}}T\right)$ is a dimensionless coupling parameter. The bulk phase behaviour of these particles is characterized by a fluid at low $\Gamma\lesssim 11$, and a hexagonally ordered solid phase at high $\Gamma\gtrsim 12$ [73, 74].

Naturally, this phase diagram is expected to change significantly in the confinement of a channel, as considered here. In particular, as the system is effectively one-dimensional, we expect only short-range ordering in the channel, and no true fluid to crystal transition. Nonetheless, at high $\Gamma$ we do expect local ordering into a hexagonal lattice, aligned with the boundaries of the channel [74].

Inside the 2d plane the particles are additionally confined in a channel geometry along the x-axis, represented by an external potential ${{V}_{\text{ext}}}(x,y)$. The lateral profile of the channel is modelled as error-function steps at the walls of the channel so the external potential is given by

$$\begin{eqnarray}&&{{V}_{\text{ext}}}(x,y)={{V}_{0}}\left[1-\frac{1}{2}\text{erf}\,\left(\frac{y+g(x)}{\sqrt{2}w}\right)+\frac{1}{2}\text{erf}\,\left(\frac{y-g(x)}{\sqrt{2}w}\right)\right],\end{eqnarray} \tag{ 5 }$$

with V₀ being the maximum potential height, $\pm g(x)$ describing the contour lines of the channel walls and w characterizing the softness of the walls. For a straight channel with width L_y and without constriction the contour functions are simply $g(x)\equiv \frac{{{L}_{y}}}{2}$. The constriction is modelled as a single cosine wave of length L_c at x₀ that is added smoothly to the channel contour. Therefore, g(x) is given by

$$\begin{eqnarray}&&g(x)=\left\{\begin{array}{*{35}{l}} \frac{{{L}_{y}}}{2}-\alpha \left[1+\cos \left(2\pi \frac{x-{{x}_{0}}}{{{L}_{c}}}\right)\right],\,\text{if}\,|x-{{x}_{0}}|<\frac{{{L}_{c}}}{2} \\ \frac{{{L}_{y}}}{2},\text{otherwise} \end{array}\right. \end{eqnarray} \tag{ 6 }$$

with amplitude $\alpha =\frac{{{L}_{y}}}{4}(1-b)$. Here, we introduced the parameter b as the ratio of constriction width over the total channel width. Consequently, $0\leqslant b\leqslant 1$, where b  =  0 refers to a completely blocked channel and b  =  1 is a channel without constriction. See figure 1(a) for an illustrative sketch of ${{V}_{\text{ext}}}(x,y)$. Note that the form of the potential in equation (5) was chosen simply to provide a steeply repulsive wall-particle interaction. Experimentally, this can be realized via e.g. optical forces or barriers which are permeable to the solvent.

Zoom In Zoom Out Reset image size

We define the accessible area as the region between the midlines of the two walls, i.e.

$$\begin{eqnarray}&&{{A}_{0}}=2{\int}_{-{{L}_{x}}/2}^{{{L}_{x}}/2}g(x)\text{d}x.\end{eqnarray} \tag{ 7 }$$

By definition, the number density in the system is given by ${{\rho}_{0}}=N/{{A}_{0}}=1/{{l}^{2}}$, with l our unit of length. In this work, we focus on channels with a width chosen such that either five or six crystalline layers reliably form within the channel, oriented such that lines of nearest-neighbours are aligned with the channel walls (see figure 1(b)). However, the number of defects in this crystal strongly depends on the commensurability between the channel width and the lattice spacing of the crystal [74]. In a perfect hexagonal lattice at density ${{\rho}_{0}}=1/{{l}^{2}}$, the distance between two crystal layers is

$$\begin{eqnarray}&&d=\sqrt{\frac{\sqrt{3}}{2}}l,\end{eqnarray} \tag{ 8 }$$

and we will adopt this definition of d for our confined system as well. In order to accommodate a crystal with a low number of defects, we therefore choose the channel width to be L_y  =  nd, with n  =  5 or 6. Both DDFT and simulations show that this indeed leads to crystals with the desired number of layers.

In order to further reduce parameter space, we fix the constriction length L_c  =  2.686l, wall softness w  =  0.25l and ${{V}_{0}}=1000{{k}_{\text{B}}}T$.

We model the dynamics of the particles in the channel via simple, overdamped Brownian dynamics, where we assume the solvent to be at rest. The equations of motions are given by:

$$\begin{eqnarray}&&{{\dot{{\mathbf{r}}}}_{i}}=\frac{D}{{{k}_{\text{B}}}T}{{\mathbf{F}}_{i}}\left({{\mathbf{r}}^{N}}\right)+\sqrt{2D}\boldsymbol{\xi}_{{i}}(t),\end{eqnarray} \tag{ 9 }$$

where ${{\mathbf{r}}_{i}}$ are the coordinates of the ith particle and ${{\mathbf{r}}^{N}}\equiv \left({{\mathbf{r}}_{1}},\ldots,{{\mathbf{r}}_{N}}\right)$ is a short-hand notation for the coordinates of all particles, D is the diffusion constant of a single particle without external forces, ${{\mathbf{F}}_{i}}\left({{\mathbf{r}}^{N}}\right)$ is the total force acting on the ith particle composed of pair interactions, external potential, and driving force:

$$\begin{eqnarray}&&{{\mathbf{F}}_{i}}=-\underset{j\ne i}{\sum}{{\nabla}_{i}}u\left(\left|{{\mathbf{r}}_{i}}-{{\mathbf{r}}_{j}}\right|\right)-{{\nabla}_{i}}{{V}_{\text{ext}}}\left({{\mathbf{r}}_{i}}\right)+f\mathbf{\hat{x}},\end{eqnarray} \tag{ 10 }$$

with u(r) given by equation (4) and ${{\nabla}_{i}}$ being the gradient operator with respect to particle coordinates ${{\mathbf{r}}_{i}}$ and the unit vector in x-direction $\mathbf{\hat{x}}$. The external force responsible for the flow of particles through the channel is modelled via a constant force f along the x-axis. Finally, $\boldsymbol{\xi}_{{i}}(t)$ is a delta-correlated Gaussian noise process modelling the thermal fluctuations. In the remainder of this work, we will fix the driving force $f=1{{k}_{\text{B}}}T/l$. As a unit of time, we will use the time it takes a particle to diffuse by a typical distance of l, i.e.

$$\begin{eqnarray}&&\tau =\frac{{{l}^{2}}}{D}.\end{eqnarray} \tag{ 11 }$$

A stochastically equivalent description of equation (9) is given by the Smoluchowski picture in which the time-dependent N-particle probability distribution $p\left({{\mathbf{r}}^{N}},t\right)$ is considered. The Smoluchowski equation is given by

$$\begin{eqnarray}&&\frac{\partial p\left({{\mathbf{r}}^{N}},t\right)}{\partial t}=D\underset{i=1}{\overset{N}{\sum}}{{\nabla}_{i}}\left[{{k}_{\text{B}}}T{{\nabla}_{i}}+{{\mathbf{F}}_{i}}\right]p\left({{\mathbf{r}}^{N}},t\right).\end{eqnarray} \tag{ 12 }$$

An integration over the probability distribution $p\left(\mathbf{r},t\right)$ with respect to all but one coordinate gives the one-particle density

$$\begin{eqnarray}&&\rho \left({{\mathbf{r}}_{1}},t\right)=N{\int}^{}\text{d}{{\mathbf{r}}_{2}}\ldots {\int}^{}\text{d}{{\mathbf{r}}_{N}}p\left({{\mathbf{r}}^{N}},t\right),\end{eqnarray} \tag{ 13 }$$

which describes the ensemble averaged particle density at time t and is the basic quantity in the DDFT.

Dynamical density functional theory (DDFT) is conveniently derived from the Smoluchowski equation (12) by projecting onto the one-particle density and invoking the additional adiabatic approximation [32]. As a result, DDFT is an approximative theory. It can be written as a continuity equation

$$\begin{eqnarray}\begin{array}{*{35}{l}} \frac{\partial \rho \left(\mathbf{r},t\right)}{\partial t} & =D\nabla \left(\rho \left(\mathbf{r},t\right)\nabla \frac{\delta \mathcal{F}\,\left[\rho \right]}{\delta \rho \left(\mathbf{r},t\right)}\right), \end{array}\end{eqnarray} \tag{ 14 }$$

which expresses the particle number concentration $\rho \left(\mathbf{r},t\right)$. The current density $\mathbf{j}\left(\mathbf{r},t\right)$ is explicitly given by a generalized Fick's law:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \mathbf{j}\left(\mathbf{r},t\right)=-D\rho \left(\mathbf{r},t\right)\nabla \frac{\delta \mathcal{F}\,\left[\rho \right]}{\delta \rho \left(\mathbf{r},t\right)}, \end{array}\end{eqnarray} \tag{ 15 }$$

with the Helmholtz free energy functional

$$\begin{eqnarray}&&\mathcal{F}\left[\rho \right]={{\mathcal{F}}_{\text{id}}}\left[\rho \right]+{{\mathcal{F}}_{\text{ext}}}\left[\rho \right]+{{\mathcal{F}}_{\text{exc}}}\left[\rho \right]\end{eqnarray} \tag{ 16 }$$

which can be split in three principal contributions. The ideal gas term

$$\begin{eqnarray}&&{{\mathcal{F}}_{\text{id}}}\left[\rho \right]={{k}_{\text{B}}}T{\int}^{}\text{d}\mathbf{r}\,\rho \left(\mathbf{r},t\right)\left(\log \left({{\Lambda}^{2}}\rho \left(\mathbf{r},t\right)\right)-1\right)\end{eqnarray} \tag{ 17 }$$

and the external potential contribution

$$\begin{eqnarray}&&{{\mathcal{F}}_{\text{ext}}}\left[\rho \right]={\int}^{}\text{d}\mathbf{r}\,\rho \left(\mathbf{r}\right)\left({{V}_{\text{ext}}}\left(\mathbf{r}\right)-fx\right)\end{eqnarray} \tag{ 18 }$$

with thermal de Broglie wavelength $\Lambda$ are known expressions. In contrast, the excess free energy functional ${{\mathcal{F}}_{\text{exc}}}\left[\rho \right]$, which describes the particle interactions, is unknown and has to be approximated. Here, we use the Ramakrishnan–Yussouff functional described in the next section. Substituting the first two terms in equation (15), the current is thus given explicitly by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \mathbf{j}\left(\mathbf{r},t\right)=-D\nabla \rho \left(\mathbf{r},t\right)+\rho \left(\mathbf{r},t\right)\nabla \left({{V}_{\text{ext}}}\left(\mathbf{r}\right)-fx+\frac{\delta {{\mathcal{F}}_{\text{exc}}}\left[\rho \right]}{\delta \rho \left(\mathbf{r},t\right)}\right). \end{array}\end{eqnarray} \tag{ 19 }$$

Since we are only interested in the flux along the channel, we define the particle flow in the x-direction, i.e.

$$\begin{eqnarray}&&{{j}_{x}}(x,t)={\int}_{-\infty}^{\infty}\text{d}y\;\mathbf{j}\left(\mathbf{r},t\right)\cdot \mathbf{\hat{x}}.\end{eqnarray} \tag{ 20 }$$

The average flow through the channel ${{\bar{j}}_{x}}$ can then simply be defined as the long-time average value of j_x(x, t):

$$\begin{eqnarray}&&{{\bar{j}}_{x}}=\underset{T\to \infty}{{\lim}}\frac{1}{T}{\int}_{0}^{T}\text{d}t\;{{j}_{x}}(x,t).\end{eqnarray} \tag{ 21 }$$

Note that as the particle density is a conserved quantity, ${{\bar{j}}_{x}}$ is independent of the position x.

We chose the Ramakrishnan–Yussouff expression [75] as an approximate excess free energy functional, which is a convenient way to model soft and long-ranged particle interactions. The functional derivative of the Ramakrishnan–Yussouff functional is given as a convolution of $\rho \left(\mathbf{r},t\right)$ and the pair (two-point) direct correlation function $c_{0}^{(2)}\left(r;{{\rho}_{0}},\Gamma\right)$ of an isotropic and homogeneous reference fluid with the prescribed density ${{\rho}_{0}}=1/{{l}^{2}}$, at interaction strength $\Gamma$:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{\delta {{\mathcal{F}}_{\text{exc}}}\left[\rho \right]}{\delta \rho \left(\mathbf{r},t\right)}=-{{k}_{\text{B}}}T{\int}^{}\text{d}{{\mathbf{r}}^{\prime}}\;\rho \left({{\mathbf{r}}^{\prime}},t\right)c_{0}^{(2)}\left(|\mathbf{r}-{{\mathbf{r}}^{\prime}}|;{{\rho}_{0}},\Gamma\right). \end{array}\end{eqnarray} \tag{ 22 }$$

We use the direct correlation functions obtained by liquid integral theory with the Rogers–Young closure which were calculated in [76], where it was shown that despite its simplicity the Ramakrishnan–Yussouff functional accounts for the freezing transition in two dimensions at $\Gamma\gtrsim 36.2$.

Since the functional derivative of the excess functional equation (22) is a convolution of $\rho \left(\mathbf{r},t\right)$ and $c_{0}^{(2)}\left(r;{{\rho}_{0}},\Gamma\right)$ we can efficiently compute its value using fast Fourier transform.

The overall length of the system is chosen as L_x  =  21.5l with periodic boundary conditions along the x-direction. As a discretisation we used ${{N}_{x}}\times {{N}_{y}}=256\times 64$ gridpoints. With prescribed density ${{\rho}_{0}}$ we have about N  =  113–120 particles in our system, depending on the constriction width b.

Starting from several initial density profiles, we solve the DDFT equation without any driving force to obtain an equilibrium density profile ${{\rho}^{(0)}}\left(\mathbf{r}\right)$. We confirmed that the equilibrium profile does not depend on the initial profile. Depending on the coupling parameter $\Gamma$ we either obtain an inhomogeneous fluid (figure 2(a)) or a crystalline profile of hexagonal order (figure 2(b)). The one-dimensional crystal in channel confinement can be observed for $\Gamma\gtrsim 30$, for both investigated channel widths L_y  =  5d and L_y  =  6d.

Zoom In Zoom Out Reset image size

For t  >  0 we switch on the driving force f, initiating the flow through the constriction. We solve equation (14) numerically using a finite volume partial differential equation solver [77].

The average flux in the system ${{\bar{j}}_{x}}$ for a range of coupling parameters $\Gamma$ and constriction widths b is shown in figures 4 and 5. In general, we observe that for stronger particle interactions the average flux is smaller, as the particles more effectively block each other from passing through the constriction. As expected, we also observe a decrease in average flux with decreasing constriction width b. We note, however, that in the simulations this trend is not always monotonic: there are regions where ${{\bar{j}}_{x}}$ decreases with increasing b.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We observe qualitative agreement between the DDFT and simulation results. The main difference occurs at high $\Gamma$, where the simulations observe complete blocking (${{\bar{j}}_{x}}=0$), while the DDFT calculations predict a finite flux. Additionally, the DDFT calculations predict only a monotonous increase in ${{\bar{j}}_{x}}$ with b for the investigated parameter range.

It should be noted here that the observed results are expected to be influenced strongly by the length of the channel: at constant number density, a longer channel implies that the driving force f is applied to a larger number of particles in front of the constriction. This results in a proportional increase in the pressure near the constriction, which is expected to enhance the flow of particles. Indeed, simulations on larger systems (N  =  400) and on systems with larger external forces confirm that doubling the channel length is approximately equivalent to doubling the external driving force on the particles.

After starting the flow, we observe four qualitatively different types of flow behavior in both our DDFT results and our simulations. First, we distinguish between systems that show a complete blockade (i.e. zero particle flow ${{\bar{j}}_{x}}$), and systems that show a finite flow of particles. In the case of a finite flow, the average flux through the constriction eventually reaches a constant value in the simulations. However, shortly after starting the flow, we often observe oscillations in the flux that decay over time. In this regime, we observe three types of decay: an almost immediate decay to a smooth flow, a brief period of transient oscillations without a clearly defined periodicity, and a long-time oscillation with a period which is independent of b and $\Gamma$. In the DDFT calculations we observe the same regimes. However, due to the lack of stochastic noise, in the long-time oscillation regime, the DDFT calculations predict periodic (i.e. non-decaying) oscillations. Below, we discuss each type of flow in detail. In figures 6 and 7, we plot the average flux through the constriction as a function of time as obtained from DDFT and simulations, respectively, for each of the four types of flow, and for channel widths L_y  =  5d and L_y  =  6d. Additionally, in figure 8, we show state diagrams for the same two channel widths from both simulations and DDFT, where we show the type of flow observed for a range of investigated values of b and $\Gamma$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.2.1. Blockade.

5.2.2. Smooth flow.

5.2.3. Three or five-particle Oscillation.

For intermediate to strong particle interactions and for intermediate to wide constrictions we observe strong oscillatory behavior in the particle flow. While in the Brownian dynamics simulations the oscillation is damped we can find for the DDFT results an undamped oscillation that is periodic after a brief transient phase. This can be understood from the fact that the damping in the simulation is due to the presence of fluctuations, which are missed in the mean-field approach of the DDFT. We expect that the fluctuations which destroy long-ranged periodic order in one dimension are also responsible for washing out the correlations in the flow dynamics. The frequency of the oscillation depends on the number of particle layers in the system. For a five layer system the frequency is lower and corresponds to five particles passing the constriction during one oscillation period. In contrast, in the six layer system we observe a higher frequency in the flow oscillation, corresponding to three particles passing the constriction. In figure 9 we illustrate the mechanism that is responsible for this qualitative difference. For both channel widths, the oscillation period represents the smallest number of particles that can pass through the constriction in such a way that the system reverts to its original configuration. In the case of an odd number of layers (i.e. five), this is simply five particles, such that the crystalline lattice shifts by one lattice spacing. For an even number of crystalline layers (i.e. L_y  =  6d), this period instead corresponds to three particles passing through the constriction, such that the lattice moves by half a lattice spacing, and then coincides with a vertically mirrored version of the initial lattice. We have confirmed with simulations that the same mechanism occurs for other (small) numbers of layers.

Zoom In Zoom Out Reset image size

In the supplementary material stacks.iop.org/JPhysCM/28/244019/mmedia we include two movies of these dynamics in a system with L_y  =  6d and b  =  0.8 as obtained from DDFT and simulations.

As can be seen in figure 8, in our DDFT findings this type of oscillatory flow is dominant in a significantly larger region of the (b, $\Gamma$) parameter space. In the simulations, this mechanism only occurs around $b\simeq 0.8$. Likely, this can be attributed to the approximative excess functional which cannot account for complex crystal configurations that are responsible for the blockade in front of the constriction.

5.2.4. Transient oscillation.