We analyze the hydrodynamic stability of force-driven parallel shear flows in nonequilibrium molecular simulations with three-dimensional periodic boundary conditions. We show that flows simulated in this way can be linearly unstable, and we derive an expression for the critical Reynolds number as a function of the geometric aspect ratio of the simulation domain. Approximate periodic extensions of Couette and Poiseuille flows are unstable at Reynolds numbers two orders of magnitude smaller than their aperiodic equivalents because the periodic boundaries impose fundamentally different constraints on the flow. This instability has important implications for simulating shear rheology and for designing nonequilibrium simulation methods that are compatible with periodic boundary conditions.

Periodic boundary conditions (PBCs) are mainstays of molecular simulations,1,2 where they are used to construct physically relevant models using material volumes that are often computationally restricted to be many orders of magnitude smaller than those in experiments. In standard PBCs, particles that exit the simulation volume reenter on the opposite face and interact with the nearest periodic images of the other particles. PBCs mitigate surface effects from the small volume, and simulated properties can agree well with experiments while modeling remarkably few particles.3,4 In addition, PBCs play an especially important role in nonequilibrium molecular simulations of shear rheology.5 Although wall-driven shear flows mimicking experimental rheometry can be modeled, surface effects unduly influence these simulations because the gap height between the surfaces that can be modeled is restricted, leading to unrealistic confinement that overemphasizes surface effects and slip compared to experiments. To overcome this limitation, simple or oscillatory shear flow of an unbounded fluid can be mimicked using Lees–Edwards PBCs,6 where the images of the fully periodic simulation cell are effectively put in relative motion at a given strain rate. However, determining properties like the shear viscosity by this method remains challenging, particularly at small strain rates, because the fluctuating stress must be carefully sampled.

Motivated by these complications, alternative simulation methods that fix the stress rather than the strain rate have been developed. Here, we focus on one class of these approaches that includes Müller-Plathe’s reverse nonequilibrium simulation (RNES) method7 and the periodic Poiseuille flow method.8 These methods do not generate flow with forcing at the boundaries but rather impose a spatially varying perturbation (i.e., a distributed body force) on the “unbounded” fluid that is compatible with the PBCs. With thoughtful choice of the implementation and form of the perturbation, rheological properties like the shear viscosity can be extracted from the measured flow field, which usually has rapidly converging statistics.7 These methods are convenient to implement and have been widely used to simulate rheological properties of fluids, such as simple liquids,9,10 water,11 ionic liquids,12,13 polymer solutions14,15 and melts,16–18 and colloidal dispersions.19–23 However, we were recently surprised to find that the RNES method can fail to produce the desired flow field in certain simulation geometries.24 RNES should generate a periodic parallel shear flow with two opposing Couette-like regions.7 We obtained this flow profile [Fig. 1(a)] in a cubic box, but undesired steady vortices [Fig. 1(b)] developed when the box was elongated in the flow direction. We found that the vortex formation depended on the presence of PBCs in the shear-gradient direction, the aspect ratio of the domain, and the shear rate. Moreover, the presence of vortices was not specific to the fluid model or using the RNES algorithm to create the flow. We speculated that this might indicate a general hydrodynamic instability underlying these nonequilibrium simulation techniques, but we were not able to establish a relationship for the conditions under which the instability occurred within the statistical accuracy of the simulations.

FIG. 1.

Streamlines for (a) stable periodic Couette flow in a cubic box at Reynolds number Re = 15 and (b) the steady secondary flows that developed when the box was elongated to aspect ratio α = 1.4 and the base flow became unstable. Both domains have PBCs in x, y, and z. The streamlines were averaged over z and colored by the flow velocity ux. All values are made dimensionless as in Eq. (3) and with x̃=x/H.

FIG. 1.

Streamlines for (a) stable periodic Couette flow in a cubic box at Reynolds number Re = 15 and (b) the steady secondary flows that developed when the box was elongated to aspect ratio α = 1.4 and the base flow became unstable. Both domains have PBCs in x, y, and z. The streamlines were averaged over z and colored by the flow velocity ux. All values are made dimensionless as in Eq. (3) and with x̃=x/H.

Close modal

Here, we analyze the hydrodynamic stability of parallel viscous shear flows that are generated by nonequilibrium (force-driven) simulation methods within three-dimensional PBCs. Hydrodynamic stability is a long-studied topic,25,26 but shear flows within a fully periodic domain27–33 have received less attention than wall-bounded flows. This may be partially due to the fact that experimentally observed flows rarely, if ever, have PBCs; however, PBCs are the norm rather than the exception in molecular simulations. We first consider the general case of a parallel shear flow in PBCs that is driven by a distributed piecewise continuous body force. We derive an expression for the critical Reynolds number for the linear instability as a function of the aspect ratio of the fully periodic domain. We then analyze the stability of three specific flows useful for simulating shear rheology: sinusoidal (Kolmogorov) flow, periodic Poiseuille flow, and a flow like the one generated by RNES that we call periodic Couette flow. We perform complementary particle-based simulations that confirm the analysis and shed valuable insight on the emergence of the instability. For a fixed Reynolds number, the results indicate that many parallel shear flows in PBCs become linearly unstable due to a geometric effect of the periodic domain, which admits unstable disturbances once it becomes sufficiently long in the flow direction. This work provides a simple framework for selecting geometries with PBCs that stabilize a flow, requiring only the desired flow’s Fourier series expansion. It also underscores the need for caution when simulating dynamic processes with PBCs, which can differ in unexpected ways from domains bounded by surfaces because the PBCs do not impose the same constraints as the surfaces.

