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.

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.

The superadiabatic-DDFT, presented in detail in Ref. 41, consists of a pair of differential equations for the coupled time-evolution of the one- and two-body densities. It is applicable to systems with pairwise interparticle interactions. The first equation of superadiabatic-DDFT is the exact expression for the one-body density,
(1)
where β=(kBT)1, D0 is the diffusion coefficient, ϕ is the interparticle pair potential, r12 = |r1r2|, and Vext is a time-dependent external potential. The second equation is an approximate equation of motion for the two-body density, given by
(2)
where the superadiabatic contribution to the two-body density is defined according to
(3)
The adiabatic two-body density, ρad(2), is obtained by evaluating the equilibrium two-body density functional at the instantaneous one-body density,
(4)
The adiabatic potential, Vad, appearing in (2) generates a fictitious external force that stabilizes the adiabatic system. This is obtained from the Yvon–Born–Green (YBG) relation of equilibrium statistical mechanics46,47
(5)
applied to the nonequilibrium system. We note that the approximate equation for the two-body density (2) becomes exact in the low density limit.
We now specialize to a bulk system, for which the external potential is set equal to zero and the one-body density becomes constant, ρ(r, t) → ρb. In this case, the adiabatic two-body density is both isotropic and translationally invariant. Consequently, the integral term in Eq. (5) vanishes, and the adiabatic potential becomes constant in space. This constant can then be set equal to zero without loss of generality. To exploit the translational invariance of the system, we now define the following relative and absolute coordinates:
(6)
that then allow us to make the replacements
(7)
(8)
where we henceforth drop the subscripts on the nabla. Due to translational invariance, ρ(2) will not depend on the absolute position coordinate.
Since we wish to investigate systems under homogeneous flow, we introduce the affine velocity field, v, which, for shear applied in the direction of the x-axis, with shear-gradient in the y-direction, is given by
(9)
where γ̇ is the shear-rate (as illustrated in Fig. 1). From Eq. (2), we thus obtain the following equation of motion for the pair distribution function:
(10)
in which we note the emergence of the pair diffusion constant, 2D0, and where the nonequilibrium pair distribution function is given by
(11)
Its superadiabatic component, which encodes the flow-induced distortion of the microstructure, is defined as
(12)
where the equilibrium radial distribution function, geq, can be obtained from the interaction pair potential using an appropriate equilibrium theory (we will use an integral equation closure).
FIG. 1.

Sketch of the geometry. We choose particle 1 as the origin of our Cartesian coordinate system (r1 is thus implicitly fixed equal to zero).

FIG. 1.

Sketch of the geometry. We choose particle 1 as the origin of our Cartesian coordinate system (r1 is thus implicitly fixed equal to zero).

Close modal
Finding a solution to Eq. (10) for arbitrary values of γ̇ and ρb is a difficult task. Even in the low density limit, considerable effort must be expended to obtain numerical solutions due to the emergence of a boundary-layer in g as the shear-rate is increased.6 We will henceforth restrict our attention to the special case of steady-shear, γ̇(t)γ̇, for which the time-independent, steady-state pair distribution function, g(r, t) → g(r), can be obtained (almost) analytically for arbitrary values of ρb. In the steady-state, Eq. (10) reduces to
(13)
In equilibrium, v(r) = 0, and the pair-current, ∇gsup(r) + gsup(r)∇βϕ(r), vanishes. This condition leads trivially to the solution gsup(r) = 0, which implies that g(r) = geq(r), as expected.
To obtain a low shear-rate solution of the steady-state equation (13) we assume gsup to be a linear function of γ̇. (This will be sufficient for our present purposes. However, we note that care should be exercised when assuming linearity since boundary-layer formation can significantly complicate the picture.19,48–50) Substituting (12) into (13) and neglecting terms quadratic and higher in γ̇ then yields the following linearized steady-state condition:
(14)
In order to correctly capture the anisotropy induced by the flow, we make the following ansatz:18 
(15)
where we have introduced the rate-of-strain tensor, E, and the isotropic radial function f(r). The rate-of-strain tensor contains all relevant information about the affine flow field, while the function f(r) depends only on the interaction potential, the bulk density, and the system dimensionality. (Additional information regarding the rate-of-strain tensor is provided in  Appendix A.)
Working through the substitution of expression (15) into Eq. (14), as described in  Appendix B, finally yields the radial balance equations required to determine f(r). These are given by
(16)
in two- and three-dimensions, respectively. The boundary conditions required to solve Eq. (16) depend on the pair potential under consideration and are discussed in detail in  Appendix C.

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.

