Analytical and Computational Studies of Correlations of Hydrodynamic Fluctuations in Shear Flow

TL;DR

This study compares analytical and DPD simulation methods to analyze hydrodynamic fluctuation correlations in shear flow, revealing good agreement at moderate shear rates and discrepancies at high shear rates, highlighting complex variable coupling.

Contribution

It introduces a combined analytical and simulation approach to study hydrodynamic fluctuations in shear flow, identifying limitations of analytical methods at high shear rates.

Findings

Analytical and DPD results agree well at moderate shear rates.

Discrepancies emerge at high shear rates due to complex variable coupling.

Results can serve as benchmarks for multiscale simulation algorithms.

Abstract

We study correlations of hydrodynamic fluctuations in shear flow analytically and also by dissipative particle dynamics~(DPD) simulations. The hydrodynamic equations are linearized around the macroscopic velocity field and then solved by a perturbation method in Fourier-transformed space. The autocorrelation functions~(ACFs) from the analytical method are compared with results obtained from DPD simulations under the same shear-flow conditions. Upto a moderate shear rate, various ACFs from the two approaches agree with each other well. At large shear rates, discrepancies between the two methods are observed, hence revealing strong additional coupling between different fluctuating variables, which is not considered in the analytical approach. In addition, the results at low and moderate shear rates can serve as benchmarks for developing multiscale algorithms for coupling of heterogeneous solvers, such as a hybrid simulation of molecular dynamics and fluctuating hydrodynamics solver, where thermal fluctuations are indispensable.

Equations

$$\left(\frac{\partial }{\partial t}+\mathbf{v}\cdot \nabla \right)\rho$$

$$\rho \left(\frac{\partial }{\partial t}+\mathbf{v}\cdot \nabla \right)\mathbf{v}$$

$$\frac{d\rho }{dt}$$

$$\rho \frac{d\mathbf{v}}{dt}$$

$${\Pi }_{\mu \sigma }$$

$${\Pi }_{\mu \sigma }^C$$

$${\Pi }_{\mu \sigma }^D$$

$$\rho \left(\frac{\partial }{\partial t}+\mathbf{v}\cdot \nabla \right)\mathbf{v}=-\nabla p+\eta {\nabla }^2\mathbf{v}+\left(\frac{\eta }{3}+\zeta \right)\nabla \left(\nabla \cdot \mathbf{v}\right).$$

$$<{\Pi }_{\mu \sigma }^R>$$

$$<{\Pi }_{\mu \sigma }^R(\mathbf{x},t){\Pi }_{ϵ\iota }^R({\mathbf{x}}^{\prime },t^{\prime })>$$

$${\Delta }_{\mu \sigma ϵ\iota }$$

$$v_{0\mu }(\mathbf{x},t)$$

$${\dot{\gamma }}_{\mu \sigma }$$

$$C_L(\mathbf{k},\tau )$$

$$C_{T_1}(\mathbf{k},\tau )$$

$$C_{T_2}(\mathbf{k},\tau )$$

$$\alpha (\mathbf{k},t)$$

$$\beta (\mathbf{k},t)$$

$${\mathbf{r}}_{ij}$$

$${\mathbf{v}}_{ij}$$

$${\mathbf{e}}_{ij}$$

$${\mathbf{F}}_{ij}^C$$

$${\mathbf{F}}_{ij}^D$$

$${\mathbf{F}}_{ij}^R$$

$$<{\theta }_{ij}(t)>$$

$$<{\theta }_{ij}(t){\theta }_{kl}(t^{\prime })>$$

$$W^D(r_{ij})$$

$$2k_{B}T\gamma$$

