Stability of force-driven shear flows in nonequilibrium molecular simulations with periodic boundaries

TL;DR

This paper investigates the stability of force-driven shear flows in nonequilibrium molecular simulations with periodic boundaries, revealing that such flows can be linearly unstable at much lower Reynolds numbers than traditional flows, impacting simulation methods.

Contribution

It derives a critical Reynolds number expression for these flows and highlights how periodic boundaries fundamentally alter flow stability compared to aperiodic conditions.

Findings

Flows in periodic domains are unstable at Reynolds numbers two orders of magnitude lower than in aperiodic flows.

The derived critical Reynolds number depends on the geometric aspect ratio of the simulation domain.

Implications for shear rheology simulations and nonequilibrium method design are discussed.

Abstract

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.

Equations

$$\nabla \cdot \mathbf{u}$$

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

$$\mathrm{Re}\left[\left(\tilde{\omega }+i\tilde{k}{\tilde{U}}_x\right)\left(\frac{d{}^2}{d{\tilde{y}}^2}-{\tilde{k}}^2\right)-i\tilde{k}\frac{d{}^2{\tilde{U}}_x}{d{\tilde{y}}^2}\right]{\tilde{v}}_y={\left(\frac{d{}^2}{d{\tilde{y}}^2}-{\tilde{k}}^2\right)}^2{\tilde{v}}_y,$$

$$-\mathrm{Re}\left[\tilde{\omega }\sum _n\left((\pi n{)}^2+{\tilde{k}}^2\right){\tilde{v}}_n{e}^{i\pi n\tilde{y}}+i\tilde{k}\sum _p\sum _n\left({\pi }^2(n^2-p^2)+{\tilde{k}}^2\right){\tilde{U}}_p{\tilde{v}}_n{e}^{i\pi (n+p)\tilde{y}}\right]$$

$$=\sum _n{\left((\pi n{)}^2+{\tilde{k}}^2\right)}^2{\tilde{v}}_n{e}^{i\pi n\tilde{y}}.$$

$$-\mathrm{Re}\left[\tilde{\omega }\left((\pi n{)}^2+{\tilde{k}}^2\right){\tilde{v}}_n+i\tilde{k}\sum _p\left({\pi }^2(n^2-2np)+{\tilde{k}}^2\right){\tilde{U}}_p{\tilde{v}}_{n-p}\right]={\left((\pi n{)}^2+{\tilde{k}}^2\right)}^2{\tilde{v}}_n.$$

$$-i\left(\frac{\pi m}{\alpha }\right)\sum _{p=-N}^N\left[1-\frac{2np}{(m/\alpha {)}^2+n^2}\right]{\tilde{U}}_p{\tilde{v}}_{n-p}-\frac{{\pi }^2[(m/\alpha {)}^2+n^2]}{\mathrm{Re}}{\tilde{v}}_n=\tilde{\omega }{\tilde{v}}_n.$$

$${\tilde{U}}_p=\frac{4}{i}\frac{\sin (p\pi /2)\sin (p\pi \tilde{d})}{\tilde{d}(1-\tilde{d})(p\pi {)}^3}.$$

$${\tilde{\omega }}_0$$

$${\tilde{\omega }}_{\pm }$$

$${\mathrm{Re}}_{\mathrm{c}}\approx \frac{\pi }{|{\tilde{U}}_1|\sqrt{2}}\frac{1+{\alpha }^{-2}}{\sqrt{1-{\alpha }^{-2}}}.$$

$$\nabla \cdot \delta \mathbf{u}$$

$$\rho \left[\frac{\partial \delta \mathbf{u}}{\partial t}+\mathbf{U}\cdot \nabla \delta \mathbf{u}+\delta \mathbf{u}\cdot \nabla \mathbf{U}\right]$$

$$ik{v}_x+v_y^{\prime }=0,$$

$$\rho (\omega v_x+ik{U}_x{v}_x+U_x^{\prime }{v}_y)$$

$$\rho (\omega v_y+ik{U}_x{v}_y)$$

$$-q=\frac{\rho }{k^2}(\omega +ik{U}_x)v_y^{\prime }+\frac{\rho }{ik}{U}_x^{\prime }{v}_y-\frac{\mu }{k^2}\left[v_y^{\prime \prime \prime }-k^2{v}_y^{\prime }\right].$$

$$U_x(y)=U\sin \left(\frac{\pi y}{H}\right),$$