For a bulk system in a steady-state, the nonequilibrium pair distribution function is related to the stress tensor according to the following exact expression:3 
(17)
where I is the unit tensor. Due to the isotropy of the equilibrium pair distribution function, only gsup contributes to the off-diagonal stress tensor elements.
The interaction part of the viscosity, η, is obtained by dividing σxy by the shear rate. In the limit of low shear-rate, γ̇0, Eq. (17) gives the low shear viscosity, η0. Substituting Eq. (15) into Eq. (17) and using the appropriate form for the rate-of-strain tensor in shear flow, given by Eq. (A5), then yields
(18)
Both the function f and the integral in (18) depend on the dimensionality of the system. We then find
(19)
in two- and three-dimensions, respectively. We have introduced the two-dimensional area fraction, Φ2D = πρb/4, and the three-dimensional volume fraction, Φ3D = πρb/6, where all lengths have been non-dimensionalized using the characteristic diameter of the particles. If the function f were known exactly, then Eq. (19) would yield the exact low-shear viscosity. Within the present approach, the adiabatic approximation employed to close the two-body equation of motion (2) yields approximate radial balance equations (16) for f and, therefore, an approximate low-shear viscosity.
For the well-studied case of low density hard-spheres in three-dimensions, Eq. (16) has the analytic solution f(r) = 1/(3r3), and Eq. (19) recovers the quadratic term in the well-known low density expansion (see p. 114 in Ref. 2),
(20)
which applies in the absence of hydrodynamic interactions between the particles. The first term in (20) is the solvent viscosity, ηs. The second term arises from the drag of the solvent on the surface of each individual sphere.51 The third term, which is exact for a system without hydrodynamic interactions, represents the influence of direct potential interactions between the particles and comes from evaluating the integral in Eq. (19) and then employing the Stokes-Einstein relation kBT/D0 = 3πηs.

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.

In Fig. 2, we show scatter plots of the rescaled superadiabatic contribution to the distorted pair correlation function, namely,
(21)
where r=x2+y2 and f(r) is obtained from the numerical integration of Eq. (16). Results are shown for four different values of the bulk density, ρb, one per panel, and are valid at low shear-rates. At the lowest considered bulk density, ρb = 0.1, the radial function f(r) obtained numerically is very close to its low-density limit value, f(r) = 1/2r2, and no packing oscillations are visible in the resulting scaled gsup. In this panel, we observe an accumulation of particles at contact, r = d, in the “compressional” quadrants, defined as the region where sign(y) = −sign(x), and an opposite depletion effect in the “extensional” quadrants, where sign(y) = sign(x). As the value of the bulk density increases, packing oscillations develop at larger values of r. These oscillations in gsup take both positive and negative values and are a nontrivial prediction of the radial balance equation (16). The packing structure of the full pair correlation function, g(r) = geq(r) + gsup(r), will thus differ from that observed in bulk as a consequence of the applied shear. However, our numerical solutions reveal that the contact value of the radial function remains unaffected by changes in ρb and is given by
(22)
in each of the panels shown in Fig. 2.
FIG. 2.

Superadiabatic-DDFT gsup. For hard-disks in two-dimensions, we show the superadiabatic contribution to the pair correlation function, gsup, in units of γ̇d2/2D0. Since we consider only low shear-rates, gsup exhibits quadripolar symmetry [see Eq. (15)]. Increasing the bulk density leads to packing oscillations.

FIG. 2.

Superadiabatic-DDFT gsup. For hard-disks in two-dimensions, we show the superadiabatic contribution to the pair correlation function, gsup, in units of γ̇d2/2D0. Since we consider only low shear-rates, gsup exhibits quadripolar symmetry [see Eq. (15)]. Increasing the bulk density leads to packing oscillations.

Close modal

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.

FIG. 3.

Low-shear viscosity. A comparison of η0 from superadiabatic-DDFT [see Eq. (23)], with a fit to simulation data taken from Ref. 54. The low-shear viscosity from simulation diverges as the system approaches random close packing, whereas the superadiabatic-DDFT retains the low density limiting form for all values of the area fraction.

FIG. 3.

Low-shear viscosity. A comparison of η0 from superadiabatic-DDFT [see Eq. (23)], with a fit to simulation data taken from Ref. 54. The low-shear viscosity from simulation diverges as the system approaches random close packing, whereas the superadiabatic-DDFT retains the low density limiting form for all values of the area fraction.

