Superadiabatic dynamical density functional theory for colloidal suspensions under homogeneous steady-shear

Tschopp, S. M.; Brader, J. M.

doi:10.1063/5.0211198

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.⁴ 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 intermediate^7–19 or high densities^20–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 exist^38–42). This is a major difference between the standard dynamical density functional theory (standard DDFT) of inhomogeneous fluids^43,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,⁴¹ 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 Gast¹⁰ 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

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,

\begin{align} \frac{1}{D_{0}} \frac{\partial ρ (r_{1}, t)}{\partial t} & = \nabla_{r_{1}} \cdot (\nabla_{r_{1}} ρ (r_{1}, t) + ρ (r_{1}, t) \nabla_{r_{1}} β V_{ext} (r_{1}, t) \\ + \int d r_{2} ρ^{(2)} (r_{1}, r_{2}, t) \nabla_{r_{1}} β ϕ (r_{12})), \end{align}

(1)

where

β = {(k_{B} T)}^{- 1}

⁠, D₀ is the diffusion coefficient, ϕ is the interparticle pair potential, r₁₂ = |r₁ − r₂|, and V_ext is a time-dependent external potential. The second equation is an approximate equation of motion for the two-body density, given by

\begin{align} \frac{1}{D_{0}} \frac{\partial ρ^{(2)} (r_{1}, r_{2}, t)}{\partial t} & = \sum_{i = 1, 2} \nabla_{r_{i}} \cdot (\nabla_{r_{i}} ρ_{\sup}^{(2)} (r_{1}, r_{2}, t) \\ + ρ_{\sup}^{(2)} (r_{1}, r_{2}, t) \nabla_{r_{i}} β ϕ (r_{12}) \\ + ρ^{(2)} (r_{1}, r_{2}, t) \nabla_{r_{i}} β V_{ext} (r_{i}) \\ - ρ_{ad}^{(2)} (r_{1}, r_{2}, t) \nabla_{r_{i}} β V_{ad} (r_{i}, t)), \end{align}

(2)

where the superadiabatic contribution to the two-body density is defined according to

ρ_{\sup}^{(2)} (r_{1}, r_{2}, t) \equiv ρ^{(2)} (r_{1}, r_{2}, t) - ρ_{ad}^{(2)} (r_{1}, r_{2}, t) .

(3)

The adiabatic two-body density,

ρ_{ad}^{(2)}

⁠, is obtained by evaluating the equilibrium two-body density functional at the instantaneous one-body density,

ρ_{ad}^{(2)} (r_{1}, r_{2}, t) \equiv ρ_{eq}^{(2)} (r_{1}, r_{2}; [ρ (r, t)]) .

(4)

The adiabatic potential, V_ad, appearing in generates a fictitious external force that stabilizes the adiabatic system. This is obtained from the Yvon–Born–Green (YBG) relation of equilibrium statistical mechanics^46,47

\begin{align} - \nabla_{r_{1}} V_{ad} (r_{1}, t) & \equiv k_{B} T \nabla_{r_{1}} \ln ρ (r_{1}, t) + \int d r_{2} \frac{ρ_{ad}^{(2)} (r_{1}, r_{2}, t)}{ρ (r_{1}, t)} \\ \times \nabla_{r_{1}} ϕ (r_{12}), \end{align}

(5)

applied to the nonequilibrium system. We note that the approximate equation for the two-body density becomes exact in the low density limit.

B. Low shear-rate solutions for the pair distribution function in bulk

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

r = r_{1} - r_{2}, R = r_{1} + r_{2},

(6)

that then allow us to make the replacements

\nabla_{r_{1}} \to \nabla_{r} = \nabla,

(7)

\nabla_{r_{2}} \to - \nabla_{r} = - \nabla,

(8)

where we henceforth drop the subscripts on the nabla. Due to translational invariance, ρ⁽²⁾ 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

v (r, t) = \dot{γ} (t) y e_{x},

(9)

where

\dot{γ}

is the shear-rate (as illustrated in Fig. 1). From Eq. , we thus obtain the following equation of motion for the pair distribution function:

\begin{align} \frac{\partial g (r, t)}{\partial t} & = - 2 D_{0} \nabla \cdot (- \nabla g_{\sup} (r, t) - g_{\sup} (r, t) \nabla β ϕ (r)) \\ - \nabla \cdot (g (r, t) v (r, t)), \end{align}

(10)

in which we note the emergence of the pair diffusion constant, 2D₀, and where the nonequilibrium pair distribution function is given by

g (r, t) = \frac{ρ^{(2)} (r, t)}{ρ_{b}^{2}} .

(11)

Its superadiabatic component, which encodes the flow-induced distortion of the microstructure, is defined as

g_{\sup} (r, t) = g (r, t) - g_{eq} (r),

(12)

where the equilibrium radial distribution function, g_eq, 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).

View large Download slide

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

Finding a solution to Eq. for arbitrary values of

\dot{γ}

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.⁶ We will henceforth restrict our attention to the special case of steady-shear,

\dot{γ} (t) \to \dot{γ}

⁠, 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. reduces to

2 D_{0} \nabla \cdot (\nabla g_{\sup} (r) + g_{\sup} (r) \nabla β ϕ (r)) - \nabla \cdot (g (r) v (r)) = 0 .

(13)

In equilibrium, v(r) = 0, and the pair-current, ∇g_sup(r) + g_sup(r)∇βϕ(r), vanishes. This condition leads trivially to the solution g_sup(r) = 0, which implies that g(r) = g_eq(r), as expected.

To obtain a low shear-rate solution of the steady-state equation we assume g_sup to be a linear function of

\dot{γ}