The rest of this article is organized as follows. We perform the stability analysis in Sec. II, and then, we describe technical details of the particle-based simulations that we used to validate it in Sec. III. We compare the analysis to the simulations in Sec. IV, showing excellent agreement between the two, before concluding in Sec. V.

In the simulation methods of interest, the fluid is modeled as discrete interacting particles within a simulation cell having PBCs in three dimensions. Shear flow is generated by imposing an external force that varies in one dimension on the individual particles. For example, the periodic Poiseuille method applies a constant-magnitude body force to all fluid particles but with a direction that depends on whether the particles are in the upper or lower half of the simulation cell,8 while RNES swaps momentum between particles in two slab regions to impose an effective average local force.7 In order to analyze the hydrodynamic stability of these particle-based methods, we formulate a continuum description of the force-driven flow that develops.

We consider the incompressible flow u of a Newtonian fluid governed by the standard continuity and Navier–Stokes equations,34 

u=0,
(1)
ρut+uu=p+μ2u+f,
(2)

with density ρ, pressure p, viscosity μ, and body force per volume f. The velocity field and stresses must be continuous through the periodic boundaries at x = ±L and y = ±H. The third dimension z is also periodic, but its details will not be required for our analysis. In keeping with common practice in force-driven simulations,7,8 we also require that there is no net acceleration (translation) of the entire domain due to f, i.e., its average value is zero.

Following the standard linear stability analysis procedure,25 we assume that the velocity u = U + δu can be separated into two parts: the steady (time-independent) base flow U and a perturbation δu. (The pressure field, p = P + δp, is expressed similarly.) Both U and u must satisfy the PBCs of the domain. In particular, we are interested in the stability of parallel shear flows U = (Ux(y), 0, 0) generated by steady body forces f = (fx(y), 0, 0), so it suffices to consider the two-dimensional disturbance,35δu = eωt+ikxv(y) with v = (vx, vy, 0) and δp = eωt+ikxq(y). The specific form of Ux and fx will depend on the simulation method. Linearization of Eq. (2)Appendix A) yields

Reω̃+ik̃Ũxd2dỹ2k̃2ik̃d2Ũxdỹ2ṽy=d2dỹ2k̃22ṽy,
(3)

which has been written in dimensionless form by defining ỹ = y/H and Ũx = Ux/U, where U is the maximum of Ux, which also naturally gives ṽy = vy/U, k̃=Hk, and ω̃=ωH/U. The Reynolds number is Re = UH/ν with ν = μ/ρ being the kinematic viscosity. The base flow is unstable when the real part of ω̃ is positive for a solution ṽy that satisfies the PBCs in ỹ.

In all but a handful of cases, Eq. (3) must be solved numerically, e.g., by discretization or by expanding ṽy in a set of orthogonal basis functions.26 Orszag solved Eq. (3) for plane Poiseuille flow using Chebyshev polynomials,36 but other basis functions can be used. Due to the PBCs, it is advantageous to expand the flow fields in complex finite Fourier series that are inherently periodic, ṽy(ỹ) = ∑nṽneiπnỹ and Ũx(ỹ) = ∑pŨpeiπpỹ. These series can be directly substituted into Eq. (3) using termwise differentiation,

Reω̃n(πn)2+k̃2ṽneiπnỹ      +ik̃pnπ2(n2p2)+k̃2Ũpṽneiπ(n+p)ỹ     =n(πn)2+k̃22ṽneiπnỹ.
(4)

We then use the orthogonality of the basis functions to obtain

Reω̃(πn)2+k̃2ṽn+ik̃pπ2(n22np)+k̃2Ũpṽnp     =(πn)2+k̃22ṽn.
(5)

Note that the sum in Eq. (5), which resulted from products between ṽy, Ũx, and their derivatives in Eq. (3), is a convolution in Fourier space, as expected.

The coefficients Ũp are determined by the Fourier representation of a given base flow, so Eq. (5) can be written as an eigenvalue problem for ṽn by truncating the series to −Nn, pN for sufficiently large N. This problem can be solved for different values of k̃ and Re to find modes having unstable ω̃. Because of the PBC in x, k̃ is restricted in the values it can take by k = πm/L with m being an integer. In dimensionless variables, this gives k̃=πm/α with α = L/H being the aspect ratio of the domain. This constraint introduces the aspect ratio into Eq. (5),

iπmαp=NN12np(m/α)2+n2Ũpṽnp     π2[(m/α)2+n2]Reṽn=ω̃ṽn.
(6)

