key: cord-0136687-vtvhu3hd
authors: Pousaneh, Faezeh; Wijn, Astrid de
title: Kinetic theory and shear viscosity of dense dipolar hard sphere liquids
date: 2020-03-25
journal: nan
DOI: nan
sha: 4e0c9a0d76919737717401480bbbf9aa2cca5bfb
doc_id: 136687
cord_uid: vtvhu3hd

Transport properties of dense fluids are fundamentally challenging, because the powerful approaches of equilibrium statistical physics cannot be applied. Polar fluids compound this problem, because the long-range interactions preclude the use of a simple effect-diameter approach based solely on hard spheres. Here, we develop a kinetic theory for dipolar hard-sphere fluids that is valid up to high density. We derive a mathematical approximation for the radial distribution function at contact directly from the equation of state, and use it to obtain the shear viscosity. We also perform molecular-dynamics simulations of this system and extract the shear viscosity numerically. The theoretical results compare favorably to the simulations.

Transport properties of dense fluids are fundamentally challenging, because it is a many-body problem out of equilibrium. A statistical approach is needed, but the powerful approaches of equilibrium statistical physics cannot be applied. Current theoretical approaches to transport in dense fluids are based on hard spheres and Enskog's heuristic extension of the Boltzmann equation and kinetic theory of gasses and liquids [1] . The only alternative to this is to resort to purely computational methods (see, for example [2] [3] [4] [5] ). Kinetic theory was heavily developed in the 60s and 70s, but little progress has been made since. In particular, there is no analytical description of high-density fluids consisting of anything more complicated than simple hard spheres (HS). This is a fundamental limitation in our current understanding, but also particularly problematic in practical applications, where kinetic theory is widely used in combination with empirical information and effective diameters to predict viscosities of some non-polar complex liquids [6, 7] .

Here, we develop kinetic theory of polar fluids, especially focusing on the viscosity. Polar fluids are a textbook example of systems where the hard-sphere approach fails, because the long-range electrostatic interactions are captured badly by instantaneous collisions. They are also ubiquitous in nature, for example in the form of water, and are increasingly important in applications in biotechnology and other fields. The physics of these systems cannot be described by simple hard spheres with an effective diameter. Moreover, Molecular Dynamics (MD) simulations involving electrostatic interactions are extremely computationally demanding and anyway cannot provide the fundamental understanding that is needed.

We choose to focus on the shear viscosity, as it is one of the most important transport properties of a fluid for practical applications. It plays a crucial role in for exam- * Electronic address: faezeh.pousaneh@ntnu.no. † Electronic address: astrid.dewijn@ntnu.no. ple lubrication and pipe flow. The viscosity of polar fluids is receiving increasing interest in practical applications, for example as the basis of environmentally-friendly lubricants, and are very promising for low-friction applications, as demonstrated for instance by the amazing effectiveness with which water-based synovial fluid lubricates our joints [8, 9] .

We derive an analytical kinetic theory for the viscosity of a simple model for a polar fluid, dipolar hard spheres (DHS). Our theoretical approach is based around Enskog's extension of the Boltzmann equation to high densities (BEk). In order to incorporate the soft and longrange electrostatic interactions between the dipoles, we extend this theory, which is originally based on simple shapes and simple interactions especially HS at low densities. We do this by explicitly including the dipole-dipole interaction into the Radial Distribution Function (RDF). We calculate the RDF using the method of [10] [11] [12] [13] from the Helmholtz free energy of DHS derived by Elfimova et al. [14] . In order to verify our theoretical results, we compare them to MD simulation of dipolar pseudo hard spheres. Our result can be used in a straight-forward manner to also calculate other transport properties such as thermal conductivity and diffusion coefficient.

BEk theory centers around the Boltzmann equation and deals with collision probabilities and collision dynamics. Boltzmann's original equation contains a crucial low-density approximation: the Stoßzahlansatz, which states that when particles collide they are uncorrelated. Solving the Boltzmann equation for transport coefficients is nontrivial, but general solutions were derived by Chapman and Enskog [15] . The general form of the zerodensity viscosity is found to be