$$F_x(y)=\frac{U\nu }{{\left(H/\pi \right)}^2}\sin \left(\frac{\pi y}{H}\right)=F\sin \left(\frac{\pi y}{H}\right).$$

$$F_x(y)=\{\begin{array}{ll}-F,& |y+H/2|\ge d\\ F,& |y-H/2|\ge d\\ 0,& \mathrm{otherwise}\end{array}.$$

$$U_x(y)/U=\{\begin{array}{ll}-g_1(H+y),& y\le -H/2-d\\ -g_2(-y),& |y+H/2|<d\\ g_1(y),& |y|\le H/2-d\\ g_2(y),& |y-H/2|<d\\ g_1(H-y),& y\ge H/2+d\end{array}$$

$$g_1(y)$$

$$g_2(y)$$

$$U=\frac{F{H}^2}{8\nu }\left[1-(1-2d/H{)}^2\right].$$

$$\frac{E_y(t)}{N_{\mathrm{p}}{k}_{\mathrm{B}}T}=\frac{1}{2}\left[(E_0+E_1)+(E_1-E_0)\tanh \left(\frac{t-{\tau }_0}{w}\right)\right].$$

Full text

I Introduction

Periodic boundary conditions (PBCs) are mainstays of molecular simulations Allen and Tildesley (1991); Frenkel and Smit (2002), 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 Alder and Wainwright (1957); Stillinger and Rahman (1973). Also, PBCs play an especially important role in nonequilibrium molecular simulations of shear rheology Hoover (1983). 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 Lees and Edwards (1972), 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) method Müller-Plathe (1999) and the periodic Poiseuille flow method Backer et al. (2005). 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 Müller-Plathe (1999). These methods are convenient to implement and have been widely used to simulate rheological properties of fluids such as simple liquids Soddemann, Dünweg, and Kremer (2003); Kelkar et al. (2007), water Mao and Zhang (2012), ionic liquids Kelkar and Maginn (2007); Zhao et al. (2008), polymer solutions Nikoubashman and Howard (2017); Moghimi et al. (2019) and melts Guo, Kremer, and Soddemann (2002); Schneider et al. (2018); Schneider and Müller (2019), and colloidal dispersions Heine, Petersen, and Grest (2010); Cerbelaud et al. (2017); Mountain, Hatch, and Shen (2017); Sambasivam, Dhakal, and Sureshkumar (2018); Olarte-Plata and Bresme (2018). However, we were recently surprised to find that the RNES method can fail to produce the desired flow field in certain simulation geometries Statt, Howard, and Panagiotopoulos (2019). RNES should generate a periodic parallel shear flow with two opposing Couette-like regions Müller-Plathe (1999). We obtained this flow profile (Fig. 1a) in a cubic box, but undesired steady vortices (Fig. 1b) 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.

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 Lin (1955); Drazin and Reid (2004), but shear flows within a fully periodic domain Meshalkin and Sinai (1961); Green (1974); Grappin, Leorat, and Londrillo (1988); Thess (1992); Bena, Malek Mansour, and Baras (1999); Sarris et al. (2007); Dullin and Worthington (2018) 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 Section II, and then we describe technical details of the particle-based simulations that we used to validate it in Section III. We compare the analysis to the simulations in Section IV, showing excellent agreement between the two, before concluding in Section V.

II Stability analysis

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 Backer et al. (2005), while RNES swaps momentum between particles in two slab regions to impose an effective average local force Müller-Plathe (1999). 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 $\mathbf{u}$ of a Newtonian fluid governed by the standard continuity and Navier–Stokes equations Deen (2012),

[TABLE]

with density $\rho$, pressure $p$, viscosity $\mu$, and body force per volume $\mathbf{f}$. The velocity field and stresses must be continuous through the periodic boundaries at $x=\pm L$ and $y=\pm 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 Müller-Plathe (1999); Backer et al. (2005), we also require that there is no net acceleration (translation) of the entire domain due to $\mathbf{f}$, i.e., its average value is zero.

Following the standard linear stability analysis procedure Lin (1955), we assume that the velocity $\mathbf{u}=\mathbf{U}+\delta\mathbf{u}$ can be separated into two parts: the steady (time-independent) base flow $\mathbf{U}$ and a perturbation $\delta\mathbf{u}$. (The pressure field, $p=P+\delta p$ is expressed similarly.) Both $\mathbf{U}$ and $\mathbf{u}$ must satisfy the PBCs of the domain. In particular, we are interested in the stability of parallel shear flows $\mathbf{U}=(U_{x}(y),0,0)$ generated by steady body forces $\mathbf{f}=(f_{x}(y),0,0)$, so it suffices to consider the two-dimensional disturbance Squire (1933), $\delta\mathbf{u}=e^{\omega t+ikx}\mathbf{v}(y)$ with $\mathbf{v}=(v_{x},v_{y},0)$ and $\delta p=e^{\omega t+ikx}q(y)$. The specific form of $U_{x}$ and $f_{x}$ will depend on the simulation method. Linearization of Eq. (2) (Appendix A) yields

[TABLE]

which has been written in dimensionless form by defining $\tilde{y}=y/H$ and $\tilde{U}_{x}=U_{x}/U$, where $U$ is the maximum of $U_{x}$, which also naturally gives $\tilde{v}_{y}=v_{y}/U$, $\tilde{k}=Hk$, and $\tilde{\omega}=\omega H/U$. The Reynolds number is ${\rm Re}=UH/\nu$ with $\nu=\mu/\rho$ being the kinematic viscosity. The base flow is unstable when the real part of $\tilde{\omega}$ is positive for a solution $\tilde{v}_{y}$ that satisfies the PBCs in $\tilde{y}$.

In all but a handful of cases, Eq. (3) must be solved numerically, e.g., by discretization or by expanding $\tilde{v}_{y}$ in a set of orthogonal basis functions Drazin and Reid (2004). Orszag solved Eq. (3) for plane Poiseuille flow using Chebyshev polynomials Orszag (1971), 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, $\tilde{v}_{y}(\tilde{y})=\sum_{n}\tilde{v}_{n}e^{i\pi n\tilde{y}}$ and $\tilde{U}_{x}(\tilde{y})=\sum_{p}\tilde{U}_{p}e^{i\pi p\tilde{y}}$. These series can be directly substituted into Eq. (3) using termwise differentiation:

[TABLE]

We then use the orthogonality of the basis functions to obtain

[TABLE]

Note that the sum in Eq. (5), which resulted from products between $\tilde{v}_{y}$, $\tilde{U}_{x}$, and their derivatives in Eq. (3), is a convolution in Fourier space, as expected.

The coefficients $\tilde{U}_{p}$ are determined by the Fourier representation of a given base flow, so Eq. (5) can be written as an eigenvalue problem for $\tilde{v}_{n}$ by truncating the series to $-N\leq n,p\leq N$ for sufficiently large $N$. This problem can be solved for different values of $\tilde{k}$ and $\mathrm{Re}$ to find modes having unstable $\tilde{\omega}$. Because of the PBC in $x$, $\tilde{k}$ is restricted in the values it can take by $k=\pi m/L$ with $m$ being an integer. In dimensionless variables, this gives $\tilde{k}=\pi m/\alpha$ with $\alpha=L/H$ being the aspect ratio of the domain. This constraint introduces the aspect ratio into Eq. (5),

[TABLE]

This approach can be considered a generalization of prior analysis for a sinusoidal flow Meshalkin and Sinai (1961) to an arbitrary (force-driven) parallel shear flow $U_{x}(y)$ in a fully periodic domain where the PBC aspect ratio $\alpha$ constrains $k$. To find the limit of stability for a given $\alpha$, we numerically determined the eigenvalues of Eq. (6) Orszag (1971); Dolph and Lewis (1958) as a function of Re and solved for the critical Reynolds number ${\rm Re}_{\rm c}$ that gave the least stable $\tilde{\omega}$ having real part $\Re(\tilde{\omega})=0$ when $m=1$. (Disturbances with larger $m$ are unstable at larger $\alpha$ because they appear as a ratio.) In doing so, we neglected the possibility of instabilities triggered by nonorthogonal eigenmodes Trefethen et al. (1993).

III Simulation methods

In order to test the stability analysis (Eq. (6)), we simulated force-driven flows using the computationally efficient multiparticle collision dynamics method Malevanets and Kapral (1999); Gompper et al. (2009); Howard, Nikoubashman, and Palmer (2019). We used the stochastic rotation dynamics collision scheme Malevanets and Kapral (1999) with cubic cells of edge length $a$, fixed $130^{\circ}$ rotation angle Allahyarov and Gompper (2002), and random grid shifting Ihle and Kroll (2001); a Maxwell–Boltzmann rescaling thermostat to maintain constant temperature $T$ Huang et al. (2015); and time $0.1\,\tau$ between collisions, where $\tau=a\sqrt{m/k_{\rm B}T}$, $m$ is the particle mass, and $k_{\rm B}$ is Boltzmann’s constant. The particle density was $\rho=5\,m/a^{3}$, giving a liquid-like Newtonian fluid with kinematic viscosity $\nu=0.79\,a^{2}/\tau$ Ihle and Kroll (2003); Ripoll et al. (2005); Padding and Louis (2006). All simulations were performed with hoomd-blue (version 2.6.0) Anderson, Lorenz, and Travesset (2008); Glaser et al. (2015); Howard, Panagiotopoulos, and Nikoubashman (2018) 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\leq\alpha\leq 3$. (The total number of particles ranged from 3 million to 15 million.) The maximum velocity of the base shear flow was restricted to $U\lesssim 0.32\,a/\tau$ to avoid artificial effects at large Mach numbers Lamura et al. (2001), giving ${\rm Re}\leq 20$.

We considered three force-driven shear flows: a sinusoidal (Kolmogorov) flow Meshalkin and Sinai (1961), a periodic Poiseuille flow Backer et al. (2005), and a periodic Couette-like flow Müller-Plathe (1999); Statt, Howard, and Panagiotopoulos (2019). The first is a historically well-studied problem Meshalkin and Sinai (1961); Green (1974); Grappin, Leorat, and Londrillo (1988); Thess (1992); Bena, Malek Mansour, and Baras (1999) 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 $F_{x}(y)$ to each particle, giving $f_{x}=\rho F_{x}$. (All particles had unit mass $m$, so $F_{x}$ is also the force per particle.) Figure 2 shows $F_{x}$ and $U_{x}$ for the three flows when ${\rm 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 $F_{x}$ 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.

The simulations were initialized by superimposing the base flow on the thermalized fluid and were run for $10^{5}\,\tau$ to reach a steady state. The velocity field $\mathbf{u}(x,y)$ was then measured by averaging the particle velocities in square bins of size $a^{2}$ every $10\,\tau$ for $0.5\times 10^{5}\,\tau$. We used the measured $\mathbf{u}$ to assess whether the targeted base flow $U_{x}$ was stable in the simulations. We adopted three empirical criteria to systematically classify the base shear flow as stable if:

The deviation of the average total kinetic energy in $y$ from equipartition ($\langle E_{y}\rangle=N_{\rm p}k_{\rm B}T/2$ for $N_{\rm p}=(\rho/m)LH^{2}$ particles) was less than 3 standard errors of the mean as computed from its own fluctuations, 2. 2.

the absolute deviation of the flow velocity from the base flow, $|u_{x}(y)-U_{x}(y)|$, was larger than $10^{-3}\,a/\tau$ in fewer than $5\%$ of the bins, and 3. 3.

the absolute value of the gradient velocity $|u_{y}(x)|$ was larger than $10^{-3}\,a/\tau$ 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 $\alpha$ that were visibly stable were properly classified. Criterion (3) tended to be very stringent and often failed even when criteria (1) and (2) held.

IV Results and discussion

We first considered the stability of the sinusoidal shear flow, $\tilde{U}_{x}(\tilde{y})=\sin(\pi\tilde{y})$ (Fig. 3a). Its Fourier coefficients are $\tilde{U}_{\pm 1}=\mp i/2$ and zero otherwise. Secondary flows like those shown in Fig. 1 readily formed in the simulations at sufficiently large Re and $\alpha$. During the accessible simulation time, the flows remained stationary, although subsequent transitions to chaotic motion have been reported Green (1974); Grappin, Leorat, and Londrillo (1988). The bottom panel of Fig. 3a 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 $\alpha\leq 1$ over this range of Re, as proved theoretically by Meshalkin and Sinai Meshalkin and Sinai (1961). However, Fig. 3a also shows that some flows in domains with $\alpha>1$ are also stable for sufficiently small Re, and over the range of $\alpha$ tested, all flows having ${\rm Re}\lesssim 5$ were stable. We confirmed that $\alpha$ 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 consistently with the boundaries in Fig. 3a when interpreted using $\alpha$ and Re.

We subsequently tested two flows commonly used in nonequilibrium molecular simulations: the periodic extension of plane Poiseuille flow (Fig. 3b) Backer et al. (2005) and an approximate periodic extension of plane Couette flow (Fig. 3c) that mimics the RNES flow Müller-Plathe (1999). Both were generated using constant forces in blocks of dimensionless half-width $\tilde{d}=d/H$, where $\tilde{d}=1/2$ for periodic Poiseuille flow and we chose $\tilde{d}=1/25$ to approximate periodic Couette flow (see Fig. 2a and Appendix B). The corresponding Fourier coefficients are

[TABLE]

Because the block forces are only piecewise continuous, $\tilde{U}_{p}\sim 1/p^{3}$ and the Fourier series for $\operatorname{d}\!{}^{2}\tilde{U}_{x}/\operatorname{d}\!{\tilde{y}}^{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 $\alpha$, in excellent agreement with the stability analysis (Figs. 3b-c). In our previous study using RNES Müller-Plathe (1999) to create periodic Couette flow, we found stable shear flows for $\alpha\lesssim 1.25$ Statt, Howard, and Panagiotopoulos (2019); 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. 3a-b). The periodic Couette flow (Fig. 3c) was stable over a larger parameter space than the other two flows (see below).

Given the similar $\alpha$ dependence of ${\rm Re}_{\rm c}$ 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$) Green (1974); Bena, Malek Mansour, and Baras (1999). This drastic simplification allows the eigenvalues of Eq. (6) to be computed analytically:

[TABLE]

$\tilde{\omega}_{0}$ and $\tilde{\omega}_{+}$ will always be negative because ${\rm Re}>0$, but $\tilde{\omega}_{-}$ can be zero or positive. We note that $\tilde{U}_{1}$ and $\tilde{U}_{-1}$ are complex conjugates, $|\tilde{U}_{1}|^{2}=\tilde{U}_{1}\tilde{U}_{-1}$, because $\tilde{U}_{x}$ is real valued. Solving for $\tilde{\omega}_{-}=0$ gives a relationship between the critical Reynolds number ${\rm Re}_{\rm c}$, $\alpha$, and $|\tilde{U}_{1}|$,

[TABLE]

There is a minimum critical Reynolds number as $\alpha\to\infty$; this result is well-known for the sinusoidal flow Meshalkin and Sinai (1961); Green (1974); Bena, Malek Mansour, and Baras (1999), where all flows having ${\rm Re}<\pi\sqrt{2}\approx 4.4$ are stable. However, ${\rm Re}_{\rm c}$ increases as $\alpha$ 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 $\alpha\gtrsim 1.5$ (Fig. 3). This explains the similar stability of the periodic Poiseuille ($|\tilde{U}_{1}|=16/\pi^{3}\approx 0.52$) and sinusoidal ($|\tilde{U}_{1}|=0.5$) flows, as their first Fourier modes are nearly identical, and the stability of the periodic Couette flow to larger Re ($|\tilde{U}_{1}|\to 4/\pi^{2}\approx 0.41$ as $d\to 0$). However, there are some discrepancies between Eqs. (6) and (10) at smaller $\alpha$. In fact, Eq. (6) predicts that periodic Poiseuille flow becomes unstable for ${\rm Re}\gtrsim 50$ when $\alpha=1$, but Eq. (10) diverges as $\alpha\to 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. 1b), although at the largest aspect ratio and Reynolds numbers simulated—$\alpha=3$ and ${\rm Re}=17.5$ and $20$—we found four stationary vortices for the sinusoidal and perioidic Poiseuille flows. We hypothesized that the structure of some of the secondary flows might be connected to the most unstable eigenmode. This eigenmode, $\mathbf{w}=e^{ikx}\mathbf{v}(y)$, can be reconstructed using the Fourier coefficients for $v_{y}$ that comprise the eigenvector having the largest $\Re(\omega)$ and computing $v_{x}=(i/k)\operatorname{d}\!{v_{y}}/\operatorname{d}\!{y}$ based on the incompressibility of $\mathbf{u}$ (Appendix A). Figure 4a shows $\mathbf{w}$ for the periodic Couette flow at ${\rm Re}=15$ and $\alpha=1.4$, which develops into the flow shown in Fig. 1b. It has two pairs of counterrotating vortices with streamlines primarily directed along the shear gradient $y$.

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, $\mathbf{u}\approx c_{1}\mathbf{U}+c_{2}\mathbf{w}$, 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 $\mathbf{u}$ for the conditions in Fig. 4, determining optimal coefficients $c_{1}=0.765$ and $c_{2}=0.373$. The fitted flow field (Fig. 4b) bears striking similarity to the simulated flow (Fig. 1b), having a dimensionless root-mean-squared error of $0.02$ per velocity component.