This approach can be considered a generalization of prior analysis for a sinusoidal flow27 to an arbitrary (force-driven) parallel shear flow Ux(y) in a fully periodic domain where the PBC aspect ratio α constrains k. To find the limit of stability for a given α, we numerically determined the eigenvalues of Eq. (6)36,37 as a function of Re and solved for the critical Reynolds number Rec that gave the least stable ω̃ having real part R(ω̃)=0 when m = 1. (Disturbances with larger m are unstable at larger α because they appear as a ratio.) In doing so, we neglected the possibility of instabilities triggered by nonorthogonal eigenmodes.38 

In order to test the stability analysis [Eq. (6)], we simulated force-driven flows using the computationally efficient multiparticle collision dynamics method.39–41 We used the stochastic rotation dynamics collision scheme39 with cubic cells of edge length a, fixed 130° rotation angle,42 and random grid shifting43; a Maxwell–Boltzmann rescaling thermostat to maintain constant temperature T;44 and time 0.1 τ between collisions, where τ=am/kBT, m is the particle mass, and kB is Boltzmann’s constant. The particle density was ρ = 5 m/a3, giving a liquid-like Newtonian fluid with kinematic viscosity ν = 0.79 a2/τ.45–47 All simulations were performed with hoomd-blue (version 2.6.0)48–50 using a three-dimensional periodic simulation box with a square cross section (H = 50 a in y and z) and L varied in x to span 0.6 ≤ α ≤ 3. (The total number of particles ranged from 3 × 106 to 15 × 106.) The maximum velocity of the base shear flow was restricted to U ≲ 0.32 a/τ to avoid artificial effects at large Mach numbers,51 giving Re ≤ 20.

We considered three force-driven shear flows: a sinusoidal (Kolmogorov) flow,27 a periodic Poiseuille flow,8 and a periodic Couette-like flow.7,24 The first is a historically well-studied problem27–31 that serves as a useful test of the simulations, while flows like the latter two are commonly used to simulate shear rheology. In our simulations, we applied a force per mass Fx(y) to each particle, giving fx = ρFx. (All particles had unit mass m, so Fx is also the force per particle.) Figure 2 shows Fx and Ux for the three flows when Re = 20, and details of their functional forms are given in  Appendix B. The periodic Poiseuille flow is the periodic extension of the classic Poiseuille flow, and it has two opposing parabolic regions in the periodic domain. The periodic Couette-like flow, which mimics the flow created by the RNES scheme, differs slightly from the true periodic extension of the classic Couette flow because it has two small parabolic regions that keep the shear stress continuous and Fx bounded, whereas the true periodic Couette flow would have a step change in the shear stress. For simplicity, we will refer to the simulated periodic Couette-like flow as periodic Couette flow.

FIG. 2.

(a) Applied force per particle mFx and (b) resulting flow Ux for the three cases simulated at Re = 20 with H = 50 aAppendix B). All quantities have units consistent with the simulations, and the maximum velocity is U ≈ 0.32 a/τ.

FIG. 2.

(a) Applied force per particle mFx and (b) resulting flow Ux for the three cases simulated at Re = 20 with H = 50 aAppendix B). All quantities have units consistent with the simulations, and the maximum velocity is U ≈ 0.32 a/τ.

Close modal

The simulations were initialized by superimposing the base flow on the thermalized fluid and were run for 105τ to reach a steady state. The velocity field u(x, y) was then measured by averaging the particle velocities in square bins of size a2 every 10 τ for 0.5 × 105τ. We used the measured u to assess whether the targeted base flow Ux was stable in the simulations. We adopted three empirical criteria to systematically classify the base shear flow as stable if:

  1. The deviation of the average total kinetic energy in y from equipartition [⟨Ey⟩ = NpkBT/2 for Np = (ρ/m)LH2 particles] was less than 3 standard errors of the mean as computed from its own fluctuations,

  2. the absolute deviation of the flow velocity from the base flow, |ux(y) − Ux(y)|, was larger than 10−3a/τ in fewer than 5% of the bins, and

  3. the absolute value of the gradient velocity |uy(x)| was larger than 10−3a/τ in fewer than 5% of the bins.

If all three criteria were not met, the base flow was classified as unstable; all other cases were unclassified. Criterion (1) is a coarse check for the development of flows directed along y, which are not present in the base parallel shear flow, while criteria (2) and (3) are more detailed tests of the measured flow against the expected flow. The thresholds for criteria (2) and (3) were chosen to be tolerant of the numerical and statistical fluctuations in equilibrium simulations (no flow) and so base flows with small α that were visibly stable were properly classified. Criterion (3) tended to be very stringent and often failed even when criteria (1) and (2) held.