Close modal
Due to the density independence of the contact value (22), the above-described route to obtaining the interaction contribution to the low-shear viscosity recovers its low-density form, namely,
(23)
Accounting for the influence of shear-flow on the bulk three-body density would give corrections to higher-order in area fraction in Eq. (23). However, this mechanism is neglected per construction within the current approximation. The fact that we recover only the leading-order contribution (23) to the low-shear viscosity is thus not surprising and is consistent with the application of the adiabatic approximation at the two-body level. It seems very likely that improved predictions for η0 could be obtained from superadiabatic-DDFT by employing the “test-particle method,” as we will discuss in Sec. IV.

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.

Integrating the many-particle Smoluchowski equation over N − 2 particle coordinates yields the following, formally exact equation of motion for the two-body density41,44
(24)
This equation can also be rewritten in the following alternative form:
(25)
which remains fully equivalent to Eq. (24), since there is no approximation involved at this stage.
Using the standard form (24) of the exact equation of motion and applying the adiabatic approximation solely to the three-body density, i.e., ρ(3)ρad(3), yields
(26)
where the adiabatic three-body density is defined as
(27)
in analogy to Eq. (4). Substitution of the second-order Yvon–Born–Green (YBG2) equation41,46
(28)
into (26) yields the second equation of superadiabatic-DDFT, namely Eq. (2). This is the procedure followed in Ref. 41.
An alternative method to implement the adiabatic closure is to start with the (reformulated) exact equation of motion (25) and to apply the approximation to the integral term,
(29)
We note that Eq. (29) is no longer equivalent to the previous Eq. (26). The essential difference between these two options is that in Eq. (26) we approximate the joint probability density, whereas in Eq. (29) we approximate the conditional probability density. Substitution of the rewritten YBG2 equation
(30)
into expression (29), then yields the following alternative equation of motion for the two-body density:
(31)
We refer to this approximation strategy as “Russel–Gast-like,” since Eq. (31) reduces to the Russel–Gast equation (see Ref. 10) for the pair distribution function in the bulk limit, as we will demonstrate in Subsection III B.
Following the same scheme as in Subsection II B, the bulk limit of Eq. (31) is given by
(32)
This provides an alternative to Eq. (10). [As a side-note, we mention that the second term on the right-hand side of Eq. (32) can be rewritten using the following rearrangement:
where we have introduced the “potential of mean force,” ϕmf(r), defined by geq(r)eβϕmf(r).] In the low density limit, geq(r) → eβϕ(r), which yields
(33)
such that Eqs. (32) and (10) then become equivalent (which is not the case in general).
In the steady-state at finite densities, Eq. (32) reduces to the condition
(34)
Using the definition (12) and again assuming that gsup is linear in the flow-rate yields the following expression to leading order:
(35)
Since Eqs. (14) and (35) are non-equivalent, the solution of Eq. (35) now requires a different ansatz than used previously. The appropriate choice in the present case is given by
(36)
where the ⋆ notation makes it explicit that the (still unknown) radial function, f(r), is different from the function, f(r), previously encountered. Substitution of ansatz (36) into the linearized steady-state equation (35) then yields the alternative radial balance equations
(37)
to determine f(r) for the cases of two- and three-dimensions. We refer the reader to  Appendix D for additional details regarding the derivation of Eq. (37) and the boundary conditions required to solve them.
Within the alternative adiabatic closure, the low-shear viscosity is given by
(38)
where we have introduced the so-called “cavity distribution function,”
(39)
familiar from liquid-state theory.46 We emphasize that f(r) appearing in Eq. (38) is given by solution of the radial balance equations (37) and is distinct from the function f(r) obtained from solution of Eq. (16). In two- and three-dimensions, Eq. (38) becomes
(40)
respectively. For the special case of low density hard-spheres in three-dimensions, the analytic solution of Eq. (37) is the same as that given by Eq. (16), namely
(41)
such that the Russel–Gast-type closure reproduces the exact low density expansion of the low-shear viscosity of hard-spheres, given by Eq. (20).
We now investigate some numerical predictions of the alternative closure, where we again consider the hard-disk system in two-dimensions and use the Percus–Yevick closure to obtain geq for input to the radial balance equation (37). In Fig. 4, we show scatter plots of the quantity
(42)
for four different values of the input bulk density. Note that the radial function f(r) is obtained from numerical integration of Eq. (37) and that r=x2+y2. The predictions of the Russel–Gast-type closure are generally very similar to those obtained from the superadiabatic-DDFT approach of Sec. II, although deviations become apparent on closer inspection.
FIG. 4.