$$W^{C,R}(r_{ij})=\left\{\begin{array}[]{l r}(1-r_{ij}/r_{c})^{k},&r_{ij}<r_{c},\\ 0,&r_{ij}\geq r_{c}.\end{array}\right.$$

$$\delta v_{\mu }(\mathbf{x},t)=v_{\mu }(\mathbf{x},t)-\dot{\gamma }{\delta }_{\mu x}{\delta }_{\sigma y}{x}_{\sigma }.$$

$$<{\widehat{\delta u}}_{\sigma }(-\mathbf{k},t){\widehat{\delta u}}_{\sigma }(\mathbf{k},t+\tau )>=\frac{1}{N_s}\sum _{s=1}^{N_s}<{\widehat{\delta u}}_{\sigma }(-\mathbf{k},t)<{\widehat{\delta u}}_{\sigma }(\mathbf{k},t+\tau ),$$

$$\widetilde{\delta \mathbf{v}}(\mathbf{k},t)$$

$${\widehat{\delta u}}_{\sigma }(\mathbf{k},t)$$

$$\alpha (\mathbf{k},\tau )=k_0^2\tau +\frac{1}{3}{\dot{\gamma }}^2{k}_x^2{\tau }^3.$$

$$\alpha (\mathbf{k},\tau )=k_0^2\tau -\dot{\gamma }{k}_x{k}_y{\tau }^2+\frac{1}{3}{\dot{\gamma }}^2{k}_x^2{\tau }^3.$$

$$\beta (\mathbf{k},\tau )=\frac{1}{2\dot{\gamma }{k}_x}\left\{k_y{k}_0-k_y(\tau )k(\tau )-k_{\perp }^{2}In\left[\frac{k_y(\tau )+k(\tau )}{k_y+k_0}\right]\right\}.$$

$$\mathbf{z}(\mathbf{x},t)=[\delta \rho (\mathbf{x},t),\delta \mathbf{v}(\mathbf{x},t)],$$

$$\delta \rho$$

Full text

1 Introduction

The complex behavior of many particles plays a significant role in atomic fluids [1, 2], chemical and biological processes [3, 4], granular materials [5], and astrophysics [6]. On the one hand, given interparticle potentials, a kinetic-theory type of description from first principles can be formulated, which, however, may be too complex to apply. On the other hand, if physical quantities vary rather slowly in space and time, a local thermodynamic equilibrium may be valid. Therefore, the system can be represented by continuous hydrodynamic fields described by the Navier-Stokes-Fourier (NSF) equations, which take into account the conservation of mass, momentum and energy. Although such partial differential equations (PDEs) are concise and practically powerful, the thermodynamic derivatives and transport coefficients must be obtained from a more fundamental theory to complete the phenomenological description. Through the seminal efforts of Einstein, Onsager, Callen, Welton, Green, Kubo and many others [7, 1], a general linear response theory has been established. Furthermore, the transport coefficients are connected to the corresponding correlation functions (CFs) of the microscopic fluctuating variables. These connections are all embraced in the fundamental Green-Kubo relations. At equilibrium or small deviations from equilibrium, a systematic connection between the CFs and the hydrodynamics equations has been established for the long wave-length and small-frequency hydrodynamic limit [8, 9, 1]. In this hydrodynamic limit, the effects of the solvent fluctuations on suspended Brownian particles have also been studied extensively [10]. A further extension for small wave-length of the fluid has been made, which connects the microscopic dynamics to the generalized hydrodynamic equations [9]. Another breakthrough is the development of the fundamental fluctuation relations at transient or stationary nonequilibrium state far from equilibrium, which was initiated by Evans et al. [11] and has later engaged many others [12, 13, 14]. Many of these works on statistical mechanics are closely related to the large deviation theory in probability theory [15].

At a stationary nonequilibrium state, it seems relatively simpler to start with the phenomenological hydrodynamic equations and work reversely to obtain various CFs of the fluctuating variables [16, 17, 18, 19, 20, 21, 22, 23, 24, 25]. This strategy has been applied by Lutsko&Dufty [26] and has been receiving continuous attention [27, 28, 29]. In the present work, we consider an isothermal shear flow at steady state as a typical setting for the nonequilibrium behavior of many particles. Following the pioneering derivations of Lutsko&Dufty [26], the equations of fluctuating hydrodynamics are linearized around the steady state by assuming small fluctuations, before they are transformed into the Fourier space. Thereafter, an equivalent generalized eigenvalue problem is solved perturbatively to provide the temporal evolution of the hydrodynamic fluctuations. Finally, various autocorrelation functions (ACFs) can be constructed in the Fourier space and transformed into the real space if needed. Some analytical ACFs have been recently compared with inelastic hard-sphere simulations and multi-particle collision dynamics at low shear rates [27, 28, 29]. As the first objective of this work, we aim to verify the analytical ACFs at low/moderate shear rates and search for deviations from the theory at large shear rates via numerical simulations. We expect that our computations will reveal a transition from decoupling to the coupling of different fluctuating variables and may shed light on the possible extension of the theory at large shear rates. To this end, we employ a mesoscopic method called dissipative particle dynamics (DPD) to quantify such deviations.

The DPD method describes the behavior of many particles at mesoscale and was invented to bridge the gap between the microscopic dynamics and the macroscopic PDEs [30]. In a DPD system, three pairwise-additive forces ${\bf F}^{C}_{ij}$, ${\bf F}^{D}_{ij}$ and ${\bf F}^{R}_{ij}$ are prescribed between neighboring particles $i$ and $j$ and they correspond to the underlying conservative, dissipative and random process, respectively [31, 32]. By postulating a steady state solution of the corresponding Fokker-Planck equations , ${\bf F}^{D}_{ij}$ and ${\bf F}^{R}_{ij}$ are found to correlate with each other so that the fluctuation-dissipation theorem is satisfied and the canonical ensemble is warranted [31]. Given molecular dynamics (MD) trajectories, the pairwise forces in DPD may be constructed via a coarse-graining procedure following the Mori-Zwanzig formalism [33, 34, 35]. Without data from a reference MD simulation, the parameters of the pairwise forces are usually tuned to achieve static and dynamic properties of a target fluid empirically [36, 37]. Although the kinetic theory for the DPD particles can qualitatively predict the transport coefficients [38, 32], a quantitative knowledge is only available via a posteriori processing of the simulation results [32, 36, 39, 40, 41, 42, 43]. DPD typically has a softer potential between particles than that of MD, therefore it allows for a larger time step. In the hydrodynamic limit with large spatial-temporal scales, it may be considered as a Lagrangian discrete counterpart of the fluctuating hydrodynamics described by the Landau-Lifshitz-Navier-Stokes equations [44, 42]. At small wave-length and high frequency, it may be considered as the representation of the generalized hydrodynamics [45, 43], especially when the pairwise forces of DPD are obtained via the coarse-graining and non-Markovian effects are nonnegligible [46, 47, 48]. In addition to simulations of simple fluid at mesoscale, DPD also finds wide applications in simulating complex fluids such as colloidal suspension, polymer solution and red blood cells under flow [49, 50, 32, 36, 4]. As the second goal of this work, for the first time we evaluate various ACFs generated by DPD simulations under shear flow by comparing with analytical solutions. This would provide a solid evidence as to when DPD may be an effective solver of the fluctuating hydrodynamics at nonequilibrium.

Recently, multiscale coupling of heterogeneous solvers (e.g., molecular dynamics and Navier-Stokes solver) has been attracting a lot of attention [51, 52, 53, 54, 55, 56, 57]. Both the accurate dynamics of a fine model and the computational efficiency of a coarse model my be exploited in such a hybrid simulation. In the course of the coupling, from a continuum perspective the thermal fluctuations are treated very often as unwanted noises to be filtered out. However, the fluctuations are unique hallmarks to micro-/meso-scopic physics within a finite volume of material and their space-time correlations encode the full dynamic information of the system [1]. Therefore, as the third goal of this work, we advocate that the various ACFs should be taken as benchmarks for a hybrid simulation whenever it is possible. For example, one transversal ACF was previously compared among different coupling algorithms for a shear flow, which led to uncovering certain artifacts of the specific coupling [58].

We shall proceed as follows. In Section 2 we shall revisit the linearized fluctuating hydrodynamics and construct the temporal ACFs of the fluctuating variables in k-space. In Section 3 we will describe the DPD method with Lees-Edwards boundary conditions, and further elaborate on some technical details on the implementation of DPD simulations. In Section 4, we compare the theory with simulations for two set of transversal autocorrelation functions (TACFs) and longitudinal autocorrelation functions (LACFs) for various wave vectors for a range of shear rates. Finally, we summarize our findings in Section 5 with discussions. Extra details on the theoretical derivations are given in Appendices A and B.

2 Theory revisited

In this section, we follow the pioneering work of Lutsko&Dufty [26, 59] to derive analytically the ACFs of fluctuating variables in shear flow. Some of the calculations have also recently been performed [27, 28, 29]. Firstly, we describe the equations of fluctuating hydrodynamics in Section 2.1. Subsequently, in Appendix A we linearize the equations around the macroscopic or averaged state by keeping only the first-order fluctuating variables. Thereafter, in Appendix B we perform a spatial Fourier transform of the linearized equations. By applying a perturbation theory, we find the approximate solutions by solving a generalized eigenvalue problem. Finally, we summarize various ACFs of the fluctuating variables in Section 2.2.

2.1 Fluctuating hydrodynamics

Due to the conservation of mass and momentum, the equations of continuity and dynamics for an isothermal fluid read as,

[TABLE]

where $\rho$ is the mass density, ${\bf v}$ is the velocity and ${\bf\Pi}$ is the stress tensor. By defining the particle derivative as $\frac{d}{dt}=(\frac{\partial}{\partial t}+{\bf v}\cdot\nabla)$, we have the hydrodynamic equations in the Lagrangian form as

[TABLE]

The Lagrangian form of the hydrodynamic equations may be particularly useful to interpret the DPD method described in Section 3.

The stress tensor $\Pi_{\mu\sigma}$ may be considered as a linear combination of three components

[TABLE]

Each of the components is according to the reversible, irreversible, and stochastic process, respectively. Assuming that the pressure is isotropic and the viscous stress depends only on the first derivatives of velocity, $\Pi^{C}$ and $\Pi^{D}$ read as [44]

[TABLE]

where $\eta$ and $\zeta$ are constant coefficients of the shear and bulk viscosity. Summation convention for Greek indices is adopted and $\delta_{\mu\sigma}$ is the Kronecker $\delta$. Furthermore, the pressure $p=p(\rho)$ is determined by an equation of state at equilibrium. This completes the description for the compressible Navier-Stokes equations

[TABLE]

For fluids at mesoscale, there are fluctuations in the state variables governed by the framework of thermodynamics, therefore local spontaneous stress does occur. By assuming an underlying Gaussian-Makovian process for the unresolved degrees of freedom, the conditions for the random stress tensor satisfying the fluctuation-dissipation theorem read as [44]

[TABLE]

Here $\delta({\bf x}-{\bf x}^{\prime})$ or $\delta(t-t^{\prime})$ is the Dirac $\delta$ function. This completes the description for the Landau-Lifshitz-Navier-Stokes equations, which sometimes are also referred to as the Navier-Stokes-Langevin equations [26].

2.2 Autocorrelation functions under shear flow

For a uniform shear flow, the macroscopic stationary state of the velocity field ${\bf v}_{0}$ reads as

[TABLE]

This corresponds to a flow along the $x$ direction, with velocity gradient $\dot{\gamma}$ in the $y$ direction, and vorticity along the $z$ direction. By assuming small deviations from the averaged fields, the fluctuating hydrodynamic equations of Eqs. (1) and (2) may be linearized as shown in Appendix A, where second-order terms in fluctuations are dropped off. Afterwards, a system of linear equations may be spatial-Fourier transformed to obtain the hydrodynamic equations (48) in k-space. The general solution to Eq. (48) is determined from a generalized eigenvalue problem of Eq. (77). The latter is solved via a perturbation theory by expanding the size of wave vector $k=|{\bf k}|$ as a small parameter to second order in the continuum limit [26, 59]. Technical details of the derivations are elaborated in Appendices A and B.

Here, we need to emphasize one important assumption of the perturbation method. The shear rate is assumed to be moderate so that $\dot{\gamma}\lesssim\nu k^{2}\ll c_{T}k$. Hence, the term of $\dot{\gamma}$ is treated as in the order of $k^{2}$ during the perturbation. With the kinematic viscosity $\nu$ fixed, we may expect the perturbation method to fail for a very large $\dot{\gamma}$ or very small $k$. For our simulations in a cubic box with length $L$ in Section 4, the size of wave vector is $|{\bf k}|=|2\pi(n_{x},n_{y},n_{z})/L|$, where $n_{i}$ is an integer and $n^{2}_{x}+n^{2}_{y}+n^{2}_{z}\geq 1$. If a proper box length $L$ is selected, there is a minimal infrared cut off as $|{\bf k}|=2\pi/L$. Therefore, the shear rate $\dot{\gamma}$ is the only free parameter in the DPD simulations that we perform to validate the theory.

Given the derivations in Appendix B, the longitudinal and two transversal ACFs are readily constructed as follows [26, 28, 29]

[TABLE]

where kinematic viscosity $\nu=\eta/\rho_{0}$, sound speed and attenuation coefficient are $c_{T}$ and $\Gamma_{T}=(2\eta/3+\zeta/2)/\rho_{0}$. Here $\widehat{\delta u}_{1}({\bf k},t)$ is the longitudinal component along wave vector, while $\widehat{\delta u}_{2}({\bf k},t)$ and $\widehat{\delta u}_{2}({\bf k},t)$ are two transversal components, as defined in Eq. (46). Therefore, $C_{L}({\bf k},\tau)$ is the normalized longitudinal autocorrelation function (LACF), while $C_{T_{1}}({\bf k},\tau)$ and $C_{T_{2}}({\bf k},\tau)$ are the normalized first and second transversal autocorrelation functions (TACFs), respectively. As a matter of fact, the normalized ACFs are just the corresponding propagator of Eq. (105). Moreover, $\alpha$ and $\beta$ are defined as

[TABLE]

where the wave vector is time dependent as ${\bf k}(t)=(k_{x},k_{y}-\dot{\gamma}tk_{x},k_{z})$ to account for the advection. It is simple to see that when $\dot{\gamma}=0$, ${\bf k}(t)\equiv{\bf k}(0)={\bf k}_{0}=(k_{x},k_{y},k_{z}$) and $\alpha({\bf k},t)=k^{2}t$. Therefore, Eqs. (15) and (16) are identical and degenerate to the solutions for equilibrium [9, 1, 42].

To have a better sense of how the shear flow alters the dissipation and frequency of sound propagation, we plot the functions $\alpha({\bf k},t)$ and $\beta({\bf k},t)$ for three representative wave vectors in Fig. 1 and 2, respectively. For each wave vector ${\bf k}$, we also represent the equilibrium behavior $\dot{\gamma}=0$ with solid line as a reference accordingly. For example, in Fig. 1(b) with both $k_{x}>0$ and $k_{y}>0$, the terms with $t^{2}$ and $t^{3}$ compete with each other as time progresses, which is indicated by Eq. (17). As a consequence, the dissipation rate of $\dot{\gamma}>0$ is slower than that of $\dot{\gamma}=0$ at short time while the reverse is true at long time. The same phenomena in Figs. 1(a) and (c) are relatively simpler, hence, indicating a monotonic increase of dissipation rate with shear rate $\dot{\gamma}$.

Similarly, from Eq. (18) we may understand how the shear flow affects the frequency of sound propagation. For example, in Fig. 2(b) with $k_{x}>0$ and $k_{y}>0$, the constant slope of $\beta({\bf k},t)$ at $\dot{\gamma}=0$ indicates a constant frequency at equilibrium, whereas $\dot{\gamma}>0$ corresponds to frequency decrease at short time and frequency increase at long time. The same phenomena in Figs. 2(a) and (c) are less complicated, that is, the shear rate always enhances the sound frequency. From Eq. (18) and Fig. 2, we may expect Doppler effects for the LACFs.

3 Dissipative particle dynamics

The method of dissipative particle dynamics (DPD) was invented two decades ago to simulate rheological properties of complex fluids at mesoscale [30]. In this section, we shall briefly revisit the classical DPD method and its implementation of boundary conditions in shear flow.

3.1 Pairwise forces

For convenience, as reference we define some simple notations

[TABLE]

where ${\bf r}_{i}$, ${\bf v}_{i}$ are position and velocity of particle $i$; ${\bf r}_{ij}$, ${\bf v}_{ij}$ are relative position and velocity of particles $i$ and $j$; $r_{ij}$ is the distance between the two and ${\bf e}_{ij}$ is the unit vector pointing $j$ to $i$. The three pairwise forces are described as follows [31, 32],

[TABLE]

where coefficients $a$, $\gamma$, and $\sigma$ reflect the strength of individual forces; $W^{C}$, $W^{D}$, and $W^{R}$ are weighting functions that monotonically decay with the relative distance $r_{ij}$; $\theta_{ij}=\theta_{ji}$ is a Gaussian white noise with

[TABLE]

The DPD version of the fluctuation-dissipation theorem reads as

[TABLE]

which warrants the canonical ensemble [31].

Given the underlying force fields of molecular dynamics (MD), the actual forms of the three pairwise forces may be constructed via the Mori-Zwanzig projection [33, 35]. Without referring to any particular MD system, a typical empirical form of the weighting kernel is suggested as [32, 36]

[TABLE]

Following [58], we take $a=25.0$, $\sigma=3.0$, $\gamma=4.5$, $r_{c}=1$ and $k_{B}T=1.0$; $k=1$ for $W^{C}$ and $k=0.25$ for $W^{R}$; particle mass $m=1$, number density $n=3.0$ and mass density $\rho=mn=3.0$. For this particular set of input parameters, the dynamic viscosity, bulk viscosity and isothermal sound speed of the fluid are $\eta=1.62$, $\zeta=2.3$ and $c_{T}=4.05$ in DPD units [42]. The velocity Verlet time integrator is employed [32] and $\delta t=0.005$ is small enough for stability.

3.2 Boundary conditions

At nonequilibrium, the equal-time correlations of fluctuations are typically long-ranged [27] and we do not wish to introduce any extra complexities due to the boundary effects from the wall [60]. Therefore, we focus on the bulk behavior of the fluid in a periodic system, which also corresponds to the condition of the theory in Section 2. Suppose there is a simple shear flow defined as in Eq. (12), that is, the flow is in the $x$ direction, the velocity gradient $\dot{\gamma}$ is along the $y$ direction, and the vorticity is along the $z$ direction. The usual periodic boundary conditions apply in the $x$ and $z$ directions while periodic boxes along the $y$ direction shift $\pm L_{y}\dot{\gamma}t$, above and below the principal box, respectively. Therefore, if a particle crosses $y=L_{y}/2$ to outside, it enters the principal box again at $y=-L_{y}/2$ with $x$ shifted $-L_{y}\dot{\gamma}t$, and $v_{x}$ shifted by $-L_{y}\dot{\gamma}$; if the particle crosses $y=-L_{y}/2$ to outside, it enters the principal box again at $y=L_{y}/2$ with $x$ shifted $L_{y}\dot{\gamma}t$, and $v_{x}$ shifted by $L_{y}\dot{\gamma}$. Furthermore, the $x$ and $z$ positions are always wrapped back into the principal box due to the periodic boundaries. This is the so-called Lees-Edwards boundary condition [12], which degenerates to the usual periodic box when $\dot{\gamma}=0$.

In practice, we utilize an implementation of the deforming triclinic box for the periodic shear flow in the LAMMPS package [61]. In nonequilibrium MD, it is the so-called SLLOD dynamics for the canonical ensemble [12]. The technical difference is that DPD has a built-in pairwise thermostat while MD relies on other classical thermostats, such as the Nosé-Hoover thermostat.

4 Results

To compare the results of simulations with that of the theory, we perform DPD simulations in a box of $[0,L_{x}]\times[-L_{y}/2,L_{y}/2]\times[0,L_{z}]$ so that the mean velocity $v_{x}=\dot{\gamma}y$ is consistent with the definition of the analytical solutions. The domain is a cube with size $L_{x}=L_{y}=L_{z}=10$. Input parameters of DPD are given in Section 33.1. In the simulations, we define the fluctuating velocity under stationary shear flow as

[TABLE]

The ACFs in k-space are calculated as

[TABLE]

where direction $\sigma=1$, $2$, and $3$, and $N_{s}$ is the number of independent simulation runs. The Fourier transform in space is defined as

[TABLE]

where $j$ is particle index and $N_{p}$ is the number of particles in each simulation. Note that fluctuating velocities are projected on unit vectors in the wave vector coordinate via Eq. (31) after transformed in Eq. (30).

At equilibrium, there is no time origin, therefore time averaging may be performed before ensemble averaging in Eq. (29) so that good statistics are obtained. At nonequilibrium, due to the time dependence of ${\bf k}(t)$, basis vectors ${\bf e}^{\sigma}$ for wave vector coordinate is also time dependent (see Appendix B), and hence it is generally much more expensive to reduce the statistical errors of ACFs. Previously we have demonstrated that when the wave wave ${\bf k}_{0}=(0,0,k_{z}\neq 0)$ is along the vorticity direction [58], the two TACFs under shear flow degenerate to be isotropic just as that of equilibrium. Therefore, in the following results, we shall consider $k_{z}=0$ and focus on the wave vectors within the shear plane.

4.1 Transversal autocorrelation functions: ${\bf k}_{0}=(2\pi/L_{x},0,0)$

If a wave vector along $x$ direction as ${\bf k}_{0}=(2\pi/L_{x},0,0)$ is selected, the dissipation of TACFs according to Eq. (17) is

[TABLE]

Therefore, it is expected that the dissipation resembles the equilibrium behavior $\sim\tau$ at short time while it is dominated by the advection behavior $\sim\tau^{3}$ at long time. The distinction of two time regimes can be observed for both TACFs as shown in Fig. 3. Equilibrium results with $\dot{\gamma}=0$ are also shown as reference; with increasing shear rate $\dot{\gamma}$, the dissipation rate is enhanced. We note that the two TACFs have different intercepts with the $x$ axis, as the extra term $k_{0}/k(\tau)$ for $C_{T_{1}}({\bf k},\tau)$ in Eq. (15) causes a further stronger decay.

For $\dot{\gamma}=1.0$, the condition of $\dot{\gamma}\lesssim\nu k^{2}\approx 0.213$ is violated and we can clearly see that $C_{T_{1}}({\bf k},\tau)$ of the simulations deviates from the theory of Eq. (15). For $\dot{\gamma}=0.5$, although the condition of $\dot{\gamma}\lesssim\nu k^{2}$ is not strictly satisfied, $\dot{\gamma}$ can still be treated as in the order of $k^{2}$ in the perturbation method. Therefore, the TACFs for $\dot{\gamma}\leq 0.5$ from the simulations agree with those from the theory very well.

4.2 Transversal autocorrelation functions: ${\bf k}_{0}=(2\pi/L_{x},\pm 2\pi/L_{y},0)$

If a wave vector is selected within the shear plane with both nonzero $x$ and $y$ components, the dissipation rate according to Eq. (17) is

[TABLE]

If $k_{y}>0$, for example, ${\bf k}_{0}=(2\pi/L_{x},2\pi/L_{y},0)$, then the negative $\tau^{2}$ term competes with positive $\tau$ and $\tau^{3}$ terms. Therefore, compared to the case of equilibrium, dissipation may decrease or increase at different time regimes. This is well depicted for both $C_{T_{1}}({\bf k},t)$ and $C_{T_{2}}({\bf k},t)$ in Fig. 4. It is noteworthy that the discrepancy between $C_{T_{1}}({\bf k},t)$ of simulations and that of theory at $\dot{\gamma}=1$ is again apparent.

If a wave vector ${\bf k}_{0}=(2\pi/L_{x},-2\pi/L_{y},0)$ is selected, the three terms with different powers of $\tau$ in Eq. (33) are all positive. Therefore, only enhancement of dissipation is expected compared to the case of equilibrium, which is confirmed in Fig. 5. A slight discrepancy is noted for $C_{T_{1}}({\bf k},\tau)$ between the simulations and the theory at $\dot{\gamma}=1$.

4.3 Longitudinal autocorrelation functions: ${\bf k}_{0}=(2\pi/L_{x},0,0)$ and ${\bf k}_{0}=(2\pi/L_{x},\pm 2\pi/L_{y},0)$

In this section, we evaluate the LACFs from both the DPD simulations and the theory. As indicated in Eq. (14), the LACF has both damping and oscillating elements. The damping is primarily determined by the dissipation rate described by $\alpha({\bf k},\tau)$ function given in Eq. (17) and much of its behavior has been seen in Sections 4.1 and 4.2. Here we shall focus on the frequency of sound propagation, which is determined by the $\beta({\bf k},\tau)$ function and it is repeated here

[TABLE]

When the wave vector is ${\bf k}_{0}=(2\pi/L_{x},0,0)$ along $x$ direction, the sound frequency is expected to increase monotonically with shear rate $\dot{\gamma}$, as already indicated in Fig. 2(a). We further show the LACFs of both DPD simulations and the theory for this wave vector in Fig. 6(a). We clearly observe an overall agreement between the results of DPD simulation and the theory, and the sound frequency indeed increases monotonically with $\dot{\gamma}$. A small discrepancy between the simulation and the theory is observed for $\dot{\gamma}=1$.

When the wave vector is ${\bf k}_{0}=(2\pi/L_{x},2\pi/L_{y},0)$, it is already suggested in Fig. 2(b) that the sound frequency may increase or decrease at different time regimes. This is explicitly confirmed by the LACFs of both DPD simulations and the theory as shown in Fig. 6(b). We may observe that simulations agree with the theory at all shear rates considered.

When the wave vector is ${\bf k}_{0}=(2\pi/L_{x},-2\pi/L_{y},0)$, it is known from Fig. 2(c) that sound frequency increases with $\dot{\gamma}$. Again this can be observed in the LACFs in Fig. 6(c), where the results of the simulations agree with the theory at all shear rates considered.

4.4 Small length scales at $\dot{\gamma}=1$

Here we revisit various ACFs at $\dot{\gamma}=1$, where some discrepancies between the simulations and the theory have been observed for the largest length scale or smallest wave number considered. For those cases, the condition $\dot{\gamma}\lesssim\nu k^{2}$ is not met. Instead, here we focus on wave vectors ${\bf k}_{0}=(2n_{w}\pi/L_{x},0,0)$ and ${\bf k}_{0}=(2n_{w}\pi/L_{x},\pm 2n_{w}\pi/L_{y},0)$, but with $n_{w}=2$ or $3$. Therefore, the condition $\dot{\gamma}\lesssim\nu k^{2}$ becomes valid again for the smaller length scales and we expect that the perturbation method in deriving the analytical solutions is accurate. This can be confirmed by the TACFs $C_{T_{1}}({\bf k},t)$ for different wave vectors with $n_{w}=2$ and $3$, as shown in Fig. 7. The same results of $n_{w}=1$ as in Fig. 3(a), 4(a) and 5(a) are plotted again in Fig. 7 for comparison.

Similarly for the LACFs at $n_{w}=2$ and $3$, results of the simulations agree well with those of the theory, since the condition $\dot{\gamma}\lesssim\nu k^{2}$ is valid, as shown in Fig. 8. The same results of $n_{w}=1$ as in Fig. 6 are shown again in Fig. 8 for comparison.

It is noteworthy that at an even smaller scale ($n_{w}>3$), another type of discrepancy between the simulations and the theory emerges. This is due to the fact that the classical fluctuating hydrodynamics may not describe well the collective behavior of the DPD simulations at this small scales [45, 43]. In this case, solutions of the generalized fluctuating hydrodynamics may be required [9], which is beyond the scope of the current work.

5 Discussion and summary

We studied the autocorrelations (ACFs) of hydrodynamic fluctuations in k-space for an isothermal fluid under shear flow, which is driven by the Lees-Edwards periodic boundary condition. We compared results of the ACFs from the dissipative particle dynamics (DPD) simulations with the theoretical approximations. DPD is a particle-based mesoscopic method with three pairwise forces between neighboring particles. After specifying the conservative, dissipative and random forces, DPD is an effective simulator for the compressible fluctuating hydrodynamics. Under the assumption of a local thermodynamic equilibrium, the dissipative and random forces act as an intrinsic thermostat and they satisfy the fluctuation-dissipation theorem. In principle, the DPD method is valid for an arbitrarily large shear rate $\dot{\gamma}$, so are its ACFs of fluctuations, as long as local thermodynamic equilibrium is valid.

In order to derive the analytical solutions for the various ACFs, a perturbation method is adopted to expand the eigenvalue and eigenvector as a power series of $k$, where the parameter $k$ is the inverse of wave-length and it is very small in the hydrodynamic limit. In such regime, $\nu k^{2}\ll c_{T}k$ holds so that the rate of hydrodynamic dissipation is much smaller than the rate of sound propagation. Here $\nu$ and $c_{T}$ are the kinematic viscosity and isothermal sound speed, respectively. In the perturbation method, it is assumed that $\dot{\gamma}\lesssim\nu k^{2}$. Therefore, terms of $\dot{\gamma}$ are treated as in the order of $k^{2}$. With such an assumption, the generalized eigenvalue problem was solved approximately and thereafter the longitudinal and the two transversal ACFs were constructed as in Eqs. (14), (15), and (16).

When the condition of $\dot{\gamma}\lesssim\nu k^{2}$ is met, that is, Reynolds number $Re=\dot{\gamma}L^{2}/\nu\lesssim 4\pi^{2}$, various ACFs of the analytical solutions are accurate and agree very well with those of DPD simulations. We observed that the two transversal ACFs in shear flow are no longer identical as in equilibrium and the dissipation rate is time dependent as well. Furthermore, depending on the individual wave vector ${\bf k}$, enhancement or attenuation of sound frequency may take place and we observed the Doppler effects. Given the increasing efforts on hybrid modeling of fluid flow, where usually a stochastic microscopic/mesoscopic solver is concurrently coupled with a continuum (Landau-Lifshitz-)Navier-Stokes solver [53, 58], the agreement between analytical and computational approaches on the results is meaningful: the temporal correlations of fluctuating variables provide a fundamental benchmark for any proposed coupling algorithm.

When $\nu k^{2}<\dot{\gamma}<c_{T}k$, some discrepancies between the ACFs of the analytical solutions and those of the DPD simulations are observed. In this regime, besides the coupling between the advection and fluctuating variables, extra couplings between different fluctuating variables at equal time are expected to be significant. In our DPD simulations these extra couplings affect significantly both $C_{T_{1}}$ and $C_{L}$, but have a negligible effect on $C_{T_{2}}$. In this regime of shear rates, the perturbation theory is inaccurate and should be modified to treat the terms of $\dot{\gamma}$ as in the order of $k$ instead of $k^{2}$. Furthermore, the contributions from the stochastic stress on the temporal correlations are also expected to emerge. However, such modifications on the theory are nontrivial to accomplish and are subjects of our future research. Nevertheless, the DPD simulations are in principle valid and provide a few guidelines on how to improve the theory. In this regime, comparison with results from other numerical methods such as the finite volume method with thermal fluctuations [62, 63] would also be helpful to confirm the results presented.

When the shear rate becomes even larger such as $\dot{\gamma}>c_{T}k$, we obtain a nonequilibrium state far from equilibrium and additional couplings between different $k$ are expected (results are not shown). The local thermodynamic equilibrium is produced and enforced by the frequent collisions between particles and the frequency is measured by the sound speed $c_{T}$ divided by the length scale, that is $c_{T}k$. Therefore, under such strong shear flow the local thermodynamic equilibrium can no longer be restored in the length scale of $1/k$. It would also be interesting to evaluate the probability of violating the Second law in such a DPD system.

Acknowledgments

This work was supported by the Computational Mathematics Program within the Department of Energy office of Advanced Scientific Computing Research as part of the Collaboratory on Mathematics for Mesoscopic Modeling of Materials (CM4) and also supported by the ARO grant W911NF-14-1-0425. Part of this research was conducted using computational resources and services at the Center for Computation and Visualization, Brown University. An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under contract DE-AC02-06CH11357. This research also used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725. X. B. acknowledges discussions with Dr. Fanhai Zeng on the generalized eigenvalue problem.

Appendix A Linearization around uniform shear flow

The fluctuations on the state variables are defined as

[TABLE]

with

[TABLE]

where the macroscopic state variables are the averaged quantities denoted by “$<\quad>$”. Therefore, the fluctuating hydrodynamic equations can be linearized as

[TABLE]

where second-order terms in fluctuations are neglected. Kinematic viscosities are defined as $\nu=\eta/\rho_{0}$ and $\kappa=\zeta/\rho_{0}$. Furthermore, the equation of state is assumed to have the property of $\delta p=c^{2}_{T}\delta\rho$, where $c_{T}$ is the isothermal sound speed.

Appendix B Hydrodynamic matrix

It proves to be convenient to solve such a linearized hydrodynamic equations of Eqs. (38) and (39) in Fourier space [8, 9, 1]. Hence, we define the fluctuations in k-space as the spatial Fourier transform of the fluctuations

[TABLE]

Suppose that the wave vector is defined as ${\bf k}=k_{x}{\bf e}_{x}+k_{y}{\bf e}_{y}+k_{z}{\bf e}_{z}=(k_{x},k_{y},k_{z})$, where ${\bf e}_{x}$, ${\bf e}_{y}$ and ${\bf e}_{z}$ are three basis vectors in the fixed Cartesian coordinate. The formulation of the problem is very succinct, if the transformed velocity is decomposed into one component parallel to ${\bf k}$ and the other two components perpendicular to ${\bf k}$ [26, 27]. Therefore, we define another three orthonormal vectors ${\bf e}^{1}$, ${\bf e}^{2}$ and ${\bf e}^{3}$ in such way that ${\bf e}^{1}$ is along ${\bf k}$ and ${\bf e}^{2,3}$ are perpendicular to ${\bf k}$,

[TABLE]

Here ${\bf e}_{y}$ is taken as reference to define the first transversal direction ${\bf e}^{2}$ so that ${\bf e}^{2}$ is in the same plane with ${\bf e}^{1}$ and ${\bf e}_{y}$, and ${\bf e}^{2}$ is perpendicular to ${\bf e}^{1}$; moreover, ${\bf e}^{3}$ is perpendicular to both ${\bf e}^{1}$ and ${\bf e}^{2}$. Also $k=|{\bf k}|$, $k_{\bot}=(k^{2}-k^{2}_{y})^{1/2}/k$. Note that superscript $1$, $2$, or $3$ refers to the pairwise orthonormal basis vectors in the wave vector coordinate or oblique coordinate, which are expressed in the fixed Cartesian coordinate.

We further define a vector of fluctuating variables in k-space as

[TABLE]

where each element is related to the Fourier-transformed variable as

[TABLE]

The fluctuating velocity $\widetilde{\delta{\bf v}}({\bf k},t)$ in k-space is the spatial Fourier transform of the fluctuating velocity $\delta{\bf v}({\bf x},t)$ in real space. The latter is the peculiar or fluctuating velocity around the macroscopic velocity field defined as

[TABLE]

It is also simple to see that $\widehat{\delta{\bf u}}({\bf k},t)$ is obtained by projecting $\widetilde{\delta{\bf v}}({\bf k},t)$ onto the basis directions of the wave vector. Hence, $\widehat{\delta u}_{1}({\bf k},t)$ is called the longitudinal fluctuating component, while $\widehat{\delta u}_{2}({\bf k},t)$ and $\widehat{\delta u}_{3}({\bf k},t)$ are the first and second transversal fluctuating components, respectively. Finally, the hydrodynamic equations in k-space are summarized in a compact form as

[TABLE]

where the matrix ${\mathcal{L}}$ is defined as

[TABLE]

with

[TABLE]

[TABLE]

and

[TABLE]

The longitudinal viscosity is defined as $\nu_{L}=4\nu/3+\kappa=(4\eta/3+\zeta)/\rho_{0}$. $\Phi_{\mu\sigma}$ is defined by

[TABLE]

Given the definitions of the basis vectors from Eqs. (41), (42), and (43) in oblique coordinate, the elements of ${\mathcal{C}}$ read

[TABLE]

The linear mode coupling is represented by the derivatives with respect to ${\bf k}$.

Finally, the stochastic terms are

[TABLE]

The general solution to the linearized hydrodynamic equation (Eq. 48) can be determined from the nonlinear eigenvalue problem,

[TABLE]

The left eigenvectors $\chi^{(\mu)}$ are defined as

[TABLE]

where $\bar{\chi}^{(\mu)}_{\epsilon}$ means conjugate of $\chi^{(\mu)}_{\epsilon}$, and $\delta_{\mu\sigma}$ is the Kronecker $\delta$

The eigenvalues $\lambda_{\mu}$ and eigenvectors $\xi^{(\mu)}$ can be obtained via the perturbation theory by the expansion of $k$ to second and first order, respectively

[TABLE]

Inserting Eqs. (79) and (80) into Eq. (48) and making use of the explicit form of $\mathcal{L}$ in Eq. (49), we have the first-order perturbation theory involving ${\mathcal{A}}$ and second-order perturbation theory involving ${\mathcal{A}}$, ${\mathcal{B}}$, and ${\mathcal{C}}$ all together. Assuming $\dot{\gamma}\lesssim\nu k^{2}\ll c_{T}k$ and treating $\dot{\gamma}$ as in the order of $k^{2}$, the two equations from the perturbation read [26, 59]

[TABLE]

From Eq. (81), we find the first set of eigenvalues as follows

[TABLE]

From Eq. (82), we find the second set of eigenvalues as follows

[TABLE]

where sound attenuation coefficient $\Gamma_{T}=\nu_{L}/2=(2\eta/3+\zeta/2)/\rho_{0}$.

Inserting the expressions for $\lambda_{\mu,1}$ and $\lambda_{\mu,2}$ back into Eq. (79), finally the eigenvalues are summarized as

[TABLE]

Furthermore, the corresponding right eigenvectors are

[TABLE]

where

[TABLE]

The corresponding left eigenvectors satisfying Eq. (78) are

[TABLE]

Given the eigenvalues and eigenvectors $\lambda_{\mu}$, $\xi^{\mu}$ and $\chi^{\mu}$, the time evolution of fluctuating variables is obtained as

[TABLE]

with a propagator defined as

[TABLE]

To account for the advection, the wave vector is time dependent ${\bf k}(t)=(k_{x},k_{y}-\dot{\gamma}tk_{x},k_{z})$. Since the random terms do not contribute in this shear rate regime [26, 29], the correlations of the fluctuating variables are readily expressed as

[TABLE]

Given the values of $\lambda_{\mu}$, $\xi^{\mu}$ and $\chi^{\mu}$,, we find that $G_{\epsilon\iota}=0$ for $\epsilon\neq\iota$ and $G_{\epsilon\iota}\neq 0$ for $\epsilon=\iota$.