We first considered the stability of the sinusoidal shear flow, Ũx(ỹ) = sin(πỹ) [Fig. 3(a)]. Its Fourier coefficients are Ũ±1 = ∓i/2 and zero otherwise. Secondary flows like those shown in Fig. 1 readily formed in the simulations at sufficiently large Re and α. During the accessible simulation time, the flows remained stationary, although subsequent transitions to chaotic motion have been reported.28,29 The bottom panel of Fig. 3(a) compares the predictions of the linear stability analysis with the simulation observations, which are in excellent agreement. Both indicate that the base flow remains stable for all α ≤ 1 over this range of Re, as proved theoretically by Meshalkin and Sinai.27 However, Fig. 3(a) also shows that some flows in domains with α > 1 are also stable for sufficiently small Re, and over the range of α tested, all flows having Re ≲ 5 were stable. We confirmed that α and Re are the appropriate dimensionless groups to use for comparison between the simulations and the stability analysis by performing selected simulations with smaller boxes and/or higher fluid viscosity. The flows were completely consistent with the boundaries in Fig. 3(a) when interpreted using α and Re.

FIG. 3.

(Top) Shear flows Ux(y) and (bottom) stability diagrams in aspect ratio α = L/H and Reynolds number Re = UH/ν for the (a) sinusoidal, (b) periodic Poiseuille, and (c) periodic Couette flows. The solid black lines indicate the critical Reynolds number Rec computed from Eq. (6), while the dashed green lines are the three-mode approximation of Eq. (10). The symbols show the observations from the simulations: stable (blue circle), unstable (red square), or unclassified (triangle).

FIG. 3.

(Top) Shear flows Ux(y) and (bottom) stability diagrams in aspect ratio α = L/H and Reynolds number Re = UH/ν for the (a) sinusoidal, (b) periodic Poiseuille, and (c) periodic Couette flows. The solid black lines indicate the critical Reynolds number Rec computed from Eq. (6), while the dashed green lines are the three-mode approximation of Eq. (10). The symbols show the observations from the simulations: stable (blue circle), unstable (red square), or unclassified (triangle).

Close modal

We subsequently tested two flows commonly used in nonequilibrium molecular simulations: the periodic extension of plane Poiseuille flow [Fig. 3(b)]8 and an approximate periodic extension of plane Couette flow [Fig. 3(c)] that mimics the RNES flow.7 Both were generated using constant forces in blocks of dimensionless half-width d̃=d/H, where d̃=1/2 for the periodic Poiseuille flow and we chose d̃=1/25 to approximate the periodic Couette flow [see Fig. 2(a) and  Appendix B]. The corresponding Fourier coefficients are

Ũp=4isin(pπ/2)sin(pπd̃)d̃(1d̃)(pπ)3.
(7)

Because the block forces are only piecewise continuous, Ũp ∼ 1/p3 and the Fourier series for d2Ũx/d ỹ2 converges pointwise rather than uniformly. We accordingly took N = 100 and confirmed that using N = 200 did not significantly affect the computed critical Reynolds numbers.

Similar to the sinusoidal flow, the simulated periodic Poiseuille and Couette flows became unstable at sufficiently large Re and α, in excellent agreement with the stability analysis [Figs. 3(b) and 3(c)]. In our previous study using RNES7 to create a periodic Couette flow, we found stable shear flows for α ≲ 1.25;24 this is in complete agreement with our new analysis for the studied parameters. We are unaware of prior reports of an instability in the periodic Poiseuille flow, but its stability was highly similar to the sinusoidal flow. This may not be surprising given that the two flows closely resemble each other [Figs. 3(a) and 3(b)]. The periodic Couette flow [Fig. 3(c)] was stable over a larger parameter space than the other two flows (see below).

Given the similar α dependence of Rec for the three flows, we posited that the stability might be controlled by only the first (long wavelength) terms in the series expansions and applied a three-mode approximation (N = 1).28,31 This drastic simplification allows the eigenvalues of Eq. (6) to be computed analytically,

ω̃0=π2(1+α2)Re,
(8)
ω̃±=π22Re1+2α2±1+8Re2π21α21+α2Ũ1Ũ1.
(9)

ω̃0 and ω̃+ will always be negative because Re > 0, but ω̃ can be zero or positive. We note that Ũ1 and Ũ−1 are complex conjugates, |Ũ1|2 = Ũ1Ũ−1, because Ũx is real valued. Solving for ω̃=0 gives a relationship between the critical Reynolds number Rec, α, and |Ũ1|,

Recπ|Ũ1|21+α21α2.
(10)

There is a minimum critical Reynolds number as α; this result is well-known for the sinusoidal flow,27,28,31 where all flows having Re<π24.4 are stable. However, Rec increases as α decreases toward one, expanding the range of Re for which flows are stable. Both are in good agreement with the simulations, and they clarify a point of uncertainty in Ref. 24 about whether there is a flow rate (Reynolds number) below which vortices should not be observed in a simulation box of fixed size.

Equation (10) approximates the solution of Eq. (6) well for all three flows when α ≳ 1.5 (Fig. 3). This explains the similar stability of the periodic Poiseuille (|Ũ1| = 16/π3 ≈ 0.52) and sinusoidal (|Ũ1| = 0.5) flows, as their first Fourier modes are nearly identical, and the stability of the periodic Couette flow to larger Re (|Ũ1| → 4/π2 ≈ 0.41 as d → 0). However, there are some discrepancies between Eqs. (6) and (10) at smaller α. In fact, Eq. (6) predicts that the periodic Poiseuille flow becomes unstable for Re ≳ 50 when α = 1, but Eq. (10) diverges as α → 1. Certain shear flows may still be unstable in cubic domains when PBCs are used. Nevertheless, Eq. (10) provides a useful estimate for either (1) choosing a Reynolds number for a fixed aspect ratio or (2) choosing an aspect ratio for a fixed Reynolds number that keeps the base flow stable. It also highlights the generality of the instability, as Eq. (10) indicates similarly shaped flows (to lowest order) become unstable under similar conditions.