Russel–Gast-type closure for gsup. Analogous plot to Fig. 2, but this time employing the alternative closure. For hard-disks, we show the superadiabatic contribution to the pair correlation function, gsup, in units of γ̇d2/2D0 [see Eq. (36)]. Since the results are qualitatively similar to those in Fig. 2 and, therefore, difficult to compare, we refer the reader to Fig. 5 for a more detailed analysis.

FIG. 4.

Russel–Gast-type closure for gsup. Analogous plot to Fig. 2, but this time employing the alternative closure. For hard-disks, we show the superadiabatic contribution to the pair correlation function, gsup, in units of γ̇d2/2D0 [see Eq. (36)]. Since the results are qualitatively similar to those in Fig. 2 and, therefore, difficult to compare, we refer the reader to Fig. 5 for a more detailed analysis.

Close modal

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.

FIG. 5.

Pair distribution function. We show the full nonequilibrium pair distribution function of hard-disks along the extensional axis, y = x, and the compressional axis, y = −x. Results are given from both the superadiabatic-DDFT and the Russel–Gast-type approximations for ρb = 0.7. The additional black curve indicates the equilibrium radial distribution function for comparison. The dimensionless shear-rate is chosen to be γ̇d2/2D0=1. The inset shows the difference between the nonequilibrium and equilibrium pair distribution functions to highlight the difference between the two approximations. Due to the symmetry of the low-shear pair distribution function for a given approximation, the curves in the compressional and extensional directions are identical (up to a sign).

FIG. 5.

Pair distribution function. We show the full nonequilibrium pair distribution function of hard-disks along the extensional axis, y = x, and the compressional axis, y = −x. Results are given from both the superadiabatic-DDFT and the Russel–Gast-type approximations for ρb = 0.7. The additional black curve indicates the equilibrium radial distribution function for comparison. The dimensionless shear-rate is chosen to be γ̇d2/2D0=1. The inset shows the difference between the nonequilibrium and equilibrium pair distribution functions to highlight the difference between the two approximations. Due to the symmetry of the low-shear pair distribution function for a given approximation, the curves in the compressional and extensional directions are identical (up to a sign).

Close modal

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 

FIG. 6.

Low-shear viscosity. Analogous plot to Fig. 3, but now including the prediction of the Russel–Gast-type approximation. This additional result predicts slightly larger values of η0 than the superadiabatic-DDFT but exhibits an unphysical curvature as the area fraction is increased.

FIG. 6.

Low-shear viscosity. Analogous plot to Fig. 3, but now including the prediction of the Russel–Gast-type approximation. This additional result predicts slightly larger values of η0 than the superadiabatic-DDFT but exhibits an unphysical curvature as the area fraction is increased.

Close modal

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.

The authors have no conflicts to disclose.

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).

The data that support the findings of this study are availablefrom the corresponding author upon reasonable request.

In a system undergoing translationally invariant, homogeneous flow, the velocity field can be conveniently expressed in the following form:
(A1)
where κ is the spatially constant velocity gradient tensor. As an example, choosing flow in the x-direction with shear-gradient in the y-direction enables κ to be represented in matrix form as
(A2)
for a three-dimensional system with shear-rate γ̇.
The right-hand side of (A1) can be decomposed into a sum of two terms,
(A3)
where
(A4)
are the (symmetric) rate-of-strain tensor and the (antisymmetric) rate-of-rotation tensor, respectively. The rate-of-strain tensor describes “pure straining motion” leading to relative motion between any two particles, whereas the rate-of-rotation tensor yields pure rotational motion, which does not affect their relative separation.65 For the aforementioned case of shear flow, we obtain
(A5)
The anisotropy of the nonequilibrium pair distribution function in a system subject to a slow translationally invariant flow is determined solely by its straining motion and does not involve its rotational component.6,18 This motivates our choice of ansatz for gsup in Eqs. (15) and (36).

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.

