## Abstract

A lattice Boltzmann (LB)-based hybrid method is developed to simulate suspensions of Brownian particles. The method uses conventional LB discretization (without fluid- level fluctuations) for suspending fluid, and treats Brownian particles as point masses with a stochastic thermal noise. LB equations are used to compute the velocity perturbations induced by the particle motion. It is shown that this method correctly reproduces the short-time and long-time diffusive behaviour of a Brownian particle. Unlike the earlier hybrid methods that use thermal fluctuations in the fluid, this method correctly reproduces the temperature of the particle and does not require an empirical rescaling of the bare friction coefficient to obtain the correct diffusive behaviour. It is observed that the present method is at least twice as fast as the earlier method. This method is best suited for flows of polymers and Brownian suspensions in microfluidic devices.

## 1. Introduction

The last decade has witnessed increased research activity in the area of design and development of lab-on-chip devices, mainly used to analyse blood and other body fluids which are suspensions of particles ranging from simple nanoparticles to complex long-chain molecules such as DNA [1,2]. Understanding the dynamical behaviour of such complex systems in these microfluidic devices is an important prerequisite to build well-designed devices. The main challenge to predict this behaviour is in the nonlinear interactions between the confined geometry, flow and the complex configurations of the suspended particles. In this paper, we present a new numerical method to correctly capture the dynamics of dilute Brownian suspension, which can be readily extended to predict the dynamics of long-chain polymer molecules in a flow in confined regions.

Hybrid methods based on the lattice Boltzmann method (LBM) are promising tools to simulate systems like lab-on-chip devices that assume complex geometries with an enormous number of particles. Such methods benefit from computational efficiency of the LBM hydrodynamic solver. Ladd [3] proposed one such hybrid method, coupling the fluctuating lattice Boltzmann method (FLBM) with molecular dynamics (MD) to simulate finite-sized spherical particle suspensions. In the FLBM–MD algorithm [3], the dynamics of the finite-sized spherical solute particles is governed by Newton’s equations of motion. No-slip boundary condition is implemented on the surface of the particles to couple the dynamics of the particle with that of the fluid. This hybrid method scales linearly with the number of particles *n* owing to the local nature of the fluid–particle coupling in space and time [3].

Ahlrichs & Dünweg [4,5] modified Ladd’s algorithm [3] to incorporate polymer systems. They modelled beads of a polymer chain as *point* particles and coupled particle–fluid dynamics using a friction model that is consistent with the global momentum conservation. Frictional coupling is appropriate here as, in polymer simulations, one is often interested only in the long-time behaviour of the system where the microscopic details of particle–fluid interaction are not important [5]. In this methodology, the solvent hydrodynamics is simulated using the FLBM algorithm given in Ladd [3], while the particle dynamics is described by the Langevin equation (FLBM–LD). Consequently, in the FLBM–LD method, the Brownian motion of particles is simulated by introducing thermal fluctuations both in the fluid and on the particles. These thermal fluctuations affect the diffusive nature of the particles and require renormalization of the diffusivity, obtained empirically [5]. In contrast, Brownian dynamics (BD) simulations treat the solvent implicitly and have a fluctuating force on the particles coupled to their frictional drag, and BD is equivalent to Langevin dynamics (LD) when the particle’s inertia can be neglected. In spite of modelling an isothermal and non-fluctuating solvent, the BD method captures the essential physics of dilute polymer suspensions [6–8] in the absence of complex flows and geometry interactions.

In the spirit of the non-fluctuating solvent of the BD method, we simplify the FLBM–LD algorithm for polymer systems. We neglect the fluctuations in the solvent and introduce an effective fluctuating force on the particle, unlike BD, which is coupled only to the local drag on it. We discuss various benchmark results from the simulation of a single particle immersed in a solvent to show that the simplified algorithm correctly reproduces its diffusion dynamics. We further show that the new algorithm is computationally efficient compared with the previous hybrid methods based on the LBM framework.

The rest of the article is organized as follows. In §2, we discuss the details of the proposed method, and present results for various validation cases in §3. In §4, we summarize the results, conclusions and next steps.