⁠. (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 into and neglecting terms quadratic and higher in

\dot{γ}

then yields the following linearized steady-state condition:

2 D_{0} \nabla \cdot (\nabla g_{\sup} (r) + g_{\sup} (r) \nabla β ϕ (r)) - \nabla \cdot (g_{eq} (r) v (r)) = 0 .

(14)

In order to correctly capture the anisotropy induced by the flow, we make the following ansatz:¹⁸

g_{\sup} (r) = - \frac{1}{2 D_{0}} (\frac{r \cdot E \cdot r}{r^{2}}) e^{- β ϕ (r)} f (r),

(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 into Eq. , as described in Appendix B, finally yields the radial balance equations required to determine f(r). These are given by

\begin{gathered} \frac{d g_{eq} (r)}{d r} \overset{2 D}{=} - \frac{1}{r^{2}} \frac{d}{d r} (r \frac{d f (r)}{d r} e^{- β ϕ (r)}) + \frac{4}{r^{3}} f (r) e^{- β ϕ (r)}, \\ \frac{d g_{eq} (r)}{d r} \overset{3 D}{=} - \frac{1}{r^{3}} \frac{d}{d r} (r^{2} \frac{d f (r)}{d r} e^{- β ϕ (r)}) + \frac{6}{r^{3}} f (r) e^{- β ϕ (r)}, \end{gathered}

(16)

in two- and three-dimensions, respectively. The boundary conditions required to solve Eq. 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.

C. Low-shear viscosity

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:³

σ = - k_{B} T ρ_{b} I + \frac{1}{2} ρ_{b}^{2} \int d r \frac{r r}{r} (\frac{d ϕ (r)}{d r}) g (r),

(17)

where

I

is the unit tensor. Due to the isotropy of the equilibrium pair distribution function, only g_sup 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,

\dot{γ} \to 0

⁠, Eq. gives the low shear viscosity, η₀. Substituting Eq. into Eq. and using the appropriate form for the rate-of-strain tensor in shear flow, given by Eq. , then yields

η_{0} = \frac{k_{B} T}{4 D_{0}} ρ_{b}^{2} \int d r \frac{x^{2} y^{2}}{r^{3}} (\frac{d}{d r} e^{- β ϕ (r)}) f (r) .

(18)

Both the function f and the integral in depend on the dimensionality of the system. We then find

\begin{matrix} η_{0} \overset{2 D}{=} \frac{k_{B} T}{π D_{0}} Φ_{2 D}^{2} \int_{0}^{\infty} d r r^{2} (\frac{d}{d r} e^{- β ϕ (r)}) f (r), \\ η_{0} \overset{3 D}{=} \frac{k_{B} T}{π D_{0}} \frac{12}{5} Φ_{3 D}^{2} \int_{0}^{\infty} d r r^{3} (\frac{d}{d r} e^{- β ϕ (r)}) f (r), \end{matrix}

(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. would yield the exact low-shear viscosity. Within the present approach, the adiabatic approximation employed to close the two-body equation of motion yields approximate radial balance equations for f and, therefore, an approximate low-shear viscosity.

For the well-studied case of low density hard-spheres in three-dimensions, Eq. has the analytic solution f(r) = 1/(3r³), and Eq. recovers the quadratic term in the well-known low density expansion (see p. 114 in Ref. 2),

η_{0} = η_{s} (1 + \frac{5}{2} Φ_{3 D} + \frac{12}{5} Φ_{3 D}^{2}),

(20)

which applies in the absence of hydrodynamic interactions between the particles. The first term in is the solvent viscosity, η_s. The second term arises from the drag of the solvent on the surface of each individual sphere.⁵¹ 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. and then employing the Stokes-Einstein relation k_BT/D₀ = 3πη_s.

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, g_eq. 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,

\frac{g_{\sup} (x, y)}{\dot{γ} d^{2} / 2 D_{0}} = - \frac{x y}{r^{2}} f (r),

(21)

where

r = \sqrt{x^{2} + y^{2}}

and f(r) is obtained from the numerical integration of Eq. . 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/2r², and no packing oscillations are visible in the resulting scaled g_sup. 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 g_sup take both positive and negative values and are a nontrivial prediction of the radial balance equation . The packing structure of the full pair correlation function, g(r) = g_eq(r) + g_sup(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

f (r = d) = 1 / (2 d^{2}), \forall ρ_{b},

(22)

in each of the panels shown in Fig. 2.

FIG. 2.

View large Download slide

Superadiabatic-DDFT g_sup. For hard-disks in two-dimensions, we show the superadiabatic contribution to the pair correlation function, g_sup, in units of $\dot{γ} d^{2} / 2 D_{0}$ ⁠. Since we consider only low shear-rates, g_sup exhibits quadripolar symmetry [see Eq. (15)]. Increasing the bulk density leads to packing oscillations.

Knowledge of the nonequilibrium pair correlation function also enables the calculation of the interaction contribution to the low-shear viscosity, η₀, 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 η₀ 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.$

View large Download slide

Low-shear viscosity. A comparison of η₀ 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.

Due to the density independence of the contact value , the above-described route to obtaining the interaction contribution to the low-shear viscosity recovers its low-density form, namely,

η_{0} = \frac{k_{B} T d^{2}}{D_{0}} \frac{Φ_{2 D}^{2}}{2 π}, \forall ρ_{b} .

(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. . However, this mechanism is neglected per construction within the current approximation. The fact that we recover only the leading-order contribution 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 η₀ could be obtained from superadiabatic-DDFT by employing the “test-particle method,” as we will discuss in Sec. IV.

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,¹⁰ 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

Integrating the many-particle Smoluchowski equation over N − 2 particle coordinates yields the following, formally exact equation of motion for the two-body density^41,44

\begin{align} \frac{1}{D_{0}} \frac{\partial ρ^{(2)} (r_{1}, r_{2}, t)}{\partial t} & = \sum_{i = 1, 2} \nabla_{r_{i}} \cdot (\nabla_{r_{i}} ρ^{(2)} (r_{1}, r_{2}, t) \\ + ρ^{(2)} (r_{1}, r_{2}, t) \nabla_{r_{i}} β (V_{ext} (r_{i}) + ϕ (r_{12})) \\ + \int d r_{3} ρ^{(3)} (r_{1}, r_{2}, r_{3}, t) \nabla_{r_{i}} β ϕ (r_{i 3})) . \end{align}

(24)

This equation can also be rewritten in the following alternative form:

\begin{align} \frac{1}{D_{0}} \frac{\partial ρ^{(2)} (r_{1}, r_{2}, t)}{\partial t} \\ = \sum_{i = 1, 2} \nabla_{r_{i}} \cdot (ρ^{(2)} (r_{1}, r_{2}, t) (\nabla_{r_{i}} \ln (ρ^{(2)} (r_{1}, r_{2}, t)) \\ + \nabla_{r_{i}} β (V_{ext} (r_{i}) + ϕ (r_{12})) \\ + \int d r_{3} \frac{ρ^{(3)} (r_{1}, r_{2}, r_{3}, t)}{ρ^{(2)} (r_{1}, r_{2}, t)} \nabla_{r_{i}} β ϕ (r_{i 3}))), \end{align}

(25)

which remains fully equivalent to Eq. , since there is no approximation involved at this stage.

Using the standard form of the exact equation of motion and applying the adiabatic approximation solely to the three-body density, i.e.,

ρ^{(3)} \to ρ_{ad}^{(3)}

⁠, yields

\begin{align} \frac{1}{D_{0}} \frac{\partial ρ^{(2)} (r_{1}, r_{2}, t)}{\partial t} & = \sum_{i = 1, 2} \nabla_{r_{i}} \cdot (\nabla_{r_{i}} ρ^{(2)} (r_{1}, r_{2}, t) \\ + ρ^{(2)} (r_{1}, r_{2}, t) \nabla_{r_{i}} β (V_{ext} (r_{i}) + ϕ (r_{12})) \\ + \int d r_{3} ρ_{ad}^{(3)} (r_{1}, r_{2}, r_{3}, t) \nabla_{r_{i}} β ϕ (r_{i 3})), \end{align}

(26)

where the adiabatic three-body density is defined as

ρ_{ad}^{(3)} (r_{1}, r_{2}, r_{3}, t) \equiv ρ_{eq}^{(3)} (r_{1}, r_{2}, r_{3}; [ρ (r, t)]),

(27)

in analogy to Eq. . Substitution of the second-order Yvon–Born–Green (YBG2) equation^41,46

\begin{align} \int d r_{3} ρ_{ad}^{(3)} (r_{1}, r_{2}, r_{3}, t) \nabla_{r_{i}} β ϕ (r_{i 3}) \\ = - \nabla_{r_{i}} ρ_{ad}^{(2)} (r_{1}, r_{2}, t) - ρ_{ad}^{(2)} (r_{1}, r_{2}, t) \nabla_{r_{i}} β \\ \times (V_{ad} (r_{i}, t) + ϕ (r_{12})), \end{align}

(28)

into yields the second equation of superadiabatic-DDFT, namely Eq. . 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 and to apply the approximation to the integral term,

\begin{align} \frac{1}{D_{0}} \frac{\partial ρ^{(2)} (r_{1}, r_{2}, t)}{\partial t} \\ = \sum_{i = 1, 2} \nabla_{r_{i}} \cdot (ρ^{(2)} (r_{1}, r_{2}, t) (\nabla_{r_{i}} \ln (ρ^{(2)} (r_{1}, r_{2}, t)) \\ + \nabla_{r_{i}} β (V_{ext} (r_{i}) + ϕ (r_{12})) \\ + \int d r_{3} \frac{ρ_{ad}^{(3)} (r_{1}, r_{2}, r_{3}, t)}{ρ_{ad}^{(2)} (r_{1}, r_{2}, t)} \nabla_{r_{i}} β ϕ (r_{i 3}))) . \end{align}

(29)

We note that Eq. is no longer equivalent to the previous Eq. . The essential difference between these two options is that in Eq. we approximate the joint probability density, whereas in Eq. we approximate the conditional probability density. Substitution of the rewritten YBG2 equation

\begin{align} \int d r_{3} \frac{ρ_{ad}^{(3)} (r_{1}, r_{2}, r_{3}, t)}{ρ_{ad}^{(2)} (r_{1}, r_{2}, t)} \nabla_{r_{i}} β ϕ (r_{i 3}) \\ = - \nabla_{r_{i}} \ln (ρ_{ad}^{(2)} (r_{1}, r_{2}, t)) - \nabla_{r_{i}} β (V_{ad} (r_{i}, t) + ϕ (r_{12})), \end{align}

(30)

into expression , then yields the following alternative equation of motion for the two-body density:

\begin{align} \frac{1}{D_{0}} \frac{\partial ρ^{(2)} (r_{1}, r_{2}, t)}{\partial t} \\ = \sum_{i = 1, 2} \nabla_{r_{i}} \cdot (ρ^{(2)} (r_{1}, r_{2}, t) (\nabla_{r_{i}} \ln (ρ^{(2)} (r_{1}, r_{2}, t)) \\ + \nabla_{r_{i}} β (V_{ext} (r_{i}) - V_{ad} (r_{i}, t)) - \nabla_{r_{i}} \ln (ρ_{ad}^{(2)} (r_{1}, r_{2}, t)))) . \end{align}

(31)

We refer to this approximation strategy as “Russel–Gast-like,” since Eq. 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.

B. Low shear-rate solutions for the pair distribution function in bulk

Following the same scheme as in Subsection II B, the bulk limit of Eq. is given by

\begin{align} \frac{\partial g (r, t)}{\partial t} & = - 2 D_{0} \nabla \cdot (- \nabla g_{\sup} (r, t) + \frac{g_{\sup} (r, t)}{g_{eq} (r)} \nabla g_{eq} (r)) \\ - \nabla \cdot (g (r, t) v (r, t)) . \end{align}

(32)

This provides an alternative to Eq. . [As a side-note, we mention that the second term on the right-hand side of Eq. can be rewritten using the following rearrangement:

\begin{array}{l} \frac{g_{\sup} (r, t)}{g_{eq} (r)} \nabla g_{eq} (r) & = - g_{\sup} (r, t) \nabla (- \ln (g_{eq} (r))) \\ = - g_{\sup} (r, t) \nabla (β ϕ_{mf} (r)), \end{array}

where we have introduced the “potential of mean force,” ϕ_mf(r), defined by

g_{eq} (r) \equiv e^{- β ϕ_{mf} (r)}

⁠.] In the low density limit, g_eq(r) → e^−βϕ(r), which yields

\frac{g_{\sup} (r, t)}{g_{eq} (r)} \nabla g_{eq} (r) \to - g_{\sup} (r, t) \nabla β ϕ (r),

(33)

such that Eqs. and then become equivalent (which is not the case in general).

In the steady-state at finite densities, Eq. reduces to the condition

2 D_{0} \nabla \cdot (\nabla g_{\sup} (r, t) - \frac{g_{\sup} (r, t)}{g_{eq} (r)} \nabla g_{eq} (r)) - \nabla \cdot (g (r, t) v (r, t)) = 0 .

(34)

Using the definition and again assuming that g_sup is linear in the flow-rate yields the following expression to leading order:

2 D_{0} \nabla \cdot (\nabla g_{\sup} (r, t) - \frac{g_{\sup} (r, t)}{g_{eq} (r)} \nabla g_{eq} (r)) - \nabla \cdot (g_{eq} (r) v (r)) = 0 .

(35)

Since Eqs. and are non-equivalent, the solution of Eq. now requires a different ansatz than used previously. The appropriate choice in the present case is given by

g_{\sup} (r) = - \frac{1}{2 D_{0}} (\frac{r \cdot E \cdot r}{r^{2}}) g_{eq} (r) f^{⋆} (r),

(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 into the linearized steady-state equation then yields the alternative radial balance equations

\begin{gathered} \frac{d g_{eq} (r)}{d r} \overset{2 D}{=} - \frac{1}{r^{2}} \frac{d}{d r} (r \frac{d f^{⋆} (r)}{d r} g_{eq} (r)) + \frac{4}{r^{3}} f^{⋆} (r) g_{eq} (r), \\ \frac{d g_{eq} (r)}{d r} \overset{3 D}{=} - \frac{1}{r^{3}} \frac{d}{d r} (r^{2} \frac{d f^{⋆} (r)}{d r} g_{eq} (r)) + \frac{6}{r^{3}} f^{⋆} (r) g_{eq} (r), \end{gathered}

(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. and the boundary conditions required to solve them.

C. Low-shear viscosity

Within the alternative adiabatic closure, the low-shear viscosity is given by

η_{0} = \frac{k_{B} T}{4 D_{0}} ρ_{b}^{2} \int d r \frac{x^{2} y^{2}}{r^{3}} (\frac{d}{d r} e^{- β ϕ (r)}) y_{eq} (r) f^{⋆} (r),

(38)

where we have introduced the so-called “cavity distribution function,”

y_{eq} (r) = g_{eq} (r) e^{β ϕ (r)},

(39)

familiar from liquid-state theory.⁴⁶ We emphasize that f^⋆(r) appearing in Eq. is given by solution of the radial balance equations and is distinct from the function f(r) obtained from solution of Eq. . In two- and three-dimensions, Eq. becomes

\begin{gathered} η_{0} \overset{2 D}{=} \frac{k_{B} T}{π D_{0}} Φ_{2 D}^{2} \int_{0}^{\infty} d r r^{2} (\frac{d}{d r} e^{- β ϕ (r)}) y_{eq} (r) f^{⋆} (r), \\ η_{0} \overset{3 D}{=} \frac{k_{B} T}{π D_{0}} \frac{12}{5} Φ_{3 D}^{2} \int_{0}^{\infty} d r r^{3} (\frac{d}{d r} e^{- β ϕ (r)}) y_{eq} (r) f^{⋆} (r), \end{gathered}

(40)

respectively. For the special case of low density hard-spheres in three-dimensions, the analytic solution of Eq. is the same as that given by Eq. , namely

f^{⋆} (r) = f (r) = 1 / (3 r^{3}),

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

D. Numerical results

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 g_eq for input to the radial balance equation . In Fig. 4, we show scatter plots of the quantity

\frac{g_{\sup} (x, y)}{\dot{γ} d^{2} / 2 D_{0}} = - \frac{x y}{r^{2}} f^{⋆} (r) g_{eq} (r),

(42)

for four different values of the input bulk density. Note that the radial function f^⋆(r) is obtained from numerical integration of Eq. and that

r = \sqrt{x^{2} + y^{2}}

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

View large Download slide

Russel–Gast-type closure for g_sup. 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, g_sup, in units of $\dot{γ} d^{2} / 2 D_{0}$ [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.

The results for g_sup at the lowest bulk density considered, ρ_b = 0.1, are very close to the known low-density limit, for which f^⋆(r) = 1/2r², 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)g_eq(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.

View large Download slide

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 $\dot{γ} d^{2} / 2 D_{0} = 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).

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, η₀ 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.¹²

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

View large Download slide

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 η₀ than the superadiabatic-DDFT but exhibits an unphysical curvature as the area fraction is increased.

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.¹⁰

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 $\dot{γ}$ ⁠. For larger shear-rates, we would need to consider the extension of our solution ansatz (15) or (36), respectively, to higher-order in $\dot{γ}$ ⁠, 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 simulation³⁷ and experiment⁵⁶ 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,⁵⁸ 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.¹² For the present Brownian dynamics, without hydrodynamic interactions, this consists of multiplying the low density limiting solution for the interaction contribution to η₀ [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,¹⁶ 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 η₀. In contrast, a rescaled Russel–Gast-type theory would continue to exhibit an unphysical curvature, as the contact value of g_eq 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 V_ext(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.⁵⁴ 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

In a system undergoing translationally invariant, homogeneous flow, the velocity field can be conveniently expressed in the following form:

v (r) = κ \cdot r,

(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

κ = (\begin{matrix} 0 & \dot{γ} & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix}),

(A2)

for a three-dimensional system with shear-rate

\dot{γ}

⁠.

The right-hand side of can be decomposed into a sum of two terms,

v (r) = E \cdot r + Ω \cdot r,

(A3)

where

E = \frac{1}{2} (κ + κ^{T}), Ω = \frac{1}{2} (κ - κ^{T}),

(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.⁶⁵ For the aforementioned case of shear flow, we obtain

E = \frac{1}{2} (\begin{matrix} 0 & \dot{γ} & 0 \\ \dot{γ} & 0 & 0 \\ 0 & 0 & 0 \end{matrix}), Ω = \frac{1}{2} (\begin{matrix} 0 & \dot{γ} & 0 \\ - \dot{γ} & 0 & 0 \\ 0 & 0 & 0 \end{matrix}) .

(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 g_sup in Eqs. and .

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.

We start by collecting a few useful identities. For a general, spatially varying, second-rank tensor M, the following holds:

\begin{align} \nabla ((\hat{r} \cdot M \cdot \hat{r}) a (r)) & = \frac{2}{r^{2}} (r \cdot M) \cdot (I - \hat{r} \hat{r}) a (r) \\ + \frac{1}{r} \frac{d a (r)}{d r} (r \cdot M) \cdot \hat{r} \hat{r}, \end{align}

(B1)

where

\hat{r} = r / r

is a unit vector,

\hat{r} \hat{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

\nabla \cdot (u \cdot M) = u \cdot (\nabla \cdot M) + M : \nabla u,

(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

\begin{aligned} \nabla \cdot (a (r) M) & = \nabla a (r) \cdot M + a (r) \nabla \cdot M \\ = (\frac{d a (r)}{d r}) \hat{r} \cdot M + a (r) \nabla \cdot M . \end{aligned}

(B3)

This identity generates the following special cases in two-dimensions:

\nabla \cdot (a (r) (I - \hat{r} \hat{r})) \overset{2 D}{=} - \frac{a (r)}{r} \hat{r},

(B4)

\nabla \cdot (a (r) \hat{r} \hat{r}) \overset{2 D}{=} \frac{1}{r} \frac{d}{d r} (r a (r)) \hat{r},

(B5)

and three-dimensions

\nabla \cdot (a (r) (I - \hat{r} \hat{r})) \overset{3 D}{=} - \frac{2 a (r)}{r} \hat{r},

(B6)

\nabla \cdot (a (r) \hat{r} \hat{r}) \overset{3 D}{=} \frac{1}{r^{2}} \frac{d}{d r} (r^{2} a (r)) \hat{r},

(B7)

respectively. In the special case where M is symmetric, it can be shown that

(\nabla (M \cdot \hat{r})) \cdot \hat{r} = 0 .

(B8)

Finally, the divergence of the scalar product

M \cdot \hat{r}

is given by the aesthetically appealing identity,

\nabla \cdot (M \cdot \hat{r}) = \frac{Tr (M)}{r} - \frac{\hat{r} \cdot M \cdot \hat{r}}{r} .

(B9)

We will now employ these identities to calculate the pair distribution function at low flow rates.

Substitution of the ansatz into the linearized steady-state equation yields

\nabla \cdot (g_{eq} (r) v (r)) = - \nabla \cdot (e^{- β ϕ (r)} \nabla ((\hat{r} \cdot E \cdot \hat{r}) f (r))),

(B10)

in which the only unknown quantity is the function f(r). Assuming incompressible flow, ∇ · v = 0, using identity and rewriting the projected velocity as

v (r) \cdot \hat{r} = r (\hat{r} \cdot E \cdot \hat{r}),

(B11)

enables us to re-express Eq. as follows:

\begin{align} r \frac{d g_{eq} (r)}{d r} (\hat{r} \cdot E \cdot \hat{r}) \\ = - \nabla \cdot ((\hat{r} \cdot E) \cdot ((I - \hat{r} \hat{r}) \frac{2 f (r)}{r} + \hat{r} \hat{r} \frac{d f (r)}{d r}) e^{- β ϕ (r)}) . \end{align}

(B12)

The advantage of this representation is the appearance of the dyadic tensors,

\hat{r} \hat{r}

and

(I - \hat{r} \hat{r})

⁠, which project either along or perpendicular to the relative position vector r.

On the right-hand side of , we have to calculate the divergence of the scalar product between the vector

\hat{r} \cdot E

and a second-rank tensor. We can thus exploit relation to obtain

\begin{align} r \frac{d g_{eq} (r)}{d r} (\hat{r} \cdot E \cdot \hat{r}) \\ = - (\hat{r} \cdot E) \cdot (\nabla \cdot (((I - \hat{r} \hat{r}) \frac{2 f (r)}{r} + \hat{r} \hat{r} \frac{d f (r)}{d r}) e^{- β ϕ (r)})) \\ - (((I - \hat{r} \hat{r}) \frac{2 f (r)}{r} + \hat{r} \hat{r} \frac{d f (r)}{d r}) e^{- β ϕ (r)}) : \nabla (E \cdot \hat{r}) . \end{align}

(B13)

We will consider separately the two terms appearing on the right-hand side of , which we henceforth refer to as (i) and (ii).

To simplify (i) in two-dimensions, we use –. This yields

- (\hat{r} \cdot E \cdot \hat{r}) (\frac{1}{r} \frac{d}{d r} (r \frac{d f (r)}{d r} e^{- β ϕ (r)}) - \frac{2 f (r) e^{- β ϕ (r)}}{r^{2}}) .

(B14)

The analogous result in three-dimensions is given by

- (\hat{r} \cdot E \cdot \hat{r}) (\frac{1}{r^{2}} \frac{d}{d r} (r^{2} \frac{d f (r)}{d r} e^{- β ϕ (r)}) - \frac{4 f (r) e^{- β ϕ (r)}}{r^{2}}),

(B15)

where we have used equations , , and .

The simplification of term (ii) does not depend on the dimensionality of the system. We first employ Eq. to re-express the factor

\nabla (E \cdot \hat{r})

as a divergence and then use identity to obtain

\begin{align} - \frac{2 f (r) e^{- β ϕ (r)}}{r} \nabla \cdot (E \cdot \hat{r}) \\ = - \frac{2 f (r) e^{- β ϕ (r)}}{r} (\frac{Tr (E)}{r} - \frac{\hat{r} \cdot E \cdot \hat{r}}{r}) \\ = \frac{2 f (r) e^{- β ϕ (r)}}{r^{2}} (\hat{r} \cdot E \cdot \hat{r}), \end{align}

(B16)

since

Tr (E) = 0

for incompressible flow.

Finally, substituting and into yields the desired two-dimensional result

\begin{align} (\hat{r} \cdot E \cdot \hat{r}) (\frac{1}{r} \frac{d}{d r} (r \frac{d f (r)}{d r} e^{- β ϕ (r)}) - \frac{4}{r^{2}} f (r) e^{- β ϕ (r)}) \\ \overset{2 D}{=} - (\hat{r} \cdot E \cdot \hat{r}) r \frac{d g_{eq} (r)}{d r}, \end{align}

(B17)

while substitution of and into gives the result in three-dimensions,

\begin{align} (\hat{r} \cdot E \cdot \hat{r}) (\frac{1}{r^{2}} \frac{d}{d r} (r^{2} \frac{d f (r)}{d r} e^{- β ϕ (r)}) - \frac{6}{r^{2}} f (r) e^{- β ϕ (r)}) \\ \overset{3 D}{=} - (\hat{r} \cdot E \cdot \hat{r}) r \frac{d g_{eq} (r)}{d r} . \end{align}

(B18)

Since these expressions remain valid for all choices of the traceless tensor E, the coefficients of the quadratic form

(\hat{r} \cdot E \cdot \hat{r})

must be equal. We thus obtain the radial balance equations stated in the main text. We note that, in addition to the explicit appearance of the pair interaction potential in , respectively , both the pair potential and the bulk density enter these equations implicitly via the equilibrium radial distribution function.

APPENDIX C: BOUNDARY CONDITIONS

The radial balance equations are second-order in spatial derivative, and their solution thus requires the specification of two boundary conditions. From the linearized steady-state equation we can identify the following expression for the pair-current at low shear-rates:

j (r) = - 2 D_{0} (\nabla g_{\sup} (r) + g_{\sup} (r) \nabla β ϕ (r)) + g_{eq} (r) v (r) .

(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

j (r) \cdot \hat{r} = 0,

(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

j (r) \cdot \hat{r} \underset{r \to 1}{=} 0,

(C3)

where we have set the particle diameter equal to unity. Using steps directly analogous to those leading from to enables the pair-current of hard-spheres to be expressed in the following form:

j (r) = (\hat{r} \cdot E) \cdot ((I - \hat{r} \hat{r}) \frac{2 f (r)}{r} + \hat{r} \hat{r} \frac{d f (r)}{d r}) + g_{eq} (r) v (r),

(C4)

for r > 1. [Note that ϕ(r > 1) = 0 for the hard-sphere potential.] Taking the scalar product of Eq. with the radial unit vector yields

j (r) \cdot \hat{r} = (\hat{r} \cdot E \cdot \hat{r}) (r g_{eq} (r) + \frac{d f (r)}{d r}),

(C5)

where we have used . The hard-sphere “zero-flux” condition at particle contact thus becomes

{\frac{d f (r)}{d r}|}_{r = 1} = - g_{eq} (1) .

(C6)

The radial balance equations , combined with and f(r → ∞) = 0, fully determine the function f(r), given that the input equilibrium radial distribution function is known.

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.

Substitution of the ansatz into the linearized steady-state equation yields

\nabla \cdot (g_{eq} (r) v (r)) = - \nabla \cdot (g_{eq} (r) \nabla ((\hat{r} \cdot E \cdot \hat{r}) f^{⋆} (r))),

(D1)

from which we can determine the function f^⋆(r). Assuming incompressible flow, ∇ · v = 0, and using the identities and , we can re-express Eq. in the following form:

\begin{align} r \frac{d g_{eq} (r)}{d r} (\hat{r} \cdot E \cdot \hat{r}) \\ = - \nabla \cdot ((\hat{r} \cdot E) \cdot ((I - \hat{r} \hat{r}) \frac{2 f^{⋆} (r)}{r} + \hat{r} \hat{r} \frac{d f^{⋆} (r)}{d r}) g_{eq} (r)) . \end{align}

(D2)

Steps analogous to those leading from Eqs. to then generate the alternative form of the radial balance equations .

From the linearized steady-state equation , we identify the pair-current at low shear-rates,

j (r, t) = - 2 D_{0} (\nabla g_{\sup} (r, t) - \frac{g_{\sup} (r, t)}{g_{eq} (r)} \nabla g_{eq} (r)) + g_{eq} (r) v (r) .

(D3)

A calculation analogous to that leading from to shows that Eq. can be re-expressed as

j (r) = (\hat{r} \cdot E) \cdot ((I - \hat{r} \hat{r}) \frac{2 f^{⋆} (r)}{r} + \hat{r} \hat{r} \frac{d f^{⋆} (r)}{d r}) g_{eq} (r) + g_{eq} (r) v (r) .

(D4)

For a general interaction potential, the boundary condition 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 at r = 1 yields

{\frac{d f^{⋆} (r)}{d r}|}_{r = 1} = - 1,

(D5)

which is different from Eq. . The boundary condition at particle contact, Eq. , together with f^⋆(r → ∞) = 0, then fully determines the solutions of the radial balance equations .

REFERENCES

1.

M. E.

Cates

and

M. R.

Evans

,

Soft and Fragile Matter

(

CRC Press

,

2000

).

Google Scholar

Crossref

2.

J.

Mewis

and

N. J.

Wagner

,

Colloidal Suspension Rheology

(

Cambridge University Press

,

2012

).

Google Scholar

3.

J. M.

Brader

, “

Nonlinear rheology of colloidal dispersions

,”

J. Phys.: Condens. Matter

22

,

363101

(

2010

).

https://doi.org/10.1088/0953-8984/22/36/363101

Google Scholar

Crossref

4.

R.

Larson

,

Structure and Rheology of Complex Fluids

,

Topics in Chemical Engineering

(

Oxford University Press

,

New York

,

1999

).

Google Scholar

5.

G. K.

Batchelor

,

An Introduction to Fluid Dynamics

,

Cambridge Mathematical Library

(

Cambridge University Press

,

2000

).

Google Scholar

Crossref

6.

J.

Bergenholtz

,

J. F.

Brady

, and

M.

Vicic

, “

The non-Newtonian rheology of dilute colloidal suspensions

,”

J. Fluid Mech.

456

,

239

–

275

(

2002

).

https://doi.org/10.1017/s0022112001007583

Google Scholar

Crossref

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

).

https://doi.org/10.1017/s0022112097006320

Google Scholar

Crossref

8.

T.

Phtsuki

, “

Dynamical properties of strongly interacting Brownian particles: 1. Dynamic shear viscosity

,”

Physica A

108

,

441

(

1981

).

https://doi.org/10.1016/0378-4371(81)90141-2

Google Scholar

Crossref

9.

D.

Ronis

, “

Theory of fluctuations in colloidal suspensions undergoing steady shear flow

,”

Phys. Rev. A

29

,

1453

(

1984

).

https://doi.org/10.1103/physreva.29.1453

Google Scholar

Crossref

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

).

https://doi.org/10.1063/1.450428

Google Scholar

Crossref

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

).

https://doi.org/10.1016/0378-4371(89)90003-4

Google Scholar

Crossref

12.

J. F.

Brady

, “

The rheological behavior of concentrated colloidal dispersions

,”

J. Chem. Phys.

99

,

567

(

1993

).

https://doi.org/10.1063/1.465782

Google Scholar

Crossref

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

).

https://doi.org/10.1122/1.550873

Google Scholar

Crossref

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

).

https://doi.org/10.1063/1.473029

Google Scholar

Crossref

15.

G.

Szamel

, “

Nonequilibrium structure and rheology of concentrated colloidal suspensions: Linear response

,”

J. Chem. Phys.

114

,

8708

(

2001

).

https://doi.org/10.1063/1.1366335

Google Scholar

Crossref

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

.

Google Scholar

Crossref

17.

J.

Morris

, “

A review of microstructure in concentrated suspensions and its implications for rheology and bulk flow

,”

Rheol. Acta

48

,

909

(

2009

).

https://doi.org/10.1007/s00397-009-0352-1

Google Scholar

Crossref

18.

W. B.

Russel

,

D. A.

Saville

, and

W. R.

Schowalter

,

Colloidal Dispersions

,

Cambridge Monographs on Mechanics

(

Cambridge University Press

,

1989

).

Google Scholar

Crossref

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

).

https://doi.org/10.1103/physreve.106.044610

Google Scholar

Crossref

PubMed

20.

M.

Fuchs

and

M. E.

Cates

, “

Theory of nonlinear rheology and yielding of dense colloidal suspensions

,”

Phys. Rev. Lett.

89

,

248304

(

2002

).

https://doi.org/10.1103/physrevlett.89.248304

Google Scholar

Crossref

PubMed

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

).

https://doi.org/10.1103/physrevlett.98.058301

Google Scholar

Crossref

PubMed

22.

J. M.

Brader

,

M. E.

Cates

, and

M.

Fuchs

, “

First-principles constitutive equation for suspension rheology

,”

Phys. Rev. Lett.

101

,

138301

(

2008

).

https://doi.org/10.1103/physrevlett.101.138301

Google Scholar

Crossref

PubMed

23.

M.

Fuchs

and

M. E.

Cates

, “

A mode coupling theory for Brownian particles in homogeneous steady shear flow

,”

J. Rheol.

53

,

957

(

2009

).

https://doi.org/10.1122/1.3119084

Google Scholar

Crossref

24.

J. M.

Brader

,

M. E.

Cates

, and

M.

Fuchs

, “

First-principles constitutive equation for suspension rheology

,”

Phys. Rev. E

86

,

021403

(

2012

).

https://doi.org/10.1103/physreve.86.021403

Google Scholar

Crossref

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

).

https://doi.org/10.1017/s0022112003006001

Google Scholar

Crossref

26.

J.

Brader

and

M.

Krüger

, “

Density profiles of a colloidal liquid at a wall under shear flow

,”

Mol. Phys.

109

,

1029

(

2011

).

https://doi.org/10.1080/00268976.2010.541889

Google Scholar

Crossref

27.

M.

Krüger

and

J. M.

Brader

, “

Controlling colloidal sedimentation using time-dependent shear

,”

Europhys. Lett.

96

,

68006

(

2011

).

https://doi.org/10.1209/0295-5075/96/68006

Google Scholar

Crossref

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

).

https://doi.org/10.1063/1.4866450

Google Scholar

Crossref

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

).

https://doi.org/10.1088/0953-8984/28/24/244023

Google Scholar

Crossref

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

).

https://doi.org/10.1039/c4sm01329h

Google Scholar

Crossref

PubMed

31.

J. C.

Crocker

and

D. G.

Grier

, “

Methods of digital video microscopy for colloidal studies

,”

J. Colloid Interface Sci.

179

,

298

(

1996

).

https://doi.org/10.1006/jcis.1996.0217

Google Scholar

Crossref

32.

V.

Prasad

,

D.

Semwogerere

, and

E. R.

Weeks

, “

Confocal microscopy of colloids

,”

J. Phys.: Condens. Matter

19

,

113102

(

2007

).

https://doi.org/10.1088/0953-8984/19/11/113102

Google Scholar

Crossref

33.

R.

Besseling

,

L.

Isa

,

E. R.

Weeks

, and

W. C.

Poon

, “

Quantitative imaging of colloidal flows

,”

Adv. Colloid Interface Sci.

146

,

1

(

2009

).

https://doi.org/10.1016/j.cis.2008.09.008

Google Scholar

Crossref

PubMed

34.

D. R.

Foss

and

J. F.

Brady

, “

Brownian dynamics simulation of hard-sphere colloidal dispersions

,”

J. Rheol.

44

,

629

(

2000

).

https://doi.org/10.1122/1.551104

Google Scholar

Crossref

35.

A.

Sierou

and

J. F.

Brady

, “

Accelerated Stokesian dynamics simulations

,”

J. Fluid Mech.

448

,

115

–

146

(

2001

).

https://doi.org/10.1017/s0022112001005912

Google Scholar

Crossref

36.

A. J.

Banchio

and

J. F.

Brady

, “

Accelerated Stokesian dynamics: Brownian motion

,”

J. Chem. Phys.

118

,

10323

(

2003

).

https://doi.org/10.1063/1.1571819

Google Scholar

Crossref

37.

Y.

Yurkovetsky

and

J. F.

Morris

, “

Triplet correlation in sheared suspensions of Brownian particles

,”

J. Chem. Phys.

124

,

204908

(

2006

).

https://doi.org/10.1063/1.2200689

Google Scholar

Crossref

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

).

https://doi.org/10.1103/physreve.102.042140

Google Scholar

Crossref

PubMed

39.

P.

Attard

, “

Spherically inhomogeneous fluids. I. Percus–Yevick hard spheres: Osmotic coefficients and triplet correlations

,”

J. Chem. Phys.

91

,

3072

(

1989

).

https://doi.org/10.1063/1.456930

Google Scholar

Crossref

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

).

https://doi.org/10.1103/physrevlett.108.037802

Google Scholar

Crossref

PubMed

41.

S. M.

Tschopp

and

J. M.

Brader

, “

First-principles superadiabatic theory for the dynamics of inhomogeneous fluids

,”

J. Chem. Phys.

157

,

234108

(

2022

).

https://doi.org/10.1063/5.0131441

Google Scholar

Crossref

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

).

https://doi.org/10.1063/5.0155856

Google Scholar

Crossref

43.

U. M. B.

Marconi

and

P.

Tarazona

, “

Dynamic density functional theory of fluids

,”

J. Chem. Phys.

110

,

8032

(

1999

).

https://doi.org/10.1063/1.478705

Google Scholar

Crossref

44.

A. J.

Archer

and

R.

Evans

, “

Dynamical density functional theory and its application to spinodal decomposition

,”

J. Chem. Phys.

121

,

4246

(

2004

).

https://doi.org/10.1063/1.1778374

Google Scholar

Crossref

45.

J.

Chakrabarti

,

J.

Dzubiella

, and

H.

Löwen

, “

Dynamical instability in driven colloids

,”

Europhys. Lett.

61

,

415

(

2003

).

https://doi.org/10.1209/epl/i2003-00193-6

Google Scholar

Crossref

46.

J.

Hansen

and

I.

McDonald

,

Theory of Simple Liquids

(

Elsevier Science

,

2006

).

Google Scholar

47.

D.

McQuarrie

,

Statistical Mechanics

(

University Science Books

,

2000

).

Google Scholar

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

).