We also noted that the waiting time $\tau_{\rm w}$ for the secondary flows to emerge in the simulations was shorter for points farther from the curve ${\rm Re}_{\rm c}(\alpha)$. 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 $\tau_{\rm w}$ for the unstable flows by empirically fitting $E_{y}$ to a hyperbolic tangent during the first $10^{5}\,\tau$ simulated; we defined $\tau_{\rm w}$ as the first time that $E_{y}$ increased by 10% of the difference between its initial and final values (Appendix C). These times can be compared to $1/\Re(\omega)$, which is the typical timescale associated with growth of the most unstable eigenmode. We found that $\tau_{\rm w}$ had a strong linear correlation with $1/\Re(\omega)$ 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\times 10^{4}\,\tau$ in dimensional units. Having this estimate of the timescale for instability is an added benefit of our analysis.

V Conclusions

The stability of force-driven shear flows simulated in fully periodic domains seems to stand in stark contrast to those in bounded ones. Plane Couette flow is linearly stable for all infinitesimal disturbances Romanov (1973), while plane Poiseuille flow is linearly stable up to Re = 5772 Orszag (1971); in practice, both become unstable in the range of ${\rm Re}\approx 10^{2}$ to $10^{3}$ Orszag and Kells (1980); Bayly, Orszag, and Herbert (1988); Lundbladh and Johansson (1991). However, their periodic extensions (Figs. 3b-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 Statt, Howard, and Panagiotopoulos (2019). 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 $\alpha<1$. The presence of a hydrodynamic instability in nonequilibrium simulation methods like 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.

Acknowledgements.

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.

Data Availability

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.Statt et al. (2019)

Appendix A Linearization

Although the derivation of Eq. (3) is widely documented Lin (1955); Drazin and Reid (2004), we include details here for completeness. Substituting $\mathbf{u}$ and $p$ into Eqs. (1) and (2), using that $\mathbf{U}$ must also be a solution of the same, and neglecting terms of $\mathcal{O}(\delta\mathbf{u}^{2})$ yields the linearized equations:

[TABLE]

Eq. (11) gives a simple relation between $v_{x}$ and $v_{y}$,

[TABLE]

where the prime denotes the ordinary derivative, $v_{y}^{\prime}=\operatorname{d}\!{v_{y}}/\operatorname{d}\!{y}$, while Eq. (12) gives

[TABLE]

Eq. (13) can be inserted in Eq. (14) and solved for $q$:

[TABLE]

Substituting Eq. (16) into Eq. (15) results in Eq. (3), which we note is the well-known Orr–Sommerfeld equation Drazin and Reid (2004) after replacing $\omega$ by $i\omega$.

Appendix B Flow fields

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

The sinusoidal flow,

[TABLE]

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

[TABLE]

We chose $U$ to obtain a certain Re, setting the force amplitude $F=(\pi\nu)^{2}{\rm{Re}}/H^{3}$.

Periodic extensions of plane Poiseuille flow and 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=\pm H/2$ of half-width $d\leq H/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:

[TABLE]

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:

[TABLE]

with

[TABLE]

$g_{1}$ is the usual linear (Couette-like) flow, and $g_{2}$ is the quadratic (Poiseuille-like) flow. The maximum velocity is

[TABLE]

For the periodic plane Poiseuille flow Backer et al. (2005), the width of the blocks was chosen to cover the full simulation box, $d=H/2$. Only $g_{2}$ contributes to this flow, establishing two opposing parabolic regions, and the force magnitude required for a given Reynolds number is $F=8\nu^{2}{\rm Re}/H^{3}$.

By making $d$ much smaller than $H/2$, the flow is dominated by $g_{1}$ 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\nu^{2}{\rm{Re}}/\left[H^{3}(1-(1-2d/H)^{2})\right]$. In the limit $d\to 0$, the shear stress has a step change at $y=\pm 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 method Müller-Plathe (1999) even though the origin of the perturbation is different.

Appendix C Timescale fitting

In order to estimate the timescale for occurrence of the instability, we empirically fit the kinetic energy $E_{y}$ to a hyperbolic tangent,

[TABLE]

We fixed $E_{0}=1/2$ because the system was initially thermalized, and so its kinetic energy must obey equipartition. We then fit the final value $E_{1}$, the midpoint $\tau_{0}$, and the width $w$, resulting in fits like those shown in Fig. 6. The waiting time was defined as $\tau_{\rm w}=\tau_{0}-1.09861w$ using the “10–90 thickness.” We found it challenging to reliably fit $E_{y}$ for a handful of conditions close to the stability curves of Fig. 2; we accordingly neglected all points having $1/\Re(\tilde{\omega})>10$ in creating and analyzing Fig. 4.