## 2. The lattice Boltzmann method–Langevin dynamics algorithm

In this section, we present a hybrid method that couples the LBM with LD (LBM–LD) to simulate the dynamics of Brownian suspensions.

### (a) The lattice Boltzmann method

The LBM is a mesoscopic method that involves solving the discrete velocity Boltzmann BhatnagarGrossKrook (BGK) equation [9], formulated such that in the long-time and large wavelength limit, the Navier–Stokes (NS) hydrodynamics is recovered. It is used as an alternative to NS methods in the incompressible limit to simulate a range of flows [9,10]. Within the LBM, along with space and time, velocity is discretized into a finite set of velocities, for which we use a D3Q19 lattice model (for explicit definition of velocity set **c**_{i} and associated weight *w*_{i}, see, for example, [9]). The resulting lattice Boltzmann equation (LBE) with the BGK collision function [11,12] reads
2.1
Here, *f*_{i} is the single particle probability distribution function for the discrete velocities in the *i*th direction, **c**_{i} is the lattice velocity vector, **r** is the spatial coordinate of a node, *δt* is the time step and *τ* is the relaxation time. The hydrodynamic quantities mass and momentum density are defined as
2.2
The functional form of the equilibrium distribution function, *f*^{eq}, used in this work is
2.3
where *c*_{s} is the speed of sound and is given in terms of lattice spacing *a* as . The relaxation time can be defined in terms of the kinematic viscosity *ν* as *τ*=*ν*/*c*^{2}_{s}. Here, we note that the polynomial expression of *f*^{eq} used in equation (2.3) is for computational convenience, and it is possible to define equilibrium as a minimizer of the entropy function [13] even for the discrete BGK model.

### (b) Langevin dynamics

The Langevin equation (LE) for a particle in a solvent is given by
2.4
where **F**^{c}_{p} is a conservative (potential) force on the particle and **F**^{d}_{p} is the drag force acting on the particle owing to its velocity relative to the surrounding fluid. **F**^{d}_{p} is modelled through a Stokes frictional model,
2.5
where *ζ* is the friction coefficient, **v**_{p} is the velocity of the point particle and **v**_{f}(**r**_{p}) is the velocity of the fluid at the location of particle **r**_{p}. The fluid velocity **v**_{f}(**r**_{p}) is a distance-weighted average velocity of the neighbouring fluid nodes. The stochastic force **F**^{r}_{p} satisfies the fluctuation–dissipation theorem (FDT),
2.6
Here, the vector **F**^{r}_{p} is written in the index notation and *k*_{B}*T* is the temperature of the fluid bath. The local velocity of a fluid node is determined by the LB equation (2.1). Apart from this, the fluid in the neighbourhood of a particle is also influenced by the motion of the particle. In order to conserve the local momentum, the net momentum imparted on the fluid by the particle is calculated as
2.7
A distance-weighted local fluid momentum is then distributed back to the neighbour nodes in the fluid [4]. The corrected momentum of the node (from the LBE propagation step and accumulated influence of particle motion steps) is then used to calculate the equilibrium distribution function *f*^{eq}.

In our simulations, we use two distinct time steps: *δt* for the LB fluid update, and *δt*_{p} for the particle integration such that *δt*≥*δt*_{p}. Since, in the present method, the fluid fluctuations are neglected, we denote this method as LBM–LD in contrast to the FLBM–LD method used in Ahlrichs & Dünweg [4].

### (c) Parameters of simulation

Parameters to simulate the Brownian motion of a single particle in a fluid are chosen from the benchmark chain model values suggested in Ahlrichs & Dünweg [5]: the solvent density *ρ*=0.864, kinematic viscosity *ν*=2.8 and temperature *k*_{B}*T*=1.20. Lattice spacing *a*=1 and lattice time step *δt*=0.01 are used. The friction coefficient *ζ*=20.8. The effective hydrodynamic diameter of the point particle, *d*_{p}=*ζ*/3*πνρ*=0.9123, is used in the calculation of the mass density *ρ*_{p} of the particle. The velocity Verlet integrator time step for the particle equation of motion is taken as *δt*_{p}=0.001, such that in one time step of the LBM the particle equations are integrated 10 times.^{1} The mass of the Brownian particle is set to match the required mass density ().