https://doi.org/10.1017/s0022112089001813

Google Scholar

Crossref

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

).

https://doi.org/10.1103/physreve.105.054606

Google Scholar

Crossref

PubMed

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

).

https://doi.org/10.1103/physreve.99.052606

Google Scholar

Crossref

PubMed

51.

J. F.

Brady

, “

The Einstein viscosity correction in n dimensions

,”

Int. J. Multiphase Flow

10

,

113

(

1983

).

https://doi.org/10.1016/0301-9322(83)90064-2

Google Scholar

Crossref

52.

F.

Lado

, “

Numerical Fourier transforms in one, two, and three dimensions for liquid state calculations

,”

J. Comput. Phys.

8

,

417

(

1971

).

https://doi.org/10.1016/0021-9991(71)90021-0

Google Scholar

Crossref

53.

F.

Lado

, “

Equation of state of the hard-disk fluid from approximate integral equations

,”

J. Chem. Phys.

49

,

3092

(

1968

).

https://doi.org/10.1063/1.1670553

Google Scholar

Crossref

54.

J.

Reinhardt

,

F.

Weysser

, and

J. M.

Brader

, “

Density functional approach to nonlinear rheology

,”

Europhys. Lett.

102

,

28011

(

2013

).

https://doi.org/10.1209/0295-5075/102/28011