We start by collecting a few useful identities. For a general, spatially varying, second-rank tensor M, the following holds:
(B1)
where r̂=r/r is a unit vector, r̂r̂ is a dyadic product, and a(r) is a scalar function of r = |r|. For an arbitrary, spatially dependent vector field, u, we have
(B2)
where the double dot notation in the second term indicates a full contraction, i.e., a scalar product followed by a trace operation. The divergence of the product a(r)M is given by
(B3)
This identity generates the following special cases in two-dimensions:
(B4)
(B5)
and three-dimensions
(B6)
(B7)
respectively. In the special case where M is symmetric, it can be shown that
(B8)
Finally, the divergence of the scalar product Mr̂ is given by the aesthetically appealing identity,
(B9)
We will now employ these identities to calculate the pair distribution function at low flow rates.
Substitution of the ansatz (15) into the linearized steady-state equation (14) yields
(B10)
in which the only unknown quantity is the function f(r). Assuming incompressible flow, ∇ · v = 0, using identity (B1) and rewriting the projected velocity as
(B11)
enables us to re-express Eq. (B10) as follows:
(B12)
The advantage of this representation is the appearance of the dyadic tensors, r̂r̂ and (Ir̂r̂), which project either along or perpendicular to the relative position vector r.
On the right-hand side of (B12), we have to calculate the divergence of the scalar product between the vector r̂E and a second-rank tensor. We can thus exploit relation (B2) to obtain
(B13)
We will consider separately the two terms appearing on the right-hand side of (B13), which we henceforth refer to as (i) and (ii).
To simplify (i) in two-dimensions, we use (B3)(B5). This yields
(B14)
The analogous result in three-dimensions is given by
(B15)
where we have used equations (B3), (B6), and (B7).
The simplification of term (ii) does not depend on the dimensionality of the system. We first employ Eq. (B8) to re-express the factor (Er̂) as a divergence and then use identity (B9) to obtain
(B16)
since TrE=0 for incompressible flow.
Finally, substituting (B14) and (B16) into (B13) yields the desired two-dimensional result
(B17)
while substitution of (B15) and (B16) into (B13) gives the result in three-dimensions,
(B18)
Since these expressions remain valid for all choices of the traceless tensor E, the coefficients of the quadratic form (r̂Er̂) must be equal. We thus obtain the radial balance equations (16) stated in the main text. We note that, in addition to the explicit appearance of the pair interaction potential in (B17), respectively (B18), both the pair potential and the bulk density enter these equations implicitly via the equilibrium radial distribution function.
The radial balance equations (16) are second-order in spatial derivative, and their solution thus requires the specification of two boundary conditions. From the linearized steady-state equation (14) we can identify the following expression for the pair-current at low shear-rates:
(C1)
Pairs of particles separated by a large distance are not spatially correlated, which implies that the radial component of the pair-current should tend to zero. Moreover, if we assume that the interparticle interaction potential has a strongly repulsive core, then we can impose that the radial component of the pair-current will also go to zero at small separations. We thus have the boundary condition
(C2)
for both r → ∞ and r → 0.
For the case of hard-disks in two-dimensions or hard-spheres in three-dimensions, the second boundary condition takes a special form since the radial pair-current must vanish when two particles touch in order to prevent unphysical overlap. The small separation boundary condition then becomes
(C3)
where we have set the particle diameter equal to unity. Using steps directly analogous to those leading from (B10) to (B12) enables the pair-current of hard-spheres to be expressed in the following form:
(C4)
for r > 1. [Note that ϕ(r > 1) = 0 for the hard-sphere potential.] Taking the scalar product of Eq. (C4) with the radial unit vector yields
(C5)
where we have used (B11). The hard-sphere “zero-flux” condition at particle contact thus becomes
(C6)
The radial balance equations (16), combined with (C6) and f(r → ∞) = 0, fully determine the function f(r), given that the input equilibrium radial distribution function is known.

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.