Most of the secondary flows that developed in our simulations had two stationary vortices [Fig. 1(b)], although at the largest aspect ratio and Reynolds numbers simulated—α = 3 and Re = 17.5 and 20—we found four stationary vortices for the sinusoidal and periodic Poiseuille flows. We hypothesized that the structure of some of the secondary flows might be connected to the most unstable eigenmode. This eigenmode, w = eikxv(y), can be reconstructed using the Fourier coefficients for vy tht comprise the eigenvector having the largest R(ω) and computing vx = (i/k)dvy/dy based on the incompressibility of uAppendix A). Figure 4(a) shows w for the periodic Couette flow at Re = 15 and α = 1.4, which develops into the flow shown in Fig. 1(b). It has two pairs of counter-rotating vortices with streamlines primarily directed along the shear gradient y.

FIG. 4.

Streamlines for (a) the most unstable eigenmode w of the periodic Couette flow of Fig. 1(b) (α = 1.4 and Re = 15) computed using Eq. (6) and (b) the best-fit linear combination of U and w to the simulated u. All quantities have been made dimensionless as in Fig. 1.

FIG. 4.

Streamlines for (a) the most unstable eigenmode w of the periodic Couette flow of Fig. 1(b) (α = 1.4 and Re = 15) computed using Eq. (6) and (b) the best-fit linear combination of U and w to the simulated u. All quantities have been made dimensionless as in Fig. 1.

Close modal

The final steady flow need not possess the same structure as either the base flow or the most unstable eigenmode. However, we considered as an ansatz that in some cases the simulated flow field might be well-approximated by a linear combination of the two, uc1U + c2w, subject to a shift of the coordinates with respect to the periodic boundaries. As an example, we performed a least-squares regression to the simulated u for the conditions in Fig. 4, determining optimal coefficients c1 = 0.765 and c2 = 0.373. The fitted flow field [Fig. 4(b)] bears striking similarity to the simulated flow [Fig. 1(b)], having a dimensionless root-mean-squared error of 0.02 per velocity component.

We also noted that the waiting time τw for the secondary flows to emerge in the simulations was shorter for points farther from the curve Rec(α). This is qualitatively expected because these conditions are less stable based on their eigenvalues and should require smaller fluctuations (and less time) to depart from the base flow. To quantify this, we computed τw for the unstable flows by empirically fitting Ey to a hyperbolic tangent during the first 105τ simulated; we defined τw as the first time that Ey increased by 10% of the difference between its initial and final values ( Appendix C). These times can be compared to 1/R(ω), which is the typical timescale associated with growth of the most unstable eigenmode. We found that τw had a strong linear correlation with 1/R(ω) for all three flows (Fig. 5). Moreover, we noted that the instability could take a surprisingly long time to emerge, up to nearly 1.5 × 104τ in dimensional units. Having this estimate of the timescale for instability is an added benefit of our analysis.

FIG. 5.

Dimensionless waiting time τ̃w=τwU/H for secondary flows to emerge in the simulations compared to the dimensionless timescale associated with the most unstable eigenmode, 1/R(ω̃), for the sinusoidal (○), periodic Poiseuille (□), and periodic Couette flows (△). The points are colored according to Re, and the line is a fit to the data with slope 2.9.

FIG. 5.

Dimensionless waiting time τ̃w=τwU/H for secondary flows to emerge in the simulations compared to the dimensionless timescale associated with the most unstable eigenmode, 1/R(ω̃), for the sinusoidal (○), periodic Poiseuille (□), and periodic Couette flows (△). The points are colored according to Re, and the line is a fit to the data with slope 2.9.

Close modal

The stability of force-driven shear flows simulated in fully periodic domains seems to stand in stark contrast to those in bounded ones. The plane Couette flow is linearly stable for all infinitesimal disturbances,52 while the plane Poiseuille flow is linearly stable up to Re = 577236; in practice, both become unstable in the range of Re ≈ 102 to Re ≈ 103.53–55 However, their periodic extensions [Figs. 3(b) and 3(c)], driven by distributed body forces, are unstable at Re up to two orders of magnitude smaller. The change of stability in PBCs is not a simple consequence of the base flow; we previously showed that an unstable periodic Couette flow could be made stable in particle-based simulations by introducing a no-penetration boundary condition at either a half period or a full period of the flow.24 Instead, the PBCs impose fundamentally different constraints: (1) the velocities and stresses are not obligated to take a certain value at a surface, expanding the types of flows that can be realized, and (2) the PBCs restrict the wavelengths of disturbances to those commensurate with the domain, introducing a strong geometric dependence once the domain admits unstable modes.