## 3. Results

The validity of the algorithm to capture the diffusive hydrodynamics can be verified by two essentially equivalent methods: (i) FDT and (ii) long-term diffusivity. We show the results of both these tests below. We further note that the simulation code is verified for consistency of equipartition of the temperature (same temperature in the three Cartesian directions equal to the set temperature *k*_{B}*T*).

### (a) Fluctuation–dissipation theorem

The FDT states that the linear response of a system to an external perturbation is the same as the response of the system to thermal fluctuations at equilibrium. This can be verified by comparing the response of a particle to an initial velocity perturbation in a quiescent fluid with the spontaneous velocity fluctuations of a Brownian particle, using the velocity autocorrelation function (VACF). The VACF of an externally perturbed particle decays to zero taking an exponential form , where the friction coefficient *ζ* is given by the Stokes drag law.

A particle immersed in a pool of liquid (without any thermal fluctuations) is initially given a velocity **v**_{p,0}. Its time evolution is simulated by Newton’s equation of motion (with **F**^{c}_{p} and **F**^{r}_{p} set to zero in equation (2.4)) along with the fluid dynamics through the LB. We apply periodic boundary conditions (PBCs) on the boundaries of the cubic simulation box of size *L*=50*a* to minimize the effect of periodicity on the particle dynamics. The lines shown in figure 1 correspond to the normalized VACF (**v**_{p}(0)⋅**v**_{p}(*t*))/(**v**_{p}(0)⋅**v**_{p}(0)) obtained in the linear response regime for externally perturbed particles of two different masses.^{2} Over short times, the VACF relaxes exponentially, and matches exactly with the theoretically expected decay (not shown). The decay in the case of a heavier particle is slower. The exponential decay is followed by an algebraic decay of approximately *t*^{−3/2} over long times owing to the unsteady flow behaviour of the surrounding fluid [14].

A particle experiencing Brownian thermal fluctuations from the surrounding fluid medium is simulated in the new algorithm by retaining the fluctuation forces governed by equation (2.6) in the LE (equation (2.4)). The normalized VACF of a Brownian particle is obtained by an ensemble average 〈**v**_{p}(0)⋅**v**_{p}(*t*)〉/〈**v**_{p}(0)⋅**v**_{p}(0)〉 over a set of independent runs. The points in figure 1 are the normalized VACF obtained for the Brownian particles with the same mass as used earlier. A good match is observed between the VACF of the Brownian particle and the VACF of the externally perturbed particle in the exponential decay region, as shown in figure 1, thus verifying the FDT (the long-time VACF regions of the Brownian particle have poor resolution because of insufficient sampling).

### (b) Long-time diffusivity

The long-time diffusivity of a Brownian particle depends on the Stokes drag coefficient *ζ* of the particle: . The long-time diffusivity can also be obtained from the mean square displacement (MSD) of the Brownian particle, using the Einstein relation, . If the correlation in the fluctuation is correctly accounted for, the long-time diffusivity obtained by the MSD should agree with that obtained from the assumed friction coefficient *ζ*. This is expected to hold good in the present case, since the diffusivity can be obtained from the VACF through a Green–Kubo formula, and we have shown the VACF to follow the correct exponential decay with the assumed *ζ*.

In figure 2, we plot scaled MSD: 〈[**r**_{p}(*t*)−**r**_{p}(0)]^{2}〉/6*t* as a function of time *t*. The steady-state value is thus the diffusivity . This agrees to a high degree of accuracy with the value of obtained through the Stokes–Einstein relation: , computed using the assumed drag coefficient *ζ* on the particles.

### (c) Comparison between the fluctuating lattice Boltzmann method–Langevin dynamics and the lattice Boltzmann method–Langevin dynamics

The two validation tests demonstrate that the proposed LBM–LD algorithm correctly reproduces the behaviour of a Brownian particle immersed in an isothermal fluid. We now compare the FLBM–LD [5] and LBM–LD algorithms.