where Ω * (2,2) is the collision integral which depends on the interactions. For HS, Ω * (2,2) = 1. With considerable effort, zero-density viscosities can also be derived for slightly more complicated interaction models, such as rough spheres [15] , spherocylinders [16] , and hard spheres with embedded point dipoles (DHS) [17] . At higher density, the equations for collision rates and dynamics become more complex, in general including correlated collisions. Enskog devised a heuristic way to incorporate some correlated collisions at higher densities, but this approach is currently limited to HS [15] and chains of hard spheres [18] , but not other types of interactions. Enskog's approach produces good agreements with simulations of HSs and experiments of very simple fluids only for low-to mid-density ranges and fails at high densities, since it still does not take into account correlated collisions. Nevertheless, Enskog's theory, though still approximate in nature, has provided a useful theoretical basis for both understanding and predicting the transport properties of hydrocarbons with short-range interactions only, including some molecules with much more complex geometry than HS [18] [19] [20] .

In order to obtain theoretical results for the viscosity of DHS, we start from the Enskog's theory for a simple dense fluid. Enskog's expression for the viscosity is [1, 15, [21] [22] [23] 

where V excl is the excluded volume of HS, V excl = (2π/3)σ 3 , ξ = πσ 3 ρ/6 is the volume fraction, and g(ξ) is the RDF at contact. The RDF in is the spherical component of the pair-distrition function. There are a number of ways to obtain good approximations for the RDF at contact of HS, such as from the Carnahan-Starling equation [24] , which gives

The zero-density limit for viscosity [see Eq. 1] for some polar interactions have been obtained. In Ref. [25] , the collision integral for the zero-density viscosity of polar gas was calculated for the Stockmayer potential. Chung et al. developed an empirical formula which works well for the viscosity of real dilute gasses [26] . Our approach for high densities is to develop the radial distribution function of DHS and apply it to the Enskog expression. The interaction between two DHSs i and j with diameter σ and dipole moments µ at distance r is given by a sum of hard sphere (U HS ij ) and dipolar (U D ij ) terms:

with the dipolar coupling constant λ = µ 2 /(k B T 4πǫ 0 σ 3 ).

In recent years, considerable effort has been made on development of the theoretical expression for the equilibrium properties of DHS [14, [27] [28] [29] [30] [31] [32] . The Helmholtz free energy of DHS can be written relative to that for a regular HS fluid F HS as

where F D is the excess free energy due to the electrostatic interaction between the dipoles. The most common approaches for dealing with DHS is thermodynamic perturbation theory with a Pade approximation and mean spherical approximation (MSA). However, because these are lower-order theories with respect to λ they do not give accurate results for low densities and virial coefficients [30, 31] .

In order to get around this problem, Elfimova et al. [30] introduced a logarithmic representation of the free energy. The result converges faster, since the logarithm of a polynomial is less sensitive to the truncation of the polynomial. The excess free energy is then written as [14] βF

The coefficients I n are obtained from the regular virial coefficients for DHS. Elfimova et al. [14] keep up to the fifth virial coefficient, corresponding to n = 4 and give explicit expressions for I 1,2,3,4 . This theory accurately captures the free energy and compares favorably with computer simulation for λ ≤ 4, even at high value of the particle volume fraction ξ ≤ 0.5. We obtain the RDF from the DHS free energy using the equation of state (EOS) [11, 33] ,

where U pot is the interaction potential. We apply the thermodynamic relations P = − ∂F ∂V | N,T and U int = ∂(βF ) ∂β to obtain the pressure and internal energy. The interaction potential is then obtained as

∂λ . Finally, we find the RDF at contact

where The viscosity is then obtained by substituting this into Eq. (2). We compare our theoretical results to MD simulations. The integration algorithms typically used for MD depend on smooth interaction, and cannot be applied to instantaneous collisions. This is worsened by the presence of long-range electrostatic interactions, which require additional techniques that also depend on smooth force fields. We circumvent this issue by employing a pseudo hard sphere model (PHS) introduced by Jover et al. [34] . The PHS potential is of Mie form where the typical powers of the LJ potential 12/6 are replaced by 50/49:

Jover et al. verified that this potential accurately captures the thermodynamics, structures, and dynamics of the HS system. It produces good results at reduced temperature T * = ǫ kB T = 2/3. This model has been shown to accurately describe the fluid-solid equilibrium [35] as well as the viscosity of HS [36] .

We use Gromacs version 5 to integrate the equations of motion and the PHS potential is implemented as a tabular form as in Ref. [36] . Our DHS consist of 5 particles on a line, as shown in Fig. 1 . The central particle has no charge or mass, but interacts with the central particles of the other DHS through a PHS potential. Two massless particles of opposite charges q and −q are at equal distance L q /2 from the center on either side and give rise to the dipole. There are also two dummy mass particles of mass m on each side at distance L m /2, controlling the moment of inertia. We use L q /σ = 0.22, since the point dipole model has been found to agree well with the extended dipole model up to L q /σ = 0.3 [5, 37, 38] . We simulate this system for different dipole moments corresponding to λ = 1, 2, 3, 4. Our simulation box contains N = 1000 DHS particles and N = 6000 for dilute cases, ρσ < 0.15. All simulations have been carried out at reduced temperature T * = ǫ/k B T = 2/3. For [14] (circle data). Solid lines are the theoretical expression in Ref. [14] .

each systems with different λ we perform simulations for a range of densities ρ * between 0 and 1. In what follows, all units are dimensionless as:

, where ρ and P denote number density and pressure respectively and η is viscosity. The reduced volume fraction is ξ * = πρ * 6 . The electrostatic interactions are treated using the Particle Mesh Ewald (PME) method with cut-off length of 2.6σ * . Time steps for simulations is δt * = 0.0011.

We first equilibrate the system and verify the equation of state (EOS), before moving on to the RDF and viscosity. Equilibration was performed in the NVT ensemble with the velocity-rescale thermostat for t * = 10 5 to t * = 6 · 10 5 depending on the system. Fig. 2 shows the EOS for DHS obtained from current simulations (blue triangle data), from previous Monte-Carlo simulations [14] (red circle data) and the theoretical expression of EOS in Ref. [14] . Our MD simulations correspond well to both.

After equilibration, we run the simulations in the NVT ensemble for an additional interval of t * = 1000 and obtain the full RDF for each density and λ. The RDF at contact is given by the maximum values of RDF. Simulation results of RDF at contact are shown in Fig. 3 (a) along with our theoretical expression Eq. (10), both with σ = 1. For comparison, Fig. 3(b) shows the results obtained by Rushbrooke et al. using the Pade approximation [39, 40] . The new theory developed here describes the simulation results significantly better and captures the trends relative to the Carnahan-Starling results. We continue to run the system in the NVT ensemble for t * = 45000 up to t * = 90000 depending on the system (for more dilute ones longer time is needed to get enough collisions). To minimize the influence of the thermostat, the temperature is controlled using a Berendsen thermostat with a slow coupling with a characteristic time of t * = 11. We obtain the shear viscosity of the DHS model using the transverse-current auto-correlation function (TCAF) method [41] . More details on our use of this method can be found in [36] .

The shear viscosities obtained from the simulations are shown in Fig. 4(a) for different λ. According to the data, the relative shear viscosity decreases upon increasing the dipole moment, whereas it shows an opposite behavior for higher densities. This is because at lower densities the dipolar particles form chain-like structure which decreases the collision rates, and consequently the viscosity. At high densities strong dipole moments cause the system to form ordered structures, which have a higher viscosity (up to infinity) than a noninteracting disordered fluid. In addition, Enskog theory for hard spheres is known to break down at higher densities even in fluids. To visualize the reason for why the viscosities deviate from the Enskog expressions for hard spheres without dipoles, we show examples of snapshots from simulations in Fig. 5 . The snapshots are for system with λ = 4 and two different densities ξ * = 0.03 and ξ * = 0.49. Similar structures are reported by simulations for dipolar fluids (and ferromagnetic particles) [42] [43] [44] [45] [46] and also by experiments [47] [48] [49] [50] .

We compare the simulation results to the theoretical results, which are given by Eq. (10) combined with (2) and (1). Since we do not have the exact low-density limit for η 0 for our system, we introduce Ω * (2,2) as a fit parameter. We estimate the range of sensible values from the results of Ref. [25] for the Stockmayer potential with a point dipole, to be in the range of 1 to 3 corresponding to λ up to λ = 5. We use an effective dipole moment µ e as a fit parameter, rather than the hard-sphere diameter. The fit parameters Ω * (2, 2) and µ e /µ are given in Table. increase by increasing the dipole moments in agreement with the trends found for the zero-density viscosity of the Stockmayer potential [25] and values in the expected range. The only previously available theory for viscosity of dense fluids is the HS Enskog theory. In order to compare our theory to this, we fit the simulations data of viscosity to the Enskog theory for HS, Eq. (2), with the HS RDF given in Eq. (3). The collision integral Ω * (2,2) is equal to unity for HS, but this is incorrect for DHS. When Enskog theory for HS is applied to real molecules this is usually taken into account by allowing Ω * (2,2) to deviate from unity and using it as a fit parameter, along with the effective diameter σ e , and we do the same here. The results are shown in Fig. 4 (b) as lines. The fit parameters Ω * (2, 2) and σ e /σ are given in Table. I. Fig. 4 clearly shows that our theory successfully describes the viscosity of dense fluids of DHS and captures qualitative behaviour that is not captured by previous HS theoretical results.

In summary, we have developed a kinetic theory for the shear viscosity of dense fluids of dipolar hard spheres (DHS). In our theory, we have included the long-range electrostatic interactions explicitly. Our theory captures the main effects of the dipole-dipole interaction on the viscosity, which were missing from previous theories. We see from our simulations that the differences between DHS and HS are mainly due to local structure. At low densities the DHS viscosity is lower due to clustering of the particles. At high densities, the DHS show orientational ordering, leading to stronger interaction and a higher viscosity. Both of these effects are captured by our theory. Our theory is in agreement with simulation results for packing fractions below about 0.35 − 0.4.

While we have focused on the viscosity, the RDF at contact is the crucial ingredient for the collision rate and consequently the density-dependence of all nonequilibrium properties of fluids. Our kinetic theory should therefore also provide for accurate descriptions of other transport properties, such as thermal conductivity and diffusion coefficient. Moreover, the approaches currently in use in applications, for viscosity as well as other transport coefficients, are all based on the simple HS results, even for much more complicated molecules. Besides the fundamental understanding of transport in polar fluids, our theory can thus also lead to significant improvements in the accuracy of calculations of transport properties in practical applications.

Viscosity of Liquids; Theory, Estimaation, Experiment, and Data

Viscosity of Dense Fluids. New York: Plenum

The Mathematicaltheory of non-uniform gases

A Concise Course on the Theory of Classical Liquids

Equations of State for Fluids and Fluid Mixtures. Elsevier

Viscosity of Dense Fluids

The work has been supported by National Infrastructure for Computational Science in Norway (UNINETT Sigma2) with computer timed for the Center for High Performance Computing (NN9573K and NN9572K). The authors acknowledge The Research Council of Norway for NFR project number 275507 and The Faculty of Engineering, Norwegian University of Science and Technology (NTNU), for financial support. FP acknowledges Prof. Ekaterina Elfimova and Prof. Philip Camp for their advise and provided data.