We have given a simple recipe for determining the stability of parallel shear flows in spatially periodic domains that requires only the Fourier series expansion of the flows, and we have applied it to understand the flows we observed in nonequilibrium molecular simulations. This recipe can be used to design well-behaved models and simulation methods and to choose appropriate simulation parameters; in the case studied, we encourage performing simulations in domains having α < 1. The presence of a hydrodynamic instability in nonequilibrium simulation methods such as RNES, which we and others had not previously appreciated, also highlights the need for caution when using PBCs to simulate certain dynamic processes, as the PBCs may fundamentally alter the underlying physics in unexpected ways.

The data that support the findings of this study are openly available in Princeton University’s DataSpace at http://arks.princeton.edu/ark:/88435/dsp01zw12z8154.56 

We thank Arash Nikoubashman and Zachary Sherman for helpful comments on this manuscript. M.P.H. and T.M.T. acknowledge support from the Welch Foundation (Grant No. F-1696). A.S. was supported by the Princeton Center for Complex Materials (PCCM), a U.S. National Science Foundation Materials Research Science and Engineering Center (Grant No. DMR-1420541). The simulations were performed using computational resources supported by the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology’s High Performance Computing Center and Visualization Laboratory at Princeton University.

Although the derivation of Eq. (3) is widely documented,25,26 we include details here for completeness. Substituting u and p into Eqs. (1) and (2), using that U must also be a solution of the same, and neglecting terms of O(δu2) yields the linearized equations,

δu=0,
(A1)
ρδut+Uδu+δuU=δp+μ2δu.
(A2)

Equation (A1) gives a simple relation between vx and vy,

ikvx+vy=0,
(A3)

where the prime denotes the ordinary derivative, vy′ = dvy/dy, while Eq. (A2) gives

ρ(ωvx+ikUxvx+Uxvy)=ikq+μvxk2vx,
(A4)
ρ(ωvy+ikUxvy)=q+μvyk2vy.
(A5)

Equation (A3) can be inserted in Eq. (A4) and solved for q,

q=ρk2(ω+ikUx)vy+ρikUxvyμk2vyk2vy.
(A6)

Substituting Eq. (A6) into Eq. (A5) results in Eq. (3), which we note is the well-known Orr–Sommerfeld equation26 after replacing ω by .

In this appendix, we give explicit functional forms for the studied periodic flow fields and the body forces that we applied to generate them. All quantities have units consistent with the multiparticle collision dynamics simulations (see Sec. III).

The sinusoidal flow,

Ux(y)=UsinπyH,
(B1)

was generated by a sinusoidal force applied to each particle based on its position,

Fx(y)=UνH/π2sinπyH=FsinπyH.
(B2)

We chose U to obtain a certain Re, setting the force amplitude F = (πν)2Re/H3.

Periodic extensions of the plane Poiseuille flow and the plane Couette flow were generated similarly using piecewise constant forces. This functional form ensured that velocities and stresses in the fluid were continuous. We defined two blocks centered at y = ±H/2 of half-width dH/2 each. Then, a positive, constant force F was applied to all particles in the upper block and a negative, constant force −F was applied to all particles in the lower block,

Fx(y)=F,|y+H/2|dF,|yH/2|d0,otherwise.
(B3)

Linear momentum was conserved on average, but there were instantaneous fluctuations due to the distribution of particles. By solving Eq. (2) piecewise using the symmetries of the problem, the flow field can be deduced,

Ux(y)/U=g1(H+y),yH/2dg2(y),|y+H/2|<dg1(y),|y|H/2dg2(y),|yH/2|<dg1(Hy),yH/2+d,
(B4)

with

g1(y)=2y/H1d/H,
(B5)
g2(y)=1(12y/H)21(12d/H)2,
(B6)

where g1 is the usual linear (Couette-like) flow, and g2 is the quadratic (Poiseuille-like) flow. The maximum velocity is

U=FH28ν1(12d/H)2.
(B7)

For the periodic plane Poiseuille flow,8 the width of the blocks was chosen to cover the full simulation box, d = H/2. Only g2 contributes to this flow, establishing two opposing parabolic regions, and the force magnitude required for a given Reynolds number is F = 8ν2Re/H3.

By making d much smaller than H/2, the flow is dominated by g1 and a periodic plane Couette-like flow can be generated having two linear regimes. In practice, d must be sufficiently large that there are enough particles in each block to reliably apply the force and achieve the targeted maximum velocity. The force magnitude required for a given Reynolds number is F=8ν2Re/H3(1(12d/H)2). In the limit d → 0, the shear stress has a step change at y = ±H/2, and the applied forces are necessarily delta functions. For finite d, however, F is bounded, and the flow has two small quadratic regions that keep the shear stress continuous. We chose d = 2 a for our boxes having H = 50 a, which gave predominantly Couette-like flow fields but also ensured the blocks contained sufficient numbers of particles. The resulting flow is highly similar to that generated by the RNES method7 even though the origin of the perturbation is different.