We compare the diffusivity estimated using the FLBM–LD algorithm^{3} with that of the LBM–LD in figure 2. The diffusivity value computed using the FLBM–LD method overestimates the theoretical Stokes–Einstein value (from assumed *ζ*). Ahlrichs & Dünweg [5] attribute this discrepancy to the back flow of the fluid surrounding the particle. In order to match the Stokes–Einstein estimate, they introduce an empirically determined effective friction coefficient *ζ*_{eff},
3.1
Here, *g*=25 is an empirically determined constant, and *η* is the viscosity of the fluid. In figure 2, the corrected theoretical estimate *k*_{B}*T*/*ζ*_{eff} is shown to match the diffusivity obtained from the MSD calculations. With the present LBM–LD algorithm, no such empirical adjustment is required, as was shown in figure 2.

### (d) Efficiency

We now compare the computational efficiencies of the FLBM–LD and the LBM–LD methods based on the time required to simulate systems for 1000 fluid time steps. We simulate systems of different sizes, ranging from *L*=25*a* to *L*=100*a*, while maintaining the number density of the particles constant across the runs. A typical C++ implementation of each of the methods was run in serial mode on an Intel Xeon quadcore, 3 GHz processor. Figure 3 shows a comparison of the relative costs of central processing unit time required for both the methods, and it can be observed that the LBM–LD method is approximately twice as fast as the FLBM–LD method. In the FLBM–LD method, the fluid-level fluctuation component of the FLBM needs six Gaussian random numbers per time step to be generated that satisfy the FDT on the stress tensor for a lattice site. In addition to this, calculation of the effect of stress fluctuations on the probability distribution function (*f*) needs projection of the stress tensor onto the *f* space. This involves projecting the stress (3×3 matrix) onto 19 lattice vectors of the D3Q19 lattice model [15]. However, in the LBM–LD method, we do not simulate fluid fluctuations explicitly. This factor makes the LBM–LD algorithm computationally efficient while simulating the correct diffusive behaviour.

## 4. Conclusions

In the present work, we have proposed a new hybrid method coupling the LBM with the LD for an efficient simulation of Brownian suspensions that is capable of handling confined geometries with ease. The proposed method (LBM–LD) correctly captures (i) the FDT and (ii) the long-time diffusivity of a Brownian particle.

The LBM–LD method obtains the correct diffusive behaviour for the Brownian particle immersed in a solvent unlike the FLBM–LD method, which needs an empirical rescaling of the friction coefficient to match the theoretical diffusivity predicted by the Stokes–Einstein relation. In the FLBM–LD method, the stress fluctuations introduced in the fluid induce extra mobility in the particles, thereby requiring renormalization.

Besides simulating the correct diffusive behaviour, the LBM–LD method is computationally efficient compared with the FLBM–LD. The LBM–LD method runs approximately twice as fast as the FLBM–LD method. In future, we intend to test the ability of the new LBM–LD method to reproduce the diffusive behaviour of linear polymer chains in a dilute system.

## Acknowledgements

We acknowledge the support of CRL, Pune, in conducting this work. M.M. thanks Ajay Nandgaonkar (CRL) for useful discussions and reviewing the article, and also thanks Rochish M. Thaokar (IITB) and Ganesh Subramanian (JNCASR) for their valuable suggestions. S.A. and P.S. acknowledge the fruitful discussions with Burkhard Dünweg, Ravi P. Jagadeeshan and Hans C. Öttinger.

## Footnotes

One contribution of 26 to a Theme Issue ‘Discrete simulation of fluid dynamics: methods’.

↵1 We also used various values for

*δt*_{p}ranging from 0.0001 to 0.01, but observe no significant difference in the results.↵2 We use perturbation velocity

*v*_{p}(0) = 0.001 and 0.01*a*/*δt*with particles of density (*ρ**) 1 and 10, respectively.↵3 Our code implementation permits treatment of the fluid fluctuations with a Boolean choice variable.

- This journal is © 2011 The Royal Society