The superadiabatic dynamical density functional theory (superadiabatic-DDFT) is a promising new method for the study of colloidal systems out-of-equilibrium. Within this approach, the viscous forces arising from interparticle interactions are accounted for in a natural way by explicitly treating the dynamics of the two-body correlations. For bulk systems subject to spatially homogeneous shear, we use the superadiabatic-DDFT framework to calculate the steady-state pair distribution function and the corresponding viscosity for low values of the shear-rate. We then consider a variant of the central approximation underlying this superadiabatic theory and obtain an inhomogeneous generalization of a rheological bulk theory due to Russel and Gast. This paper thus establishes for the first time a connection between DDFT approaches, formulated to treat inhomogeneous systems, and existing work addressing nonequilibrium microstructure and rheology in bulk colloidal suspensions.
I. INTRODUCTION
Colloidal suspensions exhibit a rich variety of rheological behavior, arising from the interplay between Brownian motion, solvent hydrodynamics, and potential interactions.1–3 For example, the phenomena of shear thinning, shear thickening, and yielding are relevant for many commercial products and industrial processes. In order to control and tune the rheological properties of a suspension for any specific application, it is necessary to have an understanding of how the microscopic interactions between the constituents influence the macroscopic response.4 The challenge for nonequilibrium statistical mechanics is to formulate robust and accurate first-principles theories based on tractable approximation schemes that capture the essential physics while remaining sufficiently simple for concrete calculations to be performed.
For the most commonly studied situation, namely bulk suspensions under homogeneous shear, there are currently a variety of microscopic approaches. Each of these aims to capture a particular aspect of cooperative particle motion within a limited range of shear-rates and thermodynamic parameters but fails to provide a unified global picture. While exact results can be obtained for low density systems at low shear-rate,5–7 systems at intermediate7–19 or high densities20–24 invoke a diverse range of approximate closure relations to account for the correlated motion of the particles. Inhomogeneous systems, for which the density and shear-rate vary in space, are more challenging to treat theoretically, and appropriate closure relations that correctly capture the coupling between gradients in shear-rate and density remain under development.25–30
The clearest path to a first-principles theory of suspension rheology is to focus on the particle correlation functions and then integrate these to obtain the macroscopic rheological properties of interest. This makes it possible to connect the macroscopic constitutive relations to the nonequilibrium microstructure of the system and, therefore, gain microscopic insight into the physical mechanisms at work.2,5 For systems with pairwise additive interparticle interactions, the two-body correlations are the primary objects of interest. In the absence of hydrodynamic interactions, these enable full calculation of the stress tensor, which is the key quantity of interest in rheology. Additional motivation to develop theories that “look inside” the flowing system is provided by developments in the visualization and tracking of particle motion in experiments (confocal microscopy),25,31–33 together with the detailed information provided by computer simulations of model systems under flow.34–37
Theoretical studies of inhomogeneous fluids, both in- and out-of-equilibrium, are primarily based on the spatially varying one-body density alone and do not usually directly involve the inhomogeneous two-body correlations (although exceptions do exist38–42). This is a major difference between the standard dynamical density functional theory (standard DDFT) of inhomogeneous fluids43,44 and the aforementioned approaches to bulk colloidal rheology. When applied to bulk systems subject to homogeneous shear flow, the standard DDFT does not present, as is, a useful framework since the one-body density is not affected by the shear and remains constant in time. There have nevertheless been attempts to supplement the one-body equation of standard DDFT with additional empirical correction terms to avert this issue.26,27,45
However, a new DDFT framework has recently been developed, the so-called superadiabatic-DDFT,41 which treats explicitly the dynamics of the two-body correlations. This scheme presents a significant improvement in describing the dynamics of inhomogeneous fluids in the presence of time-dependent external potentials by accounting, via the two-body correlations, for the structural rearrangement of the particles as the system flows.41,42 The superadiabatic-DDFT (composed of a pair of coupled equations for the one- and two-body density) per construction avoids the shortcomings of standard DDFT since its equation for the two-body density is affected by shear and can thus be directly used to study bulk homogeneous systems. By predicting the shear-induced distortion of the pair correlations, superadiabatic-DDFT allows to obtain the viscosity as an output of the theory. This is not a quantity reachable in any standard DDFT treatment of the problem, but is a very relevant bridge to the world of rheology.
In this paper, we apply superadiabatic-DDFT to bulk systems under steady-shear flow and investigate its predictions for the shear-distorted pair distribution function and the low-shear viscosity. In addition, we show that a variation of the superadiabatic-DDFT reproduces, in the bulk limit, an early theory due to Russel and Gast10 and, therefore, provides a generalization of their approach to the case of inhomogeneous fluids in external fields. This establishes a clear connection between existing microscopic theories of bulk colloidal rheology and DDFT approaches to the dynamics of inhomogeneous fluids. The numerical predictions of this Russel–Gast-type approach are compared with those of supradiabatic-DDFT.
II. SUPERADIABATIC-DDFT
A. General framework
B. Low shear-rate solutions for the pair distribution function in bulk
Although we have chosen to focus on shear, we note that the expressions presented in this subsection remain valid for any incompressible, translationally invariant flow field and so could be applied directly to, e.g., extensional flows.
C. Low-shear viscosity
D. Numerical results
To investigate the predictions of the superadiabatic-DDFT approach, we will focus on the special case of hard-disks in two-dimensions, which presents a phenomenology qualitatively similar to that of hard-spheres in three-dimensions while remaining convenient for visualization of the distorted pair correlations in the xy-plane. The solution of the radial balance equation (16) requires as input the hard-disk equilibrium radial distribution function, geq. We obtain this quantity by employing the famous Percus–Yevick closure of the homogeneous Ornstein–Zernike equation,52,53 known to be accurate for the hard-disk system at low and intermediate bulk densities.
Knowledge of the nonequilibrium pair correlation function also enables the calculation of the interaction contribution to the low-shear viscosity, η0, via Eq. (18). Results obtained by solving the radial balance equation (16) and evaluating the integral (18) for each chosen value of the bulk density are shown in Fig. 3 and compared with a fit to Brownian dynamics simulation data taken from Ref. 54. While both curves overlap for low values of ρb, the theoretical superadiabatic-DDFT curve strongly underestimates η0 relative to the simulation at higher densities. We note that the simulation curve diverges as the system approaches random close packing, located at Φ2D ≈ 0.82 for monodisperse hard-disks.
III. RUSSEL–GAST-TYPE THEORY
In Sec. II, we analyzed the properties of superadiabatic-DDFT for bulk systems subject to homogeneous shear. All predictions of the theory for the pair distribution function and viscosity arise from the adiabatic closure employed to arrive at Eq. (2). In this section, we show that there exists an alternative way to implement the adiabatic closure and that this results in a new approximate equation of motion for the inhomogeneous two-body density, different from Eq. (2). In the bulk limit, the new approximation reduces to an equation of motion for the pair distribution function first introduced in 1986 by Russel and Gast,10 which constituted one of the earliest theoretical approaches to the microstructure and rheology of colloidal suspensions at a finite volume fraction. It thus seems appropriate to call our adiabatic closure of the equation of motion for the inhomogeneous two-body density a “Russel–Gast-type” approximation.
A. Alternative adiabatic closure
B. Low shear-rate solutions for the pair distribution function in bulk
C. Low-shear viscosity
D. Numerical results
The results for gsup at the lowest bulk density considered, ρb = 0.1, are very close to the known low-density limit, for which f⋆(r) = 1/2r2, with accumulation and depletion at contact within the compressional and extensional quadrants, respectively. As the bulk density increases, we observe the emergence of packing oscillations. We recall that superadiatic-DDFT predicted that the radial function f appearing in Eq. (23) has a contact value that remains independent of the bulk density. The analogous quantity within the present approximation is the product f⋆(r)geq(r) appearing on the right-hand side of Eq. (42). We find that the contact value of this product does exhibit a nontrivial dependence on ρb and, as we will see, this causes the low-shear viscosity generated by Eq. (40) to deviate from that predicted by superadiabatic-DDFT.
In Fig. 5, we show the pair distribution function on both the extensional axis, along which the particles get pulled apart by the shear, and on the compressional axis, along which the particles get pressed together. In order to show more clearly the difference between the superadiabatic-DDFT and Russel–Gast-type approximations, we show in the inset the difference between the nonequilibrium and equilibrium pair distribution functions. On the extensional axis, both approximation schemes predict a reduction in the height of all maxima, including the first peak at particle contact. The reduction in amplitude of the peaks is slightly more pronounced within the Russel–Gast-type approximation. In both cases, the radial position of the second and higher-order peaks shifts to larger values of r, consistent with the fact that particles are being pulled away from each other by the shear flow. In contrast, on the compressional axis, the height of all maxima increases, and the radial locations of the second and higher-order peaks are shifted to smaller values of r since the particles are being pushed closer together by the shear flow. The two theories again make similar predictions, although the Russel–Gast-type approximation yields a slightly larger increase in peak height than the superadiabatic-DDFT. The general similarity of the predictions from the two approximation schemes for the nonequilibrium microstructure is reassuring, as it appears that the results are robust with respect to the details of how the adiabatic approximation is implemented.
From the nonequilibrium pair correlation function, we can calculate the zero-shear viscosity using Eq. (40), for which results are shown in Fig. 6. We find that the Russel–Gast-type approximation scheme generates a low-shear viscosity curve with values slightly larger than those predicted by the superadiabatic-DDFT. We can attribute this effect to the difference in the contact value of the pair distribution function between the two approximations (see Fig. 5). Within the Russel–Gast-type approach, the difference between the contact value on the extensional axis and the contact value on the compressional axis is larger than for superadiabatic-DDFT. This larger extensional/compressional asymmetry has the consequence that for any given area fraction, η0 from the Russel–Gast-type approach is larger than for the superadiabatic-DDFT. Although the Russel–Gast-type closure produces values of the low-shear viscosity closer to the simulation results, it has the undesirable feature of changing curvature as the area fraction increases beyond about 0.45. This questions the physicality of the alternative approximation and its rescaling potential for higher densities.12
IV. DISCUSSION
In this paper, we have investigated the predictions of superadiabatic-DDFT for the nonequilibrium pair correlation function and low-shear viscosity of bulk systems subject to the homogeneous shear flow defined by Eq. (9). This provides a conceptual link between density functional-based approaches, which are focused on the inhomogeneous one-body density,43,44,55 and microscopic theories of homogeneous bulk rheology of the type pioneered by Brady, Russel, Wagner, and others.2,3,18 A clear connection between these two substantial bodies of literature, which have coexisted for decades with little to no interaction, is established by our Russel–Gast-type closure of the equation of motion for the two-body density. As detailed in Sec. III, this new scheme generates a fully inhomogeneous DDFT, which reduces to the known Russel–Gast theory for bulk systems subject to homogeneous shear.10
The application of shear distorts the bulk pair distribution function and, for systems interacting via a pairwise additive potential, enables the interaction part of the stress tensor to be calculated. For the low shear-rates considered here, the only relevant transport coefficient is the viscosity, since the lowest-order shear-induced changes to the diagonal stress tensor elements are quadratic in . For larger shear-rates, we would need to consider the extension of our solution ansatz (15) or (36), respectively, to higher-order in , which would then generate nonvanishing normal stress differences and shear-induced modifications to the system pressure.6,29 It would be interesting to investigate the extent to which numerical solutions of the superadiabatic-DDFT can account for nonlinear rheological and microstructural changes as the shear rate is increased. For example, it is known from both simulation37 and experiment56 that the nonequilibrium bulk three-body distribution function deviates significantly from its equilibrium form, which could be difficult to capture using an adiabatic closure on the two-body level.
Focusing on the low shear-rate regime enabled a largely analytic investigation of the nonequilibrium pair distribution function, with only the solution of the radial balance equation and evaluation of the low-shear viscosity demanding (relatively simple) numerical integration. At larger values of the shear-rate, more sophisticated techniques must be employed since boundary-layer formation prevents the application of straightforward perturbation theory.6,48–50,57 An interesting extension of our study to obtain results for higher shear-rates could be to employ either multipole methods,58 specialized numerical schemes,6,59 or intermediate asymptotics.19,49 Although the theories presented here are valid for an arbitrary pair interaction potential in two- or three-dimensions, we chose to perform numerical calculations for the special case of hard-disks, as this presents the simplest continuum model for the study of shear while still retaining much of the essential phenomenology of three-dimensional systems. In Refs. 19 and 49, it is shown how the pair distribution function under shear can be calculated analytically for more general interaction potentials.
The two routes presented in this work are built on two variations of the same approximation. The numerical output of those theories is broadly similar but, as clearly pictured by the final viscosity plot, Fig. 6, is different in more demanding situations at higher bulk density. One can thus ask which scheme is more likely to succeed in further investigations. The unphysical shape of the viscosity curve obtained with the Russel–Gast-like closure suggests that this variation of the approximation may be less robust than that of superadiabatic-DDFT and could prove more difficult to improve and refine.
A natural first step toward better predictions for the viscosity at higher bulk densities would be to apply Brady’s semi-empirical rescaling method.12 For the present Brownian dynamics, without hydrodynamic interactions, this consists of multiplying the low density limiting solution for the interaction contribution to η0 [given by Eq. (23) for hard-disks] by the contact value of the equilibrium radial distribution function. The resulting low-shear viscosity has been shown to agree well with simulation data for hard-spheres,16 and we find that this is also the case for hard-disks. Since the superadiabatic-DDFT predicts that the low-density limiting result (23) holds for all values of ρb, a rescaled version of the superadiatic-DDFT would exactly reproduce Brady’s theory for η0. In contrast, a rescaled Russel–Gast-type theory would continue to exhibit an unphysical curvature, as the contact value of geq is a monotonically increasing function of bulk density. Investigating how such a rescaling factor could emerge from a systematic, first-principles extension of the superadiabatic-DDFT closure scheme will be an interesting topic for future study.
While the above considerations may provide a path to structural improvements of the superadiabatic-DDFT equations (i.e., following from a new two-body closure), there remains a great deal to be learned from the current version of the theory. In this paper, we have considered only the direct application of superadiabatic-DDFT to the bulk system, meaning that we work from the outset with a constant one-body density and focus on the properties and predictions of the resulting two-body equation of motion. Although this sheds light on the internal structure and physical content of the theory, we should keep in mind that the superadiabatic-DDFT is still a density functional theory and, as such, aims primarily to predict the dynamics of the one-body density. If we are interested solely in bulk systems under shear, as in the present work, then the full power of superadiabatic-DDFT for the one-body dynamics can be harnessed by employing a more sophisticated implementation scheme: the test-particle method. Fixing a test-particle at the coordinate origin [i.e., setting Vext(r) = ϕ(r)] induces a spatially varying steady-state one-body density, ρ(r), as particles around the test-particle accumulate in the compressional quadrants and are depleted from the extensional quadrants.54 This inhomogeneous one-body density can be related to the bulk pair distribution function according to g(r) = ρ(r)/ρb, from which the viscosity can be calculated using Eq. (17). We suspect that a test-particle implementation of superadiabatic-DDFT will yield a low-shear viscosity that greatly improves the simple quadratic expression (23), obtained by direct application of the two-body equation of motion to bulk. The proposed test-particle calculation would require explicit treatment of the inhomogeneous two-body density in the presence of a test-particle and would thus capture, to some extent, the shear-induced distortion of the bulk three-body correlations. Investigations in this direction are underway.
Finally, we mention the issue of hydrodynamic interactions. Although the majority of standard DDFT studies do not consider solvent hydrodynamics, this aspect has been incorporated into the formalism.60–64 It would thus be interesting to investigate whether a similar extension would be feasible for superadiabatic-DDFT. Many of the existing theories of homogeneous bulk rheology mentioned in the introduction, including the bulk Russel–Gast theory, were formulated to include hydrodynamic interactions to some level of approximation. These works, all of which are focused on two-body correlations, could well provide a source of inspiration for future developments.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
S. M. Tschopp: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Visualization (lead); Writing – original draft (equal); Writing – review & editing (equal). J. M. Brader: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Resources (lead); Supervision (lead); Visualization (supporting); Writing – original draft (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are availablefrom the corresponding author upon reasonable request.
APPENDIX A: RHEOLOGICAL QUANTITIES
APPENDIX B: THE RADIAL BALANCE EQUATION
Substitution of the ansatz (15) into the linearized steady-state equation (14) yields the radial balance equations (16). Since this procedure is not straightforward, we outline here the main steps of the calculation. Although we are primarily interested in the case of homogeneous shear, we formulate the problem using the rate-of-flow tensor E, which can also describe other flow types. This not only extends the generality of our results but also presents the technical advantage that certain vector/tensor identities can be employed to make the calculation as clean as possible. Unless otherwise stated, all expressions given in this appendix are valid in arbitrary dimensionality.
APPENDIX C: BOUNDARY CONDITIONS
APPENDIX D: ALTERNATIVE FORMS OF THE RADIAL BALANCE EQUATION AND BOUNDARY CONDITIONS
In the preceding two Appendixes, we provided all the details required for the derivation and solution of the radial balance equations of superadiabatic-DDFT. Since the derivation of both the alternative radial balance equations (37) and their boundary conditions are very similar, we give here only the main equations to highlight the differences between the two approximation schemes.