In order to estimate the timescale for occurrence of the instability, we empirically fit the kinetic energy Ey to a hyperbolic tangent,

Ey(t)NpkBT=12(E0+E1)+(E1E0)tanhtτ0w.
(C1)

We fixed E0 = 1/2 because the system was initially thermalized, and so its kinetic energy must obey equipartition. We then fit the final value E1, the midpoint τ0, and the width w, resulting in fits like those shown in Fig. 6. The waiting time was defined as τw = τ0 − 1.09861w using the “10–90 thickness.” We found it challenging to reliably fit Ey for a handful of conditions close to the stability curves of Fig. 2; we accordingly neglected all points having 1/R(ω̃)>10 in creating and analyzing Fig. 4.

FIG. 6.

Fit (red) to Eq. (C1) for kinetic energy in y per particle (gray) for the periodic Couette flow with α = 1.4 and Re = 15. The dashed blue line indicates the expected value of Ey/Np from equipartition. The black line marks the waiting time τw for occurrence of the instability.

FIG. 6.

Fit (red) to Eq. (C1) for kinetic energy in y per particle (gray) for the periodic Couette flow with α = 1.4 and Re = 15. The dashed blue line indicates the expected value of Ey/Np from equipartition. The black line marks the waiting time τw for occurrence of the instability.