Google Scholar

Crossref

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

).

https://doi.org/10.1063/1.4968565

Google Scholar

Crossref

56.

J.

Vermant

and

M. J.

Solomon

, “

Flow-induced structure in colloidal suspensions

,”

J. Phys.: Condens. Matter

17

,

R187

(

2005

).

https://doi.org/10.1088/0953-8984/17/4/r02

Google Scholar

Crossref

57.

J.

Bergenholtz

, “

Theory of rheology of colloidal dispersions

,”

Curr. Opin. Colloid Interface Sci.

6

,

484

(

2001

).

https://doi.org/10.1016/s1359-0294(01)00112-1

Google Scholar

Crossref

58.

J.

Bławzdziewicz

and

G.

Szamel

, “

Structure and rheology of semidilute suspension under shear

,”

Phys. Rev. E

48

,

4632

(

1993

).

https://doi.org/10.1103/physreve.48.4632

Google Scholar

Crossref

59.

E.

Nazockdast

and

J. F.

Morris

, “

Microstructural theory and the rheology of concentrated colloidal suspensions

,”

J. Fluid Mech.

713

,

420

–

452

(

2012

).

https://doi.org/10.1017/jfm.2012.467

Google Scholar

Crossref

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

).

https://doi.org/10.1103/physrevlett.101.148302