Substitution of the ansatz (36) into the linearized steady-state equation (35) yields
(D1)
from which we can determine the function f(r). Assuming incompressible flow, ∇ · v = 0, and using the identities (B1) and (B11), we can re-express Eq. (D1) in the following form:
(D2)
Steps analogous to those leading from Eqs. (B12) to (B18) then generate the alternative form of the radial balance equations (37).
From the linearized steady-state equation (35), we identify the pair-current at low shear-rates,
(D3)
A calculation analogous to that leading from (B10) to (B12) shows that Eq. (D3) can be re-expressed as
(D4)
For a general interaction potential, the boundary condition (C2) holds for both r → ∞ and r → 0.
In the special case of hard-spheres, the two-body current must be equal to zero when a pair of particles come into contact. Applying the boundary condition (C2) at r = 1 yields
(D5)
which is different from Eq. (C6). The boundary condition at particle contact, Eq. (D5), together with f(r → ∞) = 0, then fully determines the solutions of the radial balance equations (37).
1.
M. E.
Cates
and
M. R.
Evans
,
Soft and Fragile Matter
(
CRC Press
,
2000
).
2.
J.
Mewis
and
N. J.
Wagner
,
Colloidal Suspension Rheology
(
Cambridge University Press
,
2012
).
3.
J. M.
Brader
, “
Nonlinear rheology of colloidal dispersions
,”
J. Phys.: Condens. Matter
22
,
363101
(
2010
).
4.
R.
Larson
,
Structure and Rheology of Complex Fluids
,
Topics in Chemical Engineering
(
Oxford University Press
,
New York
,
1999
).
5.
G. K.
Batchelor
,
An Introduction to Fluid Dynamics
,
Cambridge Mathematical Library
(
Cambridge University Press
,
2000
).
6.
J.
Bergenholtz
,
J. F.
Brady
, and
M.
Vicic
, “
The non-Newtonian rheology of dilute colloidal suspensions
,”
J. Fluid Mech.
456
,
239
275
(
2002
).
7.
J. F.
Brady
and
J. F.
Morris
, “
Microstructure of strongly sheared suspensions and its impact on rheology and diffusion
,”
J. Fluid Mech.
348
,
103
139
(
1997
).
8.
T.
Phtsuki
, “
Dynamical properties of strongly interacting Brownian particles: 1. Dynamic shear viscosity
,”
Physica A
108
,
441
(
1981
).
9.
D.
Ronis
, “
Theory of fluctuations in colloidal suspensions undergoing steady shear flow
,”
Phys. Rev. A
29
,
1453
(
1984
).
10.
W. B.
Russel
and
A. P.
Gast
, “
Nonequilibrium statistical mechanics of concentrated colloidal dispersions: Hard spheres in weak flows
,”
J. Chem. Phys.
84
,
1815
(
1986
).
11.
N. J.
Wagner
and
W. B.
Russel
, “
Nonequilibrium statistical mechanics of concentrated colloidal dispersions: Hard spheres in weak flows with many-body thermodynamic interactions
,”
Physica A
155
,
475
(
1989
).
12.
J. F.
Brady
, “
The rheological behavior of concentrated colloidal dispersions
,”
J. Chem. Phys.
99
,
567
(
1993
).
13.
R. A.
Lionberger
and
W. B.
Russel
, “
A Smoluchowski theory with simple approximations for hydrodynamic interactions in concentrated dispersions
,”
J. Rheol.
41
,
399
(
1997
).
14.
R. A.
Lionberger
and
W. B.
Russel
, “
Effectiveness of nonequilibrium closures for the many body forces in concentrated colloidal dispersions
,”
J. Chem. Phys.
106
,
402
(
1997
).
15.
G.
Szamel
, “
Nonequilibrium structure and rheology of concentrated colloidal suspensions: Linear response
,”
J. Chem. Phys.
114
,
8708
(
2001
).
16.
R. A.
Lionberger
and
W. B.
Russel
, “
Microscopic theories of the rheology of stable colloidal dispersions
,” in
Advances in Chemical Physics
(
John Wiley & Sons, Ltd.
,
1999
), pp.
399
474
.
17.
J.
Morris
, “
A review of microstructure in concentrated suspensions and its implications for rheology and bulk flow
,”
Rheol. Acta
48
,
909
(
2009
).
18.
W. B.
Russel
,
D. A.
Saville
, and
W. R.
Schowalter
,
Colloidal Dispersions
,
Cambridge Monographs on Mechanics
(
Cambridge University Press
,
1989
).
19.
L.
Banetta
,
F.
Leone
,
C.
Anzivino
,
M. S.
Murillo
, and
A.
Zaccone
, “
Microscopic theory for the pair correlation function of liquidlike colloidal suspensions under shear flow
,”
Phys. Rev. E
106
,
044610
(
2022
).
20.
M.
Fuchs
and
M. E.
Cates
, “
Theory of nonlinear rheology and yielding of dense colloidal suspensions
,”
Phys. Rev. Lett.
89
,
248304
(
2002
).
21.
J. M.
Brader
,
T.
Voigtmann
,
M. E.
Cates
, and
M.
Fuchs
, “
Dense colloidal suspensions under time-dependent shear
,”
Phys. Rev. Lett.
98
,
058301
(
2007
).
22.
J. M.
Brader
,
M. E.
Cates
, and
M.
Fuchs
, “
First-principles constitutive equation for suspension rheology
,”
Phys. Rev. Lett.
101
,
138301
(
2008
).
23.
M.
Fuchs
and
M. E.
Cates
, “
A mode coupling theory for Brownian particles in homogeneous steady shear flow
,”
J. Rheol.
53
,
957
(
2009
).
24.
J. M.
Brader
,
M. E.
Cates
, and
M.
Fuchs
, “
First-principles constitutive equation for suspension rheology
,”
Phys. Rev. E
86
,
021403
(
2012
).
25.
M.
Frank
,
D.
Anderson
,
E. R.
Weeks
, and
J. F.
Morris
, “
Particle migration in pressure-driven flow of a Brownian suspension
,”
J. Fluid Mech.
493
,
363
378
(
2003
).
26.
J.
Brader
and
M.
Krüger
, “
Density profiles of a colloidal liquid at a wall under shear flow
,”
Mol. Phys.
109
,
1029
(
2011
).
27.
M.
Krüger
and
J. M.
Brader
, “
Controlling colloidal sedimentation using time-dependent shear
,”
Europhys. Lett.
96
,
68006
(
2011
).
28.
A. A.
Aerov
and
M.
Krüger
, “
Driven colloidal suspensions in confinement and density functional theory: Microstructure and wall-slip
,”
J. Chem. Phys.
140
,
094701
(
2014
).
29.
A.
Scacchi
,
M.
Krüger
, and
J. M.
Brader
, “
Driven colloidal fluids: Construction of dynamical density functional theories from exactly solvable limits
,”
J. Phys.: Condens. Matter
28
,
244023
(
2016
).
30.
H.
Jin
,
K.
Kang
,
K. H.
Ahn
, and
J. K. G.
Dhont
, “
Flow instability due to coupling of shear-gradients with concentration: Non-uniform flow of (hard-sphere) glasses
,”
Soft Matter
10
,
9470
(
2014
).
31.
J. C.
Crocker
and
D. G.
Grier
, “
Methods of digital video microscopy for colloidal studies
,”
J. Colloid Interface Sci.
179
,
298
(
1996
).
32.
V.
Prasad
,
D.
Semwogerere
, and
E. R.
Weeks
, “
Confocal microscopy of colloids
,”
J. Phys.: Condens. Matter
19
,
113102
(
2007
).
33.
R.
Besseling
,
L.
Isa
,
E. R.
Weeks
, and
W. C.
Poon
, “
Quantitative imaging of colloidal flows
,”
Adv. Colloid Interface Sci.
146
,
1
(
2009
).
34.
D. R.
Foss
and
J. F.
Brady
, “
Brownian dynamics simulation of hard-sphere colloidal dispersions
,”
J. Rheol.
44
,
629
(
2000
).
35.
A.
Sierou
and
J. F.
Brady
, “
Accelerated Stokesian dynamics simulations
,”
J. Fluid Mech.
448
,
115
146
(
2001
).
36.
A. J.
Banchio
and
J. F.
Brady
, “
Accelerated Stokesian dynamics: Brownian motion
,”
J. Chem. Phys.
118
,
10323
(
2003
).
37.
Y.
Yurkovetsky
and
J. F.
Morris
, “
Triplet correlation in sheared suspensions of Brownian particles
,”
J. Chem. Phys.
124
,
204908
(
2006
).
38.
S. M.
Tschopp
,
H. D.
Vuijk
,
A.
Sharma
, and
J. M.
Brader
, “
Mean-field theory of inhomogeneous fluids
,”
Phys. Rev. E
102
,
042140
(
2020
).
39.
P.
Attard
, “
Spherically inhomogeneous fluids. I. Percus–Yevick hard spheres: Osmotic coefficients and triplet correlations
,”
J. Chem. Phys.
91
,
3072
(
1989
).
40.
K.
Nygård
,
R.
Kjellander
,
S.
Sarman
,
S.
Chodankar
,
E.
Perret
,
J.
Buitenhuis
, and
J. F.
van der Veen
, “
Anisotropic pair correlations and structure factors of confined hard-sphere fluids: An experimental and theoretical study
,”
Phys. Rev. Lett.
108
,
037802
(
2012
).
41.
S. M.
Tschopp
and
J. M.
Brader
, “
First-principles superadiabatic theory for the dynamics of inhomogeneous fluids
,”
J. Chem. Phys.
157
,
234108
(
2022
).
42.
S. M.
Tschopp
,
H. D.
Vuijk
, and
J. M.
Brader
, “
Superadiabatic dynamical density functional study of Brownian hard-spheres in time-dependent external potentials
,”
J. Chem. Phys.
158
,
234904
(
2023
).
43.
U. M. B.
Marconi
and
P.
Tarazona
, “
Dynamic density functional theory of fluids
,”
J. Chem. Phys.
110
,
8032
(
1999
).
44.
A. J.
Archer
and
R.
Evans
, “
Dynamical density functional theory and its application to spinodal decomposition
,”
J. Chem. Phys.
121
,
4246
(
2004
).
45.
J.
Chakrabarti
,
J.
Dzubiella
, and
H.
Löwen
, “
Dynamical instability in driven colloids
,”
Europhys. Lett.
61
,
415
(
2003
).
46.
J.
Hansen
and
I.
McDonald
,
Theory of Simple Liquids
(
Elsevier Science
,
2006
).
47.
D.
McQuarrie
,
Statistical Mechanics
(
University Science Books
,
2000
).
48.
J. K. G.
Dhont
, “
On the distortion of the static structure factor of colloidal fluids in shear flow
,”
J. Fluid Mech.
204
,
421
431
(
1989
).
49.
S.
Riva
,
L.
Banetta
, and
A.
Zaccone
, “
Solution to the two-body Smoluchowski equation with shear flow for charge-stabilized colloids at low to moderate Péclet numbers
,”
Phys. Rev. E
105
,
054606
(
2022
).
50.
L.
Banetta
and
A.
Zaccone
, “
Radial distribution function of Lennard-Jones fluids in shear flows from intermediate asymptotics
,”
Phys. Rev. E
99
,
052606
(
2019
).
51.
J. F.
Brady
, “
The Einstein viscosity correction in n dimensions
,”
Int. J. Multiphase Flow
10
,
113
(
1983
).
52.
F.
Lado
, “
Numerical Fourier transforms in one, two, and three dimensions for liquid state calculations
,”
J. Comput. Phys.
8
,
417
(
1971
).
53.
F.
Lado
, “
Equation of state of the hard-disk fluid from approximate integral equations
,”
J. Chem. Phys.
49
,
3092
(
1968
).
54.
J.
Reinhardt
,
F.
Weysser
, and
J. M.
Brader
, “
Density functional approach to nonlinear rheology
,”
Europhys. Lett.
102
,
28011
(
2013
).
55.
B. D.
Goddard
,
A.
Nold
, and
S.
Kalliadasis
, “
Dynamical density functional theory with hydrodynamic interactions in confined geometries
,”
J. Chem. Phys.
145
,
214106
(
2016
).
56.
J.
Vermant
and
M. J.
Solomon
, “
Flow-induced structure in colloidal suspensions
,”
J. Phys.: Condens. Matter
17
,
R187
(
2005
).
57.
J.
Bergenholtz
, “
Theory of rheology of colloidal dispersions
,”
Curr. Opin. Colloid Interface Sci.
6
,
484
(
2001
).
58.
J.
Bławzdziewicz
and
G.
Szamel
, “
Structure and rheology of semidilute suspension under shear
,”
Phys. Rev. E
48
,
4632
(
1993
).
59.
E.
Nazockdast
and
J. F.
Morris
, “
Microstructural theory and the rheology of concentrated colloidal suspensions
,”
J. Fluid Mech.
713
,
420
452
(
2012
).
60.
M.
Rex
and
H.
Löwen
, “
Dynamical density functional theory with hydrodynamic interactions and colloids in unstable traps
,”
Phys. Rev. Lett.
101
,
148302
(
2008
).
61.
M.
Rex
and
H.
Löwen
, “
Dynamical density functional theory for colloidal dispersions including hydrodynamic interactions
,”
Eur. Phys. J. E
28
,
139
(
2009
).
62.
M.
Rauscher
, “
DDFT for Brownian particles and hydrodynamics
,”
J. Phys.: Condens. Matter
22
,
364109
(
2010
).
63.
B. D.
Goddard
,
A.
Nold
,
N.
Savva
,
P.
Yatsyshin
, and
S.
Kalliadasis
, “
Unification of dynamic density functional theory for colloidal fluids to include inertia and hydrodynamic interactions: Derivation and numerical experiments
,”
J. Phys.: Condens. Matter
25
,
035101
(
2012
).
64.
R. D.
Mills-Williams
,
B. D.
Goddard
, and
A. J.
Archer
, “
Dynamic density functional theory with inertia and background flow
,”
J. Chem. Phys.
160
,
174901
(
2024
).
65.
J.
Dhont
,
An Introduction to Dynamics of Colloids
(
Elsevier Science
,
1996
).