Close modal
1.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
(
Oxford University Press
,
New York
,
1991
).
2.
D.
Frenkel
and
B.
Smit
,
Understanding Molecular Simulation
, 2nd ed. (
Academic Press
,
San Diego
,
2002
).
3.
B. J.
Alder
and
T. E.
Wainwright
,
J. Chem. Phys.
27
,
1208
(
1957
).
4.
F. H.
Stillinger
and
A.
Rahman
,
J. Chem. Phys.
60
,
1545
(
1973
).
5.
W. G.
Hoover
,
Annu. Rev. Phys. Chem.
34
,
103
(
1983
).
6.
A. W.
Lees
and
S. F.
Edwards
,
J. Phys. C: Solid State Phys.
5
,
1921
(
1972
).
7.
F.
Müller-Plathe
,
Phys. Rev. E
59
,
4894
(
1999
).
8.
J. A.
Backer
,
C. P.
Lowe
,
H. C. J.
Hoefsloot
, and
P. D.
Iedema
,
J. Chem. Phys.
122
,
154503
(
2005
).
9.
T.
Soddemann
,
B.
Dünweg
, and
K.
Kremer
,
Phys. Rev. E
68
,
046702
(
2003
).
10.
M. S.
Kelkar
,
J. L.
Rafferty
,
E. J.
Maginn
, and
J. I.
Siepmann
,
Fluid Phase Equilib.
260
,
218
(
2007
).
11.
Y.
Mao
and
Y.
Zhang
,
Chem. Phys. Lett.
542
,
37
(
2012
).
12.
M. S.
Kelkar
and
E. J.
Maginn
,
J. Phys. Chem. B
111
,
4867
(
2007
).
13.
W.
Zhao
,
F.
Leroy
,
S.
Balasubramanian
, and
F.
Müller-Plathe
,
J. Phys. Chem. B
112
,
8129
(
2008
).
14.
A.
Nikoubashman
and
M. P.
Howard
,
Macromolecules
50
,
8279
(
2017
).
15.
E.
Moghimi
,
I.
Chubak
,
A.
Statt
,
M. P.
Howard
,
D.
Founta
,
G.
Polymeropoulos
,
K.
Ntetsikas
,
N.
Hadjichristidis
,
A. Z.
Panagiotopoulos
,
C. N.
Likos
, and
D.
Vlassopoulos
,
ACS Macro Lett.
8
,
766
(
2019
).
16.
H.
Guo
,
K.
Kremer
, and
T.
Soddemann
,
Phys. Rev. E
66
,
061503
(
2002
).
17.
L.
Schneider
,
M.
Heck
,
M.
Wilhelm
, and
M.
Müller
,
Macromolecules
51
,
4642
(
2018
).
18.
L.
Schneider
and
M.
Müller
,
Comput. Mater. Sci.
169
,
109107
(
2019
).
19.
D. R.
Heine
,
M. K.
Petersen
, and
G. S.
Grest
,
J. Chem. Phys.
132
,
184509
(
2010
).
20.
M.
Cerbelaud
,
A.
Maria Laganapan
,
T.
Ala-Nissila
,
R.
Ferrando
, and
A.
Videcoq
,
Soft Matter
13
,
3909
(
2017
).
21.
R. D.
Mountain
,
H. W.
Hatch
, and
V. K.
Shen
,
Fluid Phase Equilib.
440
,
87
(
2017
).
22.
A.
Sambasivam
,
S.
Dhakal
, and
R.
Sureshkumar
,
Mol. Simul.
44
,
485
(
2018
).
23.
J. D.
Olarte-Plata
and
F.
Bresme
,
Mol. Phys.
116
,
2032
(
2018
).
24.
A.
Statt
,
M. P.
Howard
, and
A. Z.
Panagiotopoulos
,
Phys. Rev. Fluids
4
,
043905
(
2019
).
25.
C. C.
Lin
,
The Theory of Hydrodynamic Stability
(
Cambridge University Press
,
Cambridge
,
1955
).
26.
P. G.
Drazin
and
W. H.
Reid
,
Hydrodynamic Stability
, 2nd ed. (
Cambridge University Press
,
Cambridge
,
2004
).
27.
L. D.
Meshalkin
and
I. G.
Sinai
,
J. Appl. Math. Mech.
25
,
1700
(
1961
).
28.
J. S. A.
Green
,
J. Fluid Mech.
62
,
273
(
1974
).
29.
R.
Grappin
,
J.
Leorat
, and
P.
Londrillo
,
J. Fluid Mech.
195
,
239
(
1988
).
30.
A.
Thess
,
Phys. Fluids A
4
,
1385
(
1992
).
31.
I.
Bena
,
M. M.
Mansour
, and
F.
Baras
,
Phys. Rev. E
59
,
5503
(
1999
).
32.
I. E.
Sarris
,
H.
Jeanmart
,
D.
Carati
, and
G.
Winckelmans
,
Phys. Fluids
19
,
095101
(
2007
).
33.
H. R.
Dullin
and
J.
Worthington
,
J. Math. Fluid Mech.
20
,
473
(
2018
).
34.
W. M.
Deen
,
Analysis of Transport Phenomena
, 2nd ed. (
Oxford University Press
,
New York
,
2012
).
35.
H. B.
Squire
,
Proc. R. Soc. A
142
,
621
(
1933
).
36.
S. A.
Orszag
,
J. Fluid Mech.
50
,
689
(
1971
).
37.
C. L.
Dolph
and
D. C.
Lewis
,
Q. Appl. Math.
16
,
97
(
1958
).
38.
L. N.
Trefethen
,
A. E.
Trefethen
,
S. C.
Reddy
, and
T. A.
Driscoll
,
Science
261
,
578
(
1993
).
39.
A.
Malevanets
and
R.
Kapral
,
J. Chem. Phys.
110
,
8605
(
1999
).
40.
G.
Gompper
,
T.
Ihle
,
D. M.
Kroll
, and
R. G.
Winkler
, in
Advanced Computer Simulation Approaches for Soft Matter Sciences III
, Advances in Polymer Science Vol. 221, edited by
C.
Holm
and
K.
Kremer
(
Springer
,
Berlin
,
2009
), pp.
1
87
.
41.
M. P.
Howard
,
A.
Nikoubashman
, and
J. C.
Palmer
,
Curr. Opin. Chem. Eng.
23
,
34
(
2019
).
42.
E.
Allahyarov
and
G.
Gompper
,
Phys. Rev. E
66
,
036702
(
2002
).
43.
T.
Ihle
and
D. M.
Kroll
,
Phys. Rev. E
63
,
020201(R)
(
2001
).
44.
C.-C.
Huang
,
A.
Varghese
,
G.
Gompper
, and
R. G.
Winkler
,
Phys. Rev. E
91
,
013310
(
2015
).
45.
T.
Ihle
and
D. M.
Kroll
,
Phys. Rev. E
67
,
066706
(
2003
).
46.
M.
Ripoll
,
K.
Mussawisade
,
R. G.
Winkler
, and
G.
Gompper
,
Phys. Rev. E
72
,
016701
(
2005
).
47.
J. T.
Padding
and
A. A.
Louis
,
Phys. Rev. E
74
,
031402
(
2006
).
48.
J. A.
Anderson
,
C. D.
Lorenz
, and
A.
Travesset
,
J. Comput. Phys.
227
,
5342
(
2008
).
49.
J.
Glaser
,
T. D.
Nguyen
,
J. A.
Anderson
,
P.
Lui
,
F.
Spiga
,
J. A.
Millan
,
D. C.
Morse
, and
S. C.
Glotzer
,
Comput. Phys. Commun.
192
,
97
(
2015
).
50.
M. P.
Howard
,
A. Z.
Panagiotopoulos
, and
A.
Nikoubashman
,
Comput. Phys. Commun.
230
,
10
(
2018
).
51.
A.
Lamura
,
G.
Gompper
,
T.
Ihle
, and
D. M.
Kroll
,
Europhys. Lett.
56
,
319
(
2001
).
52.
V. A.
Romanov
,
Funct. Anal. Appl.
7
,
137
(
1973
).
53.
S. A.
Orszag
and
L. C.
Kells
,
J. Fluid Mech.
96
,
159
(
1980
).
54.
B. J.
Bayly
,
S. A.
Orszag
, and
T.
Herbert
,
Annu. Rev. Fluid Mech.
20
,
359
(
1988
).
55.
A.
Lundbladh
and
A. V.
Johansson
,
J. Fluid Mech.
229
,
499
(
1991
).
56.
A.
Statt
,
M. P.
Howard
,
H. A.
Stone
, and
T. M.
Truskett
, Instability of Shear Flows in Spatially Periodic Domains, Princeton University DataSpace, http://arks.princeton.edu/ark:/88435/dsp01zw12z8154,
2019
.