Google Scholar

Crossref

PubMed

61.

M.

Rex

and

H.

Löwen

, “

Dynamical density functional theory for colloidal dispersions including hydrodynamic interactions

,”

Eur. Phys. J. E

28

,

139

(

2009

).

https://doi.org/10.1140/epje/i2008-10363-x

Google Scholar

Crossref

62.

M.

Rauscher

, “

DDFT for Brownian particles and hydrodynamics

,”

J. Phys.: Condens. Matter

22

,

364109

(

2010

).

https://doi.org/10.1088/0953-8984/22/36/364109

Google Scholar

Crossref

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

).

https://doi.org/10.1088/0953-8984/25/3/035101

Google Scholar

Crossref

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

).

https://doi.org/10.1063/5.0208943

Google Scholar

65.

J.

Dhont

,

An Introduction to Dynamics of Colloids

(

Elsevier Science

,

1996

).

Google Scholar

Superadiabatic dynamical density functional theory for colloidal suspensions under homogeneous steady-shear

I. INTRODUCTION

II. SUPERADIABATIC-DDFT

A. General framework

B. Low shear-rate solutions for the pair distribution function in bulk

C. Low-shear viscosity

D. Numerical results

III. RUSSEL–GAST-TYPE THEORY

A. Alternative adiabatic closure

B. Low shear-rate solutions for the pair distribution function in bulk

C. Low-shear viscosity

D. Numerical results

IV. DISCUSSION

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: RHEOLOGICAL QUANTITIES

APPENDIX B: THE RADIAL BALANCE EQUATION

APPENDIX C: BOUNDARY CONDITIONS

APPENDIX D: ALTERNATIVE FORMS OF THE RADIAL BALANCE EQUATION AND BOUNDARY CONDITIONS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Superadiabatic dynamical density functional theory for colloidal suspensions under homogeneous steady-shear

I. INTRODUCTION

II. SUPERADIABATIC-DDFT

A. General framework

B. Low shear-rate solutions for the pair distribution function in bulk

C. Low-shear viscosity

D. Numerical results

III. RUSSEL–GAST-TYPE THEORY

A. Alternative adiabatic closure

B. Low shear-rate solutions for the pair distribution function in bulk

C. Low-shear viscosity

D. Numerical results

IV. DISCUSSION

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: RHEOLOGICAL QUANTITIES

APPENDIX B: THE RADIAL BALANCE EQUATION

APPENDIX C: BOUNDARY CONDITIONS

APPENDIX D: ALTERNATIVE FORMS OF THE RADIAL BALANCE EQUATION AND BOUNDARY CONDITIONS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only